EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 615 (July 2018) A&A, 615 (2018) A93 Full HTML
Free Access
Issue
A&A
Volume 615, July 2018
Article Number A93
Number of page(s) 13
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201731892
Published online 19 July 2018

© 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 & Nesme-Ribes 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 (14C), beryllium-10 (10Be), or chlorine-36 (36Cl) 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 14C production. The isotope 10Be 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 large-scale 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 long-term (multi-millennial) 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 large-scale ocean circulation for 14C) 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 multi-proxy 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 Ice-core 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 multi-proxy 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 14C by applying the standard wiggle-matching method (Cain & Suess 1976; Muscheler et al. 2014) to the official ice-core 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.

thumbnail 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 10Be series from Greenland and Antarctica and one global 14C production series, as summarized in Table 1 and Fig. 1.

The 14C 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 Δ14C. While this series has a pseudo-annual 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.

10Be 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 10Be 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 10Be 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 Dye-3 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).

Table 1

Temporal coverage, cadence, type of data (production rate, PR; depositional flux, D; or concentration, C) and the resulting best-fit scaling factors, κ (Sect. 4), of the cosmogenic isotope series.

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ón-Carrasco 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 tree-ring 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 ice-flow 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 wiggle-matching 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 14C 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 well-suited type of wiggles is related to the so-called 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 14C 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 best-fit time adjustment d T between the analyzed 10Be and the reference 14C 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 cross-correlation coefficient R between the two series. The standard error (serr) of the correlation coefficient was calculated using the approximate formula (e.g., Cohen et al. 2003): (1)

where Rc is the maximum correlation coefficient and n is the number of the data points within the correlation window. This uncertainty serr 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 Rc = 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 = Rc − serr.

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 10Be series compared to the reference 14C 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 R2, 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 long-term series is shown in Fig. 4, calculated following the procedure described in Usoskin et al. (2009), including the significance estimate using the non-parametric random-phase 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 10Be and 14C 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.

thumbnail Fig. 2

Typical example of a wiggle in the 14C (red) series compared with the original EDML 10Be 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 14C. The black sold curve shows the EDML 10Be series after the synchronization (see Sect. 3.2).
Table 2

