EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 587 (March 2016) A&A, 587 (2016) A150 Full HTML
Free Access
Issue
A&A
Volume 587, March 2016
Article Number A150
Number of page(s) 10
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201527295
Published online 04 March 2016

© ESO, 2016

1. Introduction

Solar activity, which generically includes different manifestations of fluctuating processes in the solar convection zone on the surface and in the corona, varies on timescales from seconds to millennia. Its long-term variability, with timescales longer than several centuries, can only be studied using indirect proxies, such as cosmogenic radionuclides (e.g., Beer et al. 2012; Usoskin 2013). For this purpose, the most useful cosmogenic isotopes are radiocarbon 14C and beryllium 10Be. These isotopes allow reconstructing the past solar activity over the Holocene period, which spans the past 11 millennia, that is, since the last deglaciation. Recently, several long-term solar activity reconstructions have been published (e.g., Solanki et al. 2004; Vonmoos et al. 2006; Muscheler et al. 2007; Usoskin et al. 2007; Steinhilber et al. 2012). They revealed variability in the solar activity on centennial and millennial scales, ranging from very low (almost spotless Sun) to high activity levels. However, since these reconstructions rely on slightly different basic models and different datasets, they sometimes differ in details and overall levels. In particular, although the existence of grand minima and maxima in solar activity has been known for a long time (see the review of Usoskin 2013), it has remained a matter of debate whether grand minima and maxima are separate activity modes of the solar dynamo or simply non-Gaussian tails of its variability (e.g., Moss et al. 2008; Passos et al. 2014; Karak et al. 2015).

To settle this question, a new approach was recently developed (Usoskin et al. 2014), which takes into account the full range of uncertainties associated with a modern reconstruction of the 14C global production rate (Roth & Joos 2013), an accurate millennial-scale archeomagnetic field reconstruction (Licht et al. 2013), and a detailed 14C production model (Kovaltsov et al. 2012). Based on an analysis of the probability density function, the resulting reconstruction, limited to the past 3000 years (Licht et al. 2013), made it possible for the first time to show that grand minima of solar activity correspond to a distinct operational mode of the solar activity (Usoskin et al. 2014). Because of poor statistics, however, the result was inconclusive for the nature of the grand maxima. Constraining the solar activity over longer period is more problematic, in particular because of uncertainties acknowledged before 500–1000 BC in paleo- and archeomagnetic field reconstructions (Snowball & Muscheler 2007). However, a large set of new archeo- and paleointensity data was acquired in the past few years (see Appendix A), which, as we argue, can improve our knowledge of the geomagnetic dipole moment evolution over most of the Holocene. Here we take advantage of the new data to better constrain the long-term solar activity, as revealed from the use of a synthetic index of relative sunspot numbers. We first extend the approach of Usoskin et al. (2014) to the last 9 millennia using the reconstruction of the 14C global production rate (Roth & Joos 2013), the 10Be GRIP dataset (Yiou et al. 1997), a new dipole moment reconstruction hereafter referred to as GMAG.9k (see Appendix A), and 14C and 10Be production models (Kovaltsov et al. 2012; Kovaltsov & Usoskin 2010). The recovered reconstructions aim at scrutinizing the variability in the solar activity over millennial and centennial timescales. We also discuss the robustness of our new solar activity reconstructions using the results derived from other recent archeo- and paleomagnetic Holocene field models. Finally, we present new results providing important observational constraints on the solar dynamo.

2. Data

2.1. Cosmogenic radionuclide records

We used two sets of cosmogenic radionuclide data (14C in tree trunks and 10Be in polar ice; panels A and B in Fig. 1, respectively) as tracers of solar activity (e.g., Beer et al. 2012; Usoskin 2013).

Radiocarbon 14C is produced in the terrestrial atmosphere by cosmic rays and then takes part in the global carbon cycle (e.g., Bard et al. 1997; Beer et al. 2012; Roth & Joos 2013). The measured quantity, the relative concentration Δ14C of radiocarbon in tree rings, needs to be corrected for the apparent decay and for the carbon cycle effect to reconstruct the 14C production rate. Here we used the 14C production rate, Q(14C), as reconstructed by Roth & Joos (2013) for the Holocene, using the globally averaged IntCal09 (Reimer et al. 2009) radiocarbon database and the dynamical BERN3D-LPJ carbon cycle model, which is a new-generation carbon-cycle climate model, featuring a 3D dynamic ocean, reactive ocean sediments, and a 2D atmosphere component coupled to the Lund-Potsdam-Jena dynamic global vegetation model. The data were reduced to the decadal temporal resolution. For the decades around years 775 AD and 994 AD, the production rate was corrected to remove the modeled contribution that is due to the occurrence of two extreme solar particle events (Usoskin et al. 2013; see also the discussion in Miyake et al. 2012; Miyake et al. 2013; and Bazilevskaya et al. 2014). Finally, we only consider data before 1900 AD because of the Suess effect, which is related to extensive burning of fossil fuel, which dilutes radiocarbon in the natural reservoirs and makes the use of the 14C data after 1900 more uncertain.

thumbnail Fig. 1

Time series of cosmogenic radionuclide data. Panel A): decadal radiocarbon 14C global production rate (Roth & Joos 2013) with the 95% confidence interval plotted in gray. Panel B): quasi-decadal variability of 10Be flux in the GRIP ice core (Yiou et al. 1997). The formal 1σ error of 7% (relative to the given value) is indicated by the error bar next to the legend.

Following Usoskin et al. (2014), we used 1000 individual realizations of the Q(14C) ensemble to describe the consequences of uncertainties in the data and carbon cycle modeling (Roth & Joos 2013). The corresponding production rate is shown in Fig. 1a (mean (black) and 95% range (gray shaded area) of the 1000 realizations).

The radionuclide 10Be is produced by cosmic rays in the atmosphere through spallation reactions (Beer et al. 2012). It becomes attached to aerosols and is relatively quickly precipitated to the ground. Because of this fast precipitation, it is not completely mixed in the atmosphere and is subject to some complicated transport. We relied on the parameterization of the atmosphere transport and deposition of beryllium proposed by Heikkilä et al. (2009). Here we used a long series of 10Be depositional flux measured in central Greenland in the framework of the Greenland Ice Core Project (GRIP) for the period before 1645 AD (Yiou et al. 1997). We considered the mean data set reduced to quasi-decadal time resolution. The corresponding rate is shown in Fig. 1B, where we also plot a formal 1σ uncertainty error (estimated to be 7% in relative terms, Yiou et al. (1997)).

2.2. Axial dipole evolution over the past 9000 years

