Millennial-Scale Solar Variability in Tree Rings of Northern Fennoscandia at the End of the Holocene

ABSTRACT To investigate the possible Sun–climate connection during the Holocene, the Finnish super-long tree-ring chronology covering the period from 5634 B.C. to A.D. 2004 was analyzed. As an indicator of solar activity, we used a reconstruction of total solar irradiance (TSI) covering 9300 years, which is based on a composite using the cosmogenic radionuclide 10Be measured in polar ice cores, and also on neutron monitor data (Steinhilber et al. 2009). The Multiple Taper Method (MTM) wavelet decomposition and wavelet coherence analyses were applied to the time-series. The MTM spectral analysis identifies the main solar cycles at ca. 200 yr (de Vries or Suess), ca. 350 yr (unnamed) and ca. 900 years (Eddy). The strongest cross-wavelet correlation was discovered between the millennium-cycle components of TSI and tree-ring width variations. This Eddy cycle, which was recently discovered in solar activity, remains both strong and stable through almost the entire Holocene, and it reappears again at lower frequencies (ca. 1300 years) after ca. A.D. 200. Our results raise questions regarding the end of the Holocene and transition to the next glacial period and confirm the complex and nonlinear nature of the Sun–climate relationship during the Holocene Epoch.


INTRODUCTION
The nature of the Sun's role in current and future climate variability is now a hot topic of discussion. The increase in global surface temperature by 1.09°C over the 20 th Century is largely attributed to anthropogenic inputs of atmospheric greenhouse gases (IPCC 2021). However, some investigations demonstrate that solar variations are responsible for as much as 50% of the past century global warming (Scafetta and West 2006) and 75% of Arctic climate warming (Soon 2005). Furthermore, there is growing evidence that the Sun is one of the main drivers of the Earth's climate through the Holocene Epoch, which began about 12,000 years ago (Neff et al. 2001;Hu et al. 2003;Gupta et al. 2005;Helama et al. 2010;Kern et al. 2012;Steinhilber et al. 2012;Duan et al. 2014;Soon et al. 2014;Jiang et al. 2015;Zhao et al. 2020).
Therefore, realistic modeling of climate scenarios requires the inclusion of variable solar forcing as input.
Nowadays, the following physical mechanisms of solar forcing on the Earth's climate are considered: (1) direct heating of the Earth's atmosphere by the total solar irradiance (TSI) (Lean et al. 1995); (2) the effect of solar ultraviolet radiation (UV) on stratospheric chemistry, ozone abundance and thermal structure (Haigh 1996); and (3) the effect of galactic cosmic rays (GCR) on the cloud cover of the low atmosphere (Marsh and Svensmark 2000;Palle and Butler 2001).
Unfortunately, the directly measured data for these parameters only cover timescales of centuries. For example, the direct globally averaged surface Earth's temperature records usually cover only the last ca. 140-170 years (Zhao et al. 2020). Note that one of the most important indicators of solar variability is the sunspot number (SN), which is the solar parameter with the longest record (regular telescopic observations are available since A.D. 1610; see http://www.sidc.be/silso/). Satellite measurements of solar irradiance became available only since 1978 (Shapiro et al. 2011). Therefore, reconstructions are required to reveal a connection between solar activity and Earth's climate on centennial to millennial timescales.
Cosmogenic isotopes 10 Be in ice cores and 14 C in tree rings are known to be reliable indicators of long-term changes of GCR and solar activity Braziunas 1989, 1993;Beer 2000;Land et al. 2020;Brehm et al. 2021). The proxy records of these isotopes cover time intervals of many millennia and can be used for reconstructions of solar activity prior to the era of instrumental records (Solanki et al. 2004;Vonmoos et al. 2006;Steinhilber et al. 2009;Shapiro et al. 2011). Annually resolved tree-ring records obtained by dendrochronological methods are the best sources of paleoclimatic information (Cook and Kairiukstis 1990). These proxies can provide long, precisely dated climatic records of the past.
Comparison of super-long (5634 B.C.-A.D. 2007) temperature reconstruction based on tree rings from Northern Fennoscandia with sunspot numbers inferred from 14 C data (Solanki et al. 2004) have shown that their long-term changes are modulated by corresponding centennial-to-millennialscale solar variations (Helama et al. 2010). Sunspots are dark areas with the strong concentrations of magnetic flux at the solar surface. They are surrounded by bright regions (faculaes) that allow their use for reconstruction of TSI variations in the past (Lean et al. 1995;Solanki et al. 2004). However, there is some evidence that GCR intensity, as seen in 10 Be records from Greenland ice cores, varied with a 22-year periodicity during the Maunder minimum of solar activity (A.D. 1645-1715), when the sunspots were practically absent (Beer et al. 1998;Usoskin et al. 2001;Steinhilber et al. 2009). The Maunder minimum was a prolonged episode of low solar activity which coincided with more severe winters in the United Kingdom and Europe (Lockwood et al. 2010). Moreover, recent results identified clear periodicities in solar-tree-ring connections during and around the Grand Minima of solar activity (Kasatkina et al. 2019). These results indicate that the Sun's behaviour during Grand Minima is not just a lack of sunspots, but more complex (Steinhilber et al. 2009). The most recent 10 Be-based reconstruction of solar activity has been made by Steinhilber et al. (2009). This reconstruction represents another approach that is based on an observationally derived relationship between TSI and the open solar magnetic field (Steinhilber et al. 2009).
As noted above, the manifestations of solar cycles in climatic variations from decadal to centennial time scales have already been studied quite extensively, so in this paper we will focus on studying variations in tree growth over centennial time scales. In this paper, the super-long tree-ring (ca. 7500) chronology from Northern Fennoscandia is compared to the TSI reconstruction based on 10 Be proxies (Steinhilber et al. 2009). Here our main objective is using multi-taper spectral and wavelet-based analysis to find evidence of the possible linkages, if any, among tree growth and solar activity in Arctic during the Holocene.

