Issue 
A&A
Volume 615, July 2018



Article Number  A93  
Number of page(s)  13  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201731892  
Published online  19 July 2018 
Solar activity over nine millennia: A consistent multiproxy reconstruction^{★}
^{1}
MaxPlanckInstitut für Sonnensystemforschung,
Göttingen, Germany
^{2}
Space Climate Research Unit, University of Oulu,
Finland
email: Ilya.Usoskin@oulu.fi
^{3}
Sodankylä Geophysical Observatory, University of Oulu, Finland
^{4}
Ioffe PhysicalTechnical Institute,
194021
St. Petersburg, Russia
^{5}
CEREGE, AixMarseille Université, CNRS, NRD, INRA, Collège de France, Technopôle de l’Arbois,
AixenProvence, France
^{6}
School of Space Research, Kyung Hee University,
Yongin,
GyeonggiDo
446701, Republic of Korea
Received:
4
September
2017
Accepted:
23
March
2018
Aims. The solar activity in the past millennia can only be reconstructed from cosmogenic radionuclide proxy records in terrestrial archives. However, because of the diversity of the proxy archives, it is difficult to build a homogeneous reconstruction. All previous studies were based on individual, sometimes statistically averaged, proxy datasets. Here we aim to provide a new consistent multiproxy reconstruction of the solar activity over the last 9000 yr, using all available longspan datasets of ^{10}Be and ^{14}C in terrestrial archives.
Methods. A new method, based on a Bayesian approach, was applied for the first time to solar activity reconstruction. A Monte Carlo search (using the χ^{2} statistic) for the most probable value of the modulation potential was performed to match data from different datasets for a given time. This provides a straightforward estimate of the related uncertainties. We used six ^{10}Be series of different lengths (from 500–10 000 yr) from Greenland and Antarctica, and the global ^{14}C production series. The ^{10}Be series were resampled to match wiggles related to the grand minima in the ^{14}C reference dataset. The stability of the long data series was tested.
Results. The Greenland Icecore Project (GRIP) and the Antarctic EDML (EPICA Dronning Maud Land) ^{10}Be series diverge from each other during the second half of the Holocene, while the ^{14}C series lies in between them. A likely reason for the discrepancy is the insufficiently precise beryllium transport and deposition model for Greenland, which leads to an undercorrection of the GRIP series for the geomagnetic shielding effect. A slow 6–7 millennia variability with lows at ca. 5500 BC and 1500 AD in the longterm evolution of solar activity is found. Two components of solar activity can be statistically distinguished: the main component, corresponding to the “normal” moderate level, and a component corresponding to grand minima. A possible existence of a component representing grand maxima is indicated, but it cannot be separated from the main component in a statistically significant manner.
Conclusions. A new consistent reconstruction of solar activity over the last nine millennia is presented with the most probable values of decadal sunspot numbers and their realistic uncertainties. Independent components of solar activity corresponding to the main moderate activity and the grandminimum state are identified; they may be related to different operation modes of the dynamo.
Key words: Sun: activity / solarterrestrial relations / sunspots
A table with the reconstructed SN series is available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr(130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/615/A93
© ESO 2018
1 Introduction
The Sun is an active star whose magnetic activity varies on different timescales, from seconds to millennia. Understanding solar variability in detail is important for many reasons, ranging from applications in stellar astrophysics and dynamo theory to paleoclimatic and space weather studies. Various direct, partly multiwavelength spectroscopic solar observations cover the past few decades up to a century in the past. They provide knowledge of the solar variability that is expressed in different indices. For the preceding centuries back to 1610 AD, only visual information from simple optical observations in the form of sunspot numbers (SNs) is available (Hathaway 2015). The quality of the SN series varies, and the uncertainty increases sufficiently before the nineteenth century (Clette et al. 2014; Usoskin 2017). However, the Maunder minimum in the second half of the seventeenth century was observed sufficiently well to conclude that sunspot activity was exceptionally low (Ribes & NesmeRibes 1993; Vaquero et al. 2015; Usoskin et al. 2015). For earlier times, only indirect proxies can help assessing solar activity.^{1} Such proxies are, for example, concentrations of cosmogenic radionuclides radiocarbon (^{14}C), beryllium10 (^{10}Be), or chlorine36 (^{36}Cl) that are measured in tree trunks or polar ice cores, respectively. These archives are dated independently. The use of cosmogenic proxies for studying the solar activity in the past has been proposed long ago (e.g., Stuiver 1961; Stuiver & Quay 1980; Beer et al. 1988), and the method has been developed since then in both measurements and modeling (see reviews by Beer et al. 2012; Usoskin 2017, and references therein). These data cover timescales of up to ten millennia and more.
Cosmogenic radionuclides are produced as a byproduct of a nucleonic cascade initiated by galactic cosmic rays (GCR) in the Earth’s atmosphere. This is the only source of these nuclides in the terrestrial system (which is why they are called “cosmogenic”). Radiocarbon ^{14} C is mainly produced as a result of neutron capture (np−reaction) by nitrogen, which is responsible for > 99% of the natural ^{14}C production. The isotope ^{10}Be is produced as a result of spallation reactions of O and N nuclei caused by energetic GCR particles. After production, the two radionuclides have different processes of transport, deposition, and storage in terrestrial archives around the globe. Radiocarbon is mostly measured in dendrochronologically dated annual rings of live or dead tree trunks, while ^{10} Be is measured in glaciologically dated polar ice cores mostly from Greenland or Antarctica. In addition to the geomagnetic field, the flux of GCRs impinging on Earth is modulated by largescale heliospheric magnetic features (interplanetary magnetic fields and solar wind) so that the measured content of the nuclides may serve as a proxy of solar magnetic variability in the past (after the influence of the geomagnetic field has been removed). Different cosmogenic isotope series exhibit a high degree of similarity on timescales of a century to a millennium because they share the same production origin (Bard et al. 1997; Vonmoos et al. 2006; Beer et al. 2012; Usoskin et al. 2009, 2016; Delaygue & Bard 2011; Adolphi & Muscheler 2016). However, systematic discrepancies between longterm (multimillennial) trends in different series can be observed (Vonmoos et al. 2006; Inceoglu et al. 2015; Adolphi & Muscheler 2016; Usoskin et al. 2016), probably because of the influence of climate conditions (regional deposition pattern for ^{10} Be or largescale ocean circulation for ^{14}C) or an improper account for the geomagnetic field variation, to which the two isotopes respond differently.
Because of the differences in the isotope records, earlier reconstructions of solar activity were obtained based on individual cosmogenic series, leading to a diversity in the results. Earlier multiproxy efforts (e.g., Bard et al. 1997; McCracken et al. 2004; Vonmoos et al. 2006; Muscheler et al. 2007; Usoskin et al. 2007; Knudsen et al. 2009) were mostly based on a simple comparison of individual records. The first consistent effort to produce a merged reconstruction was made by Steinhilber et al. (2012), who used the principle component analysis (PCA) to extract the common variability signal (assumed to be solar) from the reconstructions based on three cosmogenic series (global ^{14} C, Antarctic EDML (EPICA Dronning Maud Land), and the Greenland Icecore Project, GRIP, ^{10} Be) and to remove the system effects (e.g., the deposition process, snow accumulation rate, and changes in the carbon cycle and dating uncertainties), which are different for each series. The PCA method keeps only the relative variability and looses the information on the absolute level, which needs further normalization. Moreover, this method effectively averages multiple signals without taking the accuracy of each data point and possible time lags between the signals into account (see Figs. S1–S8 in Steinhilber et al. 2012). As discussed by Adolphi & Muscheler (2016), the time mismatch between the ^{14} C and ^{10} Be series may be as large as 70 yr toward the early Holocene, however.
We here introduce a new method for a consistent multiproxy reconstruction of the solar activity that is based on the Bayesian approach to determine the most probable value (and its uncertainties) of the solar activity at any moment in time by minimizing the χ^{2} discrepancy between the modeled and the actually measured cosmogenic isotope data. This method straightforwardly accounts for error propagation and provides the most probable reconstruction and its realistic uncertainties. Since the method is sensitive to the dating accuracy of different records, we redated the ^{10} Be records to match their dating with that of ^{14}C by applying the standard wigglematching method (Cain & Suess 1976; Muscheler et al. 2014) to the official icecore chronology, as described in Sect. 3. The solar modulation potential was assessed using the Bayesian approach from all the available datasets (Sect. 4). Finally, the SN series was calculated (Sect. 5). Section 6 summarizes our results and conclusions.
Fig. 1
Data series (see Table 1) and temporal coverage. Blue shading denotes the cosmogenic isotope series, while green shading shows the geomagnetic series. 
2 Data
2.1 Cosmogenic isotope records
We here used six ^{10}Be series from Greenland and Antarctica and one global ^{14}C production series, as summarized in Table 1 and Fig. 1.
The ^{14}C series covers the entire Holocene with homogeneous resolution as given by the globally averaged ^{14} C production rate computed by Roth & Joos (2013) from the original International ^{14} C Calibration dataset INTCAL09 (Reimer et al. 2009) measurements of Δ^{14}C. While this series has a pseudoannual temporal sampling (Roth & Joos 2013), its true resolution is decadal. The series is provided as an ensemble of 1000 realizations of individual reconstructions. Each realization presents one possible reconstruction in the sense of a Monte Carlo approach (viz, one realized path of all possible paths in the parametric space), considering all known uncertainties. This enables directly assessing the error propagation throughout the entire process.
^{10}Be series have different coverage and temporal resolutions. We reduced them to two cadences. Annual series were kept as they are, while rougher resolved series were resampled to decadal cadence. For each series, we also considered 1000 realizations that were synthesized using the mean curve and the standard deviation (error bars). Two series were updated with respect to earlier studies. TheEDML series was used with the new Antarctic Ice Core Chronology, version 2012 (AICC2012; Veres et al. 2013; Bazin et al. 2013). The GRIP series was updated from its original ss09(sea) timescale (Johnsen et al. 1995, 2001) to the more recent timescale of the Greenland Ice Core Chronology, version 2005 (GICC05; Vinther et al. 2006).
All ^{10}Be series were converted into units of production/flux [atoms cm^{−2} s^{−1}] when possible, which is natural for the isotope production by cosmic rays. Originally, ^{10} Be measurements are given in units of concentration [atoms g^{−1}]. To convert them into the flux, the independently measured or obtained snow accumulation rate at each site is needed. The accumulation rate was considered individually for each ice core, using the same chronology as for the ^{10}Be data (Veres et al. 2013; Bazin et al. 2013 for EDML on AICC2012 and Rasmussen et al. 2014; Seierstad et al. 2014 for GRIP on the GICC05 timescale). For the Dye3 and South Pole series, accumulation data are not provided, and in these cases, we used concentration data assuming that the concentration is proportional to the depositional flux, with the scaling factor being a free parameter (Table 1).
2.2 Geomagnetic data
As the geomagnetic data for the last millennia, we used the recent archeo/paleomagnetic model GMAG.9k of the virtual axial dipole moment (VADM) as published by Usoskin et al. (2016) for the period since ca. 7000 BC. This model provides an ensemble of 1000 individual VADM reconstructions that include all the uncertainties. The range of the ensemble reconstruction covers other archeo and paleomagnetic models (e.g., Genevey et al. 2008; Knudsen et al. 2008; Licht et al. 2013; Nilsson et al. 2014; PavónCarrasco et al. 2014), as is shown in Fig. 2 of Usoskin et al. (2016), which means that the related uncertainties are covered. The length of the geomagnetic series limits our study to the period since 6760 BC.
3 Temporal synchronization of the records: wiggle matching
While the radiocarbon is absolutely dated via dendrochronology by treering counting, dating of ice cores is less precise. In principle, dating of ice cores can be done, especially on short timescales, by counting annual layers of sulfate or sodium when the accumulation and the characteristics of the site allow it (e.g., Sigl et al. 2015). In practice, however, dating of long series is performed by applying iceflow models between tie points that are related to known events, such as volcano eruptions that leave clear markers in ice. While the accuracy of dating is quite good around the tie points, the exact ages of the samples between the tie points may be rather uncertain. This ranges from several years during the last millennium to up to 70–100 yr in the earlier part of the Holocene (Muscheler et al. 2014; Sigl et al. 2015; Adolphi & Muscheler 2016).
Since our method is based on the minimization of the χ^{2} statistics of all series for a given moment in time, it is crucial that these series are well synchronized. Accordingly, we performed a formal synchronization of the series based on the wigglematching procedure, which is a standard method for synchronizing time series (e.g., Cain & Suess 1976; Hoek & Bohncke 2001; Muscheler et al. 2014). We took the ^{14}C chronology as the reference and adjusted the timing of all other series to it.
3.1 Choice of wiggles
Since cosmogenic isotope series exhibit strong fluctuations of various possible origins, it is important to match only those wiggles that can presumably be assigned to the production changes, that is, cosmic ray (solar) variability, and exclude possible regional climate spells. One wellsuited type of wiggles is related to the socalled grand minima of solar activity, which are characterized by a fast and pronounced drop in solar activity for several decades. The most famous example of a grand minimum is the Maunder minimum, which occurred between 1645–1715 (Eddy 1976; Usoskin et al. 2015). Grand minima can be clearly identified in the cosmogenic isotope records as sharp spikes (Usoskin et al. 2007; Inceoglu et al. 2015). A typical (although not the most pronounced) example of such wiggles (grand minima) is shown in Fig. 2 for the (scaled) original EDML ^{10} Be series in a dashed curve along with the ^{14}C production data in red. Although the overall variability of the two series looks similar, the mismatch in the timing between them is clear and is roughly a few decades.
For the further analysis, we selected periods with clear wiggles (spikes, corresponding to grand minima, as listed in earlier works – see references in the notes to Table 2) in the ^{14} C series as listed in Table 2, along with their centers and time span.
3.2 Synchronization of the wiggles
For all the selected wiggles (Table 2), we found the bestfit time adjustment d T between the analyzed ^{10}Be and the reference ^{14}C production series by maximizing the cross correlation between the series, calculated within a time window centered at the middle of the wiggle. The data were annually interpolated within the time windows so that the time step in defining d T was one year. The length of the correlation window was chosen as twice the length of the wiggle (see Table 2). For each wiggle, we repeatedly calculated the Pearson linear correlation coefficients between the ^{14} C and ^{10} Be series, and we selected the value of dT that maximized the crosscorrelation coefficient R between the two series. The standard error (s_{err}) of the correlation coefficient was calculated using the approximate formula (e.g., Cohen et al. 2003): (1)
where R_{c} is the maximum correlation coefficient and n is the number of the data points within the correlation window. This uncertainty s_{err} was translated into the 1σ confidence interval for dT, as illustrated in Fig. 3, which shows the correlation coefficient, R (black curve), as a function of the time shift dT. It reaches its maximum R_{c} = 0.86 at d T = −25 yr, as indicated by the vertical solid line. The dotted lines bound the 68% confidence interval for d T defined as R = R_{c} − s_{err}.
The “momentary” time adjustments dT were considered as “tie points” (listed in Table 2) for the ^{10} Be series, with a linear interpolation used between them. Adjustments for the EDML series, based here on a new Antarctic Ice Core Chronology (AICC2012; Veres et al. 2013), lie within +20/−40 yr, while for the GRIP series, they vary within +20/−50 yr. The adjustment range for GRIP is concordant with the Greenland chronology correction function (e.g., Muscheler et al. 2014) within the uncertainties of the GICC05 timescale (Seierstad et al. 2014). Some discrepancies may be caused by the differences in the datasets and applied method. There are no earlier results for the EDML synchronization chronology to be compared with. We emphasize here that we do not pretend to perform a full chronological scale update, but only to match wiggles between a single beryllium series and ^{14} C data, which is sufficient forthis work. In particular, earlier synchronization studies produced smooth correction curves (e.g., Knudsen et al. 2009; Muscheler et al. 2014), where individual wiggles may still be slightly mismatched, while we are focused here on matching each wiggle in each series separately.
An example of the resulting redated ^{10}Be series compared to the reference ^{14}C record is shown with the solid black curve in Fig. 2. The synchronization obviously improves the cross correlation between the series, as shown in Table 3. The improvement (in terms of the ratio of R^{2}, which is a measure of the power of covariability between the original and synchronized series) is significant for the long (1.1 and 2.14 for the GRIP and EDML series, respectively) and shorter Greenland series (NGRIP and Dye3), but small for the short Antarctic series.
The pairwise wavelet coherence between longterm series is shown in Fig. 4, calculated following the procedure described in Usoskin et al. (2009), including the significance estimate using the nonparametric randomphase method (Ebisuzaki 1997). The coherence between EDML ^{10} Be and ^{14} C series (panel a) is good on timescales shorter than 1000 yr and insignificant on timescales longer than 2000–3000 yr, with no coherence in between. The coherence between the GRIP ^{10}Be and ^{14}C series is good at all timescales longer than 400–500 yr. The coherence between the GRIP and EDML series is intermittent on timescales shorter than 1000 yr and insignificant on longer timescales.
Fig. 2
Typical example of a wiggle in the ^{14}C (red) series compared with the original EDML ^{10}Be series (dashed black, scaled up by a factor of 172). Solid and dashed vertical lines denote the middle and the span of the wiggle considered for ^{14}C. The black sold curve shows the EDML ^{10}Be series after the synchronization (see Sect. 3.2). 
Wiggles (central date and the length in years) used for synchronization of the various ^{10} Be time series to ^{14}C (Eddy 1976; Stuiver & Quay 1980; Stuiver & Braziunas 1989; Goslar 2003; Inceoglu et al. 2015; Usoskin et al. 2016), and the synchronization time (in years) for the GRIP and EDML series.
Fig. 3
Example of the calculation of the bestfit time adjustment dT = −25 yr (solid vertical line) and its 68% confidence interval (dashed lines, −31 and −19 yr) for the wiggle case shown in Fig. 2. The Pearson linear correlation was calculated between the ^{14}C and EDML ^{10}Be series for the ± 100yr time window around the center of the wiggle. 
Squared correlation coefficients between the six ^{10}Be series and the ^{14}C series for the originally dated (R_{o}) and synchronized (R_{s}) series.
4 Reconstruction of the solar modulation potential
Since cosmogenic isotopes are produced by cosmic rays in the Earth’s atmosphere (Beer et al. 2012), their measured production/depositional flux reflects changes in the cosmic ray flux in the past. In turn, cosmic rays are modulated by solar magnetic activity, which is often quantified in terms of the modulation potential ϕ. The latter is a useful parameter to describe the solar modulation of GCRs using the socalled forcefield parametrization formalism (e.g., CaballeroLopez & Moraal 2004; Usoskin et al. 2005). In an ideal case, when both the production rate of cosmogenic isotopes and the geomagnetic field at a given time are known, the corresponding modulation parameter can be calculated for the given isotope using a production model that considers in great detail all the processes of the nucleonicmuonelectromagnetic cascade that are triggered by energetic cosmic rays in the atmosphere. Here we used the production model by Poluianov et al. (2016), which is a recent update of the widely used cosmic ray atmosphericcascade model (CRAC; Kovaltsov & Usoskin 2010; Kovaltsov et al. 2012). This model provides absolute production rates and is in full agreement with other modern models (Pavlov et al. 2017). The modulation potential ϕ is defined here (see the full formalism in Usoskin et al. 2005) for the local interstellar spectrum according to Burger et al. (2000). Since the modulation potential is a modeldependent parameter, our result cannot be directly compared to the ϕ−values based on different assumptions (e.g., Steinhilber et al. 2012) without a recalibration (Herbst et al. 2010; Asvestari et al. 2017). Thus, the ϕ−series is an intermediate result that is further converted into an physical index of the open magnetic flux and subsequently into the SNs.
The relation between the isotope production rate and its measured content (Δ ^{14}C for ^{14} C and depositional flux or concentration for ^{10}Be) depends on the corresponding atmospheric or terrestrial cycle of the isotope. Radiocarbon is involved, as carbondioxide gas, in the global carbon cycle and is almost completely mixed and homogenized over the global hemisphere. Here we used the globally averaged ^{14} C production rate as computed by Roth & Joos (2013) from the INTCAL09 standard Δ ^{14} C dataset (Reimer et al. 2009) using a newgeneration dynamic carbon cycle model that includes coupling with the diffusive ocean. However, the effect of extensive fossil fuel burning (Suess effect) makes it difficult to use ^{14} C data after the midninteenth century because of large and poorly constrained uncertainties (Roth & Joos 2013). Radiocarbon data cannot be used after the 1950s because of manmade nuclear explosions that led to massive production of ^{14} C. Accordingly, we did not extend our analysis to the twentieth century (cf., e.g., Knudsen et al. 2009).
In contrast to ^{14}C, ^{10} Be is not globally mixed, and its transport or deposition in the atmosphere is quite complicated and subject to local and regional conditions. Here we applied the ^{10}Be production model by Poluianov et al. (2016), and atmospheric transport and deposition were considered via the parametrization by Heikkilä et al. (2009, 2013), who performed a full 3D simulation of the beryllium transport and deposition in the Earth’s atmosphere. However, the existing models consider only the largescale atmospheric transport and do not address in full detail the deposition at each specific location; this may differ significantly from site to site. This remains an unknown factor (up to 1.5 in either direction) between the modeled and actually measured deposition flux of ^{10} Be at any given location (e.g., Sukhodolov et al. 2017). On the other hand, a free conversion factor exists if the ^{10} Be data are provided in concentration units rather than depositional flux. Therefore, we considered a single scaling factor that enters the production rate Q so that it matches the measured values. This factor was adjusted for each ^{10} Be series separately.
Sporadic solar energetic particle (SEP) events are sometimes produced by the Sun, with strong fluxes of energetic particles impinging on the Earth’s atmosphere. Although these solar particle storms usually have a short duration and soft energy spectrum, the SEP flux can produce additional cosmogenic isotopes in the atmosphere. For extreme events, the enhancement of the isotope production may greatly exceed the annual yield from GCR (Usoskin & Kovaltsov 2012). If not properly accounted for, these events may mimic periods of reduced solar activity (McCracken & Beer 2015), since the enhanced isotope production is erroneously interpreted in terms of the enhanced GCR flux, and consequently, in terms of reduced solar activity. Two extreme SEP events are known that can lead to such an erroneous interpretation (Bazilevskaya et al. 2014): the strongest event occurred around 775 AD (Miyake et al. 2012), and a weaker event was reported in 994 AD (Miyake et al. 2013). The energy spectra of these events have been assessed elsewhere (Usoskin et al. 2013; Mekhaldi et al. 2015). The production effect of the event on the ^{10} Be data in polar ice was calculated by Sukhodolov et al. (2017) and removed from the original data. One potential candidate around 5480 BC studied by Miyake et al. (2017) appears to be an unusual solar minimum rather than an SEP event. Accordingly, we kept it as a wiggle and did not correct for the possible SEP effect.
Here we first calculated the modulation potential ϕ in the past for each series individually, and then for all series together. The reconstruction process is described in detail below.
Fig. 4
Wavelet coherence between long series, redated using the wiggle matching, considered here: panel a: ^{10} Be EDML vs. ^{14}C, panel b: ^{10}Be GRIP vs. ^{14}C, and panel c: ^{10}Be GRIP vs. ^{10}Be EDML. The color scale ranges from 0 (deep blue) to 1 (dark red). Arrows denote the relative phasing between the series: arrows pointing right denote phase matching, while arrows pointing left show the antiphase. The white curves denote the cone of influence (COI) beyond which the result is unreliable. 
Fig. 5
Panel a: original (black) mean ^{14}C production rate (Roth & Joos 2013) and the one reduced to the standard geomagnetic conditions (red). Panel b: mean VADM reconstruction (Usoskin et al. 2016). 
4.1 Reducing the series to the reference geomagnetic conditions
The temporal variability of cosmogenic isotope production contains two signals: solar modulation, and changes in geomagneticfield. These signals are independent of each other and can thus be separated. Since here we are interested in the solar variability, we removed the geomagnetic signal by reducing all production rates to the reference geomagnetic conditions, defined as follows: the geomagnetic field is dipolelike, aligned with the geographical axis, and its virtual axial dipole moment (VADM) is 8 × 10^{22} A m^{2}. The exact VADM value is not important for the procedure, therefore we chose a rounded value close to the mean value for the twentieth century. The reduction was made in two steps:
 1.
From the given isotope production rate Q(t) and the geomagnetic VADM M(t), we calculated the value of ϕ(t) using the production model by Poluianov et al. (2016).
 2.
From the value of ϕ(t) computed in step 1, we calculated the reduced production rate Q^{*}(t) using thesame production model, but now fixing the VADM at M = 8 × 10^{22} A m^{2}. This value roughly corresponds to conditions in the twentieth century. This yields the isotope production rate as it would have been during time t if the geomagnetic field had been kept constant at this value of M.
The new Q^{*} series is now free (in the framework of the adopted model) of the geomagnetic changes and is used in the subsequent reconstructions. An example of the original series and the series corrected for variations in the geomagnetic field is shown in Fig. 5.
Fig. 6
Example of the χ^{2} vs. ϕ dependence for 805 AD for the ^{14}C series. The dashed lines represent the 68% confidence interval for ϕ. 
4.2 Reconstruction based solely on ^{14}C
First, we calculated the solar modulation potential based solely on the radiocarbon data using a Monte Carlo method similar to the method developed by Usoskin et al. (2014, 2016). The reconstruction includes the following steps:
 1.
For each moment t in time, we used 1000 realizations Q_{i}(t) of the full ensemble provided by Roth & Joos (2013). These realizations include uncertainties of the carbon cycle and of the measurement errors. At the same time, we also used 1000 realizations of the VADM M_{j}(t) from Usoskin et al. (2016), which include uncertainties and cover the range of available archeomagnetic reconstructions to calculate 10^{6} values of , as described in Sect. 4.1. The whole ensemble provides a natural way to represent the range and uncertainties of the reconstructed quantities (Q or VADM). For this ensemble, we calculated the mean ⟨Q^{*}(t)⟩ and the standard deviation σ_{Q}(t).
 2.
Using the statistics of the Q^{*}(t) ensemble, we defined for each time t, the bestfit ϕ(t) that minimizes the value of χ^{2}:
where Q′(ϕ) is the value of Q^{*} computed for a given value of ϕ, which was scanned over the range 0–2000 MeV. The bestfit value of ϕ_{0} is defined as the value corresponding to the minimum (→0 for a single series used). The 68% confidence interval of ϕ is defined as the interval bounded by the values of . An example of the χ^{2}(ϕ) dependence and definition of the bestfit ϕvalues and its uncertainties is shown in Fig. 6.
 3.
The series of reconstructed ϕ based on ^{14}C (ϕ_{14C}) was then computed, along with the uncertainties, as shown in Fig. 7.
4.3 Comparison between the long ^{10}Be and ^{14}C series
Next we compared the longterm behavior, in the sense of the scaling factors, of the longrunning ^{10} Be series vs. the reference ^{14}C series. Firstwe considered a 1000yr window and calculated the mean value of ⟨ϕ_{14C}⟩ within this period, as described in Sect. 4.2. We then scaled each ^{10} Be series individually with a scaling factor κ and reconstructed the value of ϕ_{κ} for the rescaled ^{10} Be series in this time window. The value ofκ was defined such that the mean value of ϕ_{10Be} agrees with the value for ^{14}C for the same period, that is, ⟨ϕ_{10Be}⟩ = ⟨ϕ_{14C}⟩. Then, the 1000 yr window was moved by 100 yr and the procedure was repeated. The resulting ^{10} Be scaling factors κ are shown as a function of time in Fig. 8.
The κfactor for the GRIP series is relatively stable during the period 6760 BC–3000 BC and depicts a steady monotonous decrease around 3000 BC and reaching −12% with respect to the final κ_{GRIP} around 1000 AD. The Pearson squared correlation between the two curves is significant, R^{2} = 0.78 (p−value 0.02)^{2}, implying a possible residual effect of the geomagnetic field in the reconstruction.
It is interesting to note that the two ^{10}Be series show very different trends during the second half of the Holocene. The κfactor for the EDML series shows a weak growing trend against the ^{14}C series over the entire period of their overlap, with κ varying from −5% to +10%. A shallow wavy variability can be noticed, with a quasiperiod of approximately 2400 yr, which is probably related to the Hallstatt cycle (Damon & Sonett 1991; Usoskin et al. 2016). The correlation between the κ−factor for EDML and the VADM is insignificant R^{2} = 0.52 (p = 0.12). This suggests that the undercorrected geomagnetic field effect is most likely only related to the GRIP series.
The discrepancy between GRIP and ^{14}C series is wellknown (e.g., Vonmoos et al. 2006; Inceoglu et al. 2015), but has typically been ascribed to the early part of the Holocene because both series are normalized to the modern period. With this normalization, the records agree with each other over the last millennia but diverge before ca. 2000 BC (Inceoglu et al. 2015). This discrepancies has sometimes (e.g., Usoskin 2017, and references therein) also been explained as a possible delayed effect (not perfectly stable thermohaline circulation) of the deglaciation in the carbon cycle (e.g., Muscheler et al. 2004). However, as we show here, this explanation is unlikely for two reasons.
More details of the longterm relation between the series are shown in Fig. 9, where we plot the pairwise wavelet coherence between the modulation potential (shown in Fig. 12) from the three long series. Panel a shows that the EDML series is coherent with the ^{14} C series at all timescales shorter than ≈2000 yr and again (but insignificantly) longer than 4000 yr, while no coherence exists between 2000 and 4000 yr. The GRIPbased series (panel b) is well coherent with the ^{14}C one on timescales longer than 100–200 yr, and the coherence disappears at timescales longer than 2000–3000 yr. The two ^{10} Be are hardly coherent with each other on timescales longer than 1000 yr (panel c), with a small insignificant isle of coherence beyond the cone of influence. It is noteworthy that each of the ^{10} Be series exhibits a higher coherence with the ^{14}C one than with each other. A similar conclusion was made by Usoskin et al. (2009) using different ^{10} Be data series. A reason for this discrepancy is yet to be discovered.
Since the reason for this discrepancy is unknown, we did not correct for it. In the following we apply fixed coefficients for ^{10} Be series, thus keeping the longterm variability as it exists in the original data.
It is interesting that while the coherence between EDML and ^{14} Cbased reconstructions has been improved compared to the coherence between the raw data series (Fig. 4), the coherence is degraded for the GRIP series (panels b) in timescales longer than 2000 yr. This suggests that beryllium may be better mixed for Greenland than was derived through the applied beryllium transport parametrization based on Heikkilä et al. (2009). The reconstruction based on the ^{14} C series has no residual correlation with VADM (R^{2} < 0.01). A new fullsize model of the atmospheric transport and deposition of ^{10} Be in polar regions, particularly in Greenland, is needed to resolve the issue.
Fig. 7
Series of the modulation potential ϕ computed based only on the ^{14}C data. Shading denotes the 68% confidence interval. 
Fig. 8
The ^{10}Be scaling factor κ, with 68%uncertainties indicated by the grey shading, in the sliding 1000yr window as a function of time. Panels a and b are for the GRIP and EDML series, respectively. Dashed lines depict the κ factor defined for the entire period (Sect. 4.5). 
4.4 Reconstructions from individual series
Two beryllium series, GRIP and EDML, were used here for a preliminary assessment of this method for their overlap. The mean modulation potential for the period 6800 BC–700 AD based on radiocarbon is ⟨ϕ_{14C}⟩ = 0.468 GV.
First, we scaled each ^{10}Be series individually with a free scaling factor κ_{0} and reconstructed the series of ϕ. The values of κ_{0} were defined such that the mean value ⟨ϕ_{10Be}⟩ is equal to that for ^{14}C for the same period. The scaling factor we found for the GRIP data is close to unity, κ_{0} = 1.1. This implies that ^{14}C and ^{10}Be data are fully consistent in the framework of the model, and that our modeling of the cosmogenic isotope production and transport or deposition is appropriate. The bestfit scaling factor for EDML is κ_{0} = 0.8, which is about 20% lower than unity, which is reasonable considering specific features of the deposition at this site. These values of κ_{0} were considered as initial guesses for a more precise search for the scaling factors, as described below.
Fig. 9
Wavelet coherence between individual longterm series of the modulation potential ϕ reconstructed from individual longterm cosmogenic series (shown in Fig. 12). Panel a: ^{10} Be EDML vs. ^{14}C; panel b: ^{10}Be GRIP vs. ^{14}C; panel c: ^{10}Be GRIP vs. ^{10}Be EDML. Notations are similar to Fig. 4. 
4.5 Combinedrecord reconstruction of ϕ
In the preceding section we calculated three series of ϕ from different cosmogenic isotope datasets, using their intercalibration to the mean ϕ value obtained from the ^{14}C data over the overlap period. This is similar to what has been done previously (Vonmoos et al. 2006; Steinhilber et al. 2012; Usoskin et al. 2014, 2016). Here we proceed and perform a consistent Bayesianbased reconstruction of the modulation potential ϕ. Using the measured data and knowledge of the other complementary parameters, we determine for each moment in time the most probable value of ϕ and its uncertainty.
 1.
First, we fixed the scaling factors for the GRIP and EDML series at their initial guess values κ_{0} as described above.
 2.
We then calculated, as described in Sect. 4.2 (step 1), 10^{6} realizations of the isotope production rates Q^{*}(t) reduced to the standard geomagnetic conditions.
 3.
The mean ⟨Q^{*}(t)⟩ and the standard deviation were calculated over the 10^{6} ensemble members for each time point t.
 4.
For each t, the value of χ^{2}(ϕ) was calculated as
where the index i takes values 1 through 3 for the ^{14}C, GRIP, and EDML series, respectively. An example of χ^{2} as a function of ϕ at one pointin time for the three individual series as well as for their sum is shown in Fig. 10. Each of the individual datasets (panels a–c, each similar to Fig. 6) yields a very sharp and welldefined dip in the χ^{2} value; this dip approaches zero. The reason is that for a given single value of Q, the corresponding value of ϕ can be defined precisely. However, the obtainedindividual ϕ−values are not identical for different datasets, which leads to a smooth overall χ^{2} vsϕ dependence (panel d), the minimum χ^{2} of which is about 1.7 or 0.85 per degree of freedom (DoF). This implies that the same value of ϕ =0.58 GV satisfies all three isotope data records within statistical confidence.
 5.
We calculated the sum of individual χ^{2} values (Eq. (3)) as over all 750 time points during 6760 BC–730 AD. The corresponding sum is or ≈ 2.9 per DoF, indicating a likely systematic difference between the series. The DoF number is defined as 750 × 3 (numberof points in the three series) minus 752 (the number of fitted parameters), thus 1498.
 6.
Since thevalues of the scaling factors κ_{0} were initially defined by normalizing the mean values of ϕ, which is not optimal, we redefined them via χ^{2} as follows. We repeated steps (1)–(5) above, now scanning the values of κ in the range of 0.85–1.25 and 0.6–1.1, with steps 0.01 and 0.015, for the GRIP and EDML series, respectively. The corresponding were calculated. The distribution of as a function of κ_{GRIP} and κ_{EDML} is shown in Fig. 11. The distribution has a clear minimum (), and the values of κ_{GRIP} and κ_{EDML} are 1.028 and 0.815, respectively, which is close to the initial guess (Sect. 4.3).
Reconstructions of ϕ based on the bestfit scaling for individual series are shown in Fig. 12 as colored curves. The differences between the data points from the different series are normally distributed, with a mean value of about zero and the standard deviation of 100–200 MV. This serves as an estimate of the accuracy of individual reconstructions.
Fig. 11
Distribution of vs. κ_{GRIP} and κ_{EDML}. The gray scale is on the right. The minimum of the distribution () occurs at κ_{GRIP} = 1.028 and κ_{EDML} = 0.815, as marked by the crosshairs. 
Fig. 10
Dependence of χ^{2} (Eq. (3)) on ϕ for the decade centered at 2095 BC: panel a: ^{14}C with no scaling, panel b: ^{10}Be GRIP with scaling (κ = 1.028), panel c: EDML with scaling κ = 0.815, and panel d: sum of the three χ^{2} components. 
4.6 Full reconstruction
Next, we analyzed all the available data series, also including four shorter ^{10} Be series. This extends the ϕ reconstruction to 1900 AD. The four shorter series are (see Table 1) NGRIP and Dye3 from Greenland, and Dome Fuji (DF) and South Pole (SP) from Antarctica.
We first estimated the bestguess scaling κ−factors for each short series similar to what is described in Sect. 4.3, that is, we equalized the mean values of ⟨ϕ⟩ for the series in question to the value of the reference ^{14}C series for the period of their overlap. We then slightly varied the κ−value around these bestguessed values to calculate the corresponding χ^{2} for each individual series, as shown in Fig. 13. The vertical solid lines mark the bestfit κvalues that minimize the fivepoint smoothed χ^{2}. The bestfitscaling factors are listed in Table 1. The scaling factor for ^{10} Be series with depositional flux data (GRIP, EDML, NGRIP, and Dome Fuji, see Table 1 and Fig. 8) are close to unity (within 20%), which again implies that our model is quite realistic. The 68% uncertainties are defined as corresponding to .
Next, we performed the full reconstruction using the χ^{2} method described in Sect. 4.5, but now for all the series (see Fig. 12). The number of different series used to reconstruct individual data points varied in time between two and five (see Fig. 1). The mean χ^{2} per DoF is 0.57, for 84% of data points χ^{2} < 1 per DoF, implying that the agreement between different series is good.
Fig. 12
Reconstruction of the modulation potential ϕ using only individual cosmogenic isotope series (color curves as denoted in the legend) with the bestfit scaling (see Table 1) and the final composite series (thick black curve). The top and bottom panels depict the two halves of the entire interval. Only mean values are shown without uncertainties. 
5 Reconstruction of the sunspot number
The modulation potential series alone are not very useful as a solar activity proxy since the modulation potential is a relative index whoseabsolute value is model dependent (Usoskin et al. 2005; Herbst et al. 2010, 2017). Therefore, we converted the modulation potential, reconstructed in Sect. 4, into a more definitive index, the SN. This was done via the open solar magnetic flux F_{o}, following an established procedure (e.g., Usoskin et al. 2003, 2016; Solanki et al. 2004). Applying the updated SATIREM model (Vieira & Solanki 2010; Vieira et al. 2011; Wu et al. 2018), we can write the relation (Usoskin et al. 2007, see the appendix therein) between the two indices as (4)
where the SN during the ith decade, SN_{i}, is defined by the modulation potential (expressed in GV) during the contemporary and the following decades. The negative offset term (−16) reflects the fact that zero SN during grand minima does not imply the absence of GCR modulation, so that even when SN = 0, the value of ϕ is about 0.14 GV (e.g., Owens et al. 2012). Since SNs cannot be negative, the SN was assigned zero values when the values of ϕ dropped below the sunspot formation threshold. The uncertainties of the reconstructed SN values were defined by converting the low and upperbound (68%) ϕ values (defined as described above) for each decade into SNs. The reconstructed SN series^{3} is shown in Fig. 14 along with its 68% confidence interval and is compared with the international SN series (SILSO ISN version 2, Clette et al. 2014)^{4}. Since we used SNs in their “classical” definition, the ISN v.2 data were scaled down by a factor 0.6 (as described in Clette et al. 2014).
The reconstructed solar activity varies at different timescales from decades to millennia. In particular, a long period with relatively high activity occurred between roughly 4000 BC and 1500 BC, and periods of lower activity occurred at ca. 5500 BC and 1500 AD. The origin of this is unclear. For example, Usoskin et al. (2016) suggested, based on the fact that the GRIP and ^{14} C series behave differently at this timescale, that it is a global climate effect rather than solar. As a result, they removed this wave from the data in an ad hoc manner. However, as we show here, it is more likely related to the undercorrection of the GRIPbased series and therefore may be related to solar activity, so that we retained it in the final dataset.
The new series is generally consistent with previous reconstructions (e.g., Usoskin et al. 2016), but it also has some new features. In particular, it implies a lower activity during the sixth millennium BC. We note that a possible overestimation of the solar activity for this period has been suspected previously (Usoskin et al. 2007, 2016).
Although the overall level of solar activity may vary significantly, the periods of grand minima may correspond to a special state of the solar dynamo (Schmitt et al. 1996; Küker et al. 1999; Moss et al. 2008; Choudhuri & Karak 2012; Käpylä et al. 2016)and thus are expected to provide roughly the same low level of activity, corresponding to a virtual absence of sunspots. Thus, the level of the reconstructed activity during clearly distinguishable grand minima may serve as a rough estimate of the “stability” of the reconstruction when we assume that activity always drops to the same nearly zero level during each grand minimum (Sokoloff & NesmeRibes 1994; Usoskin et al. 2014). The expected level of solar activity during grand minima (SN → 0) is indicated by the horizontal dashed line in Fig. 14. The level of the reconstructed grand minima (observed as sharp dips) is roughly consistent with SN → 0 throughout almost the entire period, considering an uncertainty on the order of 10 in SN units. However, periods of 3500 BC–2500 BC and before 5500 BC are characterized by a slightly higher SN level during the grand minima, suggesting a possible overestimation of activity during these times. Interestingly, no clear grand minima occurred during 2500 BC–1000 BC, suggesting that it was a long period of stable operation of the solar activity main mode.
The overall level of solar activity remained roughly constant, around 45–50, during most of the time, but appeared somewhat lower (around 40) before 5000 BC, suggesting a possible slow variability with a timescale of about 6–7 millennia (cf., Usoskin et al. 2016). The reason for this is unknown, but it might be due to (1) climate influence, although this is expected to affect ^{14} C and ^{10} Be isotopes differently; (2) solar activity, or (3) large systematic uncertainty in the geomagnetic field reconstruction, which is poorly known before about 3000 BC.
Figure 15 shows the kernel density estimation of the probability density function (KDF) of the reconstructed decadal SNsfor the entire period (877 decades). The high peak of the distribution with values around 40 is clearly visible; this corresponds to moderate activity. A lowactivity component with values below 15 is clearly visible as well; this corresponds to a grand minimum. In addition, a bump is faintly visible at the highactivity tail (SN greater than 60, visible only as a small excess above the main Gaussian curve), suggesting a grand maximum component. In order to illustrate this, we applied a formal multipeak (Gaussian) fit to the KDF, which is shown as the blue dotted curves. The main peak is centered at decadal SN = 42 (σ = 10) and represents the main component. Another peak can be found as a small bump around SN = 27 (σ = 6), but it is not sufficiently separated from the main component to justify calling it a separate component. The grand minimum component (Gaussian with σ = 6 centered at SN = 12) is significantly separated from the main component (statistical significance p < 0.05). The separation of the grand maximum component, while visible by eye as a deviation from the Gaussian shape at high values, is not statistically significant. The statistical separation of the special grand minimum component was first shown by Usoskin et al. (2014) for the last three millennia, while the result for the grand maximum component was inconclusive.
Here we fully confirm this result over nine millennia, which implies that grand minimum and normal activity components form a robust feature of solar variability. This also suggests that the new reconstruction is more robust and less noisy than previous reconstructions (e.g., Steinhilber et al. 2012; Usoskin et al. 2016) and allows statistically identifying the grand minimum component over 9000 yr.
Fig. 13
Values of χ^{2} vs. the scaling factors κ for four short ^{10}Be series, as marked in each panel. The fivepoint running mean curves are shown in red. The locations of the bestfit κ−values are shown as straight black lines. 
Fig. 14
Reconstructed SN along with its 68% confidence interval (gray shading). This series is available in the ancillary data. The red line depicts the decadally resampled international SN (version 2, scaled by 0.6) from Clette et al. (2014). The dashed line denotes the level of SN = 10. 
Fig. 15
Probability kernel density function estimate (black dots) of the reconstructed decadal SNs (Gaussian kernel, width = 3) and its multipeak fit (red). 
6 Conclusions
We have provided a new fully consistent multiproxy reconstruction of the solar activity over about nine millennia, based for the first time on a Bayesian approach. We used all the available datasets of cosmogenic radioisotopes with sufficient length and quality in terrestrial archives anduptodate models of isotope production and transport or deposition as well as a recent archeomagnetic model. We used six ^{10}Be series of different lengths from Greenland and Antarctica, and the official INTCAL global ^{14} C series. Earlier reconstructions were based on either individual datasets or on a statistical superposition of them (e.g., Steinhilber et al. 2012). Our new method is based on finding the most probable value of the solar modulation potential that matches all the data for a given point in time, providing also a straightforward estimate of the uncertainties. Prior to the analysis, long ^{10} Be series were formally redated to match wiggles in the ^{14}C data. All employed cosmogenic isotope series were reduced to the reference geomagnetic field conditions.
We have also tested the stability of the two long ^{10}Be series and found that they appear to diverge from each other during the second half of the Holocene, while the ^{14} C series lies in between them. The GRIPbased series appears to anticorrelate with the longterm geomagnetic VADM series, implying a possible undercorrection for the geomagnetic shielding effect for the GRIP location, while EDML and ^{14} Cbased reconstructions do not show any significant residual correlation with VADM. This suggests that the applied model of beryllium transport and deposition does not work properly for the Greenland site but is reasonably good for the Antarctic site. A fullsize transport and deposition model (e.g., Sukhodolov et al. 2017) needs to be applied in the future to resolve this issue.
The reconstructed series (Fig. 14) shows variability on different timescales. Grand minima of activity are of particular interest. They are visible as strong dips in the time series, in which the level of the decadal SNs sinks below 10–15. Another feature of the longterm evolution of the solar activity is a slow variability on the 6–7 millennia timescale with lows occurring in ca. 5500 BC and 1500 AD. This behavior has been interpreted by Usoskin et al. (2016) as a possible effect of climate influence on the carbon cycle, but the ^{14} C series lies between the two diverging ^{10}Be series on this timescale, which makes this explanation unlikely. The cause of the feature remains unknown.
Most of the time, the solar activity varies slightly around the moderate level of SN ≈ 40, which corresponds to the main component of the solar activity. Grand minima form a statistically distinguishable component, however. The existence of this component was first determined for the last three millennia by Usoskin et al. (2014), who interpreted it as special mode of the solar dynamo. Its confirmation here for nine millennia implies that it is a robust feature of the solar activity. At the same time, the possible existence of a component representing grand maxima is indicated (cf. Usoskin et al. 2014), although it cannot be separated from the main component in a statistically significant manner. We speculate that the different components of the activity distribution may be related to different modes of the dynamo operation.
Finally, a new consistent reconstruction of the solar activity (in the form of decadal SNs) was presented that offers a more reliable estimate of the longterm evolution of the solar variability and poses robust constraints on the development of solar and stellar dynamo models as well as solarterrestrial studies.
Acknowledgements
We acknowledge Mads Knudsen for consultations on the GRIP series in GICC05 timescale and Raimund Muscheler for discussions about the ice core chronologies. Support by the Academy of Finland to the ReSoLVE Center of Excellence (project no. 272157) is acknowledged. C.J.W. acknowledges postgraduate fellowship of the International Max Planck Research School for Solar System Science. This work has been partly supported by the BK21 plus program through the National Research Foundation (NRF) funded by the Ministry of Education of Korea.
References
 Adolphi, F., & Muscheler, R. 2016, Clim. Past, 12, 15 [CrossRef] [Google Scholar]
 Asvestari, E., Gil, A., Kovaltsov, G. A., & Usoskin, I. G. 2017, J. Geophys. Res.: Space Phys., 122, 9790 [Google Scholar]
 Bard, E., Raisbeck, G., Yiou, F., & Jouzel, J. 1997, Earth Planet. Sci. Lett., 150, 453 [NASA ADS] [CrossRef] [Google Scholar]
 Bazilevskaya, G. A., Cliver, E. W., Kovaltsov, G. A., et al. 2014, Space Sci. Rev., 186, 409 [Google Scholar]
 Bazin, L., Landais, A., LemieuxDudon, B., et al. 2013, Clim. Past, 9, 1715 [NASA ADS] [CrossRef] [Google Scholar]
 Beer, J., Siegenthaler, U., Oeschger, H., Bonani, G., & Finkel, R. 1988, Nature, 331, 675 [NASA ADS] [CrossRef] [Google Scholar]
 Beer, J., Blinov, A., Bonani, G., Hofmann, H., & Finkel, R. 1990, Nature, 347, 164 [NASA ADS] [CrossRef] [Google Scholar]
 Beer, J., McCracken, K., & von Steiger, R. 2012, Cosmogenic Radionuclides: Theory and Applications in the Terrestrial and Space Environments, Physics of Earth and Space Environments (Berlin: Springer) [Google Scholar]
 Berggren, A.M., Beer, J., Possnert, G., et al. 2009, Geophys. Res. Lett., 36, L11801 [Google Scholar]
 Burger, R., Potgieter, M., & Heber, B. 2000, J. Geophys. Res., 105, 27447 [Google Scholar]
 CaballeroLopez, R., & Moraal, H. 2004, J. Geophys. Res., 109, A01101 [Google Scholar]
 Cain, W. F., & Suess, H. E. 1976, J. Geophys. Res., 81, 3688 [NASA ADS] [CrossRef] [Google Scholar]
 Choudhuri, A. R., & Karak, B. B. 2012, Phys. Rev. Lett., 109, 171103 [NASA ADS] [CrossRef] [Google Scholar]
 Clette, F., Svalgaard, L., Vaquero, J., & Cliver, E. 2014, Space Sci. Rev., 186, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Cohen, J., Cohen, P., West, S., & Aiken, L. 2003, Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences (New York: Routledge) [Google Scholar]
 Damon, P., & Sonett, C. 1991, in The Sun in Time, eds. C. Sonett, M. Giampapa, & M. Matthews (Tucson: University of Arizona Press), 360 [Google Scholar]
 Delaygue, G., & Bard, E. 2011, Clim. Dyn., 36, 2201 [Google Scholar]
 Ebisuzaki, W. 1997, J. Clim., 10, 2147 [NASA ADS] [CrossRef] [Google Scholar]
 Eddy, J. 1976, Science, 192, 1189 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Genevey, A., Gallet, Y., Constable, C. G., Korte, M., & Hulot, G. 2008, Geochem. Geophys. Geosyst., 9, Q04038 [NASA ADS] [CrossRef] [Google Scholar]
 Goslar, T. 2003, PAGES News, 11, 12 [CrossRef] [Google Scholar]
 Hathaway, D. H. 2015, Liv. Rev. Sol. Phys., 12, 4 [Google Scholar]
 Heikkilä, U., Beer, J., & Feichter, J. 2009, Atmos. Chem. Phys., 9, 515 [Google Scholar]
 Heikkilä, U., Beer, J., Abreu, J. A., & Steinhilber, F. 2013, Space Sci. Rev., 176, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Herbst, K., Kopp, A., Heber, B., et al. 2010, J. Geophys. Res., 115, D00I20 [Google Scholar]
 Herbst, K., Muscheler, R., & Heber, B. 2017, J. Geophys. Res.: Space Phys., 122, 23 [Google Scholar]
 Hoek, W. Z., & Bohncke, S. J. P. 2001, Quat. Sci. Rev., 20, 1251 [NASA ADS] [CrossRef] [Google Scholar]
 Horiuchi, K., Uchida, T., Sakamoto, Y., et al. 2008, Quat. Geochronol., 3, 253 [CrossRef] [Google Scholar]
 Inceoglu, F., Simoniello, R., Knudsen, V. F., et al. 2015, A&A, 577, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johnsen, S. J., DahlJensen, D., Dansgaard, W., & Gundestrup, N. 1995, Tellus Ser. B, 47, 624 [NASA ADS] [CrossRef] [Google Scholar]
 Johnsen, S. J., DahlJensen, D., Gundestrup, N., et al. 2001, J. Quat. Sci., 16, 299 [CrossRef] [Google Scholar]
 Käpylä, M. J., Käpylä, P. J., Olspert, N., et al. 2016, A&A, 589, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Knudsen, M. F., Riisager, P., Donadini, F., et al. 2008, Earth Planet. Sci. Lett., 272, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Knudsen, M. F., Riisager, P., Jacobsen, B. H., et al. 2009, Geophys. Res. Lett., 36, L16701 [NASA ADS] [CrossRef] [Google Scholar]
 Kovaltsov, G., & Usoskin, I. 2010, Earth Planet. Sci. Lett., 291, 182 [NASA ADS] [CrossRef] [Google Scholar]
 Kovaltsov, G., Mishev, A., & Usoskin, I. 2012, Earth Planet. Sci. Lett., 337, 114 [Google Scholar]
 Küker, M., Arlt, R., & Rüdiger, G. 1999, A&A, 343, 977 [NASA ADS] [Google Scholar]
 Licht, A., Hulot, G., Gallet, Y., & Thébault, E. 2013, Phys. Earth Planet. Inter., 224, 38 [NASA ADS] [CrossRef] [Google Scholar]
 McCracken, K. G., & Beer, J. 2015, Sol. Phys., 290, 3051 [Google Scholar]
 McCracken, K., McDonald, F., Beer, J., Raisbeck, G., & Yiou, F. 2004, J. Geophys. Res., 109, 12103 [NASA ADS] [CrossRef] [Google Scholar]
 Mekhaldi, F., Muscheler, R., Adolphi, F., et al. 2015, Nat. Commun., 6, 8611 [NASA ADS] [CrossRef] [Google Scholar]
 Miyake, F., Nagaya, K., Masuda, K., & Nakamura, T. 2012, Nature, 486, 240 [NASA ADS] [CrossRef] [Google Scholar]
 Miyake, F., Masuda, K., & Nakamura, T. 2013, Nat. Commun., 4, 1748 [NASA ADS] [CrossRef] [Google Scholar]
 Miyake, F., Jull, A. J. T., Panyushkina, I. P., et al. 2017, Proc. Natl. Acad. Sci., 114, 881 [Google Scholar]
 Moss, D., Sokoloff, D., Usoskin, I., & Tutubalin, V. 2008, Sol. Phys., 250, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Muscheler, R., Beer, J., Wagner, G., et al. 2004, Earth Planet. Sci. Lett., 219, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Muscheler, R., Joos, F., Beer, J., et al. 2007, Quat. Sci. Rev., 26, 82 [Google Scholar]
 Muscheler, R., Adolphi, F., & Knudsen, M. F. 2014, Quat. Sci. Rev., 106, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Nilsson, A., Holme, R., Korte, M., Suttie, N., & Hill, M. 2014, Geophys. J. Int., 198, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Owens, M. J., Usoskin, I., & Lockwood, M. 2012, Geophys. Res. Lett., 39, L19102 [Google Scholar]
 Pavlov, A. K., Blinov, A. V., Frolov, D. A., et al. 2017, J. Atmos. Sol.Terr. Phys., 164, 308 [NASA ADS] [CrossRef] [Google Scholar]
 PavónCarrasco, F. J., Osete, M. L., Torta, J. M., & De Santis, A. 2014, Earth Planet. Sci. Lett., 388, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Poluianov, S. V., Kovaltsov, G. A., Mishev, A. L., & Usoskin, I. G. 2016, J. Geophys. Res. (Atm.), 121, 8125 [NASA ADS] [Google Scholar]
 Raisbeck, G., Yiou, F., Jouzel, J., & Petit, J. 1990, Roy. Soc. London Philos. Trans. Ser. A, 330, 463 [NASA ADS] [CrossRef] [Google Scholar]
 Rasmussen, S. O., Bigler, M., Blockley, S. P., et al. 2014, Quat. Sci. Rev., 106, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Reimer, P. J., Baillie, M. G. L., Bard, E., et al. 2009, Radiocarbon, 51, 1111 [Google Scholar]
 Ribes, J., & NesmeRibes, E. 1993, A&A, 276, 549 [NASA ADS] [Google Scholar]
 Roth, R., & Joos, F. 2013, Clim. Past, 9, 1879 [NASA ADS] [CrossRef] [Google Scholar]
 Schmitt, D., Schüssler, M., & FerrizMas, A. 1996, A&A, 311, L1 [NASA ADS] [Google Scholar]
 Seierstad, I. K., Abbott, P. M., Bigler, M., et al. 2014, Quat. Sci. Rev., 106, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Sigl, M., Winstrup, M., McConnell, J. R., et al. 2015, Nature, 523, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Sokoloff, D., & NesmeRibes, E. 1994, A&A, 288, 293 [NASA ADS] [Google Scholar]
 Solanki, S. K., Usoskin, I. G., Kromer, B., Schüssler, M., & Beer, J. 2004, Nature, 431, 1084 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Steinhilber, F., Abreu, J., Beer, J., et al. 2012, Proc. Natl. Acad. Sci. USA, 109, 5967 [Google Scholar]
 Stuiver, M. 1961, J. Geophys. Res., 66, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Stuiver, M., & Braziunas, T. 1989, Nature, 338, 405 [NASA ADS] [CrossRef] [Google Scholar]
 Stuiver, M., & Quay, P. 1980, Science, 207, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Sukhodolov, T., Usoskin, I., Rozanov, E., et al. 2017, Sci. Rep., 7, 45257 [NASA ADS] [CrossRef] [Google Scholar]
 Usoskin, I. G. 2017, Liv. Rev. Sol. Phys., 14, 3 [Google Scholar]
 Usoskin, I. G., & Kovaltsov, G. A. 2012, ApJ, 757, 92 [Google Scholar]
 Usoskin, I. G., Solanki, S. K., Schüssler, M., Mursula, K., & Alanko, K. 2003, Phys. Rev. Lett., 91, 211101 [Google Scholar]
 Usoskin, I. G., AlankoHuotari, K., Kovaltsov, G. A., & Mursula, K. 2005, J. Geophys. Res., 110, A12108 [Google Scholar]
 Usoskin, I. G., Solanki, S. K., & Kovaltsov, G. A. 2007, A&A, 471, 301 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Usoskin, I. G., Horiuchi, K., Solanki, S., Kovaltsov, G. A., & Bard, E. 2009, J. Geophys. Res., 114, A03112 [Google Scholar]
 Usoskin, I. G., Kromer, B., Ludlow, F., et al. 2013, A&A, 552, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Usoskin, I. G., Hulot, G., Gallet, Y., et al. 2014, A&A,562, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Usoskin, I. G., Arlt, R., Asvestari, E., et al. 2015, A&A, 581, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Usoskin, I. G., Gallet, Y., Lopes, F., Kovaltsov, G. A., & Hulot, G. 2016, A&A, 587, A150 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vaquero, J. M., Kovaltsov, G. A., Usoskin, I. G., Carrasco, V. M. S., & Gallego, M. C. 2015, A&A, 577, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Veres, D., Bazin, L., Landais, A., et al. 2013, Clim. Past, 9, 1733 [CrossRef] [Google Scholar]
 Vieira, L. E. A., & Solanki, S. K. 2010, A&A, 509, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vieira, L. E. A., Solanki, S. K., Krivova, N. A., & Usoskin, I. 2011, A&A, 531, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vinther, B. M., Clausen, H. B., Johnsen, S. J., et al. 2006, J. Geophys. Res.: Atm., 111, D13102 [NASA ADS] [CrossRef] [Google Scholar]
 Vonmoos, M., Beer, J., & Muscheler, R. 2006, J. Geophys. Res., 111, A10105 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, C.J., Krivova, N. A., Solanki, S. K., & Usoskin, I. G. 2018, A&A, submitted [Google Scholar]
 Yiou, F., Raisbeck, G., Baumgartner, S., et al. 1997, J. Geophys. Res., 102, 26783 [NASA ADS] [CrossRef] [Google Scholar]
Although nakedeye observations of sunspots are available, they do not provide quantitative assessments (Usoskin 2017).
The significance is estimated using the nonparametric randomphase method by Ebisuzaki (1997).
Available as a table at the CDS and at the MPS sunclimate webpage http://www2.mps.mpg.de/projects/sunclimate/data/SN_composite.txt.
File SN_y_tot_V2.0.txt available at http://www.sidc.be/silso/infosnytot.
All Tables
Wiggles (central date and the length in years) used for synchronization of the various ^{10} Be time series to ^{14}C (Eddy 1976; Stuiver & Quay 1980; Stuiver & Braziunas 1989; Goslar 2003; Inceoglu et al. 2015; Usoskin et al. 2016), and the synchronization time (in years) for the GRIP and EDML series.
Squared correlation coefficients between the six ^{10}Be series and the ^{14}C series for the originally dated (R_{o}) and synchronized (R_{s}) series.
All Figures
Fig. 1
Data series (see Table 1) and temporal coverage. Blue shading denotes the cosmogenic isotope series, while green shading shows the geomagnetic series. 

In the text 
Fig. 2
Typical example of a wiggle in the ^{14}C (red) series compared with the original EDML ^{10}Be series (dashed black, scaled up by a factor of 172). Solid and dashed vertical lines denote the middle and the span of the wiggle considered for ^{14}C. The black sold curve shows the EDML ^{10}Be series after the synchronization (see Sect. 3.2). 

In the text 
Fig. 3
Example of the calculation of the bestfit time adjustment dT = −25 yr (solid vertical line) and its 68% confidence interval (dashed lines, −31 and −19 yr) for the wiggle case shown in Fig. 2. The Pearson linear correlation was calculated between the ^{14}C and EDML ^{10}Be series for the ± 100yr time window around the center of the wiggle. 

In the text 
Fig. 4
Wavelet coherence between long series, redated using the wiggle matching, considered here: panel a: ^{10} Be EDML vs. ^{14}C, panel b: ^{10}Be GRIP vs. ^{14}C, and panel c: ^{10}Be GRIP vs. ^{10}Be EDML. The color scale ranges from 0 (deep blue) to 1 (dark red). Arrows denote the relative phasing between the series: arrows pointing right denote phase matching, while arrows pointing left show the antiphase. The white curves denote the cone of influence (COI) beyond which the result is unreliable. 

In the text 
Fig. 5
Panel a: original (black) mean ^{14}C production rate (Roth & Joos 2013) and the one reduced to the standard geomagnetic conditions (red). Panel b: mean VADM reconstruction (Usoskin et al. 2016). 

In the text 
Fig. 6
Example of the χ^{2} vs. ϕ dependence for 805 AD for the ^{14}C series. The dashed lines represent the 68% confidence interval for ϕ. 

In the text 
Fig. 7
Series of the modulation potential ϕ computed based only on the ^{14}C data. Shading denotes the 68% confidence interval. 

In the text 
Fig. 8
The ^{10}Be scaling factor κ, with 68%uncertainties indicated by the grey shading, in the sliding 1000yr window as a function of time. Panels a and b are for the GRIP and EDML series, respectively. Dashed lines depict the κ factor defined for the entire period (Sect. 4.5). 

In the text 
Fig. 9
Wavelet coherence between individual longterm series of the modulation potential ϕ reconstructed from individual longterm cosmogenic series (shown in Fig. 12). Panel a: ^{10} Be EDML vs. ^{14}C; panel b: ^{10}Be GRIP vs. ^{14}C; panel c: ^{10}Be GRIP vs. ^{10}Be EDML. Notations are similar to Fig. 4. 

In the text 
Fig. 11
Distribution of vs. κ_{GRIP} and κ_{EDML}. The gray scale is on the right. The minimum of the distribution () occurs at κ_{GRIP} = 1.028 and κ_{EDML} = 0.815, as marked by the crosshairs. 

In the text 
Fig. 10
Dependence of χ^{2} (Eq. (3)) on ϕ for the decade centered at 2095 BC: panel a: ^{14}C with no scaling, panel b: ^{10}Be GRIP with scaling (κ = 1.028), panel c: EDML with scaling κ = 0.815, and panel d: sum of the three χ^{2} components. 

In the text 
Fig. 12
Reconstruction of the modulation potential ϕ using only individual cosmogenic isotope series (color curves as denoted in the legend) with the bestfit scaling (see Table 1) and the final composite series (thick black curve). The top and bottom panels depict the two halves of the entire interval. Only mean values are shown without uncertainties. 

In the text 
Fig. 13
Values of χ^{2} vs. the scaling factors κ for four short ^{10}Be series, as marked in each panel. The fivepoint running mean curves are shown in red. The locations of the bestfit κ−values are shown as straight black lines. 

In the text 
Fig. 14
Reconstructed SN along with its 68% confidence interval (gray shading). This series is available in the ancillary data. The red line depicts the decadally resampled international SN (version 2, scaled by 0.6) from Clette et al. (2014). The dashed line denotes the level of SN = 10. 

In the text 
Fig. 15
Probability kernel density function estimate (black dots) of the reconstructed decadal SNs (Gaussian kernel, width = 3) and its multipeak fit (red). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.