Two approaches can be used to constrain the axial dipole moment evolution over the past few millennia. The first consists in constructing global geomagnetic field models in the form of time-varying series of Gauss coefficients by taking advantage of all available archeo- and paleomagnetic data (see, e.g., Korte & Constable 2005, 2011; Korte et al. 2011; Licht et al. 2013; Pavón-Carrasco et al. 2014; Nilsson et al. 2014) and using the corresponding axial dipole component . Differences among these models mainly come from the treatment applied to the data, in particular the way experimental and dating uncertainties are being handled (see discussion and details in the references above). In the present study, we considered three recent models, referred to as A_FM (Licht et al. 2013), SHA.DIF.14k (Pavón-Carrasco et al. 2014), and pfm9k.1a/b (Nilsson et al. 2014). We note that A_FM and SHA.DIF.14k were built using archeomagnetic and volcanic data sets, whereas pfm9k.1a/b also took sedimentary data into account. We assume that pfm9k.1a/b supersedes the slightly older CALS10k.1b field model constructed by Korte et al. (2011) using practically the same dataset.

The second approach is based on archeo- and paleointensity data collected worldwide, using archeological artifacts and volcanic rocks. The data are transformed into virtual axial dipole moments (VADM), which are then carefully weighted to produce a worldwide average VADM. To be valid, this “paleomagnetic” approach requires a dual averaging of the data, both in time and space, to best smooth out non-dipole field components (for a discussion, see, e.g., Korte & Constable 2005; Genevey et al. 2008). Here we consider the two most recent mean VADM curves built in this way, one by Genevey et al. (2008), which encompasses the past 3000 years, and a second one by Knudsen et al. (2008), which covers the entire Holocene. In addition, and because quite a large number of additional intensity data have recently been collected, an updated mean VADM curve was also produced for the purpose of the present study. For this we used the GEOMAGIA50.v3 data base (Brown et al. 2015), to which we added or modified about 390 individual intensity data points (see more details in Appendix A). The new data compilation contains 4764 intensity values dated to between 7000 BC and 2000 AD.

To build this new VADM curve, we first carried out a series of computations to explore the effects of changing the width of the temporal averaging sliding windows (200, 500, and 1000 years) and the size of the region of spatial weighting (over regions of 10°, 20°, and 30° width). We also used a bootstrap technique to simulate the effect of noise in the intensity data within their age uncertainties and within their 2σ experimental error bars (see, e.g., Korte et al. 2009; Thébault & Gallet 2010). For each set of parameters, an ensemble of 1000 individual curves was computed, allowing us to obtain at each epoch the mean VADM and its standard deviation, together with the maximum and minimum VADM from the 1000 possible values. Results from different computations are shown in Appendix A. This analysis revealed that VADMs derived using sliding windows with widths of 500 years and 1000 years are very similar. VADM values also appear to be relatively insensitive to the size of the area chosen for the regional averaging. Some differences, but still quite limited, are observed when the width of the sliding window is reduced to 200 years, which reveals enhanced variations compared to that obtained when using sliding windows of larger widths, as expected. Averaging over such a narrow window, however, is reasonable only for the most recent time interval (here, the past 3500 years), which is documented by a rather large number of data points (3756 among the 4764 available data) with a relatively wide but still uneven geographical distribution (see for instance Fig. 4 in Knudsen et al. 2008; Genevey et al. 2008). Because of this, we finally decided to build a composite VADM variation curve, which we hereafter refer to as GMAG.9k. This curve was computed for every 10 years, using sliding windows of widths 200 years between 1500 BC and 2000 AD and 500 years between 7000 BC and 1500 BC, with a spatial weighting over regions of 10° in size for both time intervals (numerical values for this composite curve are provided in the CDS, Tables C.1 and C.2).

Figure 2 shows this new GMAG.9k curve. This updated VADM curve does not markedly differ from previous dipole moment curves. Its behavior over the past 3000 years is very similar to that of the VADM curve obtained by Genevey et al. (2008), who also used a spatial weighting, but with a smaller number of different and distant regions (30° in size) and a smaller dataset selected based on specific quality criteria. We note, however, that the new VADM curve tends to lie slightly below that of Genevey et al. (2008). Differences with the VADM variation curve of Knudsen et al. (2008) are larger. However, the latter was computed using sliding windows of width 500 years with no geographical weighting (these authors concluded that they had no significant bias despite the poor spatial data distribution). Comparison with the A_FM, SHA-DIF.14k and pfm9k.1a/b curves (Licht et al. 2013; Pavón-Carrasco et al. 2014; Nilsson et al. 2014) reveals a fairly similar evolution except for a small offset between ≈500 AD and ≈1500 AD. This offset also persists when a sliding-window duration of 500 years is used for the VADM computations. However, the A_FM, SHA-DIF.14k and pfm9k.1a/b curves generally lie within the envelope of possible VADM values (Fig. 2). The same observation holds for all the periods before 1000 BC (Fig. 2A). This encouragingly suggests that relying on the GMAG.9k ensemble of 1000 individual VADM curves to reconstruct the solar activity as done in the present study, can be considered as a conservative procedure from a geomagnetic point of view.

thumbnail Fig. 2

Time series of the axial dipole moment reconstructions spanning the past 9000 years (panel A)), with a zoom for the last 3200 years (panel B)). The black solid line depicts GMAG.9k (the reconstruction presented and used in this work) with ± 1σ and the full range variability presented by the gray shading and the hatching, respectively. Other reconstructions shown are (Licht et al. 2013) denoted as AF_M, (Genevey et al. 2008) denoted as G08, (Nilsson et al. 2014) denoted as pfm9k.1b and pfm9k.1a, (Knudsen et al. 2008) denoted as Kn08, and (Pavón-Carrasco et al. 2014) denoted as SHA-DIF.14k. For better readability, error bars were omitted for these curves, but this does not affect the discussion of the results (see text).

3. Reconstructing the solar activity

Since details of the reconstruction of solar activity from cosmogenic nuclides are described elsewhere (e.g., Beer et al. 2012; Usoskin 2013), we only briefly describe this reconstruction and recall important relevant information. Cosmogenic isotopes are produced by cosmic rays in the terrestrial atmosphere. Since cosmic rays are modulated by solar magnetic activity, the variability of cosmogenic isotope production reflects the latter. However, two terrestrial processes may disturb this relation. One process is additional shielding of Earth from cosmic rays by the geomagnetic field, whose changes must be known independently. For this purpose, we relied on the GMAG.9k axial dipole evolution constructed as described in Sect. 2.2. As discussed in Sect. 2.1, another important process is transport and deposition of the nuclides in the terrestrial system. Because of the poorly known details of climate variability in the past, the related transport models are commonly adjusted to modern conditions, which may lead to some uncertainties in the older part of the time interval.

We converted the cosmogenic isotope production rate to the Galactic cosmic ray (GCR) flux variability using recent production models. The global production of 14C was modeled using the model of Kovaltsov et al. (2012), while the production of 10Be was modeled using an updated version of the model of Kovaltsov & Usoskin (2010). Cosmic ray variability was calculated in terms of the heliospheric modulation potential (see definitions and formalism in Usoskin et al. 2005), also considering α-particles and heavier species of cosmic rays (Webber & Higbie 2009). This modulation potential was furthermore converted into decadal (solar-cycle averaged) sunspot number through the open solar magnetic flux model (Solanki et al. 2000; Krivova et al. 2007), which relates the solar surface magnetic cycle to the emerging magnetic flux (Cameron & Schüssler 2015). We note that the overall reconstruction method used here is similar to that previously used by Usoskin et al. (2014).