MATERIAL AND METHODS
To investigate the possible Sun-climate relationship during the Holocene, the super-long Scots pine (Pinus sylvestris L.) tree-ring chronology covering the period of ca. 7600 years from 5634 B.C. to A.D. 2004A.D. (Helama et al. 2008) was analyzed. The chronology consists of 1249 tree-ring-width series including samples collected from living and dead trees, and from subfossil wood from small lakes in the area of subarctic timberline in northern Fennoscandia (68°-70°N, 20°-30°E). All series were dendrochronologically processed and dated to calendar years (Helama et al. 2008).
As an indicator of solar activity, we've used a reconstruction of total solar irradiance (TSI) covering 9300 years, which is based on a composite using the cosmogenic radionuclide 10 Be measured in ice cores GRIP (Greenland) and South Pole (Antarctic), and also on neutron monitor data (Steinhilber et al. 2009). This reconstruction is based on the observationally derived relationship between TSI and the open solar magnetic field that can be obtained from 10 Be (Steinhilber et al. 2009). The 10 Be abundance in polar ice depends on the intensity of the GCR flux modulated by solar magnetic and geomagnetic fields, which can enter the Earth's atmosphere (Beer 2000). In atmosphere, GCR particles interact with oxygen and nitrogen atoms producing 10 Be, and therefore 10 Be records in polar ice are reliable indicators of GCR flux and solar irradiance changes over the Holocene (Beer 2000;Steinhilber et al. 2009). The data were 40-year running means and were resampled to a 5-year time resolution (Steinhilber et al. 2009).
For further analysis, these two records were homogenized ( Figure 1). In the present study, we applied the Multiple Taper Method (MTM) of spectral analysis to identify the main periodicities in these time-series. The MTM is a non-parametric method that avoids some limitations associated with the use of an a priori model of the process and applies multiple orthogonal windows (or tapers) to minimize the spectral leakage caused by the finite length of the data set (Thomson 1982). The MTM spectral analysis was performed with the help of the SSA-MTM software toolkit (Ghil et al. 2002). Significance levels were tested against a red-noise background (Ghil et al. 2002).
Wavelet analysis was performed to analyze the main periodicities and their evolution in timefrequency space. This method is a powerful tool to study non-stationary signals that can vary both in amplitude and frequency with time (Torrence and Compo 1998). The Morlet function was applied in this work as being the most fruitful for geophysical time-series analysis (Grinsted et al. 2004). The Morlet wavelet was defined as plane wave modulated by a Gaussian function (Torrence and Compo 1998): where s is the wavelet scale, w 0 was dimensionless frequency, here taken to be 6 to provide a good balance between time and frequency localization (Torrence and Compo 1998). Cross-wavelet transform was also performed between TSI reconstruction and tree-ring data (Torrence and Compo 1998;Grinsted et al. 2004): where W 1 and W 2 denote continuous wavelet transforms of tree-ring and TSI time-series, respectively; * denotes the complex conjugate.
The wavelet coherence (WTC) was calculated to evaluate a localized correlation coefficient in time frequency space between TSI and tree-ring timeseries (Grinsted et al. 2004): where S denotes smoothing in both time and scale. The statistical significance levels of wavelet transform and WTC were estimated using the Monte Carlo method against a red noise model (Torrence and Compo 1998;Grinsted et al. 2004).
To evaluate fluctuations in power over a frequency band s 1 to s 2 , the scale-averaged wavelet power was calculated (Torrence and Compo 1998): where δt is the sampling interval, δj denotes the factor for scale averaging, C = 0.776 for the Morlet wavelet. This scale-averaged wavelet power can be used to examine modulation of one frequency by another within the same time-series.
To assess the periodicities in low-frequency bands, the discrete wavelet decomposition (Farge 1992) of the original time-series was performed. It decomposes the original signal in approximations (A) and in details (D) containing the frequency subbands limited by numbers 2 n and 2 n+1 (where n related to the time-series step). Using the orthogonal discrete Meyer wavelet, we performed an n = 9-level decomposition. Because the variable characteristics of long-term fluctuations for solar activity and tree growth over centennial time scales are studied in this work, the wavelet decomposition was carried out starting from a period of 160 years. Therefore, three frequency bands were selected for further analysis: D 5 (160-320 years), D 6 (320-640 years) and D 7 (640-1280 years). We computed the non-linear Spearman correlation coefficients between tree-ring widths and TSI time-series in these frequency bands. All the wavelet results in this paper were done for the standardized (zero mean, unit standard deviation) time-series.