Wiggles (central date and the length in years) used for synchronization of the various 10 Be time series to 14C (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.

thumbnail Fig. 3

Example of the calculation of the best-fit 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 14C and EDML 10Be series for the ± 100-yr time window around the center of the wiggle.
Table 3

Squared correlation coefficients between the six 10Be series and the 14C series for the originally dated (Ro) and synchronized (Rs) 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 so-called force-field parametrization formalism (e.g., Caballero-Lopez & 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 nucleonic-muon-electromagnetic 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 model-dependent 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 (Δ 14C for 14 C and depositional flux or concentration for 10Be) depends on the corresponding atmospheric or terrestrial cycle of the isotope. Radiocarbon is involved, as carbon-dioxide 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 new-generation 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 mid-ninteenth century because of large and poorly constrained uncertainties (Roth & Joos 2013). Radiocarbon data cannot be used after the 1950s because of man-made 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 14C, 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 10Be 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 large-scale 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.

thumbnail Fig. 4

Wavelet coherence between long series, redated using the wiggle matching, considered here: panel a: 10 Be EDML vs. 14C, panel b: 10Be GRIP vs. 14C, and panel c: 10Be GRIP vs. 10Be 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.
thumbnail Fig. 5

Panel a: original (black) mean 14C 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 dipole-like, aligned with the geographical axis, and its virtual axial dipole moment (VADM) is 8 × 1022 A m2. 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 × 1022 A m2. 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.

thumbnail Fig. 6

Example of the χ2 vs. ϕ dependence for 805 AD for the 14C series. The dashed lines represent the 68% confidence interval for ϕ.

4.2 Reconstruction based solely on 14C

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 Qi(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 Mj(t) from Usoskin et al. (2016), which include uncertainties and cover the range of available archeomagnetic reconstructions to calculate 106 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 best-fit ϕ(t) that minimizes the value of χ2:

(2)

where Q′(ϕ) is the value of Q* computed for a given value of ϕ, which was scanned over the range 0–2000 MeV. The best-fit 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 best-fit ϕvalues and its uncertainties is shown in Fig. 6.

  • 3.

    The series of reconstructed ϕ based on 14C (ϕ14C) was then computed, along with the uncertainties, as shown in Fig. 7.

4.3 Comparison between the long 10Be and 14C series

Next we compared the long-term behavior, in the sense of the scaling factors, of the long-running 10 Be series vs. the reference 14C series. Firstwe considered a 1000-yr 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 14C 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, R2 = 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 10Be 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 14C series over the entire period of their overlap, with κ varying from −5% to +10%. A shallow wavy variability can be noticed, with a quasi-period 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 R2 = 0.52 (p = 0.12). This suggests that the under-corrected geomagnetic field effect is most likely only related to the GRIP series.

The discrepancy between GRIP and 14C 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 long-term 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 GRIP-based series (panel b) is well coherent with the 14C 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 14C 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 long-term variability as it exists in the original data.

It is interesting that while the coherence between EDML- and 14 C-based 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 (R2 < 0.01). A new full-size model of the atmospheric transport and deposition of 10 Be in polar regions, particularly in Greenland, is needed to resolve the issue.

thumbnail Fig. 7

Series of the modulation potential ϕ computed based only on the 14C data. Shading denotes the 68% confidence interval.
thumbnail Fig. 8

The 10Be scaling factor κ, with 68%uncertainties indicated by the grey shading, in the sliding 1000-yr 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 10Be 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 14C for the same period. The scaling factor we found for the GRIP data is close to unity, κ0 = 1.1. This implies that 14C and 10Be 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 best-fit 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.

thumbnail Fig. 9

Wavelet coherence between individual long-term series of the modulation potential ϕ reconstructed from individual long-term cosmogenic series (shown in Fig. 12). Panel a: 10 Be EDML vs. 14C; panel b: 10Be GRIP vs. 14C; panel c: 10Be GRIP vs. 10Be EDML. Notations are similar to Fig. 4.

4.5 Combined-record reconstruction of ϕ

In the preceding section we calculated three series of ϕ from different cosmogenic isotope datasets, using their inter-calibration to the mean ϕ value obtained from the 14C 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 Bayesian-based 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), 106 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 106 ensemble members for each time point t.

  • 4.

    For each t, the value of χ2(ϕ) was calculated as

(3)

where the index i takes values 1 through 3 for the 14C, 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 well-defined 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 best-fit 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.

thumbnail 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 cross-hairs.
thumbnail Fig. 10

Dependence of χ2 (Eq. (3)) on ϕ for the decade centered at 2095 BC: panel a: 14C with no scaling, panel b: 10Be 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 best-guess 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 14C series for the period of their overlap. We then slightly varied the κ−value around these best-guessed values to calculate the corresponding χ2 for each individual series, as shown in Fig. 13. The vertical solid lines mark the best-fit κ-values that minimize the five-point smoothed χ2. The best-fitscaling 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.

thumbnail Fig. 12

Reconstruction of the modulation potential ϕ using only individual cosmogenic isotope series (color curves as denoted in the legend) with the best-fit 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 Fo, following an established procedure (e.g., Usoskin et al. 2003, 2016; Solanki et al. 2004). Applying the updated SATIRE-M 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, SNi, 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 upper-bound (68%) ϕ values (defined as described above) for each decade into SNs. The reconstructed SN series3 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 GRIP-based 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 & Nesme-Ribes 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 low-activity 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 high-activity 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 multi-peak (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.

thumbnail Fig. 13

Values of χ2 vs. the scaling factors κ for four short 10Be series, as marked in each panel. The five-point running mean curves are shown in red. The locations of the best-fit κ−values are shown as straight black lines.
thumbnail 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.
thumbnail Fig. 15

Probability kernel density function estimate (black dots) of the reconstructed decadal SNs (Gaussian kernel, width = 3) and its multi-peak fit (red).

6 Conclusions

We have provided a new fully consistent multi-proxy 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 andup-to-date models of isotope production and transport or deposition as well as a recent archeomagnetic model. We used six 10Be 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 14C data. All employed cosmogenic isotope series were reduced to the reference geomagnetic field conditions.

We have also tested the stability of the two long 10Be 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 GRIP-based series appears to anticorrelate with the long-term geomagnetic VADM series, implying a possible undercorrection for the geomagnetic shielding effect for the GRIP location, while EDML and 14 C-based 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 full-size 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 long-term 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 10Be 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 long-term evolution of the solar variability and poses robust constraints on the development of solar and stellar dynamo models as well as solar-terrestrial 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

  1. Adolphi, F., & Muscheler, R. 2016, Clim. Past, 12, 15 [CrossRef] [Google Scholar]
  2. Asvestari, E., Gil, A., Kovaltsov, G. A., & Usoskin, I. G. 2017, J. Geophys. Res.: Space Phys., 122, 9790 [Google Scholar]
  3. Bard, E., Raisbeck, G., Yiou, F., & Jouzel, J. 1997, Earth Planet. Sci. Lett., 150, 453 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bazilevskaya, G. A., Cliver, E. W., Kovaltsov, G. A., et al. 2014, Space Sci. Rev., 186, 409 [Google Scholar]
  5. Bazin, L., Landais, A., Lemieux-Dudon, B., et al. 2013, Clim. Past, 9, 1715 [NASA ADS] [CrossRef] [Google Scholar]
  6. Beer, J., Siegenthaler, U., Oeschger, H., Bonani, G., & Finkel, R. 1988, Nature, 331, 675 [NASA ADS] [CrossRef] [Google Scholar]
  7. Beer, J., Blinov, A., Bonani, G., Hofmann, H., & Finkel, R. 1990, Nature, 347, 164 [NASA ADS] [CrossRef] [Google Scholar]
  8. 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]
  9. Berggren, A.-M., Beer, J., Possnert, G., et al. 2009, Geophys. Res. Lett., 36, L11801 [Google Scholar]
  10. Burger, R., Potgieter, M., & Heber, B. 2000, J. Geophys. Res., 105, 27447 [Google Scholar]
  11. Caballero-Lopez, R., & Moraal, H. 2004, J. Geophys. Res., 109, A01101 [Google Scholar]
  12. Cain, W. F., & Suess, H. E. 1976, J. Geophys. Res., 81, 3688 [NASA ADS] [CrossRef] [Google Scholar]
  13. Choudhuri, A. R., & Karak, B. B. 2012, Phys. Rev. Lett., 109, 171103 [NASA ADS] [CrossRef] [Google Scholar]
  14. Clette, F., Svalgaard, L., Vaquero, J., & Cliver, E. 2014, Space Sci. Rev., 186, 35 [NASA ADS] [CrossRef] [Google Scholar]
  15. Cohen, J., Cohen, P., West, S., & Aiken, L. 2003, Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences (New York: Routledge) [Google Scholar]
  16. 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]
  17. Delaygue, G., & Bard, E. 2011, Clim. Dyn., 36, 2201 [Google Scholar]
  18. Ebisuzaki, W. 1997, J. Clim., 10, 2147 [NASA ADS] [CrossRef] [Google Scholar]
  19. Eddy, J. 1976, Science, 192, 1189 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  20. Genevey, A., Gallet, Y., Constable, C. G., Korte, M., & Hulot, G. 2008, Geochem. Geophys. Geosyst., 9, Q04038 [NASA ADS] [CrossRef] [Google Scholar]
  21. Goslar, T. 2003, PAGES News, 11, 12 [CrossRef] [Google Scholar]
  22. Hathaway, D. H. 2015, Liv. Rev. Sol. Phys., 12, 4 [Google Scholar]
  23. Heikkilä, U., Beer, J., & Feichter, J. 2009, Atmos. Chem. Phys., 9, 515 [Google Scholar]
  24. Heikkilä, U., Beer, J., Abreu, J. A., & Steinhilber, F. 2013, Space Sci. Rev., 176, 321 [NASA ADS] [CrossRef] [Google Scholar]
  25. Herbst, K., Kopp, A., Heber, B., et al. 2010, J. Geophys. Res., 115, D00I20 [Google Scholar]
  26. Herbst, K., Muscheler, R., & Heber, B. 2017, J. Geophys. Res.: Space Phys., 122, 23 [Google Scholar]
  27. Hoek, W. Z., & Bohncke, S. J. P. 2001, Quat. Sci. Rev., 20, 1251 [NASA ADS] [CrossRef] [Google Scholar]
  28. Horiuchi, K., Uchida, T., Sakamoto, Y., et al. 2008, Quat. Geochronol., 3, 253 [CrossRef] [Google Scholar]
  29. Inceoglu, F., Simoniello, R., Knudsen, V. F., et al. 2015, A&A, 577, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Johnsen, S. J., Dahl-Jensen, D., Dansgaard, W., & Gundestrup, N. 1995, Tellus Ser. B, 47, 624 [NASA ADS] [CrossRef] [Google Scholar]
  31. Johnsen, S. J., Dahl-Jensen, D., Gundestrup, N., et al. 2001, J. Quat. Sci., 16, 299 [CrossRef] [Google Scholar]
  32. Käpylä, M. J., Käpylä, P. J., Olspert, N., et al. 2016, A&A, 589, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Knudsen, M. F., Riisager, P., Donadini, F., et al. 2008, Earth Planet. Sci. Lett., 272, 319 [NASA ADS] [CrossRef] [Google Scholar]
  34. Knudsen, M. F., Riisager, P., Jacobsen, B. H., et al. 2009, Geophys. Res. Lett., 36, L16701 [NASA ADS] [CrossRef] [Google Scholar]
  35. Kovaltsov, G., & Usoskin, I. 2010, Earth Planet. Sci. Lett., 291, 182 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kovaltsov, G., Mishev, A., & Usoskin, I. 2012, Earth Planet. Sci. Lett., 337, 114 [Google Scholar]
  37. Küker, M., Arlt, R., & Rüdiger, G. 1999, A&A, 343, 977 [NASA ADS] [Google Scholar]
  38. Licht, A., Hulot, G., Gallet, Y., & Thébault, E. 2013, Phys. Earth Planet. Inter., 224, 38 [NASA ADS] [CrossRef] [Google Scholar]
  39. McCracken, K. G., & Beer, J. 2015, Sol. Phys., 290, 3051 [Google Scholar]
  40. McCracken, K., McDonald, F., Beer, J., Raisbeck, G., & Yiou, F. 2004, J. Geophys. Res., 109, 12103 [NASA ADS] [CrossRef] [Google Scholar]
  41. Mekhaldi, F., Muscheler, R., Adolphi, F., et al. 2015, Nat. Commun., 6, 8611 [NASA ADS] [CrossRef] [Google Scholar]
  42. Miyake, F., Nagaya, K., Masuda, K., & Nakamura, T. 2012, Nature, 486, 240 [NASA ADS] [CrossRef] [Google Scholar]
  43. Miyake, F., Masuda, K., & Nakamura, T. 2013, Nat. Commun., 4, 1748 [NASA ADS] [CrossRef] [Google Scholar]
  44. Miyake, F., Jull, A. J. T., Panyushkina, I. P., et al. 2017, Proc. Natl. Acad. Sci., 114, 881 [Google Scholar]
  45. Moss, D., Sokoloff, D., Usoskin, I., & Tutubalin, V. 2008, Sol. Phys., 250, 221 [NASA ADS] [CrossRef] [Google Scholar]
  46. Muscheler, R., Beer, J., Wagner, G., et al. 2004, Earth Planet. Sci. Lett., 219, 325 [NASA ADS] [CrossRef] [Google Scholar]
  47. Muscheler, R., Joos, F., Beer, J., et al. 2007, Quat. Sci. Rev., 26, 82 [Google Scholar]
  48. Muscheler, R., Adolphi, F., & Knudsen, M. F. 2014, Quat. Sci. Rev., 106, 81 [NASA ADS] [CrossRef] [Google Scholar]
  49. Nilsson, A., Holme, R., Korte, M., Suttie, N., & Hill, M. 2014, Geophys. J. Int., 198, 229 [NASA ADS] [CrossRef] [Google Scholar]
  50. Owens, M. J., Usoskin, I., & Lockwood, M. 2012, Geophys. Res. Lett., 39, L19102 [Google Scholar]
  51. Pavlov, A. K., Blinov, A. V., Frolov, D. A., et al. 2017, J. Atmos. Sol.-Terr. Phys., 164, 308 [NASA ADS] [CrossRef] [Google Scholar]
  52. Pavón-Carrasco, F. J., Osete, M. L., Torta, J. M., & De Santis, A. 2014, Earth Planet. Sci. Lett., 388, 98 [NASA ADS] [CrossRef] [Google Scholar]
  53. Poluianov, S. V., Kovaltsov, G. A., Mishev, A. L., & Usoskin, I. G. 2016, J. Geophys. Res. (Atm.), 121, 8125 [NASA ADS] [Google Scholar]
  54. Raisbeck, G., Yiou, F., Jouzel, J., & Petit, J. 1990, Roy. Soc. London Philos. Trans. Ser. A, 330, 463 [NASA ADS] [CrossRef] [Google Scholar]
  55. Rasmussen, S. O., Bigler, M., Blockley, S. P., et al. 2014, Quat. Sci. Rev., 106, 14 [NASA ADS] [CrossRef] [Google Scholar]
  56. Reimer, P. J., Baillie, M. G. L., Bard, E., et al. 2009, Radiocarbon, 51, 1111 [Google Scholar]
  57. Ribes, J., & Nesme-Ribes, E. 1993, A&A, 276, 549 [NASA ADS] [Google Scholar]
  58. Roth, R., & Joos, F. 2013, Clim. Past, 9, 1879 [NASA ADS] [CrossRef] [Google Scholar]
  59. Schmitt, D., Schüssler, M., & Ferriz-Mas, A. 1996, A&A, 311, L1 [NASA ADS] [Google Scholar]
  60. Seierstad, I. K., Abbott, P. M., Bigler, M., et al. 2014, Quat. Sci. Rev., 106, 29 [NASA ADS] [CrossRef] [Google Scholar]
  61. Sigl, M., Winstrup, M., McConnell, J. R., et al. 2015, Nature, 523, 543 [NASA ADS] [CrossRef] [Google Scholar]
  62. Sokoloff, D., & Nesme-Ribes, E. 1994, A&A, 288, 293 [NASA ADS] [Google Scholar]
  63. Solanki, S. K., Usoskin, I. G., Kromer, B., Schüssler, M., & Beer, J. 2004, Nature, 431, 1084 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  64. Steinhilber, F., Abreu, J., Beer, J., et al. 2012, Proc. Natl. Acad. Sci. USA, 109, 5967 [Google Scholar]
  65. Stuiver, M. 1961, J. Geophys. Res., 66, 273 [NASA ADS] [CrossRef] [Google Scholar]
  66. Stuiver, M., & Braziunas, T. 1989, Nature, 338, 405 [NASA ADS] [CrossRef] [Google Scholar]
  67. Stuiver, M., & Quay, P. 1980, Science, 207, 11 [NASA ADS] [CrossRef] [Google Scholar]
  68. Sukhodolov, T., Usoskin, I., Rozanov, E., et al. 2017, Sci. Rep., 7, 45257 [NASA ADS] [CrossRef] [Google Scholar]
  69. Usoskin, I. G. 2017, Liv. Rev. Sol. Phys., 14, 3 [Google Scholar]
  70. Usoskin, I. G., & Kovaltsov, G. A. 2012, ApJ, 757, 92 [Google Scholar]
  71. Usoskin, I. G., Solanki, S. K., Schüssler, M., Mursula, K., & Alanko, K. 2003, Phys. Rev. Lett., 91, 211101 [Google Scholar]
  72. Usoskin, I. G., Alanko-Huotari, K., Kovaltsov, G. A., & Mursula, K. 2005, J. Geophys. Res., 110, A12108 [Google Scholar]
  73. Usoskin, I. G., Solanki, S. K., & Kovaltsov, G. A. 2007, A&A, 471, 301 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Usoskin, I. G., Horiuchi, K., Solanki, S., Kovaltsov, G. A., & Bard, E. 2009, J. Geophys. Res., 114, A03112 [Google Scholar]
  75. Usoskin, I. G., Kromer, B., Ludlow, F., et al. 2013, A&A, 552, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Usoskin, I. G., Hulot, G., Gallet, Y., et al. 2014, A&A,562, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Usoskin, I. G., Arlt, R., Asvestari, E., et al. 2015, A&A, 581, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Usoskin, I. G., Gallet, Y., Lopes, F., Kovaltsov, G. A., & Hulot, G. 2016, A&A, 587, A150 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. 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]
  80. Veres, D., Bazin, L., Landais, A., et al. 2013, Clim. Past, 9, 1733 [CrossRef] [Google Scholar]
  81. Vieira, L. E. A., & Solanki, S. K. 2010, A&A, 509, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Vieira, L. E. A., Solanki, S. K., Krivova, N. A., & Usoskin, I. 2011, A&A, 531, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Vinther, B. M., Clausen, H. B., Johnsen, S. J., et al. 2006, J. Geophys. Res.: Atm., 111, D13102 [NASA ADS] [CrossRef] [Google Scholar]
  84. Vonmoos, M., Beer, J., & Muscheler, R. 2006, J. Geophys. Res., 111, A10105 [NASA ADS] [CrossRef] [Google Scholar]
  85. Wu, C.-J., Krivova, N. A., Solanki, S. K., & Usoskin, I. G. 2018, A&A, submitted [Google Scholar]
  86. Yiou, F., Raisbeck, G., Baumgartner, S., et al. 1997, J. Geophys. Res., 102, 26783 [NASA ADS] [CrossRef] [Google Scholar]