Uncertainties were assessed straightforwardly by computing a large ensemble of individual reconstructions. We used the set of 1000 time-varying individual archeomagnetic reconstructions of GMAG.9k (see Sect. 2.2), which account in particular for experimental and age uncertainties. This ensemble was cross-used with a similar ensemble of 1000 production rates of 14C to account for measurement and compilation uncertainties in the IntCal09 and SHCal04 data, in the air-sea gas exchange rate, in the terrestrial primary production, and in the closure of the atmospheric CO2 budget (Roth & Joos 2013). GMAG.9k was also cross-used in the same way with a set of 1000 10Be series. In that case, however, 10Be series of decadal values were generated around the mean provided by GRIP using normally distributed random numbers (with a standard deviation equal to 7% of the mean value, Yiou et al. 1997) to reflect known errors. In both cases, all possible combinations of the ensembles yielded 106 series of the reconstructed heliospheric modulation potential, next converted into 106 series of sunspot numbers. These series reflect the error propagation through all the intermediate steps. An additional random error with σ = 0.5 was finally added to each computed sunspot number to account for the small possible error related to the conversion between the modulation potential and the solar open magnetic flux (Solanki et al. 2004).

Decadal sunspot numbers reconstructed in this way from the 14C and 10Be data are henceforth denoted as SN-14C and SN-10Be, respectively. Ensemble means of these SN-14C and SN-10Be are shown in Fig. 3A. These reconstructions agree well with the earlier reconstruction (Usoskin et al. 2014) that cover the last 3000 years.

thumbnail Fig. 3

Panel A): raw reconstructions of the sunspot numbers (mean curves) SN-14C (blue) and SN-10Be (red), compared to the recent 3 kyr reconstruction (Usoskin et al. 2014,– green curve). Panel B): first component of the singular spectrum analysis (SSA – see Appendix B) for the SN-14C (blue) and SN-10Be (red) series. The shaded areas depict the uncertainties related to the parameter L of the SSA. Panel C): same as in panel B), but for the second SSA components of the SN-14C (blue) and SN-10Be (red) series. The large dots and red stars denote times of the grand minima (see Table 1) and grand maxima (Table 2), respectively.

We checked the influence of the choice of axial dipole moment reconstruction on the SN-14C reconstruction by also considering alternative geomagnetic field models. Results are shown in Fig. 4. This figure clearly shows that all 14C-based SN reconstructions lie close to each other and reveal a common general pattern. In particular, Figs. 3A and 4 both display a long-term trend in the 14C-based SN reconstructions over the past 9000 years. This trend, however, is different from that of the 10Be-based SN reconstruction. What causes these long-term trends is unclear: they may reflect various combinations of climate effects or improperly corrected geomagnetic field effects. We discuss this important point in the next section.

thumbnail Fig. 4

Comparison of alternative SN-14 sunspot number reconstructions when relying on different axial dipole reconstructions, (same notations as in Fig. 2). Only the mean curves of the corresponding ensembles are shown.

Table 1

Grand minima with their centers, approximate duration, and comments (1 – listed in Usoskin et al. 2007; 2 – listed in Inceoglu et al. 2015).

Table 2

Grand maxima with their centers, approximate duration, and comments (1 – listed in Usoskin et al. 2007; 2 – listed in Inceoglu et al. 2015).

4. Long-term behavior

To investigate the long-term behavior of the reconstructed solar activity, we relied on the singular spectral analysis (SSA) method (Vautard & Ghil 1989; Vautard et al. 1992) as described in Appendix B. This SSA was applied to both the SN-14C and SN-10Be series, and the corresponding two first SSA components are shown in Figs. 3B and C. The robustness of these SSA results has also been assessed by considering a wide range of values for the embedding dimension L. The resulting uncertainties are indicated by means of shaded areas in these figures.

thumbnail Fig. 5

Corrected sunspot number reconstructions SN-14C-C (panel A)) and SN-10Be-C (panel B)), after removing long-term trends (see Sect. 4.1). The black curves and the gray shading depict the mean and the 95% range (over 106 ensemble members) of the reconstructions, respectively, and the green curve represents the 3 kyr reconstruction by Usoskin et al. (2014). Stars and circles denote grand maxima and minima, respectively, as in Fig. 3. Tables for this plot are available at the CDS, including the mean values and the uncertainties of the sunspot numbers reconstructed here as shown in Fig. 5.

4.1. Multimillennial trend: possible climate influence

Here, we first consider the long-term primary SSA components of the SN-14C and SN-10Be series. These are shown in Fig. 3B. The shaded areas show the full range of computations for L-values ranging between 150 and 200 for 14C and between 120 and 170 for 10Be. These primary components are well identified.

Just as clearly, it also appears that these primary components are different for the two series: SN-14C yields a single, nearly symmetric wave along the entire time interval of 9 millennia with a range of about 20 in sunspot number, while SN-10Be yields a nearly monotonous trend within the same range of about 20 in sunspot number. The fact that these trends differ so much implies that they can hardly be related to a common process. This makes terrestrial processes, in particular transport and deposition, a much more likely cause. Differences in the very long term Holocene trends between the two isotopes have previously been noted and ascribed to such terrestrial processes (Vonmoos et al. 2006; Usoskin et al. 2009; Steinhilber et al. 2012; Inceoglu et al. 2015). Climate change, in particular, is a likely cause because it affects the two isotopes in very different ways (Beer et al. 2012), with 14C being sensitive to long-term changes in the ocean circulation (e.g., Hua et al. 2015), while 10Be is mainly sensitive to large-scale atmospheric dynamics. In any case, it is quite clear that these long-term trends are unlikely to be of solar origin. For this reason, we decided to remove them from our original SN-14C and SN-10Be series to produce what we hereafter refer to as the SN-14C-C and SN-10Be-C series, where the last “C” stands for “corrected”. The corresponding solar activity reconstructions are shown in Fig. 5. (We note that since these reconstructions were corrected for long-term trends, they only reflect relative changes within the solar activity, strictly speaking).

thumbnail Fig. 6

Decadally averaged sunspot numbers for the period 1600–1900 AD. The thick gray curve represents this work (uncertainties are not shown). Other curves correspond to sunspot number series: L14 (Lockwood et al. 2014), H98 (Hoyt & Schatten 1998), SS15 (Svalgaard & Schatten 2015), and the C14 international sunspot number (v.2) scaled with a factor 0.6 (Clette et al. 2014).