RESULTS
It must be stressed that the data consist of averaged values, where the averaging was done over a 40-year interval. Therefore, the results are only indicative at longer-term variations and do not include and represent solar cycles from decadal to centennial time scales. The MTM spectra of super-long tree-ring (FTR) and TSI reconstruction records are shown in Figure 2. Spectral analysis of FTR shows periodicities of ca. 200, 350-450 years, and ca. 900 years above the 90% confidence level (see Figure 2b), suggesting cyclic changes at multicentennial scales.
The WTC analysis between tree-ring width and TSI clearly illustrates a relationship between these time-series. It shows a statistically significant (or close to the significant level) high power in frequency bands indicating solar cycles: D 5 (160-320 years), D 6 (320-640 years) and D 7 (640-1280 years). However, it should be noted that the periods in D 5 and D 6 frequency bands are intermittent, and the phase angles are very variable in time, indicating a non-linear connection between the variables. Figure  3 displays the level D 5 (160-320 years) of decomposition ( Figure 3a) and the wavelet coherence plot of the FTR-TSI data for this frequency band ( Figure  3b). This level of decomposition encompasses the ca. 200-year de Vries solar cycle. The same plots are presented in Figures 4 and 5 for the levels D 6 (320-640 years) and D 7 (640-1280 years). These levels encompass the ca. 350-year and ca. 900-yr solar cycles, respectively. The correlation is sufficient and significant (r = 0.6, p < 0.001) only for the millennial-cycle components.
Variability of tree rings within period of 160-320 years demonstrates intermittent oscillatory characteristics throughout the large portion of record (see Figure 3a). The de Vries solar cycle with intensity above the 95% confidence level intermittently appears in the FRT-TSI connection around 3000 B.C. and after ca. A.D. 1000 until now (see Figure 3b). Black arrows in Figure 3b indicate the phase angles between the time-series (with in-phase pointing right, anti-phase pointing left). The signal is weak (below the 95% confidence level) around 5500-4500 B.C. and A.D. 500-1000. Arrows indicate that the two series have a non-stationary relative phase connection during these time intervals (see Figure 3b). The FTR-TSI connection in the D 6 frequency band is also irregular and shows a statistically significant (near/or above 95% confidence level) power for the time intervals around 3500-2000 B.C., 1000-100 B.C. and after ca. A.D. 1000, but almost outside of the COI (see Figure 4b). The signal demonstrates a non-stationary phase behavior. The D 7 band that encompasses the ca. 900-year Eddy cycle demonstrates a high (above the 95% confidence level) power and almost in-phase relationship over a period of 5634 B.C.-A.D. 1000 (see Figure 5b). Note the interval around of 5634-4800 B.C. is out of the cone of influence (COI, see Figure 4b). The wavelet coherence in this frequency band (ca. 900 years) demonstrates also a high (above 95% confidence level) power through 5634-1500 B.C. with an in-phase relationship, and it reappears again at lower frequencies (ca. 1280 years) after ca. A.D. 200 (see Figure 5b). Figure 5a shows that the millennial cycle in tree rings steadily decreases in amplitude throughout the whole period of record. The calculation indicates that periodic length of the ca. 900-year period in tree rings mainly changes through 800 to 1000 years, and just exceeds 1000 years in a short-time span before ca. A.D. 1000 (see Figure 5a). To quantitatively estimate the power of the millennial-year oscillation signals in TSI and tree rings at each time epoch, scale-averaged wavelet power values over all scales between 800 and 1300 years are constructed, and the results are plotted in Figure 6. Both time-series show a similar evolution with a significant (above 95% confidence level) maximum around 5000 B.C. and a subsequent gradual decrease in power, and they are in a state with a low wavelet power (<95% confidence level) ca. 3500 B.C.-A.D. 1000 (see

DISCUSSION
Spectral, wavelet decomposition and wavelet coherence analysis of the super-long tree-ring chronology from northern Fennoscandia revealed pronounced periodicities in tree-ring growth, corresponding to solar variations at multi-centennial to millennial scales. All these periods may be attributed to variations of solar activity (Figure 2a), and several of them coincide with well known solar cycles such as the Suess or de Vries cycle (ca. 200 years) and the Eddy cycle (ca. 1000 years) (Ma 2007;Vasiliev and Dergachev 2008;Kern et al. 2012;Steinhilber et al. 2012;Hanslmeier et al. 2013;Soon et al. 2014;Zhao et al. 2020). The unnamed periodicities of ca. 350 years have been documented in solar proxy records (Vasiliev and Dergachev 2008;Kern et al. 2012;Steinhilber et al. 2012).
These results agreed with previous findings indicating a solar-climate connection at multicentennial time scales at ca. 200 and ca. 350 years (Wang and Zhang 2011;Breitenmoser et al. 2012;Kern et al. 2012;Steinhilber et al. 2012;Czymzik et al. 2016;Rimbu et al. 2020). The strongest crosswavelet correlation was discovered between the millennium-cycle components of TSI and tree-ring width variations. That seems to indicate that the ca. 900-year solar cycle has a greater impact on tree growth and climate in Northern Fennoscandia than the multidecadal to multi-centennial solar cycles.  (Kasatkina et al. 2019). There is some evidence of non-linearity of the solar-climate connections at decadal to centennial time scales with complex feedback processes in the ocean-atmosphere system, which can also amplify the response to solar irradiance variations (Kasatkina et al. 2006Scafetta and West 2006;White and Liu 2008;Gusev and Martin 2012;Hanslmeier et al. 2013, Mares et al. 2021). In addition, volcanic activity may reinforce or mask the solar signal in space and time (Shumilov et al. 2011;Breitenmoser et al. 2012;Kasatkina et al. 2018). Otherwise, another interpretation has been given by Hanslmeier et al. (2013) who indicated that solar proxies seem to exhibit a more regular and predictable behavior at larger time scales (well above 100 years).
The ca. 900-year cycle identified in the solartree-ring connection is similar to the mean 1370 ± 500-year periodicity associated with ice-rafted debris coinciding with rapid climate changes during the Holocene (Bond et al. 2001). A millennial periodicity was previously detected also by Helama et al. (2010), who performed a comparison of solar activity and summer temperature records derived from these tree-ring data. The authors used Figure 6. Scale-averaged wavelet power over the 800-1300-year band for the tree-ring width (dashed) and TSI (solid) time-series. The thin solid line is the 95% confidence level for tree-ring series, and the thin dashed line is the 95% level for the TSI (assuming red noise). the sunspot reconstruction of Solanki et al. (2004) as a proxy of solar activity (Helama et al. 2010). To achieve a better understanding of this solar pattern, we used the TSI reconstruction based on 10 Be data by Steinhilber et al. (2009). This more recent reconstruction, compared to that of Solanki et al. (2004), represents another approach that is based on modulation of cosmic rays by solar magnetic field and on an observationally derived relationship between TSI and the open solar magnetic field (Steinhilber et al. 2009). This Eddy cycle, which was recently discovered in solar activity (Ma 2007;Yin et al. 2007;Dima and Lohmann 2009;Vasiliev and Dergachev 2009;Hanslmeier et al. 2012;Steinhilber et al. 2012;Soon et al. 2014), remains both strong and stable during the mid-Holocene and almost dissipated during the past two millennia. Figures 5a and 6 demonstrate that the thousand-year cycles in tree rings and TSI reach maximum power at ca. 5000 B.C. with a subsequent constant decrease, thus coinciding with a gradual temperature decline occurred after the Holocene Climatic Optimum (ca. 9000 to 5500 B.P.) (Kalis et al. 2003). Furthermore, the periodic length of this cycle in tree rings is shorter in early times and is slightly longer in the present (see Figure 5a). Both series show similar behavior, indicating a possible modulation with a ca. 6500-year period (see Figure 6). A ca. 6500-year quasi period has also been previously found in solar sunspot number proxy record (Yin et al. 2007;Dima and Lohmann 2009). However, the time-series presented may be not long enough to show a ca. 6500year period of solar or geomagnetic origin.
As mentioned above, the following solar factors acting on atmosphere and climate are considered: total solar irradiance, UV-radiation and GCR (Lean et al. 1995;Haigh 1996;Marsh and Svensmark 2000;Palle and Butler 2001). Cosmic rays influence climatic patterns through cloud formation in the low atmosphere, and consequently, changing the radiation balance and hence higher cosmic ray flux can therefore be associated with lower solar irradiance, and vice versa. In turn, solar radiation reaching the Earth's surface influences the tree's growth either directly (sunlight) or indirectly through temperature variations. Temperature and sunlight seem to be main growth-limiting factors for trees living at high latitudes in extreme cold conditions (e.g. close to northern timberline). At high latitudes, these factors are more essential for photosynthesis (Muraki et al. 2015). Photosynthetically active solar radiation contains wavelengths between 400 nm (visible light) and 700 nm (nearinfrared range) (Fitter and Hay 2002). Recently, Kasatkina et al. (2019) suggested photosynthetically active solar spectral solar irradiance is an important agent of solar activity responsible for tree growth beyond the Arctic Circle, especially during and around the Grand minima of solar activity. This can also be indirectly confirmed by the recent Solar Radiation and Climate Experiment (SORCE) satellite data according to which the observed solar radiation flux in the visible and near-infrared spectral ranges increased while the solar activity level decreased (Harder et al. 2009). That is, the trees may continue an increase in growth even in the case of future decline in the Sun's activity.
According to some results, the intensity of low-frequency spectral peaks in solar proxy records exceeds the spectral intensities for the 11-year solar cycle by a factor more than 20 (Shindell et al. 2001;Vasiliev and Dergachev 2009). For example, a 0.1% decrease in TSI related to the 11-year solar cycle could lead to changes in stratospheric ozone and temperature, producing changes in surface climate (Haigh 1996). On centennial scales, ca. 0.25% reduction (or about 3-4 W/m 2 less than at present) of TSI during the Maunder minimum of solar activity resulted in cold snaps in the Northern Hemisphere (Shindell et al. 2001;Lockwood et al. 2010). Moreover, the increase in solar irradiance in separate spectral bands from the Maunder Minimum to the present was much larger than 0.25% for TSI (Shapiro et al. 2011). The highest wavelet coherence between TSI and tree rings at ca. 900 years seems to indicate that the solar irradiance has a greater impact on tree growth and climate in Northern Fennoscandia on the millennial scale than on the multidecadal to multi-centennial scales. According to the solar dynamo theory, these long-term quasiregular oscillations in the Sun's magnetic and radiative outputs can result from the nonlinear interplay between the magnetic field, differential rotation and convective motions (Soon et al. 2014).
The results presented support an interpretation that solar forcing seems to regulate, at least in part, multi-centennial and millennial cyclicity in the polar tree-ring growth during the Holocene.

CONCLUSIONS
We have used MTM, wavelet decomposition and wavelet coherence analysis of tree-ring and solar activity proxies to examine the relationship between these time-series at multi-centennial-tomillennial time scales throughout the Holocene. The proxies used represent the super-long (ca. 7600 years) tree-ring chronology from northern Fennoscandia and the total solar irradiance reconstruction that is based upon the variations of 10 Be in ice core, Greenland.
Summarizing the obtained results, we conclude that: 1. Our analysis identifies a significant solar-treering (climatic) connection at multi-centennial and millennial timescales around 200, 350 and 900 years. All these periodicities coincide with well known solar cycles: ca. 200 years (the de Vries cycle), ca. 350 years (the unnamed cycle) and ca. 900 years (the Eddy cycle). 2. Periodicities at ca. 200 and ca. 350 years are intermittent in time, indicating a non-linear connection between the variables. Compared with other periodicities, the millennial-year cycle in the solar-tree-ring connection is stable and significant (above 95% confidence level) during most of the Holocene (throughout 5634-1500 B.C.), and it reappears again at lower (ca. 1300 years) frequencies after ca. A.D. 200. Could this imply the end of the Holocene and transition to the next glacial period? The answer to this question requires additional studies of millennialscale variations in the Sun-climate relationships using different reconstructions of solar activity and climate changes.
Our new data provide additional and robust evidence that solar forcing seems to regulate, at least in part, multi-centennial and millennial cyclicity in the polar tree-ring growth during the Holocene. In any case, the results presented confirm the evident, but complex and nonlinear nature of the Sunclimate relationship during the Holocene epoch.

ACKNOWLEDGMENTS
The authors would like to thank the anonymous reviewers and the Associate Editor for their critical comments and valuable suggestions that led to a significant improvement of this paper.