1

Although naked-eye observations of sunspots are available, they do not provide quantitative assessments (Usoskin 2017).

2

The significance is estimated using the non-parametric random-phase method by Ebisuzaki (1997).

3

Available as a table at the CDS and at the MPS sun-climate web-page http://www2.mps.mpg.de/projects/sun-climate/data/SN_composite.txt.

4

File SN_y_tot_V2.0.txt available at http://www.sidc.be/silso/infosnytot.

All Tables

Table 1

Temporal coverage, cadence, type of data (production rate, PR; depositional flux, D; or concentration, C) and the resulting best-fit scaling factors, κ (Sect. 4), of the cosmogenic isotope series.

In the text
Table 2

Wiggles (central date and the length in years) used for synchronization of the various 10 Be time series to 14C (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.

In the text
Table 3

Squared correlation coefficients between the six 10Be series and the 14C series for the originally dated (Ro) and synchronized (Rs) series.

In the text

All Figures

thumbnail 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
thumbnail Fig. 2

Typical example of a wiggle in the 14C (red) series compared with the original EDML 10Be 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 14C. The black sold curve shows the EDML 10Be series after the synchronization (see Sect. 3.2).
In the text
thumbnail Fig. 3

Example of the calculation of the best-fit 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 14C and EDML 10Be series for the ± 100-yr time window around the center of the wiggle.
In the text
thumbnail Fig. 4

Wavelet coherence between long series, redated using the wiggle matching, considered here: panel a: 10 Be EDML vs. 14C, panel b: 10Be GRIP vs. 14C, and panel c: 10Be GRIP vs. 10Be 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
thumbnail Fig. 5

Panel a: original (black) mean 14C 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
thumbnail Fig. 6

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

Series of the modulation potential ϕ computed based only on the 14C data. Shading denotes the 68% confidence interval.
In the text
thumbnail Fig. 8

The 10Be scaling factor κ, with 68%uncertainties indicated by the grey shading, in the sliding 1000-yr 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
thumbnail Fig. 9

Wavelet coherence between individual long-term series of the modulation potential ϕ reconstructed from individual long-term cosmogenic series (shown in Fig. 12). Panel a: 10 Be EDML vs. 14C; panel b: 10Be GRIP vs. 14C; panel c: 10Be GRIP vs. 10Be EDML. Notations are similar to Fig. 4.
In the text
thumbnail 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 cross-hairs.
In the text
thumbnail Fig. 10

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

Reconstruction of the modulation potential ϕ using only individual cosmogenic isotope series (color curves as denoted in the legend) with the best-fit 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
thumbnail Fig. 13

Values of χ2 vs. the scaling factors κ for four short 10Be series, as marked in each panel. The five-point running mean curves are shown in red. The locations of the best-fit κ−values are shown as straight black lines.
In the text
thumbnail 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
thumbnail Fig. 15

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

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.