We also show in Fig. 6 the comparison of the sunspot number series for the period 1600–1900 AD reconstructed here from 14C with other series based on sunspot counts and drawings (Hoyt & Schatten 1998; Lockwood et al. 2014; Clette et al. 2014; Svalgaard & Schatten 2015). All series exhibit the same centennial-scale evolution, while at shorter timescales the variations differ from one series to the next. However, the reconstruction clearly lies within the range of the other sunspot series.

4.2. Common signature of the Hallstatt cycle

We now consider the second SSA components as shown in Fig. 3C. These components are dominated by a ≈2400 yr periodicity that is remarkably coherent between the SN-14C and SN-10Be series. The formal Pearson correlation coefficient between the two curves shown in Fig. 3C is 0.77 ± 0.01, which is highly significant with p< 10-5 estimated using the non-parametric random phase method (Ebisuzaki 1997; Usoskin et al. 2006).

To confirm the significance of this observation, we also carried out a wavelet coherence analysis of the SN-14C and SM-10Be series. The wavelet coherence is a normalized cross-spectrum of the two series and provides a measure of their covariance in time-frequency domain. It was calculated using the Morlet basis and a code originally provided by Grinsted et al. (2004), but modified to adopt the non-parametric random phase method for assessing confidence level (Ebisuzaki 1997; Usoskin et al. 2006). The corresponding wavelet coherence is displayed in Fig. 7. It shows that while the coherence is strong but intermittent at shorter timescales and nearly absent at the longest timescales (cf. Usoskin et al. 2009), there is a wide band of very high coherence that is in phase along the entire time interval for the periods of 2000–3000 years, consistent with the periods seen in the second SSA component plotted in Fig. 3C.

It is therefore clear that the ≈2400 yr quasi-periodicity is common to both series and that it dominates their super-millennial timescale variability. It is related to the so-called Hallstatt cycle that is known in Δ14C (e.g., Damon & Sonett 1991; Vasiliev & Dergachev 2002; Ma 2007), but has been poorly documented until now in the 10Be data (McCracken et al. 2011; Hanslmeier et al. 2013). We also note that the principal component analysis applied by Steinhilber et al. (2012) to a composite solar activity reconstruction likewise revealed a Hallstatt cycle in the heliospheric modulation potential that is synchronous with the reconstruction shown in Fig. 3C, although it was not explicitly characterized.

This Hallstatt cycle has so far either been ascribed to climate variability (Vasiliev & Dergachev 2002) or to geomagnetic fluctuations, particularly geomagnetic pole migration (Vasiliev et al. 2012). However, the fact that the signal we found is in phase and of the same magnitude in the two cosmogenic isotope reconstruction implies that it can hardly be of climatic origin. As already pointed out, 14C and 10Be respond differently to climate changes. In particular, 14C is mostly affected by the ocean ventilation and mixing, while 10Be (in particular its deposition on central Greenland) is mainly affected by the large-scale atmospheric circulation, particularily in the North Atlantic region (Field et al. 2006; Heikkilä et al. 2009). It can also hardly be of geomagnetic origin and related to geomagnetic pole migration, since 14C is globally (hemispherically) mixed in the terrestrial system and insensitive to the migration of these poles. To reproduce the observed Hallstatt cycle, the dipole moment (VADM) would have to vary, with the corresponding period, in a range of ≈2 × 1022 A m2 (i.e., by about 20%). This is not supported by any geomagnetic field reconstruction (Fig. 2). The SSA analysis of the geomagnetic series (not shown) does not yield the Hallstatt cycle.

We thus conclude that the ≈2400-yr Hallstatt cycle is most likely a property of the long-term solar activity.

thumbnail Fig. 7

Wavelet coherence between the SN-14C and SN-10Be series. The color code gives the value of the coherence from 0 (blue) to 1 (red). The arrows denote the relative phase between the series so that the right-pointing arrows correspond to an exact in-phase and the left-pointing arrows to an exact anti-phase relation. Black contours delimit the areas of high coherence (95% confidence level). White curves delimit the cone of influence where results can be influenced by the edges of the time series (beyond which the analysis is possibly biased).

thumbnail Fig. 8

Probability density function (pdf) of the time of occurrence of grand minima (panel A)) and grand maxima (panel B)) relative to the time of occurrence of the nearest low and high, respectively, of the Hallstatt cycle, using the superposed epoch analysis. The times of occurrence of highs and lows of the Hallstatt cycle are defined by considering the average of the two curves (second SSA components of the SN-14C and SN-10Be series) shown in Fig. 3C (leading to 5530 BC, 3650 BC, 810 BC, and 1510 AD for the lows and to 6670 BC, 4600 BC, 1970 BC, and 160 AD for the highs).

5. New constraints on the temporal distribution of grand minimum and grand maximum events

Using the reconstructed SN-14C-C and SN-10Be-C time series shown in Fig. 5, we provide a list of grand minima and maxima (Tables 1 and 2, respectively). In contrast to earlier work (Usoskin et al. 2007; Steinhilber et al. 2008), we propose a conservative list based on both the SN-14C-C and SN-10Be-C reconstructions, meaning that we only list events that are simultaneously seen in both reconstructions (a time adjustment of the 10Be-based series for ±40 years was allowed owing to the dating uncertainties (Muscheler et al. 2014)). To identify grand minima, the following criterion was used (with one exception, see below): the event in both reconstructions (using the mean of the ensemble) must correspond to a SN value below a threshold value of SN = 20 for at least 30 years. Although the event ca. 6385 BC lies slightly above this threshold, we also considered it as a grand minimum because it clearly has a Spörer-type (i.e., prolonged) shape and occurs in both series. It is possible that the level of activity was slightly overestimated for this event because of uncertainties in dipole moment evolution during the older part of the time interval. To identify grand maxima, we similarly requested the events to have a SN value exceeding the threshold of SN = 55 for at least 30 years in both reconstructions. We thus defined 20 grand minima with a total duration of 1460 years (≈17% of time) and 14 grand maxima with a total duration of 750 years (≈8%). These numbers are similar to those estimated earlier by Usoskin et al. (2007), although we note that more grand maxima are now identified. In contrast, these numbers are significantly lower than those recently estimated by Inceoglu et al. (2015), who relied on a different type of analysis and set of criteria (somewhat less restrictive, leading to the identification of a larger set of events, which is essentially inclusive of the set of events we identified here, see Tables 1 and 2).

Times of grand minima and maxima identified in this way are shown in Figs. 3C and 5 as circles and stars, respectively. An important observation from this figure is that the occurrence of grand minima and maxima appears intermittent in time. Grand minima seem to be closely related to the Hallstatt cycle, being clearly more numerous during the lows of the Hallstatt cycle (see Fig. 3C). This feature was only hinted at in passing by Steinhilber et al. (2012). Figure 8 shows the pdf (built using the superposed epoch analysis) of the grand minima and maxima times of occurrence relative to the time of occurrence of the nearest Hallstatt cycle low and high. A tendency to cluster is clearly observed. We checked that a similar tendency also appears when considering the lists of grand maxima and minima provided by Usoskin et al. (2007) and Inceoglu et al. (2015) over the same time period. We speculate that this clustering might mean that the probability of a switch of the solar dynamo from the normal mode to the grand minimum mode (resp. grand maximum mode), according to Usoskin et al. (2014), is modulated by the Hallstatt cycle. Although this interpretation relies on (necessarily) arbitrary choices for defining grand minimum and maximum events, it clearly deserves more investigation to better constrain the behavior of the solar dynamo on centennial and millennial timescales.

6. Conclusions

Here we presented new reconstructions of solar activity (quantified in terms of sunspot numbers) spanning the past 9000 yr (tables are available at the CDS) and assessed their accuracy using different geomagnetic field reconstructions and current cosmogenic isotope production models.

We found that the primary SSA components of the reconstructions that are based on 14C and 10Be are significantly different. These primary components probably reflect long-term changes in the terrestrial system that affect the 14C and 10Be isotopes in different ways (ocean circulation and/or large-scale atmospheric transport). These components were therefore removed to produce meaningful corrected sunspot number series. In contrast, the secondary SSA components of the 14C and 10Be based reconstructions revealed a common remarkably synchronous ≈2400-year quasi-periodicity. We therefore concluded that this so-called Hallstatt periodicity most likely reflects some periodicity in the solar activity.

From the two cosmogenic isotope records, we finally defined a conservative list of grand minima and grand maxima covering the past 9 millennia. An important finding is that the grand minima and maxima occurred intermittently over the studied period, with clustering near minima and maxima of the Hallstatt cycle, respectively. The Hallstatt cycle thus appears to be a long-term feature of solar activity that needs to be taken into account in models of solar dynamo.

Acknowledgments

We are thankful to Raphael Roth and Fortunate Joos for providing data on 14C production rates and for useful discussions on the carbon cycle. Support by the Academy of Finland to the ReSoLVE Center of Excellence (project no. 272157) and by IPGP through its invitation program is acknowledged. Y.G. and G.K. were partly financed by grant N 14.Z50.31.0017 of the Russian Ministry of Science and Education. This is IPGP contribution No. 3699.

References

  1. Bard, E., Raisbeck, G., Yiou, F., & Jouzel, J. 1997, Earth Planet. Sci. Lett., 150, 453 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bazilevskaya, G. A., Cliver, E. W., Kovaltsov, G. A., et al. 2014, Space Sci. Rev., 186, 409 [NASA ADS] [CrossRef] [Google Scholar]
  3. Beer, J., McCracken, K., & von Steiger, R. 2012, Cosmogenic Radionuclides: Theory and Applications in the Terrestrial and Space Environments (Berlin: Springer) [Google Scholar]
  4. Brown, M. C., Donadini, F., Korte, M., et al. 2015, Earth Plan. Space, 67, 83 [NASA ADS] [CrossRef] [Google Scholar]
  5. Cai, S., Tauxe, L., Deng, C., et al. 2014, Earth Planet. Sci. Lett., 392, 217 [NASA ADS] [CrossRef] [Google Scholar]
  6. Cai, S., Chen, W., Tauxe, L., et al. 2015, J. Geophys. Res., 120, 2056 [NASA ADS] [CrossRef] [Google Scholar]
  7. Cameron, R., & Schüssler, M. 2015, Science, 347, 1333 [NASA ADS] [CrossRef] [Google Scholar]
  8. Clette, F., Svalgaard, L., Vaquero, J., & Cliver, E. 2014, Space Sci. Rev., 186, 35 [NASA ADS] [CrossRef] [Google Scholar]
  9. Cromwell, G., Tauxe, L., Staudigel, H., & Ron, H. 2015, Phys. Earth Planet. Inter., 241, 44 [NASA ADS] [CrossRef] [Google Scholar]
  10. Damon, P., & Sonett, C. 1991, in The Sun in Time, eds. C. Sonett, M. Giampapa, & M. Matthews (Tucson, USA: University of Arizona Press), 360 [Google Scholar]
  11. de Groot, L. V., Béguin, A., Kosters, M. E., et al. 2015, Earth Planet. Sci. Lett., 419, 154 [NASA ADS] [CrossRef] [Google Scholar]
  12. Di Chiara, A., Tauxe, L., & Speranza, F. 2014, Phys. Earth Planet. Inter., 234, 1 [NASA ADS] [CrossRef] [Google Scholar]
  13. Donadini, F., Korhonen, K., Riisager, P., & Pesonen, L. J. 2006, EOS Transactions, 87, 137 [NASA ADS] [CrossRef] [Google Scholar]
  14. Ebisuzaki, W. 1997, J. Climate, 10, 2147 [Google Scholar]
  15. Field, C., Schmidt, G., Koch, D., & Salyk, C. 2006, J. Geophys. Res., 111, D15107 [NASA ADS] [CrossRef] [Google Scholar]
  16. Gallet, Y., & Al Maqdissi, M. 2010, Akkadica, 131, 29 [Google Scholar]
  17. Gallet, Y., & Butterlin, P. 2015, Archaeometry, 57, S1, 263 [CrossRef] [Google Scholar]
  18. Gallet, Y., Le Goff, M., Genevey, A., Margueron, J., & Matthiae, P. 2008, Geophys. Res. Lett., 35, 2307 [NASA ADS] [CrossRef] [Google Scholar]
  19. Gallet, Y., Genevey, A., Le Goff, M., et al. 2009, Comptes Rendus Physique, 10, 630 [NASA ADS] [CrossRef] [Google Scholar]
  20. Gallet, Y., MolistMontaña, M., Genevey, A., et al. 2015, Phys. Earth Planet. Inter., 238, 89 [NASA ADS] [CrossRef] [Google Scholar]
  21. Genevey, A., Gallet, Y., Constable, C. G., Korte, M., & Hulot, G. 2008, Geochem., Geophys., Geosyst., 9, 04038 [Google Scholar]
  22. Genevey, A., Gallet, Y., Thébault, E., Jesset, S., & Le Goff, M. 2013, Geochem., Geophys., Geosys., 14, 2858 [CrossRef] [Google Scholar]
  23. Golub, G., & Kahan, W. 1965, J. Soc. Indust. Appl. Math. B: Num. Anal., 2, 205 [CrossRef] [Google Scholar]
  24. Golyandina, N., Nekrutkin, V., & Zhigljavsky, A. 2001, Analysis of Time Series Structure: SSA and related techniques (Boca Raton, London, New York, Washington DC: Chapman and Hall/CRC) [Google Scholar]
  25. Grinsted, A., Moore, J. C., & Jevrejeva, S. 2004, Nonlin. Processes Geophys., 11, 561 [NASA ADS] [CrossRef] [Google Scholar]
  26. Hanslmeier, A., Brajša, R., Čalogović, J., et al. 2013, A&A, 550, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Hassani, H. 2007, Journal of Data Science, 5, 239 [Google Scholar]
  28. Heikkilä, U., Beer, J., & Feichter, J. 2009, Atmos. Chem. Phys., 9, 515 [Google Scholar]
  29. Hong, H., Yu, Y., Lee, C. H., et al. 2013, Earth Planet. Sci. Lett., 383, 142 [NASA ADS] [CrossRef] [Google Scholar]
  30. Hoyt, D. V., & Schatten, K. H. 1998, Sol. Phys., 179, 189 [NASA ADS] [CrossRef] [Google Scholar]
  31. Hua, Q., Webb, G. E., Zhao, J., et al. 2015, Earth Planet. Sci. Lett., 422, 33 [NASA ADS] [CrossRef] [Google Scholar]
  32. Inceoglu, F., Simoniello, R., Knudsen, V. F., et al. 2015, A&A, 577, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Kapper, K. L., Donadini, F., & Hirt, A. M. 2015, Phys. Earth Planet. Inter., 241, 21 [NASA ADS] [CrossRef] [Google Scholar]
  34. Karak, B. B., Kitchatinov, L. L., & Brandenburg, A. 2015, ApJ, 803, 95 [NASA ADS] [CrossRef] [Google Scholar]
  35. Kissel, C., Laj, C., Rodriguez-Gonzalez, A., et al. 2015, Earth Planet. Sci. Lett., 430, 178 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kittler, J., & Young, P. C. 1973, Pattern Recognition, 5, 335 [CrossRef] [Google Scholar]
  37. Knudsen, M. F., Riisager, P., Donadini, F., et al. 2008, Earth Planet. Sci. Lett., 272, 319 [NASA ADS] [CrossRef] [Google Scholar]
  38. Korhonen, K., Donadini, F., Riisager, P., & Pesonen, L. J. 2008, Geochem., Geophys., Geosys., 9, Q04029 [Google Scholar]
  39. Korte, M., & Constable, C. 2005, Earth Planet. Sci. Lett., 236, 348 [CrossRef] [Google Scholar]
  40. Korte, M., & Constable, C. 2011, Phys. Earth Planet. Inter., 188, 247 [NASA ADS] [CrossRef] [Google Scholar]
  41. Korte, M., Donadini, F., & Constable, C. G. 2009, Geochem. Geophys. Geosys., 10, Q06008 [Google Scholar]
  42. Korte, M., Constable, C., Donadini, F., & Holme, R. 2011, Earth Planet. Sci. Lett., 312, 497 [NASA ADS] [CrossRef] [Google Scholar]
  43. Kovaltsov, G., & Usoskin, I. 2010, Earth Planet. Sci. Lett., 291, 182 [NASA ADS] [CrossRef] [Google Scholar]
  44. Kovaltsov, G., Mishev, A., & Usoskin, I. 2012, Earth Planet. Sci. Lett., 337, 114 [Google Scholar]
  45. Krivova, N., Balmaceda, L., & Solanki, S. 2007, A&A, 467, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Licht, A., Hulot, G., Gallet, Y., & Thébault, E. 2013, Phys. Earth Planet. Inter., 224, 38 [NASA ADS] [CrossRef] [Google Scholar]
  47. Lockwood, M., Owens, M. J., & Barnard, L. 2014, J. Geophys. Res., Space Phys., 119, 5172 [Google Scholar]
  48. Ma, L. H. 2007, Sol. Phys., 245, 411 [NASA ADS] [CrossRef] [Google Scholar]
  49. Mane, R. 1981, in Dynamical Systems and Turbulence, eds. D. Rand, & L. Young (Springer-Verlag), Lect. Notes Math., 898, 230 [Google Scholar]
  50. McCracken, K., Beer, J., Steinhilber, F., & Abreu, J. 2011, Space Sci. Rev., 324 [Google Scholar]
  51. Miyake, F., Nagaya, K., Masuda, K., & Nakamura, T. 2012, Nature, 486, 240 [NASA ADS] [CrossRef] [Google Scholar]
  52. Miyake, F., Masuda, K., & Nakamura, T. 2013, Nat. Comm., 4, 1748 [Google Scholar]
  53. Moss, D., Sokoloff, D., Usoskin, I., & Tutubalin, V. 2008, Sol. Phys., 250, 221 [NASA ADS] [CrossRef] [Google Scholar]
  54. Muscheler, R., Joos, F., Beer, J., et al. 2007, Quat. Sci. Rev., 26, 82 [Google Scholar]
  55. Muscheler, R., Adolphi, F., & Knudsen, M. F. 2014, Quat. Sci. Rev., 106, 81 [NASA ADS] [CrossRef] [Google Scholar]
  56. Nachasova, I., & Burakov, K. 1995, Geomagn. Aeron., 35 [Google Scholar]
  57. Nachasova, I., & Burakov, K. 1998, Geomagn. Aeron., 38 [Google Scholar]
  58. Nilsson, A., Holme, R., Korte, M., Suttie, N., & Hill, M. 2014, Geophys. J. Int., 198, 229 [NASA ADS] [CrossRef] [Google Scholar]
  59. Osete, M. L., Catanzariti, G., Chauvin, A., et al. 2015, Phys. Earth Planet. Inter., 242, 24 [NASA ADS] [CrossRef] [Google Scholar]
  60. Passos, D., Nandy, D., Hazra, S., & Lopes, I. 2014, A&A, 563, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. 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]
  62. Reimer, P. J., Baillie, M. G. L., Bard, E., et al. 2009, Radiocarbon, 51, 1111 [Google Scholar]
  63. Roth, R., & Joos, F. 2013, Clim. Past, 9, 1879 [NASA ADS] [CrossRef] [Google Scholar]
  64. Shaar, R., Tauxe, L., Ben-Yosef, E., et al. 2015, Geochem., Geophys., Geosyst., 16, 195 [CrossRef] [Google Scholar]
  65. Snowball, I., & Muscheler, R. 2007, Holocene, 17, 851 [NASA ADS] [CrossRef] [Google Scholar]
  66. Solanki, S., Schüssler, M., & Fligge, M. 2000, Nature, 408, 445 [Google Scholar]
  67. Solanki, S., Usoskin, I., Kromer, B., Schüssler, M., & Beer, J. 2004, Nature, 431, 1084 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  68. Steinhilber, F., Abreu, J. A., & Beer, J. 2008, Astrophys. Space Sci. Trans., 4, 1 [NASA ADS] [CrossRef] [Google Scholar]
  69. Steinhilber, F., Abreu, J., Beer, J., et al. 2012, Proc. Nat. Acad. Sci. USA, 109, 5967 [Google Scholar]
  70. Stillinger, M. D., Feinberg, J. M., & Frahm, E. 2015, J. Archaeol. Sci., 53, 345 [CrossRef] [Google Scholar]
  71. Svalgaard, L., & Schatten, K. H. 2015, ArXiv e-prints [arXiv:1506.00755] [Google Scholar]
  72. Takens, F. 1981, in Dynamical Systems and Turbulence, eds. D. Rand, & L. Young (Springer-Verlag), Lect. Notes Math., 898, 366 [Google Scholar]
  73. Thébault, E., & Gallet, Y. 2010, Geophys. Res. Lett., 37, L22303 [NASA ADS] [Google Scholar]
  74. Usoskin, I. G. 2013, Liv. Rev. Sol. Phys., 10, 1 [Google Scholar]
  75. Usoskin, I. G., Alanko-Huotari, K., Kovaltsov, G. A., & Mursula, K. 2005, J. Geophys. Res., 110, A12108 [NASA ADS] [CrossRef] [Google Scholar]
  76. Usoskin, I. G., Voiculescu, M., Kovaltsov, G. A., & Mursula, K. 2006, J. Atmosph. Sol.-Terr. Phys., 68, 2164 [NASA ADS] [CrossRef] [Google Scholar]
  77. Usoskin, I. G., Solanki, S. K., & Kovaltsov, G. A. 2007, A&A, 471, 301 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Usoskin, I. G., Horiuchi, K., Solanki, S., Kovaltsov, G. A., & Bard, E. 2009, J. Geophys. Res., 114, A03112 [NASA ADS] [Google Scholar]
  79. Usoskin, I. G., Kromer, B., Ludlow, F., et al. 2013, A&A, 552, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Usoskin, I. G., Hulot, G., Gallet, Y., et al. 2014, A&A, 562, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Vasiliev, S., & Dergachev, V. 2002, Annales Geophys., 20, 115 [Google Scholar]
  82. Vasiliev, S. S., Dergachev, V. A., Raspopov, O. M., & Jungner, H. 2012, Geomagn. Aeron., 52, 121 [NASA ADS] [CrossRef] [Google Scholar]
  83. Vautard, R., & Ghil, M. 1989, Physica D Nonlin. Phenom., 35, 395 [Google Scholar]
  84. Vautard, R., Yiou, P., & Ghil, M. 1992, Physica D Nonlin. Phen., 58, 95 [Google Scholar]
  85. Vonmoos, M., Beer, J., & Muscheler, R. 2006, J. Geophys. Res., 111, A10105 [NASA ADS] [CrossRef] [Google Scholar]
  86. Webber, W. R., & Higbie, P. R. 2009, J. Geophys. Res., 114, A02103 [NASA ADS] [CrossRef] [Google Scholar]
  87. Yiou, F., Raisbeck, G., Baumgartner, S., et al. 1997, J. Geophys. Res., 102, 26783 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: GMAG.9k axial dipole evolution

Here we used the GEOMAGIA50.v3 database (Donadini et al. 2006; Korhonen et al. 2008; Brown et al. 2015) to which we added recent archeo- and paleointensity results (Cai et al. 2014, 2015; Cromwell et al. 2015; de Groot et al. 2015; Di Chiara et al. 2014; Gallet et al. 2008, 2009, 2015; Gallet & Al Maqdissi 2010; Hong et al. 2013; Kapper et al. 2015; Kissel et al. 2015; Osete et al. 2015; Shaar et al. 2015; Stillinger et al. 2015). Concerning the data compiled in this version of GEOMAGIA, the Mesopotamian data from Nachasova & Burakov (1995, 1998) were modified according to Gallet et al. (2015). Further revisions were also discussed in Genevey et al. (2013; also A. Genevey, pers. comm.) and Gallet & Butterlin (2015).

Axial dipole evolution over the past 9000 years, referred to as GMAG.9k, was constrained by virtual axial dipole moments (VADM) averaged over sliding windows of 200 years between 1500 BC and 2000 AD and of 500 years between 7000 BC and 1500 BC shifted by 10 years and using a regional weighting scheme over regions of 10° width. Weights of one third or two thirds were assigned to the regions that contained at a given time interval only one or two individual intensity (VADM) data points, respectively. For those intensity data with no age uncertainties provided in the GEOMAGIA50.v3 database, we used the same approach as in Licht et al. (2013). We computed the means of the known age uncertainties over 500 yr long time intervals between 1000 BC and 2000 AD, over 2000 yr interval between 3000 BC and 1000 BC and over 4000 yr interval between 7000 BC and 3000 BC. After multiplication by a factor of 1.5, the corresponding values were assigned to the intensity data with unknown age uncertainties within the periods of concern. Following Knudsen et al. (2008), when no experimental errors were provided on the data, we assigned errors amounting 25% of the corresponding intensity values. We also note that for the archeomagnetic data for which no attempt was made to take the cooling rate effect on thermoremanent magnetization acquisition into account, we systematically implemented a cooling rate correction of 5% decrease (see for instance in Genevey et al. 2008). Mean VADM estimates were then derived using a bootstrap technique to account for the noise in the available paleo- and archeointensity data within their age uncertainties and within their 2σ experimental uncertainties. 1000 VADM curves, also derived using different randomly attributed locations of the weighting regions, were hence determined, whose statistics are summarized in Tables C.1 and C.2 (available in CDS) for the periods 1500 BC-2000 AD and 7000 BC-1500 BC, respectively. These tables provide for each epoch (first column) a mean VADM (second column), a standard deviation (third column), and the maximum and minimum values defining the envelope of possible VADM results (fourth and fifth columns). The variability of VADM is shown in Figs. A.1 (between 1500 BC and 2000 AD) and A.2 (between 7000 BC and 2000 AD).

thumbnail Fig. A.1

Comparison between VADM curves computed between 1500 BC and 2000 AD using sliding windows of ΔT = 200 years (panels A), C)) and 500 years (B), D)) shifted by 50 years, and using a weighting over regions of ΔS = 10° × 10° (A), B)) and ΔS = 30° × 30° (C), D) width. The thick black line exhibits the averaged VADM computed using a bootstrap scheme (see main text and legend of Tables C.1 and C.2 available at CDS), with its 1σ uncertainties (dotted lines) and the envelope of possible VADM values (gray lines).

thumbnail Fig. A.2

Comparison between VADM curves computed between 7000 BC and 2000 AD using sliding windows of ΔT = 500 years (panels A), C)) and 1000 years (B), D)) shifted by 50 years, and using a weighting over regions of ΔS = 10° × 10° (A), B)) and ΔS = 30° × 30° (C), D) width. The thick black line exhibits the averaged VADM computed using a bootstrap scheme (see main text and legend of Tables C.1 and C.2 available at CDS), with its 1σ uncertainties (dotted lines) and the envelope of possible VADM values (gray lines).

Appendix B: Identification of long-term trends by singular spectrum analysis

The non-parametric singular spectrum analysis (SSA) of time series is based on the Karhunen-Loeve spectral decomposition theorem (Kittler & Young 1973) and the Mané-Takens embedded theorem (Mane 1981; Takens 1981). It allows a time series to be decomposed into several components with distinct temporal behaviors and is very convenient to identify long-term trends and quasi-periodic oscillations.

The basic version of SSA that we used consists of four straightforward steps (see, e.g., Golyandina et al. 2001; Hassani 2007): embedding, singular value decomposition, grouping, and reconstructions.

When considering a real-value time series x (x1, x2, ..., xN), the first step of this SSA consists of embedding this series into an L-dimensional vector space, using lagged copies of x to form the so-called trajectory (Hankel) matrix (where K = N−L + 1 and L is a parameter to be chosen), (B.1)The second step consists of performing a singular value decomposition (Golub & Kahan 1965) of the trajectory matrix. This provides a set of L eigenvalues λi (arranged in decreasing order λ1 ≥ λ2 ≥ ··· ≥ λL ≥ 0) and eigenvectors Ui (often called

“empirical orthogonal functions”) of the matrix D = XXT. If we denote with d the number of nonzero eigenvalues, we may next define (i = 1,...,d). Then, the trajectory matrix can be written as a sum of elementary matrices X = X1 + ··· + Xd, where .

Once this decomposition has been completed, the third step consists of the construction of groups of components by rearranging X into X = XG1 + XG2 + ..., where each XG is the sum (group) of a number of Xi. The choice of the components to be considered in each group is made empirically by grouping eigentriples () with similar eigenvalues. Finally, a diagonal averaging is applied to each XG to make it take the form of a trajectory matrix, from which the associated time series component of length N can be recovered (for details, see, e.g., Golyandina et al. 2001; Hassani 2007).

All Tables

Table 1

Grand minima with their centers, approximate duration, and comments (1 – listed in Usoskin et al. 2007; 2 – listed in Inceoglu et al. 2015).

In the text
Table 2

Grand maxima with their centers, approximate duration, and comments (1 – listed in Usoskin et al. 2007; 2 – listed in Inceoglu et al. 2015).

In the text

All Figures

thumbnail Fig. 1

Time series of cosmogenic radionuclide data. Panel A): decadal radiocarbon 14C global production rate (Roth & Joos 2013) with the 95% confidence interval plotted in gray. Panel B): quasi-decadal variability of 10Be flux in the GRIP ice core (Yiou et al. 1997). The formal 1σ error of 7% (relative to the given value) is indicated by the error bar next to the legend.

In the text
thumbnail Fig. 2

Time series of the axial dipole moment reconstructions spanning the past 9000 years (panel A)), with a zoom for the last 3200 years (panel B)). The black solid line depicts GMAG.9k (the reconstruction presented and used in this work) with ± 1σ and the full range variability presented by the gray shading and the hatching, respectively. Other reconstructions shown are (Licht et al. 2013) denoted as AF_M, (Genevey et al. 2008) denoted as G08, (Nilsson et al. 2014) denoted as pfm9k.1b and pfm9k.1a, (Knudsen et al. 2008) denoted as Kn08, and (Pavón-Carrasco et al. 2014) denoted as SHA-DIF.14k. For better readability, error bars were omitted for these curves, but this does not affect the discussion of the results (see text).

In the text
thumbnail Fig. 3

Panel A): raw reconstructions of the sunspot numbers (mean curves) SN-14C (blue) and SN-10Be (red), compared to the recent 3 kyr reconstruction (Usoskin et al. 2014,– green curve). Panel B): first component of the singular spectrum analysis (SSA – see Appendix B) for the SN-14C (blue) and SN-10Be (red) series. The shaded areas depict the uncertainties related to the parameter L of the SSA. Panel C): same as in panel B), but for the second SSA components of the SN-14C (blue) and SN-10Be (red) series. The large dots and red stars denote times of the grand minima (see Table 1) and grand maxima (Table 2), respectively.

In the text
thumbnail Fig. 4

Comparison of alternative SN-14 sunspot number reconstructions when relying on different axial dipole reconstructions, (same notations as in Fig. 2). Only the mean curves of the corresponding ensembles are shown.

In the text
thumbnail Fig. 5

Corrected sunspot number reconstructions SN-14C-C (panel A)) and SN-10Be-C (panel B)), after removing long-term trends (see Sect. 4.1). The black curves and the gray shading depict the mean and the 95% range (over 106 ensemble members) of the reconstructions, respectively, and the green curve represents the 3 kyr reconstruction by Usoskin et al. (2014). Stars and circles denote grand maxima and minima, respectively, as in Fig. 3. Tables for this plot are available at the CDS, including the mean values and the uncertainties of the sunspot numbers reconstructed here as shown in Fig. 5.

In the text
thumbnail Fig. 6

Decadally averaged sunspot numbers for the period 1600–1900 AD. The thick gray curve represents this work (uncertainties are not shown). Other curves correspond to sunspot number series: L14 (Lockwood et al. 2014), H98 (Hoyt & Schatten 1998), SS15 (Svalgaard & Schatten 2015), and the C14 international sunspot number (v.2) scaled with a factor 0.6 (Clette et al. 2014).

In the text
thumbnail Fig. 7

Wavelet coherence between the SN-14C and SN-10Be series. The color code gives the value of the coherence from 0 (blue) to 1 (red). The arrows denote the relative phase between the series so that the right-pointing arrows correspond to an exact in-phase and the left-pointing arrows to an exact anti-phase relation. Black contours delimit the areas of high coherence (95% confidence level). White curves delimit the cone of influence where results can be influenced by the edges of the time series (beyond which the analysis is possibly biased).
In the text
thumbnail Fig. 8

Probability density function (pdf) of the time of occurrence of grand minima (panel A)) and grand maxima (panel B)) relative to the time of occurrence of the nearest low and high, respectively, of the Hallstatt cycle, using the superposed epoch analysis. The times of occurrence of highs and lows of the Hallstatt cycle are defined by considering the average of the two curves (second SSA components of the SN-14C and SN-10Be series) shown in Fig. 3C (leading to 5530 BC, 3650 BC, 810 BC, and 1510 AD for the lows and to 6670 BC, 4600 BC, 1970 BC, and 160 AD for the highs).

In the text
thumbnail Fig. A.1

Comparison between VADM curves computed between 1500 BC and 2000 AD using sliding windows of ΔT = 200 years (panels A), C)) and 500 years (B), D)) shifted by 50 years, and using a weighting over regions of ΔS = 10° × 10° (A), B)) and ΔS = 30° × 30° (C), D) width. The thick black line exhibits the averaged VADM computed using a bootstrap scheme (see main text and legend of Tables C.1 and C.2 available at CDS), with its 1σ uncertainties (dotted lines) and the envelope of possible VADM values (gray lines).

In the text
thumbnail Fig. A.2

Comparison between VADM curves computed between 7000 BC and 2000 AD using sliding windows of ΔT = 500 years (panels A), C)) and 1000 years (B), D)) shifted by 50 years, and using a weighting over regions of ΔS = 10° × 10° (A), B)) and ΔS = 30° × 30° (C), D) width. The thick black line exhibits the averaged VADM computed using a bootstrap scheme (see main text and legend of Tables C.1 and C.2 available at CDS), with its 1σ uncertainties (dotted lines) and the envelope of possible VADM values (gray lines).

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.