A History of Solar Activity over Millennia

Presented here is a review of present knowledge of the long-term behavior of solar activity on a multi-millennial timescale, as reconstructed using the indirect proxy method. The concept of solar activity is discussed along with an overview of the special indices used to quantify different aspects of variable solar activity, with special emphasis upon sunspot number. Over long timescales, quantitative information about past solar activity can only be obtained using a method based upon indirect proxies, such as the cosmogenic isotopes \super{14}C and \super{10}Be in natural stratified archives (e.g., tree rings or ice cores). We give an historical overview of the development of the proxy-based method for past solar-activity reconstruction over millennia, as well as a description of the modern state. Special attention is paid to the verification and cross-calibration of reconstructions. It is argued that this method of cosmogenic isotopes makes a solid basis for studies of solar variability in the past on a long timescale (centuries to millennia) during the Holocene. A separate section is devoted to reconstructions of strong solar energetic-particle (SEP) events in the past, that suggest that the present-day average SEP flux is broadly consistent with estimates on longer timescales, and that the occurrence of extra-strong events is unlikely. Finally, the main features of the long-term evolution of solar magnetic activity, including the statistics of grand minima and maxima occurrence, are summarized and their possible implications, especially for solar/stellar dynamo theory, are discussed.


Introduction
The concept of the perfectness and constancy of the sun, postulated by Aristotle, was a strong belief for centuries and an official doctrine of Christian and Muslim countries. However, as people had noticed even before the time of Aristotle, some slight transient changes of the sun can be observed even with the naked eye. Although scientists knew about the existence of "imperfect" spots on the sun since the early 17th century, it was only in the 19th century that the scientific community recognized that solar activity varies in the course of an 11-year solar cycle. Solar variability was later found to have many different manifestations, including the fact that the "solar constant" (the amount of total incoming solar electromagnetic radiation in all wavelengths per unit area) is not a constant. The sun appears much more complicated and active than a static hot plasma ball, with a great variety of nonstationary active processes going beyond the adiabatic equilibrium foreseen in the basic theory of sun-as-star. Such transient nonstationary (often eruptive) processes can be broadly regarded as solar activity, in contrast to the so-called "quiet" sun. Solar activity includes active transient and long-lived phenomena on the solar surface, such as spectacular solar flares, sunspots, prominences, coronal mass ejections (CMEs), etc.
The very fact of the existence of solar activity poses an enigma for solar physics, leading to the development of sophisticated models of an upper layer known as the convection zone and the solar corona. The sun is the only star, which can be studied in great detail and thus can be considered as a proxy for cool stars. Quite a number of dedicated ground-based and space-borne experiments are being carried out to learn more about solar variability. The use of the sun as a paradigm for cool stars, leads to a better understanding of the processes driving the broader population of cool sun-like stars. Therefore, studying and modelling solar activity can increase the level of our understanding of nature.
On the other hand, the study of variable solar activity is not of purely academic interest, as it directly affects the terrestrial environment. Although changes in the sun are barely visible without the aid of precise scientific instruments, these changes have great impact on many aspects of our lives. In particular, the heliosphere (a spatial region of about 100 astronomical units) is mainly controlled by the solar magnetic field. This leads to the modulation of galactic cosmic rays (GCRs). Additionally, eruptive and transient phenomena in the sun/corona and in the interplanetary medium can lead to the acceleration of energetic particles with greatly enhanced flux. Such processes can modify the radiation environment on Earth and need to be taken into account for planning and maintaining space missions and even transpolar jet flights. Solar activity can cause, through coupling of solar wind and the Earth's magnetosphere, strong geomagnetic storms in the magnetosphere and ionosphere, which may disturb radio-wave propagation and navigation-system stability, or induce dangerous spurious currents in long pipes or power lines. Another important aspect is the link between solar-activity variations and the Earth's climate (see, e.g., the review by Haigh, 2007).
It is important to study solar variability on different timescales. The primary basis for such studies is observational (or reconstructed) data. The sun's activity is systematically explored in different ways (solar, heliospheric, interplanetary, magnetospheric, terrestrial), including groundbased and space-borne experiments and dedicated missions during the last few decades, thus covering 3 -4 solar cycles. However, it should be noted that the modern epoch is characterized by unusually-high solar activity dominated by an 11-year cyclicity, and it is not straightforward to extrapolate present knowledge (especially empirical and semi-empirical relationships and models) to a longer timescale. Therefore, the behavior of solar activity in the past, before the era of direct measurements, is of great importance for a variety of reasons. For example, it allows an improved knowledge of the statistical behavior of the solar-dynamo process, which generates the cyclicallyvarying solar-magnetic field, making it possible to estimate the fractions of time the sun spends in states of very-low activity, what are called grand minima. Such studies require a long time series of solar-activity data. The longest direct series of solar activity is the 400-year-long sunspot-number series, which depicts the dramatic contrast between the (almost spotless) Maunder minimum and the modern period of very high activity. Thanks to the recent development of precise technologies, including accelerator mass spectrometry, solar activity can be reconstructed over multiple millennia from concentrations of cosmogenic isotopes 14 C and 10 Be in terrestrial archives. This allows one to study the temporal evolution of solar magnetic activity, and thus of the solar dynamo, on much longer timescales than are available from direct measurements.
This paper gives an overview of the present status of our knowledge of long-term solar activity, covering the Holocene (the last 11 millennia). A description of the concept of solar activity and a discussion of observational methods and indices are presented in Section 2. The proxy method of solar-activity reconstruction is described in some detail in Section 3. Section 4 gives an overview of what is known about past solar activity. The long-term averaged flux of solar energetic particles is discussed in Section 5. Finally, conclusions are summarized in Section 6.
2 Solar Activity: Concept and Observations

The concept of solar activity
The sun is known to be far from a static state, the so-called "quiet" sun described by simple stellar-evolution theories, but instead goes through various nonstationary active processes. Such nonstationary and nonequilibrium (often eruptive) processes can be broadly regarded as solar activity. Whereas the concept of solar activity is quite a common term nowadays, it is not straightforwardly interpreted nor unambiguously defined. For instance, solar-surface magnetic variability, eruption phenomena, coronal activity, radiation of the sun as a star or even interplanetary transients and geomagnetic disturbances can be related to the concept of solar activity. A variety of indices quantifying solar activity have been proposed in order to represent different observables and caused effects. Most of the indices are highly correlated to each other due to the dominant 11-year cycle, but may differ in fine details and/or long-term trends. The indices can be roughly divided into physical indices (i.e., those representing real physical quantities measurable for the sun) and synthetic indices. In addition to the solar indices, indirect proxy data is often used to quantify solar activity via its presumably known effect on the magnetosphere or heliosphere. The indices of solar activity that are often used for long-term studies are reviewed below.

Indices of solar activity
Solar (as well as other) indices can be divided into physical and synthetic according to the way they are obtained/calculated. Physical indices quantify the directly-measurable values of a real physical observable, such as the radioflux, and thus have clear physical meaning. Physical indices quantify physical features of different aspects of solar activity and their effects. Synthetic indices (the most common being sunspot number) are calculated (or synthesized) using a special algorithm from observed (often not measurable in physical units) data or phenomena. Additionally, solar activity indices can be either direct (i.e., directly relating to the sun) or indirect (relating to indirect effects caused by solar activity), as discussed in subsequent subsections 2.2.1, 2.2.2.

Direct solar indices
The most commonly used index of solar activity is based on sunspot number. Sunspots are dark areas on the solar disc (of size up to tens of thousands of km, lifetime up to half-a-year), characterized by a strong magnetic field, which leads to a lower temperature (about 4000 K compared to 5800 K in the photosphere) and observed darkening.
Sunspot number is a synthetic, rather than a physical, index, but it has still become quite a useful parameter in quantifying the level of solar activity. This index presents the weighted number of individual sunspots and/or sunspot groups, calculated in a prescribed manner from simple visual solar observations. The use of the sunspot number makes it possible to combine together thousands and thousands of regular and fragmentary solar observations made by earlier professional and amateur astronomers. The technique, initially developed by Rudolf Wolf, yielded the longest series of directly and regularly-observed scientific quantities. Therefore, it is common to quantify solar magnetic activity via sunspot numbers. For details see the review on sunspot numbers and solar cycles (Hathaway and Wilson, 2004;Hathaway, 2008).

Wolf sunspot number (WSN) series
The concept of the sunspot number was developed by Rudolf Wolf of the Zürich observatory in the middle of the 19th century. The sunspot series, initiated by him, is called the Zürich or Wolf sunspot number (WSN) series. The relative sunspot number R z is defined as where G is the number of sunspot groups, N is the number of individual sunspots in all groups visible on the solar disc and k denotes the individual correction factor, which compensates for differences in observational techniques and instruments used by different observers, and is used to normalize different observations to each other. The value of R z (see Figure 1a) is calculated for each day using only one observation made by the "primary" observer (judged as the most reliable observer during a given time) for the day. The primary observers were Staudacher (1749 -1787), Flaugergues (1788 -1825), Schwabe (1826 -1847), Wolf (1848Wolf ( -1893, Wolfer (1893-1928), Brunner (1929-1944, Waldmeier (1945Waldmeier ( -1980 and Koeckelenbergh (since 1980). If observations by the primary observer are not available for a certain day, the secondary, tertiary, etc. observers are used (see the hierarchy of observers in Waldmeier, 1961). The use of only one observer for each day aims to make R z a homogeneous time series. As a drawback, such an approach ignores all other observations available for the day, which constitute a large fraction of the existing information. Moreover, possible errors of the primary observer cannot be caught or estimated. If no sunspot observations are available for some period, the data gap is filled, without note in the final WSN series, using an interpolation between the available data and by employing some proxy data. The observational uncertainties in the monthly R z can be up to 25% (e.g., Vitinsky et al., 1986). The WSN series is based on observations performed at the Zürich Observatory during 1849 -1981 using almost the same technique. This part of the series is fairly stable and homogeneous. However, prior to that there have been many gaps in the data that were interpolated. Therefore, the WSN series is a combination of direct observations and interpolations for the period before 1849, leading to possible errors and inhomogeneities as discussed, e.g, by Vitinsky et al. (1986); Wilson (1998); Letfus (1999). The quality of the Wolf series before 1749 is rather poor and hardly reliable (Hoyt et al., 1994;Hoyt and Schatten, 1998;Hathaway and Wilson, 2004).
Note that the sun has been routinely photographed since 1876 so that full information on daily sunspot activity is available (the Greenwich series) with observational uncertainties being negligible for the last 140 years. The routine production of the WSN series was terminated in Zürich in 1982. Since then, the sunspot number series is routinely updated as the International sunspot number R i , provided by the Solar Influences Data Analysis Center in Belgium (Clette et al., 2007). The international sunspot number series is computed using the same definition (Equation 1) as WSN but it has a significant distinction from the WSN; it is based not on a single primary solar observation for each day but instead uses a weighted average of more than 20 approved observers.
In addition to the standard sunspot number R i , there is also a series of hemispheric sunspot numbers R N and R S , which account for spots only in the northern and southern solar hemispheres, respectively (note that R i = R N + R S ). These series are used to study the N-S asymmetry of solar activity (Temmer et al., 2002).

Group sunspot number (GSN) series
Since the WSN series is of lower quality before the 1850's and is hardly reliable before 1750, there was a need to re-evaluate early sunspot data. This tremendous work has been done by Hoyt and Schatten (1998), who performed an extensive archive search and nearly doubled the amount of original information compared to the Wolf series. They have produced a new series of sunspot activity called the group sunspot numbers (GSN -see Figure 1b), including all available archival records. The daily group sunspot number R g is defined as follows: where G i is the number of sunspot groups recorded by the i-th observer, k is the observer's individual correction factor, n is the number of observers for the particular day, and 12.08 is a normalization number scaling R g to R z values for the period of 1874 -1976. R g is more robust than R z since it is based on more easily identified sunspot groups and does not include the number of individual spots. The GSN series includes not only one "primary" observation, but all available observations, and covers the period since 1610, being, thus, 140 years longer than the original WSN series. It is particularly interesting that the period of the Maunder minimum (1645 -1715) was surprisingly well covered with daily observations allowing for a detailed analysis of sunspot activity during this grand minimum (see also Section 4.3). Systematic uncertainties of the R g values are estimated to be about 10% before 1640, less than 5% from 1640 -1728 and from 1800 -1849, 15 -20% from 1728 -1799, and about 1% since 1849 (Hoyt and Schatten, 1998). The GSN series is more reliable and homogeneous than the WSN series before 1849. The two series are nearly identical after the 1870's (Hoyt and Schatten, 1998;Letfus, 1999;Hathaway and Wilson, 2004). However, the GSN series still contains some lacunas, uncertainties and possible inhomogeneities (see, e.g., Letfus, 2000;Usoskin et al., 2003a). The search for other lost or missing records of past solar instrumental observations has not ended even since the extensive work by Hoyt and Schatten. Archival searches still give new interesting findings of forgotten sunspot observations, often outside major observatories (e.g., Casas et al., 2006;Vaquero et al., 2005Vaquero et al., , 2007.

Other indices
An example of a synthetic index of solar activity is the flare index, representing solar flare activity (e.g., Özgüç et al., 2003;Kleczek, 1952). The flare index quantifies daily flare activity in the following manner; it is computed as a product of the flare's relative importance I in the H α -range and duration t, Q = I t, thus being a rough measure of the total energy emitted by the flare. The daily flare index is produced by Bogazici University (Özgüç et al., 2003) and is available since 1936.
A traditional physical index of solar activity is related to the radioflux of the sun in the wavelength range of 10.7 cm and is called the F10.7 index (e.g., Tapping and Charrois, 1994). This index represents the flux (in solar flux units, 1sfu = 10 −22 W m −2 Hz −1 ) of solar radio emission at a centimetric wavelength. This emission is close to the peak of solar radio emission and is produced as a result of the nonradiative heating of coronal plasma over active regions. It is a good quantitative measure of the level of solar activity, which is not directly related to sunspots. Close correlation between the F10.7 index and sunspot number indicates that the latter is a good index of general solar activity, including coronal activity. The solar F10.7 cm record has been measured continuously since 1947.
Another physical index is the coronal index (e.g., Rybanský et al., 2005), which is a measure of the irradiance of the sun as a star in the coronal green line. Computation of the coronal index is based on observations of green corona intensities (Fe XIV emission line at 530.3 nm wavelength) from coronal stations all over the world, the data being transformed to the Lomnický Stit photometric scale. This index is considered a basic optical index of solar activity. A synthesized homogeneous database of the Fe XIV 530.3 nm coronal-emission line intensities has existed since 1943 and covers six solar cycles.
Often sunspot area is considered as a physical index representing solar activity (e.g., Baranyi et al., 2001;Balmaceda et al., 2005). This index gives the total area of visible spots on the solar disc in units of millionths of the sun's visible hemisphere, corrected for apparent distortion due to the curvature of the solar surface. The area of individual groups may vary between tens of millionths (for small groups) up to several thousands of millionths for huge groups. This index has a physical meaning related to the solar magnetic flux emerging at sunspots. Sunspot areas are available since 1874 in the Greenwich series obtained from daily photographic images of the sun. In addition, some fragmentary data of sunspot areas, obtained from solar drawings, are available for earlier periods (Vaquero et al., 2004;Arlt, 2008).
An important quantity is solar irradiance, total and spectral. Irradiance variations are physically related to solar magnetic variability (e.g., Solanki et al., 2000), and are often considered manifestations of solar activity, which is of primary importance for solar-terrestrial relations. For details see the Living Review on solar irradiance (Krivova and Schmutz).
All the above indices are closely correlated to sunspot numbers on the solar-cycle scale, but may depict quite different behavior on short or long timescales.

Indirect indices
Sometimes quantitative measures of solar-variability effects are also considered as indices of solar activity. These are related not to solar activity per se, but rather to its effect on different environments. Accordingly, such indices are called indirect, and can be roughly divided into terrestrial/geomagnetic and heliospheric/interplanetary. Geomagnetic indices quantify different effects of geomagnetic activity ultimately caused by solar variability, mostly by variations of solar-wind properties and the interplanetary magnetic field. For example, the aa-index, which provides a global index of magnetic activity relative to a quiet-day curve for a pair of antipodal magnetic observatories (in England and Australia), is available from 1868 (Mayaud, 1972). An extension of the geomagnetic series is available from the 1840s using the Helsinki Ak(H) index (Nevanlinna, 2004a,b). A review of the geomagnetic effects of solar activity can be found, e.g., in Pulkkinen (2007). It is noteworthy that geomagnetic indices, in particular low-latitude aurorae (Silverman, 2006), are associated with coronal/interplanetary activity (high-speed solar-wind streams, interplanetary transients, etc.) that may not be directly related to the sunspot-cycle phase and amplitude, and therefore serve only as an approximate index of solar activity. One of the earliest instrumental geomagnetic indices is related to the daily magnetic declination range, the range of diurnal variation of magnetic needle readings at a fixed location, and is available from the 1780s (Nevanlinna, 1995). However, this data exists as several fragmentary sets, which are difficult to combine into a homogeneous data series.
Heliospheric indices are related to features of the solar wind or the interplanetary magnetic field measured (or estimated) in the interplanetary space. For example, the time evolution of the total (or open) solar magnetic flux is extensively debated (e.g., Lockwood et al., 1999;Wang et al., 2005;Krivova et al., 2007). Some details of the derivation and use of these indirect indices for long-term solar-activity studies can be found, e.g., in the Living Review by Lockwood. A special case of heliospheric indices is related to the galactic cosmic-ray intensity recorded in natural terrestrial archives. Since this indirect proxy is based on data recorded naturally throughout the ages and revealed now, it makes possible the reconstruction of solar activity changes on long timescales, as discussed in Section 3.

Solar activity observations in the pre-telescopic epoch
Instrumental solar data is based on regular observation (drawings or counting of spots) of the sun using optical instruments, e.g., the telescope invented by Galileo in the early 17th century. These observations have mostly been made by professional astronomers whose qualifications and scientific thoroughness were doubtless. They form the basis of the Group sunspot-number series (Hoyt and Schatten, 1998), which can be more-or-less reliably extended back to 1610 (see discussion in Section 2.2.1). However, some fragmentary records of qualitative solar and geomagnetic observations exist even for earlier times, as discussed below (Sections 2.3.1 -2.3.2).

Instrumental observations: Camera obscura
The invention of the telescope revolutionized astronomy. However, another solar astronomical instrument, the camera obscura, also made it possible to provide relatively good solar images and was still in use until the late 18th century. Small camera obscuras were known from early times, and they have been used in major cathedrals to define the sun's position (see the review by Vaquero, 2007). The earliest known drawing of the solar disc was made by Frisius, who observed the solar eclipse in 1544 using a camera obscura. That observation was performed during the Spörer minimum and no spots were observed on the sun. The first known observation of a sunspot using a camera obscura was done by Kepler in May 1607, who erroneously ascribed the spot on the sun to a transit of Mercury. Although such observations were sparse and related to other phenomena (solar eclipses or transits of planets), there were also regular solar observations by camera obscura. For example, about 300 pages of logs of solar observations made in the cathedral of San Petronio in Bologna from 1655 -1736 were published by Eustachio Manfredi in 1736 (see the full story in Vaquero, 2007).
Therefore, observations and drawings made using camera obscura can be regarded as instrumental observations.

Naked-eye observations
Even before regular professional observations performed with the aid of specially-developed instruments (what we now regard as scientific observations) people were interested in unusual phenomena. Several historical records exist based on naked-eye observations of transient phenomena on the sun or in the sky.
From even before the telescopic era, a large amount of evidence of spots being observed on the solar disc can be traced back as far as to the middle of the 4th century BC (Theophrastus of Athens). The earliest known drawing of sunspots is dated to December 8, 1128 AD as published in "The Chronicle of John of Worcester" (Willis and Stephenson, 2001). However, such evidence from occidental and Moslem sources is scarce and mostly related to observations of transits of inner planets over the sun's disc, probably because of the dominance of the dogma on the perfectness of the sun's body, which dates back to Aristotle's doctrine (Bray and Loughhead, 1964). Oriental sources are much richer for naked-eye sunspot records, but that data is also fragmentary and irregular (see, e.g., Clark and Stephenson, 1978;Wittmann and Xu, 1987;Yau and Stephenson, 1988). Spots on the sun are mentioned in official Chinese and Korean chronicles from 165 BC to 1918 AD. While these chronicles are fairly reliable, the data is not straightforward to interpret since it can be influenced by meteorological phenomena, e.g., dust loading in the atmosphere due to dust storms (Willis et al., 1980) or volcanic eruptions (Scuderi, 1990) can facilitate sunspots observations. Direct comparison of Oriental naked-eye sunspot observations and European telescopic data shows that naked-eye observations can serve only as a qualitative indicator of sunspot activity, but can hardly be quantitatively interpreted (see, e.g., Willis et al., 1996, and references therein). Moreover, as a modern experiment of naked-eye observations (Mossman, 1989) shows, Oriental chronicles contain only a tiny ( 1 / 200 -1 / 1000 ) fraction of the number of sunspots potentially visible with the naked eye (Eddy et al., 1989). This indicates that records of sunspot observations in the official chronicles were highly irregular (Eddy, 1983) and probably dependent on dominating traditions during specific historical periods (Clark and Stephenson, 1978). Although naked-eye observations tend to qualitatively follow the general trend in solar activity according to a posteriori information (Vaquero et al., 2002, e.g.,), extraction of any independent quantitative information from these records seems impossible.
Visual observations of aurorae borealis at middle latitudes form another proxy for solar activity (e.g., Siscoe, 1980;Schove, 1983;Křivský, 1984;Silverman, 1992;Schröder, 1992;Lee et al., 2004;Basurah, 2004). Fragmentary records of aurorae can be found in both occidental and oriental sources since antiquity. The first known dated notation of an aurora is from March 12, 567 BC from Babylon (Stephenson et al., 2004). Aurorae may appear at middle latitudes as a result of enhanced geomagnetic activity due to transient interplanetary phenomena. Although auroral activity reflects coronal and interplanetary features rather than magnetic fields on the solar surface, there is a strong correlation between long-term variations of sunspot numbers and the frequency of aurora occurrences. Because of the phenomenon's short duration and low brightness, the probability of seeing aurora is severely affected by other factors such as the weather, the Moon's phase, season (night duration), etc. The fact that these observations were not systematic in early times (before the beginning of the 18th century) makes it difficult to produce a homogeneous data set. Moreover, the geomagnetic latitude of the same geographical location may change quite dramatically over centuries, due to the migration of the geomagnetic axis, which also affects the probability of watching aurorae (Siscoe and Verosub, 1983;Oguti and Egeland, 1995). For example, the geomagnetic latitude of Seoul (37.5 • N 127 • E), which is currently less than 30 • , was about 40 • a millennium ago . This dramatic change alone can explain the enhanced frequency of aurorae observations recorded in oriental chronicles.

Mathematical/statistical extrapolations
Due to the lack of reliable information regarding solar activity in the pre-instrumental era, it seems natural to try to extend the sunspot series back in time, before 1610, by means of extrapolating its statistical properties. Indeed, numerous attempts of this kind have been made even recently (e.g., Nagovitsyn, 1997;de Meyer, 1998;Rigozo et al., 2001). Such models aim to find the main feature of the actually-observed sunspot series, e.g., a modulated carrier frequency or a multi-harmonic representation, which is then extrapolated backwards in time. The main disadvantage of this approach is that it is not a reconstruction based upon measured or observed quantities, but rather a "post-diction" based on extrapolation. This method is often used for short-term predictions, but it can hardly be used for the reliable long-term reconstruction of solar activity. In particular, it assumes that the sunspot time series is stationary, i.e., a limited time realization contains full information on its future and past. Clearly such models cannot include periods exceeding the time span of observations upon which the extrapolation is based. Hence, the pre or post-diction becomes increasingly unreliable with growing extrapolation time and its accuracy is hard to estimate.
Sometimes a combination of the above approaches is used, i.e., a fit of the mathematical model to indirect qualitative proxy data. In such models a mathematical extrapolation of the sunspot series is slightly tuned and fitted to some proxy data for earlier times. For example, Schove (1955Schove ( , 1979 fitted the slightly variable but phase-locked carrier frequency (about 11 years) to fragmentary data from naked-eye sunspot observations and auroral sightings. The phase locking is achieved by assuming exactly nine solar cycles per calendar century. This series, known as Schove series, reflects qualitative long-term variations of the solar activity, including some grand minima, but cannot pretend to be a quantitative representation in solar activity level. The Schove series played an important historical role in the 1960s. In particular, a comparison of the ∆ 14 C data with this series succeeded in convincing the scientific community that secular variations of 14 C in tree rings have solar and not climatic origins (Stuiver, 1961). This formed a cornerstone of the precise method of solar-activity reconstruction, which uses cosmogenic isotopes from terrestrial archives. However, attempts to reconstruct the phase and amplitude of the 11-year cycle, using this method, were unsuccessful. For example, Schove (1955) made predictions of forthcoming solar cycles up to 2005, which failed. We note that all these works are not able to reproduce, for example, the Maunder minimum (which cannot be represented as a result of the superposition of different harmonic oscillations), yielding too high sunspot activity compared to that observed. From the modern point of view, the Schove series can be regarded as archaic, but it is still in use in some studies. This approach has been modified using nonlinear relations (Nagovitsyn, 1997), but its shortcomings are still limiting the reliability of reconstructions.

Quasi-periodicities
The main feature of solar activity is its pronounced quasi-periodicity with a period of about 11 years, known as the Schwabe cycle. However, the cycle varies in both amplitude and duration. The first observation of a possible regular variability in sunspot numbers was made by the Danish astronomer Christian Horrebow in the 1770s on the basis of his sunspot observations from 1761 -1769 (see details in Gleissberg, 1952;Vitinsky, 1965), but the results were forgotten. It took over 70 years before the amateur astronomer Schwabe announced in 1844 that sunspot activity varies cyclically with a period of about 10 years. This cycle, called the 11-year or Schwabe cycle, is the most prominent variability in the sunspot-number series. It is recognized now as a fundamental feature of solar activity originating from the solar-dynamo process. This 11-year cyclicity is prominent in many other parameters including solar, heliospheric, geomagnetic, space weather, climate and others. The background for the 11-year Schwabe cycle is the 22-year Hale magnetic polarity cycle. Hale found that the polarity of sunspot magnetic fields changes in both hemispheres when a new 11-year cycle starts (Hale et al., 1919). This relates to the reversal of the global magnetic field of the sun with the period of 22 years. It is often considered that the 11-year Schwabe cycle is the modulo of the sign-alternating Hale cycle (e.g., Sonett, 1983;Bracewell, 1986;Kurths and Ruzmaikin, 1990;de Meyer, 1998;Mininni et al., 2001). A detailed review of solar-cyclic variability can be found in (Hathaway, 2008).
Sometimes the regular time evolution of solar activity is broken up by periods of greatly depressed activity called grand minima. The last grand minimum (and the only one covered by direct solar observations) was the famous Maunder minimum from 1645 -1715 (Eddy, 1976(Eddy, , 1983. Other grand minima in the past, known from cosmogenic isotope data, include, e.g., the Spörer minimum around 1450 -1550 and the Wolf minimum around the 14th century (see the detailed discussion in Section 4.3). Sometimes the Dalton minimum (ca. 1790 -1820) is also considered to be a grand minimum. However, sunspot activity was not completely suppressed and still showed Schwabe cyclicity during the Dalton minimum. As suggested by Schüssler et al. (1997), this can be a separate, intermediate state of the dynamo between the grand minimum and normal activity, or an unsuccessful attempt of the sun to switch to the grand minimum state (Frick et al., 1997;Sokoloff, 2004). This is observed as the phase catastrophe of solar-activity evolution (e.g., Vitinsky et al., 1986;Kremliovsky, 1994). A peculiarity in the phase evolution of sunspot activity around 1800 was also noted by Sonett (1983) who ascribed it to a possible error in Wolf sunspot data and by Wilson (1988a), who reported on a possible misplacement of sunspot minima for cycles 4 -6 in the WSN series. It has been also suggested that the phase catastrophe can be related to a tiny cycle, which might have been lost at the end of the 18th century because of very sparse observations (Usoskin et al., 2001b(Usoskin et al., , 2002a(Usoskin et al., , 2003bZolotova and Ponyavin, 2007). Another period associated with a phase catastrophe is the period of sunspot-activity recovery after the Maunder minimum. The recovery of the 11-year cycle passed through the time of distorted phase evolution (Usoskin et al., 2001a) with the concurrent loss of the 88-year cycle phase (Feynman and Gabriel, 1990).
The long-term change (trend) in the Schwabe cycle amplitude is known as the secular Gleissberg cycle (Gleissberg, 1939). However, the Gleissberg cycle is not a cycle in the strict periodic sense but rather a modulation of the cycle envelope with a varying timescale of 60 -120 years (e.g., Gleissberg, 1971;Kuklin, 1976;Ogurtsov et al., 2002). This secular cycle has also been reported, using a spectral analysis of radiocarbon data as a proxy for solar activity (see Section 3), to exist on long timescales (Feynman and Gabriel, 1990;Peristykh and Damon, 2003), but the question of its phase locking and persistency/intermittency still remain open.
Longer (super-secular) cycles cannot be studied using direct solar observations, but several such cycles have been found in cosmogenic isotope data. A cycle with a period of 205 -210 years, called the de Vries or Suess cycle in different sources, is a prominent feature, observed in varying cosmogenic data (e.g., Suess, 1980;Sonett and Finney, 1990;Zhentao, 1990;. Sometimes variations with a characteristic time of 600 -700 years or 1000 -1200 years are discussed (e.g., Vitinsky et al., 1986;Sonett and Finney, 1990;Vasiliev and Dergachev, 2002), but they are intermittent and can hardly be regarded as a typical feature of solar activity. A 2000 -2400-year cycle is also noticeable in radiocarbon data series (see, e.g., Vitinsky et al., 1986;Damon and Sonett, 1991;Vasiliev and Dergachev, 2002). However, the nonsolar origin of these super-secular cycles (e.g., geomagnetic or climatic variability) cannot be excluded.

Randomness vs. regularity
The short-term (days -months) variability of sunspot numbers is greater than the observational uncertainties indicating the presence of random fluctuations (noise). As typical for most real signals, this noise is not uniform (white), but rather red or correlated noise (e.g., Ostryakov and Usoskin, 1990;Oliver and Ballester, 1996;Frick et al., 1997), namely, its variance depends on the level of the signal. While the existence of regularity and randomness in sunspot series is apparent, their relationship is not clear (e.g., Wilson, 1994) -are they mutually independent or intrinsically tied together? Moreover, the question of whether randomness in sunspot data is due to chaotic or stochastic processes is still open.
Earlier it was common to describe sunspot activity as a multi-harmonic process with several basic harmonics (e.g., Vitinsky, 1965;Sonett, 1983;Vitinsky et al., 1986) with an addition of random noise, which plays no role in the solar-cycle evolution. However, it has been shown (e.g., Rozelot, 1994;Weiss and Tobias, 2000;Charbonneau, 2001;Mininni et al., 2002) that such an oversimplified approach depends on the chosen reference time interval and does not adequately describe the long-term evolution of solar activity. A multi-harmonic representation is based on an assumption of the stationarity of the benchmark series, but this assumption is broadly invalid for solar activity (e.g., Kremliovsky, 1994;Sello, 2000;Polygiannakis et al., 2003). Moreover, a multi-harmonic representation cannot, for an apparent reason, be extrapolated to a timescale larger than that covered by the benchmark series. The fact that purely mathematical/statistical models cannot give good predictions of solar activity (as will be discussed later) implies that the nature of the solar cycle is not a multi-periodic or other purely deterministic process, but random (chaotic or stochastic) processes play an essential role in sunspot formation. Different numeric tests, such as an analysis of the Lyapunov exponents (Ostriakov and Usoskin, 1990;Mundt et al., 1991;Kremliovsky, 1995;Sello, 2000), Kolmogorov entropy (Carbonell et al., 1994;Sello, 2000) and Hurst exponent (Ruzmaikin et al., 1994;Oliver and Ballester, 1998), confirm the chaotic/stochastic nature of the solar-activity time evolution (see, e.g., the recent review by Panchev and Tsekov, 2007).
It was suggested quite a while ago that the variability of the solar cycle may be a temporal realization of a low-dimensional chaotic system (Ruzmaikin, 1981, e.g.,). This concept became popular in the early 1990s, when many authors considered solar activity as an example of lowdimensional deterministic chaos, described by the strange attractor (Kurths and Ruzmaikin, 1990;Ostriakov and Usoskin, 1990;Morfill et al., 1991;Mundt et al., 1991;Rozelot, 1995;Salakhutdinova, 1999;Serre and Nesme-Ribes, 2000, e.g.,). Such a process naturally contains randomness, which is an intrinsic feature of the system rather than an independent additive or multiplicative noise. However, although this approach easily produces features seemingly similar to those of solar activity, quantitative parameters of the low-dimensional attractor have varied greatly as obtained by different authors. Later it was realized that the analyzed data set was too short (Carbonell et al., 1993(Carbonell et al., , 1994, and the results were strongly dependent on the choice of filtering methods (Price et al., 1992). Developing this approach, Mininni et al. (2000Mininni et al. ( , 2001 suggest that one consider sunspot activity as an example of a 2D Van der Pol relaxation oscillator with an intrinsic stochastic component. Such phenomenological or basic principles models, while succeeding in reproducing (to some extent) the observed features of solar-activity variability, do not provide insight into the nature of regular and random components of solar variability. In this sense efforts to understand the nature of randomness in sunspot activity in the framework of dynamo theory are more advanced. Corresponding theoretical dynamo models have been developed (see reviews by Ossendrijver, 2003;Charbonneau, 2005), which include stochastic processes (Weiss et al., 1984;Feynman and Gabriel, 1990;Schmalz and Stix, 1991;Moss et al., 1992;Brooke and Moss, 1994;Lawrence et al., 1995;Schmitt et al., 1996;Charbonneau and Dikpati, 2000, e.g.,). For example, Feynman and Gabriel (1990) suggest that the transition from a regular to a chaotic dynamo passes through bifurcation. Charbonneau and Dikpati (2000) studied stochastic fluctuations in a Babcock-Leighton dynamo model and succeeded in the qualitative reproduction of the anti-correlation between cycle amplitude and length (Waldmeier rule). Their model also predicts a phase-lock of the Schwabe cycle, i.e., that the 11-year cycle is an internal "clock" of the sun. Note that a significant fluctuating component (with the amplitude more than 100% of the regular component) is essential in all these model.

A note on solar activity predictions
Randomness (see Section 2.4.2) in the SN series is directly related to the predictability of solar activity. Forecasting solar activity has been a subject of intense study for many years (e.g., Yule, 1927;Newton, 1928;Gleissberg, 1948;Vitinsky, 1965) and has greatly intensified recently with a hundred journal articles being published before 2008 (see, e.g., the review by Kane, 2007) following the boosting of space-technology development and increasing debates on solar-terrestrial relations.
All prediction methods can be generically divided into precursor and statistical techniques or their combinations (Hathaway et al., 1999). The precursor methods are usually based on phe-nomenological, but sometimes physical, links between the poloidal solar-magnetic field, estimated, e.g., from geomagnetic activity in the declining phase of the preceding cycle or in the minimum time (Hathaway et al., 1999, e.g.,), with the toroidal field responsible for sunspot formation. These methods usually yield better short-term predictions of a forthcoming cycle maximum than the statistical methods, but cannot be applied to timescales longer than one solar cycle. Statistical methods, including a low-dimensional solar-attractor representation (Kurths and Ruzmaikin, 1990), are based solely on the statistical properties of sunspot activity and may give a reasonable result for short-term forecasting, but yield very poor results for long-term predictions (see reviews by Conway, 1998;Hathaway et al., 1999;Li et al., 2001;Usoskin and Mursula, 2003;Kane, 2007, e.g.,) because of chaotic/stochastic behavior (see Section 2.4.2).
Some models, mostly based on precursor method, succeed in reasonable predictions of a forthcoming solar cycle (i.e., several years ahead), but they do not pretend to extend further in time.
On the other hand, many claims of the solar activity forecast for 40 -50 years ahead and even beyond have been made recently, often without sensible argumentation. However, so far there is no evidence of any method giving a reasonable prediction of solar activity beyond the solar-cycle scale (see, e.g., Section 2.3.3), probably because of the intrinsic limit of solar-activity predictability due to its stochastic/chaotic nature (Kremliovsky, 1995). Accordingly, such attempts can be regarded as speculative, unless they are verified by the actual behavior of solar activity. Note that even an exact prediction of the amplitude of one solar cycle can be just a random coincidence and cannot serve as a proof of the method's veracity. Only a sequence of successful predictions can form a basis for confidence, which requires several decades.

Summary
In this section, the concept of solar activity and quantifying indices is discussed, as well as the main features of solar-activity temporal behavior.
The concept of solar activity is quite broad and covers nonstationary and nonequilibrium (often eruptive) processes, in contrast to the "quiet" sun concept, and their effects upon the terrestrial and heliospheric environment. Many indices are used to quantify different aspects of variable solar activity. Quantitative indices include direct (i.e., related directly to solar variability) and indirect (i.e., related to terrestrial and interplanetary effects caused by solar activity), they can be physical or synthetic. While all indices depict the dominant 11-year cyclic variability, their relationships on other timescales (short scale or long-term trends) may vary to a great extent.
The most common and the longest available index of solar activity is the sunspot number, which is a synthetic index and is very useful for the quantitative representation of overall solar activity outside the grand minimum. The sunspot number series is available for the period from 1610 AD, after the invention of the telescope, and covers, in particular, the Maunder minimum in the late 17th century. Fragmentary noninstrumental observations of the sun before 1610, while giving a possible hint of relative changes in solar activity, cannot be interpreted in a quantitative manner.
Solar activity in all its manifestations is dominated by the 11-year Schwabe cycle, which has, in fact, a variable length of 9 -14 years for individual cycles. The amplitude of the Schwabe cycle varies greatly -from the almost spotless Maunder minimum to the very high cycle 19, possibly in relation to the Gleissberg or secular cycle. Longer super-secular characteristic times can also be found in various proxies of solar activity, as discussed in Section 4. Solar activity contains essential chaotic/stochastic components, that lead to irregular variations and make the prediction of solar activity for a timescale exceeding one solar cycle impossible.

The Proxy Method of Past Solar-Activity Reconstruction
In addition to direct solar observations, described in Section 2.2.1, there are also indirect solar proxies, which are used to study solar activity in the pre-telescopic era. Unfortunately, we do not have any reliable data that could give a direct index of solar variability before the beginning of the sunspot-number series. Therefore, one must use indirect proxies, i.e., quantitative parameters, which can be measured nowadays but represent different effects of solar magnetic activity in the past. It is common to use, for this purpose, signatures of terrestrial indirect effects induced by variable solar-magnetic activity, that is stored in natural archives. Such traceable signatures can be related to nuclear (used in the cosmogenic-isotope method) or chemical (used in the nitrate method -Section 5.2) effects caused by cosmic rays (CRs) in the Earth's atmosphere, lunar rocks or meteorites.
The most common proxy of solar activity is formed by the data on cosmogenic radionuclides (e.g., 10 Be and 14 C), which are produced by cosmic rays in the Earth's atmosphere (e.g, Stuiver and Quay, 1980;Beer et al., 1990;Bard et al., 1997;Beer, 2000). Other cosmogenic nuclides, which are used in geological and paleomagnetic dating, are less suitable for studies of solar activity (see e.g., Beer, 2000). Cosmic rays are the main source of cosmogenic nuclides in the atmosphere (excluding anthropogenic factors during the last decades) with the maximum production being in the upper troposphere/stratosphere. After a complicated transport in the atmosphere, the cosmogenic isotopes are stored in natural archives such as polar ice, trees, marine sediments, etc. This process is also affected by changes in the geomagnetic field and climate. Cosmic rays experience heliospheric modulation due to solar wind and the frozen-in solar magnetic field. The intensity of modulation depends on solar activity and, therefore, cosmic-ray flux and the ensuing cosmogenic isotope intensity depends inversely on solar activity. An important advantage of the cosmogenic data is that primary archiving is done naturally in a similar manner throughout the ages, and these archives are measured nowadays in laboratories using modern techniques. If necessary, all measurements can be repeated and improved, as has been done for some radiocarbon samples. In contrast to fixed historical archival data (such as sunspot or auroral observations) this approach makes it possible to obtain homogeneous data sets of stable quality and to improve the quality of data with the invention of new methods (such as accelerator mass spectrometry). Cosmogenic isotope data is the only regular indicator of solar activity on the very long-term scale but it cannot always resolve the details of individual solar cycles. The redistribution of nuclides in terrestrial reservoirs and archiving may be affected by local and global climate/circulation processes, which are, to a large extent, unknown for the past. However, a combined study of different nuclides data, whose responses to terrestrial effects are very different, may allow for disentangling external and terrestrial signals.

Heliospheric modulation of cosmic rays
Before reaching the vicinity of Earth, galactic cosmic rays experience complicated transport in the heliosphere that leads to modulation of their flux. Heliospheric transport of GCR is described by Parker's theory (Parker, 1965;Toptygin, 1985) and includes four basic processes: the diffusion of particles due to their scattering on magnetic inhomogeneities, the convection of particles by out-blowing solar wind, adiabatic energy losses in expanding solar wind, drifts of particles in the magnetic field, including the gradient-curvature drift in the regular heliospheric magnetic field, and the drift along the heliospheric current sheet, which is a thin magnetic interface between the two heliomagnetic hemispheres. Because of variable solar-magnetic activity, CR flux in the vicinity of Earth is strongly modulated (see Figure 2). The most prominent feature in CR modulation is the 11-year cycle, which is in inverse relation to solar activity. The 11-year cycle in CR is delayed (from a month up to two years) with respect to the sunspots (Usoskin et al., 1998). The time profile of cosmic-ray flux, as measured by the longest-running neutron-monitor climax, is shown in Figure 2 (panel b) together with the sunspot numbers (panel a). Besides the inverse relation between them, some other features can also be noted. A 22-year cyclicity manifests itself in cosmic-ray modulation through the alteration of sharp and flat maxima in cosmic-ray data, originated from the chargedependent drift mechanism. One may also note short-term fluctuations, which are not directly related to sunspot numbers but are driven by interplanetary transients caused by solar eruptive events, e.g., flares or CMEs. For the last 50 years of high and roughly-stable solar activity, no trends have been observed in CR data; however, as will be discussed later, the overall level of CR has changed significantly on the centurial-millennial timescales. Full solution of the CR transport problems is a complicated task and requires sophisticated 3D time-dependent self-consistent modelling. However, the problem can be essentially simplified for applications at a long-timescale. An assumption on the azimuthal symmetry (requires times longer that the solar-rotation period) and quasi-steady changes reduces it to a 2D quasi-steady problem. Further assumption of the spherical symmetry of the heliosphere reduces the problem to a 1D case. This approximation can be used only for rough estimates, since it neglects the drift effect, but it is useful for long-term studies, when the heliospheric parameters cannot be evaluated independently. Further, but still reasonable, assumptions (constant solar-wind speed, roughly power-law CR energy spectrum, slow spatial changes of the CR density) lead to the forcefield approximation (Gleeson and Axford, 1968), which can be solved analytically. The differential intensity J i of the cosmic-ray nuclei of type i with kinetic energy T at 1 AU is given in this case as where Φ i = (Z i e/A i )φ for a cosmic nuclei of i-th type (charge and mass numbers are Z i and A i ), T and φ are expressed in MeV/nucleon and in MV, respectively, T r = 938 MeV. T is the CR particle's kinetic energy, and φ is the modulation potential. The local interstellar spectrum (LIS) J LIS forms the boundary condition for the heliospheric transport problem. Since LIS is not measured directly, i.e., outside the heliosphere, it is not well known in the energy range affected by CR modulation (below 100 GeV). Presently-used approximations for LIS (e.g., Garcia-Muñoz et al., 1975;Burger et al., 2000;Webber and Higbie, 2003) agree with each other for energies above 20 GeV but may contain uncertainties of up to a factor of 1.5 around 1 GeV. These uncertainties in the boundary conditions make the results of the modulation theory slightly model-dependent (see discussion in Usoskin et al., 2005a) and require the LIS model to be explicitly cited. This approach gives results, which are at least dimensionally consistent with the full theory and can be used for long-term studies 1 (Usoskin et al., 2002b;Caballero-Lopez and Moraal, 2004). Differential CR intensity is described by the only time-variable parameter, called the modulation potential φ, which is mathematically interpreted as the averaged rigidity (i.e., the particle's momentum per unit of charge) loss of a CR particle in the heliosphere. However, it is only a formal spectral index whose physical interpretation is not straightforward, especially on short timescales and during active periods of the sun (Caballero-Lopez and Moraal, 2004). Despite its cloudy physical meaning, this force-field approach provides a very useful and simple single-parametric approximation for the differential spectrum of GCR, since the spectrum of different GCR species directly measured near the Earth can be perfectly fitted by Equation 3 using only the parameter φ in a wide range of solar activity level (Usoskin et al., 2005a). Therefore, changes in the whole energy spectrum (in the energy range from 100 MeV/nucleon to 100 GeV/nucleon) of cosmic rays due to the solar modulation can be described by this single number within the framework of the adopted LIS. The concept of modulation potential is a key concept for the method of solar-activity reconstruction by cosmogenic isotope proxy.

Geomagnetic shielding
Cosmic rays are charged particles and therefore are affected by the Earth's magnetic field. Thus the geomagnetic field puts an additional shielding on the incoming flux of cosmic rays. The shielding effect of the geomagnetic field is usually expressed in terms of the cutoff rigidity P c , which is the minimum rigidity a CR particle must posses (on average) in order to reach the ground at a given location and time (Cooke et al., 1991). Neglecting such effects as the East-West asymmetry, which is roughly averaged out for the isotropic particle flux, or nondipole magnetic momenta, which decay rapidly with distance, one can come to a simple approximation, called the Störmer's equation, that describes the vertical geomagnetic cutoff rigidity P c : where M is the geomagnetic dipole moment (in 10 25 G cm 3 ), R o is the Earth's mean radius, R is the distance from the given location to the dipole center, and λ G is the geomagnetic latitude. This approximation provides a good compromise between simplicity and reality, especially when using the eccentric dipole description of the geomagnetic field (Fraser-Smith, 1987). The eccentric dipole has the same dipole moment and orientation as the centered dipole, but the dipole's center and consequently the poles, defined as crossings of the axis with the surface, are shifted with respect to geographical ones. The shielding effect is the strongest at the geomagnetic equator, where the present-day value of P c may reach up to 17 GV in the region of India. There is almost no cutoff in the geomagnetic polar regions (λ G ≥ 60 • ). However, even in the latter case the atmospheric cutoff becomes important, i.e., particles must have rigidity above 0.5 GV in order to initiate the atmospheric cascade (see subsection 3.1.3).
The geomagnetic field is seemingly stable on the short-term scale, but it changes essentially on centurial-to-millennial timescales (e.g., Korte and Constable, 2006). Such past changes can be evaluated based on measurements of the residual magnetization of independently-dated samples. These can be paleo (i.e., natural stratified archives such as lake or marine sediments) or archaeological (e.g., clay bricks that preserve magnetization upon baking) samples. Most paleo-magnetic data preserve not only the magnetic field intensity but also the direction of the local field, while archeo-magnetic samples provide information on the intensity only. Using a large database of such samples, it is possible to reconstruct (under reasonable assumptions) the large-scale magnetic field of the Earth. Data available provides good global coverage for the last 3 millennia, allowing for a reliable paleomagnetic reconstruction of the true dipole moment (DM) or virtual dipole moment 2 (VDM) and its orientation (the ArcheoInt collaboration - Genevey et al., 2008). Less precise, but still reliable reconstructions of the DM and its orientation are possible for the last seven millennia (the CALS7K.2 model by Korte and Constable, 2005), however they may slightly underestimate the dipole moment, especially in the earlier part of the period . No reliable directional paleomagnetic reconstruction is presently possible on a longer timescale, because of the spatial sparseness of the paleo/archeo-magnetic samples in the earlier part of the Holocene. Therefore, only the virtual axial dipole moment 3 (VADM) can be estimated before ca. 5000 BC (Y00 model by Yang et al., 2000). Note that the strong assumption of the coincidence of magnetic and geographical axes may lead to an overestimate of the true dipole moment, if the magnetic axis is actually tilted. Accordingly, we assume that the two models of Y00 and CALS7K.2 bind the true dipole moment variations. These paleomagnetic reconstructions are shown in Figure 3. All paleomagnetic models depict a similar long-term trend -an enhanced intensity during the period between 1500 BC and 500 AD and a significantly lower field before that.
Changes in the dipole moment M inversely modulate the flux of CR at Earth, with strong effects in tropical regions and globally. The migration of the geomagnetic axis, which changes the geomagnetic latitude λ G of a given geographical location is also important; while not affecting the global flux of CR, it can dramatically change the CR effect regionally, especially at middle and high latitude . These changes affect the flux of CR impinging on the Earth's atmosphere both locally and globally and must be taken into account when reconstructing solar activity from terrestrial proxy data. Accounting for these effects is quite straightforward provided the geomagnetic changes in the past are known independently, e.g., from archeo and paleomagnetic studies. However, because of progressively increasing uncertainties of paleomagnetic reconstructions back in time, it presently forms the main difficulty for the proxy method on the long-term scale , especially in the early part of the Holocene. On the other hand, the geomagnetic field variations are relatively well known for the last few millennia (Genevey et al., 2008;Korte and Constable, 2008).

Cosmic-ray-induced atmospheric cascade
When an energetic CR particle enters the atmosphere, it first moves straight in the upper layers, suffering mostly from ionization energy losses that lead to the ionization of the ambient rarefied air. However, after traversing some amount of matter (the nuclear interaction mean-free path is on the order of 100 g/cm 2 for a proton in the air) the CR particle may collide with a nucleus in the atmosphere, producing a number of secondaries. These secondaries have their own fate in the atmosphere, in particular they may suffer further collisions and interactions forming an atmospheric cascade (Dorman, 2004, e.g.,). Because of the thickness of the Earth's atmosphere (1033 g/cm 2 at sea level) the number of subsequent interactions can be large, leading to a fullydeveloped cascade (also called an air shower) consisting of secondary rather than primary particles. Three main components can be separated in the cascade: • The "hadronic" nucleonic component is formed by the products of nuclear collisions of primary cosmic rays and their secondaries with the atmospheric nuclei, and consists mostly of superthermal protons and neutrons.
• The "soft" or electromagnetic component consists of electrons, positrons and photons.
• The "hard" or muon component consists mostly of muons; pions are short lived and decay almost immediately upon production, feeding muons and the "soft" component.
The development of the cascade depends mostly on the amount of matter traversed and is usually linked to residual atmospheric depth rather than to the actual altitude, that may vary depending on the exact atmospheric density profile.
Cosmogenic isotopes are a by-product of the hadronic branch of the cascade (details are given below). Accordingly, in order to evaluate cosmic-ray flux from the cosmogenic isotope data, one needs to know the physics of cascade development. Several models have been developed for this cascade, in particular its hadronic branch with emphasis on the generation of cosmogenic isotope production. The first models were simplified quasi-analytical (e.g., Lingenfelter, 1963;O'Brien and Burke, 1973) or semi-empirical models (e.g., Castagnoli and Lal, 1980). With the fast advance of computing facilities it became possible to exploit the best numerical method suitable for such problems -Monte-Carlo (e.g., Masarik and Beer, 1999;Webber and Higbie, 2003;Usoskin and Kovaltsov, 2008). The fact that models, based on different independent Monte-Carlo packages, namely, a general GEANT tool and a specific CORSIKA code, yield similar results provides additional verification of the approach.

Transport and deposition
A scheme for the transport and redistribution of the two most useful cosmogenic isotopes, 14 C and 10 Be, is shown in Figure 4. After a more-or-less similar production, the two isotopes follow quite different fates, as discussed in detail in the subsections 3.2.3 and 3.3.3. Therefore, expected terrestrial effects are quite different for the isotopes and comparing them with each other can help in disentangling solar and climatic effects (see Section 3.6.3).

Radioisotope 14 C
The most commonly used cosmogenic isotope is radiocarbon 14 C. This radionuclide is an unstable isotope of carbon with a half-life T 1/2 of about 5730 years. Since the radiocarbon method is extensively used in other science disciplines where accurate dating is a key issue (e.g., archeology, paleoclimatology, quaternary geology), it was developed primarily for this task. The solar-activityreconstruction method, based on radiocarbon data, was initially developed as a by-product of the dating techniques used in archeology and Quaternary geology, in an effort to improve the quality of the dating by means of better information on the 14 C variable source function. The present-day radiocarbon calibration curve, based on a dendrochronological scale, uninterruptedly covers the whole Holocene (and extending to 26,000 BP) and provides a solid quantitative basis for studying solar activity variations on the multi-millennial time scale.

Measurements
Radiocarbon is usually measured in tree rings, which allows an absolute dating of the samples. Using a complicated technique, the 14 C activity 4 A is measured in an independently dated sample, which is then corrected for age as where t and T 1/2 are the age of the sample and the half-life of the isotope, respectively. Then the relative deviation from the standard activity A o of oxalic acid (the National Bureau of Standards) is calculated: After correction for the carbon isotope fractionating (account for the 13 C isotope) of the sample, the radiocarbon value of ∆ 14 C is calculated (see details in Stuiver and Pollach, 1977).
where δ 13 C is the per mille deviation of the 13 C content in the sample from that in the standard belemnite sample calculated similarly to Equation 6. The value of ∆ 14 C (measured in per mille ) is further used as the index of radiocarbon relative activity. The series of ∆ 14 C for the Holocene is presented in Figure 5A as  (Stuiver et al., 1998;Reimer et al., 2004). The long-term trend is caused by the geomagnetic field variations and the slow response of the oceans. Lower panel: Production rate of 14 C in the atmosphere, reconstructed from the measured ∆ 14 C .
series -Stuiver et al., 1998) collaboration of 21 dating laboratories as a result of systematic precise measurements of dated samples from around the world . A potentially interesting approach has been recently made by Lal et al. (2005), who measured the amount of 14 C directly produced by CR in polar ice. Although this method is free of the carbon-cycle influence, the first results, while being in general agreement with other methods, are not precise.

Production
The main source of radioisotope 14 C (except anthropogenic sources during the last decades) is cosmic rays in the atmosphere. It is produced as a result of the capture of a thermal neutron by atmospheric nitrogen 14 N + n → 14 C + p.
Neutrons are always present in the atmosphere as a product of the cosmic-ray-induced cascade (see Section 3.1.3) but their flux varies in time along with the modulation of cosmic-ray flux. This provides continuous source of the isotope in the atmosphere, while the sinks are isotope decay and transport into other reservoirs as described below (the carbon cycle). The connection between the cosmogenic-isotope-production rate, Q, at a given location (quantified via the geomagnetic latitude λ G ) and the cosmic-ray flux is given by Energy (GeV) Figure 6: Differential production rate for cosmogenic isotopes and ground-based neutron monitors as a function of cosmic-ray energy. Panel A): Yield functions of the globally averaged and polar 10 Be production (Webber and Higbie, 2003), global 14 C production (Castagnoli and Lal, 1980), polar neutron monitor (Clem and Dorman, 2000) as well as the energy spectrum of galactic cosmic protons for medium modulation (φ = 550 MV). Panel B): The differential production rate for global and polar 10 Be production, global 14 C production, and polar neutron monitor.
where P c is the local cosmic-ray-rigidity cutoff (see Section 3.1.2), S(P, φ) is the differential energy spectrum of CR (see Section 3.1.1) and Y (P ) is the differential yield function of cosmogenic isotope production (Castagnoli and Lal, 1980;Masarik and Beer, 1999). Because of the global nature of the carbon cycle and its long attenuation time, the radiocarbon is globally mixed before the final deposition, and Equation 9 should be integrated over the globe. The yield function Y (P ) of the 14 C production is shown in Figure 6A together with those for 10 Be (see Section 3.3.2) and for a ground-based neutron monitor (NM), which is the main instrument for studying cosmic-ray variability during the modern epoch. One can see that the yield function increases with the energy of CR. On the other hand, the energy spectrum of CR decreases with energy. Accordingly, the differential production rate (i.e., the product of the yield function and the spectrum, F = Y · Sthe integrand of Equation 9), shown in Figure 6B, is more informative. The differential production rate reflects the sensitivity to cosmic rays, and the total production rate is simply an integral of F over energy above the geomagnetic threshold. One can see that the sensitivity of 14 C to CR peaks at a few GeV and is quite close to that of a neutron monitor (Beer, 2000). : Globally-averaged production rate of 14 C as a function of the modulation potential φ and geomagnetic dipole moment M , computed using the yield function by Castagnoli and Lal (1980), LIS by Burger et al. (2000) and cosmic-ray-modulation model by Usoskin et al. (2005a). Another often used model (Masarik and Beer, 1999) yields a similar result.
Thanks to the development of atmospheric cascade models (Section 3.1.3), there are numerical models that allow one to compute the radiocarbon production rate as a function of the modulation potential φ and the geomagnetic dipole moment M . The overall production of 14 C is shown in Figure 7.
The production rate of radiocarbon, Q14 C , can vary as affected by different factors (see, e.g., Damon and Sonett, 1991): • Variations of the cosmic-ray flux on a geological timescale due to the changing galactic background (e.g., a nearby supernova explosion or crossing the dense galactic arm).
• Secular-to-millennial variations are caused by the slowly-changing geomagnetic field. This is an important component of the variability, which needs to be independently evaluated from paleo and archeo-magnetic studies.
• Modulation of cosmic rays in the heliosphere by solar magnetic activity. This variation is the primary aim of the present method.
• Short-term variability of CR on a daily scale (suppression due to interplanetary transients or enhancement due to solar energetic-particle events) cannot be resolved in radiocarbon data.
Therefore, the production rate of 14 C in the atmosphere can be modelled for a given time (namely, the modulation potential and geomagnetic dipole moment) and location. The global production rate Q is then obtained as a result of global averaging.
There is still a small unresolved discrepancy in the absolute value of the modeled 14 C production rate. Different models yield the global-average production rate of 1.75 -2.3 atoms/cm 2 /sec (see, e.g., O'Brien, 1979;Masarik and Beer, 1999;Goslar, 2001;Usoskin et al., 2006b, and references therein), which is consistent with a direct estimate of the radiocarbon reservoir, based on analyses of the specific 14 C activity on the ground, 1.76 -2.0 (Suess, 1965;Damon et al., 1978;O'Brien, 1979). On the other hand, the steady-state production calculated from the 14 C inventory in the carbon-cycle model (see Section 3.2.3) typically yields 1.6 -1.7 atoms/cm 2 /sec for the preindustrial period (e.g., Goslar, 2001, and references therein). Thus, results obtained from the carbon cycle models and production models agree only marginally in the absolute values, and this needs further detailed studies. In 14 C-based reconstructions, the pre-industrial steady-state production is commonly used.

Transport and deposition
Upon production cosmogenic radiocarbon gets quickly oxidized to carbon dioxide CO 2 and takes part in the regular carbon cycle of interrelated systems: atmosphere-biosphere-ocean (Figure 4). Because of the long residence time, radiocarbon becomes globally mixed in the atmosphere and involved in an exchange with the ocean. It is common to distinguish between an upper layer of the ocean, which can directly exchange CO 2 with the air and deeper layers. The measured ∆ 14 C comes from the biosphere (trees), which receives radiocarbon from the atmosphere. Therefore, the processes involved in the carbon cycle are quite complicated. The carbon cycle is usually described using a box model (Oeschger et al., 1974;Siegenthaler et al., 1980), where it is represented by fluxes between different carbon reservoirs and mixing within the ocean reservoir(s), as shown in Figure 8. Production and radioactive decay are also included in box models. Free parameters in a typical box model are the 14 C production rate Q, the air-sea exchange rate (expressed as turnover rate κ), and the vertical-eddy-diffusion coefficient K, which quantifies ocean ventilation. Starting from the original representation (Oeschger et al., 1974), a variety of box models have been developed, which take into account subdivisions of the ocean reservoir and direct exchange between the deep ocean and the atmosphere at high latitudes. More complex models, including a diffusive approach, are able to simulate more realistic scenarios, but they require knowledge of a large number of model parameters. These parameters can be evaluated for the present time using the bomb test -studying the transport and distribution of the radiocarbon produced during the atmospheric nuclear tests. However, for long-term studies, only the production rate is considered variable, while the gas-exchange rate and ocean mixing are kept constant. Under such assumptions, there is no sense in subdividing reservoirs or processes, and a simple carbon box model is sufficient.
Using the carbon cycle model and assuming that all its parameters are constant in time, one can evaluate the production rate Q from the measured ∆ 14 C data. This assumption is well validated for the the Holocene (Damon et al., 1978;Stuiver et al., 1991) as there is no evidence of considerable oceanic change or other natural variability of the carbon cycle (Gerber et al., 2002), and accordingly all variations of ∆ 14 C predominantly reflect the production rate. This is supported by the strong similarity of the fluctuations of the 10 Be data in polar ice cores (Section 3.3) compared to 14 C, despite their completely different geochemical fate (Bard et al., 1997). However, the changes in the carbon cycle during the last glaciation and deglaciation were dramatic, especially regarding ocean ventilation; this and the lack of independent information about the carbon cycle parameters, make it hardly possible to qualitatively estimate solar activity from 14 C before the Holocene.
First attempts to extract information on production-rate variations from measured ∆ 14 C were based on simple frequency separations of the signals. All slow changes were ascribed to climatic and geomagnetic variations, while short-term fluctuations were believed to be of solar origin. This was done by removing the long-term trend from the ∆ 14 C series and claiming the residual as being a series of solar variability (e.g., Peristykh and Damon, 2003). This oversimplified approach was natural at earlier times, before the development of carbon cycle models, but later it was replaced by the inversion of the carbon cycle (i.e., the reconstruction of the production rate from  (Broeker and Peng, 1986;Siegenthaler et al., 1980).The number on each individual box is the steady-state ∆ 14 C of this particular reservoir expressed in per mil. (After Bard et al., 1997) the measured 14 C concentration). Although mathematically this problem can be solved correctly as a system of linear differential equations, the presence of fluctuating noise with large magnitude makes it not straightforward, since the time derivative cannot be reliably identified. One traditional approach (e.g., Stuiver and Quay, 1980) is based on an iterative procedure, first assuming a constant production rate, and then fitting the calculated ∆ 14 C variations to the actual measurements using a feedback scheme. A concurrent approach based on the presentation of the carbon cycle as a Fourier filter  produces similar results. The production rate Q14 C for the Holocene is shown in Figure 5 and depicts both short-term fluctuations as well as slower variations, mostly due to geomagnetic field changes (see Section 3.2.5).

The Suess effect and nuclear bomb tests
Unfortunately, cosmogenic 14 C data cannot be easily used for the last century, primarily because of the extensive burning of fossil fuels. Since fossil fuels do not contain 14 C, the produced CO 2 dilutes the atmospheric 14 CO 2 concentration with respect to the pre-industrial epoch. Therefore, the measured ∆ 14 C cannot be straightforwardly translated into the production rate Q after the late 19th century, and a special correction for fossil fuel burning is needed. This effect, known as the Suess effect (e.g., Suess, 1955) can be up to −25 in ∆ 14 C in 1950 (Tans et al., 1979), which is an order of magnitude larger than the amplitude of the 11-year cycle of a few per mil. Moreover, while the cosmogenic production of 14 C is roughly homogeneous over the globe and time, the use of fossil fuels is highly nonuniform (e.g., de Jong and Mook, 1982) both spatially (developed countries, in the northern hemisphere) and temporarily (World Wars, Great Depression, industrialization, etc.). This makes it very difficult to perform an absolute normalization of the radiocarbon production to the direct measurements. Sophisticated numerical models (e.g., Sabine et al., 2004;Mikaloff Fletcher et al., 2006) aim to account for the Suess effect and make good progress. However, the results obtained indicate that the determination of the Suess effect does not yet reach the accuracy required for the precise modelling and reconstruction of the 14 C production for the industrial epoch. As noted by Matsumoto et al. (2004), "...Not all is well with the current generation of ocean carbon cycle models. At the same time, this highlights the danger in simply using the available models to represent state-of-the-art modeling without considering the credibility of each model." Note that the atmospheric concentration of another carbon isotope 13 C is partly affected by land use, which has also been modified during the last century.
Another anthropogenic activity greatly disturbing the natural variability of 14 C is related to the atmospheric nuclear bomb tests actively performed in the 1960s. For example, the radiocarbon concentration nearly doubled in the early 1960s in the northern hemisphere after nuclear tests performed by the USSR and the USA in 1961 (Damon et al., 1978). On one hand, such sources of momentary spot injections of radioactive tracers (including 14 C) provide a good opportunity to verify and calibrate the exchange parameters for different carbon -cycle reservoirs and circulation models (e.g., Bard et al., 1987;Sweeney et al., 2007). Thus, the present-day carbon cycle is moreor-less known. On the other hand, the extensive additional production of isotopes during nuclear tests makes it hardly possible to use the 14 C as a proxy for solar activity after the 1950s (Joos, 1994).
These anthropogenic effects do not allow one to make a straightforward link between preindustrial data and direct experiments performed during more recent decades. Therefore, the question of the absolute normalization of 14 C model is still open (see, e.g., the discussion in Solanki et al., , 2005Muscheler et al., 2005).

The effect of the geomagnetic field
As discussed in Section 3.1.2, knowledge of geomagnetic shielding is an important aspect of the cosmogenic isotope method. Since radiocarbon is globally mixed in the atmosphere before deposition, its production is affected by changes in the geomagnetic dipole moment M , while magnetic-axis migration plays hardly any role in 14 C data.
The crucial role of paleomagnetic reconstructions has long been known (e.g., Elsasser et al., 1956;Kigoshi and Hasegawa, 1966). Many earlier corrections for possible geomagnetic-field changes were performed by detrending the measured ∆ 14 C abundance or production rate Q (Stuiver and Quay, 1980;Voss et al., 1996;Peristykh and Damon, 2003), under the assumption that geomagnetic and solar signals can be disentangled in the frequency domain. Accordingly, the temporal series of either measured ∆ 14 C or its production rate Q is decomposed into the slow changing trend and faster oscillations. The trend is supposed to be entirely due to geomagnetic changes, while the oscillations are ascribed to solar variability. Such a method, however, obliterates all information on possible long-term variations of solar activity. Simplified empirical correction factors were also often used (e.g., Stuiver and Quay, 1980;Stuiver et al., 1991). The modern approach is based on a physics-based model (e.g., Vonmoos et al., 2006) and allows the quantitative reconstruction of solar activity, explicitly using independent reconstructions of the geomagnetic field. In this case the major source of errors in solar activity reconstructions is related to uncertainties in the paleomagnetic data . These errors are insignificant for the last millennium (Usoskin et al., 2006a), but become increasingly important for earlier times.

Measurements
The cosmogenic isotope 10 Be is useful for long-term studies of solar activity because of its long halflife of around 1.5 × 10 6 years. Its concentration is usually measured in stratified ice cores allowing for independent dating. Because of its long life, the beryllium concentration is difficult to measure by the decay rate (Beer, 2000). Accordingly, the 10 Be/ 9 Be ratio needs to be precisely measured at an accuracy better than 10 −13 . This can be done using AMS (Accelerator Mass Spectrometry) technique, which make the measurements complicated and expensive. Correction for the decay is straightforward and does not include isotope fractionating. From the measured samples, first the 10 Be concentration is defined, usually in units of 10 4 atoms/g. Sometimes, a correction for the snow precipitation amount is considered leading to the observable 10 Be flux, which is the number of atoms per cm 2 per second.
There exist different 10 Be series suitable for studies of long-term solar activity, coming from ice cores in Greenland and Antarctica. They have been obtained from different cores with different resolutions, and include data from Milcent, Greenland (Beer et al., 1983), Camp Century, Greenland (Beer et al., 1988), Dye 3, Greenland (Beer et al., 1990), Dome Concordia and South Pole, Antarctica (Raisbeck et al., 1990), GRIP, Greenland , GISP2, Greenland (Finkel and Nishiizumi, 1997), Dome Fuji, Antarctica (Horiuchi et al., 2007(Horiuchi et al., , 2008. We note that data on 10 Be in other archives, e.g., lake sediments, is usually more complicated to interpret because of the potential influence of the climate (Horiuchi et al., 1999;Belmaker et al., 2008).
Details of the 10 Be series and their comparison with each other can be found in Beer (2000); Muscheler et al. (2007).

Production
The isotope 10 Be is produced as a result of spallation of atmospheric nitrogen and oxygen (carbon is less abundant by far in the atmosphere and makes a negligible contribution) by the nucleonic component of the cosmic-ray-induced atmospheric cascade (Section 3.1.3). The cross section (a few mb) of the spallation reactions is almost independent of the energy of impacting particles and has a threshold of about 15 MeV. Thus, the production of 10 Be is defined solely by the multiplicity of the nucleonic component, which increases with the energy of primary cosmic rays (see Figure 6). Maximum production occurs at an altitude of 10 -15 km due to a balance between the total energy of the cascade (which increases with altitude) and the number of secondaries (decreasing with altitude). Most of the global 10 Be is produced in the stratosphere (55 -70%) and the rest in the troposphere (Lal and Peters, 1967;Masarik and Beer, 1999;Usoskin and Kovaltsov, 2008).
Computation of 10 Be isotope production is straightforward, provided a model of the atmospheric cascade is available. The first consistent model was developed by D. Lal et al. (Bhandari et al., 1966;Lal and Peters, 1967;Lal and Suess, 1968), using an empirical approach based on fitting simplified model calculations to measurements of the isotope concentrations and "star" (inelastic nuclear collisions) formations in the atmosphere. Next was an analytical model by O'Brien (1979), who solved the problem of the GCR-induced cascade in the atmosphere using an analytical stationary approximation in the form of the Boltzmann equation, which has also been normalized per "star" formation. Those models were based on calculating the rate of inelastic collisions or "stars" and then applying the mean spallation yield per "star". A new step in the modelling of isotope production was made by Masarik and Beer (1999), who performed a full Monte-Carlo simulation of a GCR-initiated cascade in the atmosphere and used cross sections of spallation reactions directly instead of the average "star" efficiency. A recent model by Webber and Higbie (2003) and Webber et al. (2007) is also based on a full Monte-Carlo simulation of the atmospheric cascade, using improved cross sections. The global production rate of 10 Be is about 0.02 -0.03 atoms cm −2 sec −1 (Masarik and Beer, 1999;Webber et al., 2007;Usoskin and Kovaltsov, 2008), which is lower than that for 14 C by two orders of magnitude (about 2 atoms cm −2 sec −1 ; see Section 3.2.2). The yield function of 10 Be production is shown in Figure 6A and the differential production rate in Figure 6B. One can see that the peak of 10 Be sensitivity, especially in polar regions, is shifted towards lower energies (below 1 GeV) compared with both 14 C and neutron monitor sensitivities. This implies that the 10 Be isotope is relatively more sensitive to less energetic CR and is, therefore, more affected by solar energetic particles .
Although the production of 10 Be can be more or less precisely modelled, a simple normalization "surface", similar to that shown in Figure 7 for 14 C, is not easy to produce because of partial mixing in the atmosphere (see Section 3.3.3). Simplified models, assuming either only global (e.g., Beer, 2000) or polar production (Bard et al., 1997;, have been used until recently. However, it has been recognized that a more realistic model of the limited atmospheric mixing should be used. Without detailed knowledge of 10 Be transport in the atmosphere, it is impossible to relate the quantitatively-measured concentration to the production (as done for 14 C using the carbon cycle), and one has to assume that the measured abundance is proportional (with an unknown coefficient) to the production rate in a specific geographical region (see Section 3.3.3).

Transport
After production, the 10 Be isotope has a seemingly simple (Figure 4) but difficult-to-accountfor fate in the atmosphere. Its atmospheric residence time depends on scavenging, stratospheretroposphere exchange and inter-tropospheric mixing (e.g., McHargue and Damon, 1991). Soon after production, the isotope becomes attached to atmospheric aerosols and follows their fate. In addition, it may be removed from the lower troposphere by wet deposition (rain and snow). The mean residence time of the aerosol-bound radionuclide in the atmosphere is quite different for the troposphere, being a few weeks, and stratosphere, where it is one to two years (Raisbeck et al., 1981). Accordingly, 10 Be produced in the troposphere is deposited mostly locally, i.e., in the polar regions, while stratospheric 10 Be can be partly or totally mixed. In addition, because of the seasonal (usually Spring) intrusion of stratospheric air into the troposphere at mid-latitudes, there is an additional contribution of stratospheric 10 Be. Therefore, the measured 10 Be concentration (or flux) in polar ice is modulated not only by production but also by climate/precipitation effects (e.g., Steig et al., 1996;Bard et al., 1997). This led Lal (1987) to the extreme conclusion that variations of polar 10 Be reflect a meteorological, rather than solar, signal. However, comparison between Greenland and Antarctic 10 Be series and between 10 Be and 14 C data (e.g., Bard et al., 1997;Horiuchi et al., 2008) suggests that the beryllium data mostly depicts production variations (i.e., solar signal) on top of which some meteorological effects can be superposed (see also Section 3.6.3).
Since both assumptions of the global and purely-local polar production of 10 Be archived in polar ice are over-simplified, several attempts have been made to overcome this problem. For instance, McCracken (2004) proposed several simple mathematical models of partial atmospheric mixing (without division in the troposphere and stratosphere) and compared them with observed data. From this semi-empirical approach McCracken concluded that M2 (full mixing above 60 • latitude and a limited mixing between 40 • and 60 • latitude) is a reasonable model for Antarctica. Vonmoos et al. (2006) assumed that the production of 10 Be recorded in Greenland is related to the entire hemisphere in the stratosphere (i.e, global stratospheric mixing) but is limited to latitudes above 40 • latitude in the troposphere (partial tropospheric mixing). This approach uses either semi-empirical or indirect arguments in choosing the unknown degree of mixing.
Recent efforts in employing modern atmospheric 3D circulation models for simulations of 10 Be transport and deposition, including realistic air-mass transport and dry-vs-wet deposition (Field et al., 2006;Heikkilä et al., 2007), look more promising. An example of 10 Be deposition computed on the world grid using the NASA GISS model (Field et al., 2006) is shown in Figure 9. Precision of the models allows one to distinguish local effects, e.g., for Greenland (Heikkilä et al., 2007). A simulation performed by combining a detailed 10 Be-production model with an air-dynamics model can result in an absolute model relating production and deposition of the radionuclide. We may expect this break through to occur in the near future.

Effect of the geomagnetic field
In order to properly account for geomagnetic changes (Section 3.1.2), one needs to know the effective region in which the radionuclide is produced before being stored in the archive analyzed. For instance, if the concentration of 10 Be measured in polar ice reflects mainly the isotope's production in the polar atmosphere (as, e.g., assumed by Usoskin et al., 2003c), no strong geomagnetic signal is expected to be observed, since the geographical poles are mostly related to high geomagnetic latitudes. On the other hand, assuming global mixing of atmospheric 10 Be before deposition in polar ice (e.g., Masarik and Beer, 1999), one expects that only changes in the geomagnetic dipole moment affect will the signal. However, because of partial mixing, which can be different in the stratosphere and troposphere, taking into account migration and displacement of the geomagnetic dipole axis may be essential for a reliable reconstruction of solar variability from 10 Be data (Mc-Cracken, 2004). Therefore, only a full combination of the transport and production models, the latter explicitly including geomagnetic effects estimated from paleomagnetic reconstructions, can adequately account for geomagnetic changes and separate the solar signal. These will form the next generation of physics-based models for the cosmogenic-isotope proxy method. We note that paleomagnetic data should ideally not only provide the dipole moment (VADM or VDM) but should also provide estimates of the geomagnetic axis attitude and displacement of the dipole center.

Towards a quantitative physical model
Several methods have been developed historically to convert measured cosmogenic-isotope data into a solar activity index, ranging from very simple regressions to physics-based models. A new step in long-term solar-activity reconstruction has been made recently, which is the development of the proxy method in which physics-based models are used, instead of a phenomenological regression, to link SN with cosmogenic-isotope production (Usoskin et al., 2003c;Vonmoos et al., 2006;Muscheler et al., 2007). Due to recent theoretical developments, it is now possible to construct a chain of physical models to simulate the entire relationship between solar activity and cosmogenic data.
The physics-based reconstruction of solar activity (in terms of sunspot numbers) from cosmogenic proxy data includes several steps: • Computation of the isotope's production rate in the atmosphere from the measured concentration in the archive (Sections 3.2.2 and 3.3.2); • Computation, considering independently-known secular geomagnetic changes (see Section 3.2.5) and a model of the CR-induced atmospheric cascade, of the GCR spectrum parameter quantified via the modulation potential φ (Section 3.4.2), some reconstructions being terminated at this point; • Computation of a heliospheric index, whether of the open solar magnetic flux or of the average HMF intensity at the Earth's orbit (Section 3.4.2) • Computation of a solar index (sunspot number series), corresponding to the above-derived heliospheric parameter (Section 3.4.3).
Presently, all these steps can be completed using appropriate models. Although the uncertainties of the models may be considerable, the models allow a full basic quantitative reconstruction of solar activity in the past. However, much needs to be done, both theoretically and experimentally, to obtain an improved reconstruction.

Regression models
Mathematical regression is the most apparent and often used (even recently) method of solaractivity reconstruction from proxy data (see, e.g., Stuiver and Quay, 1980;Ogurtsov, 2004). The reconstruction of solar activity is performed in two consecutive steps. First, a phenomenological regression (either linear or nonlinear) is built between proxy data set and a direct solar-activity index for the available "training" period (e.g., since 1750 for WSN or since 1610 for GSN). Then this regression is extended backwards to evaluate SN from the proxy data. The main shortcoming of the regression method is that it depends on the time resolution and choice of the "training" period. The former is illustrated by Figure 10, which shows the scatter plot of the 10 Be concentration vs. GSN for the annual and 11-year smoothed data.
One can see that the slope of the 10 Be-vs-GSN relation (about −500 g/atom) within individual cycles is significantly different from the slope of the long-term relation (about −100 g/atom, i.e., individual cycles do not lie on the line of the 11-year averaged cycles. Therefore, a formal regression built using the annual data for 1610 -1985 yields a much stronger GSN-vs-10 Be dependence than for the cycle-averaged data (see Figure 10b), leading to a potentially-erroneous evaluation of the sunspot number from the 10 Be proxy data.
It is equally dangerous to evaluate other solar/heliospheric/terrestrial indices from sunspot numbers, by extrapolating an empirical relation obtained for the last few decades back in time. This is because the last few decades (after the 1950s), which are well covered by direct observations of solar, terrestrial and heliospheric parameters, correspond to a very high level of solar activity. After a steep rise in activity level between the late 19th and mid 20th centuries, the activity remained at a roughly constant high level, being totally dominated by the 11-year cycle without an indication of long-term trends. Accordingly, all empirical relations built based on data for this period are focused on the 11-year variability and can overlook possible long-term trends . This may affect all regression-based reconstructions, whose results cannot be independently (directly or indirectly) tested. In particular, this may be related to solar irradiance reconstructions, which are often based on regression-like models, built and verified using data for the last three solar cycles, when there was no strong trend in solar activity.
As an example let us consider an attempt (Belov et al., 2006) to reconstruct cosmic-ray intensity since 1610 from sunspot numbers using a (nonlinear) regression. The regression between the count rate of a neutron monitor and sunspot numbers, established for the last 30 years, yields an agreement at a 95% level for the period 1976 -2003. Based on that, Belov et al. (2006) extrapolated the regression back in time to produce a reconstruction of cosmic-ray intensity (quantified in NM count rate) to 1560 (see Figure 11). One can see that there is no notable long-term trend in the reconstruction, and the fact that all CR maxima essentially lie at the same level, from the Maunder minimum to modern times, is noteworthy. It would be difficult to dispute such a result if there was Figure 11: An unsuccessful attempt at the reconstruction of cosmic-ray intensity in the past using a regression with sunspot numbers (After Belov et al., 2006). Dots represent the observed cosmic-ray intensity since 1951. Note the absence of a long-term trend. no direct test for CR levels in the past. Independent reconstructions based on cosmogenic isotopes or theoretical considerations (Usoskin et al., 2002b;Scherer and Fichtner, 2004, e.g.,) provide clear evidence that cosmic-ray intensity was essentially higher during the Maunder minimum than nowadays. This example shows how easy it is to overlook an essential feature in a reconstruction based on a regression extrapolated far beyond the period it is based on. Fortunately, for this particular case we do have independent information that can prevent us from making big errors. In many other cases, however, such information does not exist (e.g., for total or spectral solar irradiance), and those who make such unverifiable reconstructions should be careful about the validity of their models beyond the range of the established relations.

Reconstruction of heliospheric parameters
The modulation potential φ (see Section 3.1.1) is most directly related to cosmogenic isotope production in the atmosphere. It is a parameter describing the spectrum of galactic cosmic rays (see the definition and full description of this index in Usoskin et al. (2005a)) and is sometimes used as a stand-alone index of solar (or, actually, heliospheric) activity. We note that, provided the isotope production rate Q is estimated and geomagnetic changes can be properly accounted for, it is straightforward to obtain a time series of the modulation potential, using, e.g., the relation shown in Figure 7. Several reconstructions of modulation potential for the last few centuries are shown in Figure 12. While being quite consistent in the relative changes, they differ in the absolute level and fine details. Reconstructions of solar activity often end at this point, representing solar activity by the modulation potential, as some authors (e.g., Beer et al., 2003;Vonmoos et al., 2006;Muscheler et al., 2007) believe that further steps (see subsection 3.4.3) may introduce additional uncertainties. However, since φ is a heliospheric, rather than solar, index, the same uncertainties remain when using it as an index of solar activity. Moreover, the modulation potential is a model-dependent quantity (see discussion in Section 3.1.1) and therefore does not provide an unambiguous measure of heliospheric activity. In addition, the modulation potential is not a physical index but rather a formal fitting parameter to describe the GCR spectrum near Earth (Usoskin et al., 2005a) and, thus, does not seem to be a universal solar-activity index.
Modulation of GCR in the heliosphere (see Section 3.1.1) is mostly defined by the turbulent heliospheric magnetic field (HMF), which ultimately originates from the sun and is thus related to   Figure 13: An example of reconstruction of the heliospheric magnetic field at Earth orbit for the last 600 years from 10 Be data (McCracken, 2007). solar activity. It has been shown, using a theoretical model of the heliospheric transport of cosmic rays (e.g., Usoskin et al., 2002b), that on the long-term scale (beyond the 11-year solar cycle) the modulation potential φ is closely related to the open solar magnetic flux F o , which is a physical quantity describing the solar magnetic variability (e.g., Solanki et al., 2000;Krivova et al., 2007).
Sometimes, instead of the open magnetic flux, the mean HMF intensity at Earth orbit, B, is used as a heliospheric index (Caballero-Lopez and Moraal, 2004;McCracken, 2007). Note that B is linearly related to F o assuming constant solar-wind speed, which is valid on long-term scales. An example of HMF reconstruction for the last 600 years is shown in Figure 13. In addition, the count rate of a "pseudo" neutron monitor (i.e., a count rate of a neutron monitor if it was operated in the past) is considered as a solar/heliospheric index (e.g., Beer, 2000;McCracken and Beer, 2007).

A link to sunspot numbers
The open solar magnetic flux F o described above is related to the solar surface magnetic phenomena such as sunspots or faculae. Modern physics-based models allows one to calculate the open solar magnetic flux from data of solar observation, in particular sunspots (Solanki et al., 2000(Solanki et al., , 2002Krivova et al., 2007). Besides the solar-active regions, the model includes ephemeral regions. Although this model is based on physical principals, it contains one adjustable parameter, the decay time of the open flux, which cannot be measured or theoretically calculated and has to be found by means of fitting the model to data. This free parameter has been determined by requiring the model output to reproduce the best available data sets for the last 30 years with the help of a genetic algorithm. Inversion of the model, i.e., the computation of sunspot numbers for given F o values is formally a straightforward solution of a system of linear differential equations, however, the presence of noise in the real data makes it only possible in a numerical-statistical way (see, e.g., Usoskin et al., , 2007. By inverting this model one can compute the sunspot-number series corresponding to the reconstructed open flux, thus forging the final link in a chain quantitatively connecting solar activity to the measured cosmogenic isotope abundance. A sunspot-number series reconstructed for the Holocene using 14 C isotope data is shown in Figure 14. While the definition of the grand minima (Section 4.3) is virtually insensitive to the uncertainties of paleomagnetic data, : Long-term sunspot-number reconstruction from 14 C data (after Usoskin et al., 2007). All data are decade averages. Solid (denoted as 'Y00') and grey ('K05') curves are based on the paleo-geomagnetic reconstructions of Yang et al. (2000) and Korte and Constable (2005), respectively. Observed group sunspot numbers (Hoyt and Schatten, 1998) are shown after 1610.
the definition of grand maxima depends on the paleomagnetic model used .
Since the Y00 paleomagnetic model forms an upper bound for the true geomagnetic strength (Section 3.1.2), the corresponding solar-activity reconstructions may underestimate the solar-activity level. Accordingly, the grand maxima defined using the Y00 model are robust and can be regarded as "maximum maximorum" (see Section 4.4).

Solar activity reconstructions
Detailed computational models of cosmogenic isotope production in the atmosphere (e.g., Masarik and Beer, 1999) have opened up a new possibility for long-term solar-activity reconstruction (e.g., Beer, 2000). The first quantitative reconstructions of solar activity from cosmogenic proxy appeared in the early 2000s based on 10 Be deposited in polar ice (Beer et al., 2003;Usoskin et al., 2003c). Beer et al. (2003) reconstructed the modulation potential on a multi-millennial timescale using the model computations by Masarik and Beer (1999) and the 10 Be data from the GISP2 core in Greenland. This result has been extended, even including the 14 C data set, and presently covers the whole Holocene (Vonmoos et al., 2006). Usoskin et al. (2003c) presented the reconstruction of sunspot activity over the last millennium, based on 10 Be data from both Greenland and Antarctica, using a physics-based model described in detail in . This result reproduces the four known grand minima of solar activity -Maunder, Spörer, Wolf and Oort minima (see Section 4.3). Later  reconstructed 10-year-averaged sunspot numbers from the 14 C content in tree rings throughout the Holocene and estimated its uncertainties (see Figure 15). This result was disputed by Muscheler et al. (2005), whose concurrent model, however, rested on an erroneous normalization, as argued in Solanki et al. (2005). The reconstruction of  has been recently updated by Usoskin et al. (2006a), using a newer paleomagnetic reconstruction by Korte and Constable (2005), and was later slightly revised (Usoskin et Krivova et al. (2007). Reconstruction of the HMF from 10 Be data has been performed by Caballero-Lopez and Moraal (2004), using a model of CR modulation in the heliosphere and a 10 Be production model by Webber and Higbie (2003). Recently it was revised (McCracken, 2007) to present a detailed reconstruction of HMF intensity since 1428. The obtained results are discussed in Section 4.

Verification of reconstructions
Because of the diversity of the methods and results of solar-activity reconstruction, it is vitally important to verify them. Even though a full verification is not possible, there are different means of indirect or partial verification, as discussed below. Several solar-activity reconstructions on the millennium timescale, which differ from each other to some degree and are based on terrestrial cosmogenic -isotope data, have been published recently by various groups. . Also, they may suffer from systematic effects. Therefore, there is a need for an independent method to verify/calibrate these results in order to provide a reliable quantitative estimate of the level of solar activity in the past, prior to the era of direct observations.

Comparison with direct data
The most direct verification of solar-activity reconstruction is a comparison with the actual GSN data for the last few centuries. However, regression-based models (see Section 3.4.1) cannot be tested in this way, since it would require a long set of independent direct data outside the "training" interval. It is usual to include all available data into the "training" period to increase the statistics of the regression, which rules out the possibility of testing the model. On the other hand, such a comparison to the actual GSN since 1610 can be regarded as a direct test for a physics-based model since it does not include phenomenological links over the same time interval. The period of the last four centuries is pretty good for testing purposes since it includes the whole range of solar activity levels from the nearly spotless Maunder minimum to the modern period of a very active sun. As an example, a comparison between the measured GSN and the 14 C-based  and 10 Be-based (Usoskin et al., 2003c) reconstructions is shown in Figure 15. The agreement between the actual and reconstructed sunspot numbers is quite good, the correlation coefficient for the 14 C-based series is r = 0.93 with the RMS deviation between the two series being six for the period of 1610 -1900 . We want to stress that this reconstruction is fully physics based and does not include any fitting to the whole GSN data series; thus this comparison verifies the approach in both absolute level and relative variations. The agreement between GSN and 10 Be-based reconstructions (Usoskin et al., 2003c) is also good (r = 0.78, RMS = 10 for 1700 -1985). In this case, however, the comparison can only test the relative variation because of the unknown proportionality coefficient between the measured concentration of 10 Be and the production rate (Section 3.3.3), which is fitted to match the overall level of the reconstructed solar activity. One can see that the reconstructed sunspot series generally follows the real GSN series, depicting the same main features, namely, the Maunder minimum, the tiny Dalton minimum, a slight decrease of activity around 1900 (known as the modern minimum) as well as a steep rise in the first half of the 20th century. Note, however, that individual 11-year cycles are poorly resolved in these reconstructions. This validates the reliability of the physics-based reconstruction of sunspot numbers. Models focused on the reconstruction of heliospheric parameters (HMF or the modulation potential φ) cannot be verified in this manner since no heliospheric data exists before the 1950s. Comparison to direct cosmic-ray data after the 1950s (or, with caveats, after the 1930s -McCracken and  is less conclusive, since the latter are of shorter length and correspond to a period of high solar activity, leading to larger uncertainties during grand minima.

Meteorites and lunar rocks: A direct probe of the galactic cosmic-ray flux
Another more-or-less direct test of solar/heliospheric activity in the past comes from cosmogenic isotopes measured in lunar rock or meteorites. Cosmogenic isotopes, produced in meteoritic or lunar rocks during their exposure to CR in interplanetary space, provide a direct measure of cosmicray flux. Uncertainties due to imprecisely known terrestrial processes, including the geomagnetic shielding and redistribution process, are naturally avoided in this case, since the nuclides are directly produced by cosmic rays in the body, where they remain until they are measured, without any transport or redistribution. The activity of a cosmogenic isotope in meteorite/lunar rock represents an integral of the balance between the isotope's production and decay, thus representing the time-integrated CR flux over a period determined by the mean life of the radioisotope. The results of different analyses of measurements of cosmogenic isotopes in meteoritic and lunar rocks show that the average GCR flux remained roughly constant -within 10% over the last million years and within a factor of 1.5 for longer periods of up to 10 9 years (e.g., Vogt et al., 1990;Grieder, 2001).
By means of measuring the abundance of relatively short-lived cosmogenic isotopes in meteorites, which fell through the ages, one can evaluate the variability of the CR flux, since the production of cosmogenic isotopes ceases after the fall of the meteorite. A nearly ideal isotope for studying centurial-scale variability is 44 Ti with a half-life of 59.2±0.6 yr (a lifetime of about 85 years). The isotope is produced in nuclear interactions of energetic CR with nuclei of iron and nickel in the body of a meteorite (Bonino et al., 1995;Taricco et al., 2006). Because of its mean life, 44 Ti is relatively insensitive to variations in cosmic-ray flux on decade (11-year Schwabe cycle) or shorter timescales, but is very sensitive to the level of CR flux and its variations on a centurial scale. Using a full model of 44 Ti production in a stony meteorite (Michel and Neumann, 1998) and data on the measured activity of cosmogenic isotope 44 Ti in meteorites, which fell during the past 235 years (Taricco et al., 2006), Usoskin et al. (2006c) tested, in a straightforward manner, several recent reconstructions of heliospheric activity after the Maunder minimum. First, the expected 44 Ti activity has been calculated from the reconstructed series of the modulation potential, and then compared with the results of actual measurements (see Figure 16). It has been shown that 44 Ti data can distinguish between various reconstructions of past solar activity, allowing unrealistic models to be ruled out. Presently, the Torino group (Taricco et al., 2008) is working hard on improving the quality of 44 Ti measurements in meteorites, reducing the error bars in Figure 16, which will allow for more precise estimates in the near future.
Because of the long life time of the 44 Ti nuclide (about 85 years), this method does not allow for the reconstruction of solar/heliospheric activity, but it serves as a direct way to test existing reconstructions independently. Most of the reconstructions appear consistent with the measured 44 Ti activity in meteorites, including the last decades, thus validating their veracity. The only apparently-inconsistent model is the one by Muscheler et al. (2005), which is based on erroneous normalization (as discussed in Solanki et al., 2005). In particular, the 44 Ti data confirms significant secular variations of the solar magnetic flux during the last century (cf. Lockwood et al., 1999;Solanki et al., 2000;Wang et al., 2005).

Comparison between isotopes
As an indirect test of the solar-activity reconstruction, one can compare different isotopes. The idea behind this test is that two isotopes, 14 C and 10 Be, have essentially different terrestrial fates, so that only the production signal, namely, solar modulation of cosmic rays, can be regarded as common in the two series. Processes of transport/deposition are completely different (moreover, the 14 C series is obtained as an average of the world-wide-distributed samples). The effect of changing geomagnetic fields is also different (although not completely) for the two isotopes, since radiocarbon is globally mixed, while 10 Be is only partly mixed before being stored in an archive. Even comparison between data of the same 10 Be isotope, but measured in far-spaced ice cores (e.g., Greenland and Antarctica), may help in separating climatic and extraterrestrial factors, since meteorology in the two opposite polar areas is quite different.
The first thorough consistent comparison between 10 Be and 14 C records for the last millennium was performed by Bard et al. (1997). They assumed that the measured 10 Be concentration in Antarctica is directly related to CR variations. Accordingly, 14 C production was considered as proportional to 10 Be data. Then, applying a 12-box carbon-cycle model, Bard et al. (1997) computed the expected ∆ 14 C synthetic record. Finally, these 10 Be-based ∆ 14 C variations were compared with the actual measurements of ∆ 14 C in tree rings, which depicted a close agreement in the profile of temporal variation (coefficient of linear correlation r =0.81 with exact phasing). Despite some fine discrepancies, which can indicate periods of climatic influence, that result has clearly proven the dominance of solar modulation of cosmogenic nuclide production variations during the last millennium. This conclusion has been confirmed (Usoskin et al., 2003c;Muscheler et al., 2007, e.g.,) in the sense that quantitative solar-activity reconstructions, based on 10 Be and 14 C data series for the last millennium, yield very similar results, which differ only in small details. However, a longer comparison over the entire Holocene timescale suggests that, while centennial variations of solar activity reconstructed from the two isotopes are very close to each other, there might be a discrepancy in the very long-term trend (Vonmoos et al., 2006;Muscheler et al., 2007), whose nature is not clear (climate changes, geomagnetic effects or model uncertainties).
Thus, comparison of the results obtained from different sources implies that the variations of cosmogenic nuclides on the long-term scale (centuries to millennia) during the Holocene are primarily defined by the solar modulation of CR.

Summary
In this section, a proxy method of past-solar-activity reconstruction is described in detail.  (Taricco et al., 2006). Curves correspond to the theoretically expected 44 Ti activity, computed using the method of Usoskin et al. (2006c) and different reconstructions of φ shown in Figure 12.
This method is based on the use of indirect proxies of solar activity, i.e., quantitative parameters, which can be measured now, but represent signatures, stored in natural archives, of the different effects of solar magnetic activity in the past. Such traceable signatures can be related to nuclear or chemical effects caused by cosmic rays in the Earth's atmosphere, lunar rocks or meteorites. This approach allows one to obtain homogeneous data sets with stable quality and to improve the quality of data when new measurement techniques become available. It provides the only possible regular indicator of solar activity on a very long-term scale.
The most common proxy of solar activity is formed by data of the cosmogenic radionuclides, 10 Be and 14 C, produced by cosmic rays in the Earth's atmosphere. After a complicated transport in the atmosphere, these cosmogenic isotopes are stored in natural archives such as polar ice, trees, marine sediments, from where they can now be measured. This process is also affected by changes in the geomagnetic field and the climate.
Radioisotope 14 C, measured in independently-dated tree rings, forms a very useful proxy for long-term solar-activity variability. It participates in the complicated carbon cycle, which smoothes out spatial and short-term variability of isotope production. For the Holocene period, with its stable climate, it provides a useful tool for studying solar activity in the past. Existing models allow the quantitative conversion between the measured relative abundance of 14 C and the production rate in the atmosphere. The use of radiocarbon for earlier periods, the glacial and deglaciation epochs, is limited by severe climate and ocean ventilation changes. Radiocarbon data cannot be used after the end of the 19th century because of the Suess effect and atmospheric nuclear tests.
Another solar activity proxy is the cosmogenic 10 Be isotope measured in stratified polar ice cores. Atmospheric transport of 10 Be is relatively straightforward, but its details are as of yet unresolved, leading to the lack of a reliable quantitative model relating the measured isotope concentration in ice to the atmospheric production. Presently, it is common to assume that the production rate is proportional, with an unknown coefficient, to the measured concentration. However, a newly-developed generation of models, which include 3D atmospheric-circulation models, will hopefully solve this problem soon.
Modern physics-based models make it possible to build a chain, which quantitatively connects isotope production rate and sunspot activity, including subsequently the GCR flux quantified via the modulation potential, the heliospheric index, quantified via the open solar magnetic flux or the average HMF intensity at the Earth's orbit, and finally the sunspot-number series. Presently, all these steps can be made using appropriate models allowing for a full basic quantitative reconstruction of solar activity in the past. The main uncertainties in the solar-activity reconstruction arise from paleo-magnetic models and the overall normalization.
An independent verification of the reconstructions, including direct comparison with sunspot numbers, cosmogenic isotopes in meteorites and the comparison of different models with each other, confirms their veracity in both relative variations and absolute level. It also implies that the variations in cosmogenic nuclides on the long-term scale (centuries to millennia) during the Holocene are primarily defined by the solar modulation of CR.

Variability of Solar Activity Over Millennia
Several reconstructions of solar activity on multi-millennial timescales have been performed recently using physics-based models (see Section 3) from measurements of 14 C in tree rings and 10 Be in polar ice. The validity of these models for the last few centuries was discussed in Section 3.6. In this section we discuss the temporal variability of thus-reconstructed solar activity on a longer scale.

The overall activity level
Here we consider the 14 C-based decade reconstruction of sunspot numbers (shown in Figure 17). It is identical to that shown in Figure 14, but includes also a Gleissberg 1-2-2-2-1 filter in order to suppress noise and short-term fluctuations. This series forms the basis for the forthcoming analysis, while differences related to the use of other reconstructions are discussed.
First, we analyze the distribution of the occurrence frequency of sunspot-number values. The histogram for sunspot-number distribution (for the series shown in Figure 17) is shown in Figure 18. The over-decades filtered sunspot numbers range between 0 and 95. The bulk of the distribution can be roughly approximated by a normal Gaussian distribution with a mean of 31 and standard deviation of σ = 30. However, there are apparent excesses both at very low (SN<10) and high (SN>70) sunspot numbers. These excesses are very unlikely to be a result of random fluctuations or noise in the data and, as argued in Section 4.3, correspond to special states of the solar dynamo, namely, the grand minima and grand maxima. It is important that the entire distribution is moreor-less consistent with the directly-observed sunspot series after 1610, suggesting that the latter can serve as a representative sample for sunspot-activity statistics, including a grand minimum (the Maunder minimum) and the modern maximum.

Quasi-periodicities and characteristic times
In order to discuss spectral features of long-term solar-activity dynamics, we show in Figure 19 a wavelet spectral decomposition of the sunspot number reconstruction throughout the Holocene shown in Figure 17. The left-hand panels show the conventional wavelet decomposition in the timefrequency domain, while the right-hand panels depict the global spectrum, namely, an integral over the time domain, which is comparable to a Fourier spectrum. The peak in the global spectrum at about an 80-year period corresponds to the Gleissberg periodicity, known from a simple Fourier analysis of the ∆ 14 C series (Peristykh and Damon, 2003). The peak at an approximately 150 year period does not correspond to a persistent periodicity, but is formed by a few time intervals (mostly 6000 -4000 BC) and can be related to another "branch" of the secular cycle, according to Ogurtsov et al. (2002). The de Vries/Suess cycle, with a period of about 210 years, is prominent in the global spectrum, but it is intermittent and tends to become strong with around 2400 clustering time (Usoskin and Kovaltsov, 2004). Another variation with a period of around 350 years can be observed after 6000 BC. Variations with a characteristic time of 600 -700 years are intermittent and can be hardly regarded as a typical feature of solar activity. Of special interest is the 2000 -2400 year Hallstatt cycle (see, e.g., Vitinsky et al., 1986;Damon and Sonett, 1991;Vasiliev and Dergachev, 2002), which is relatively stable and mostly manifests itself as a modulation of long-term solar activity, leading to the clustering of grand minima .
On the other hand, an analysis of the occurrence of grand minima (see Section 4.3) shows no clear periodicity except for a marginal 2400 year clustering, implying that the occurrence of grand minima and maxima is not a result of long-term cyclic variability but is defined by stochastic/chaotic processes.  Usoskin et al. (2007) using geomagnetic data by Yang et al. (2000). Blue and red areas denote grand minima and maxima, respectively.

Grand minima of solar activity
A very particular type of solar activity is the grand minimum, when solar activity is greatly reduced. The most famous is the Maunder minimum in the late 17th century, which is discussed below in some detail (for a detailed review see the book by Soon and Yaskell, 2003). grand minima are believed to correspond to a special state of the dynamo (Sokoloff, 2004;Miyahara et al., 2006b), and its very existence poses a challenge for the solar-dynamo theory. It is noteworthy that dynamo models do not agree on how often such episodes occur in the sun's history and whether their appearance is regular or random. For example, the commonly used mean-field dynamo yields a fairly-regular 11-year cycle (Charbonneau, 2005), while dynamo models including a stochastic driver predict the intermittency of solar magnetic activity (Choudhuri, 1992;Schüssler et al., 1994;Schmitt et al., 1996;Ossendrijver, 2000;Weiss and Tobias, 2000;Mininni et al., 2001;Charbonneau, 2001).

The Maunder minimum
The Maunder minimum is a representative of grand minima in solar activity (e.g., Eddy, 1976), when sunspots have almost completely vanished from the solar surface, while the solar wind keeps blowing, although at a reduced pace (Cliver et al., 1998;Usoskin et al., 2001a). There is some uncertainty in the definition of its duration; the "formal" duration is 1645 -1715 (Eddy, 1976), while its deep phase with the absence of apparent sunspot cyclic activity is often considered as 1645 -1700, with the low, but very clear, solar cycle of 1700 -1712 being ascribed to a recovery or transition phase (Usoskin et al., 2000). The Maunder minimum was amazingly well covered (more  Figure 17. Hatched areas correspond to directly-observed sunspots after 1610. The curve represents the best fit normal distribution. than 95% of days) by direct sunspot observations (Hoyt and Schatten, 1996), especially in its late phase (Ribes and Nesme-Ribes, 1993). On the other hand, sunspots appeared rarely (during ∼ 2% of the days) and seemingly sporadically, without an indication of the 11-year cycle . This makes it almost impossible to apply standard methods of time-series analysis to sunspot data during the Maunder minimum (e.g., Frick et al., 1997)). Therefore, special methods such as the distribution of spotless days vs. days with sunspots (e.g., Harvey and White, 1999) or an analysis of sparsely-occurring events (Usoskin et al., 2000) should be applied in this case. Using these methods, Usoskin et al. (2001a) have shown that sunspot occurrence during the Maunder minimum was gathered into two large clusters (1652 -1662 and 1672 -1689), with the mass centers of these clusters being in 1658 and 1679 -1680. Together with the sunspot maxima before (1640) and after (1705) the deep Maunder minimum, this implies a dominant 22-year periodicity in sunspot activity throughout the Maunder minimum (Mursula et al., 2001), with a subdominant 11-year cycle emerging towards the end of the Maunder minimum (Ribes and Nesme-Ribes, 1993;Mendoza, 1997;Usoskin et al., 2000) and becoming dominant again after 1700. Similar behavior of a dominant 22-year cycle and a weak subdominant Schwabe cycle during the Maunder minimum has been found in other indirect solar proxy data: auroral occurrence (Křivský and Pejml, 1988;Schlamminger, 1990;Silverman, 1992) and 14 C data (Stuiver and Braziunas, 1993;Kocharov et al., 1995;Peristykh and Damon, 1998;Miyahara et al., 2006b). This is in general agreement with the concept of "immersion" of 11-year cycles during the Maunder minimum (Vitinsky et al., 1986, and references therein). This concept means that full cycles cannot be resolved and sunspot activity only appears as pulses around cycle-maximum times. The time behavior of sunspot activity during the Maunder minimum yields the following general scenario (Vitinsky et al., 1986;Ribes and Nesme-Ribes, 1993;Sokoloff and Nesme-Ribes, 1994;Usoskin et al., 2000Usoskin et al., , 2001aMiyahara et al., 2006b). Transition from the normal high activity to the deep minimum was sudden (within a few years) without any apparent precursor. A 22year cycle was dominant in sunspot occurrence during the deep minimum (1645 -1700), with the subdominant 11-year cycle, which became visible only in the late phase of the Maunder minimum. The 11-year Schwabe cycle started dominating solar activity after 1700. Recovery of sunspot activity from the deep minimum to normal activity was gradual, passing through a period of nearly-linear amplification of the 11-year cycle. It is interesting to note that such a qualitative evolution of a grand minimum is consistent with predictions of the stochastically-forced return map (Charbonneau, 2001).
Although the Maunder minimum is the only one with available direct sunspot observations, its predecessor, the Spörer minimum from 1450 -1550, is covered by precise bi-annual measurements of 14 C (Miyahara et al., 2006a). An analysis of this data (Miyahara et al., 2006a,b) reveals a similar pattern with the dominant 22-year cycle and suppressed 11-year cycle, thus supporting the idea that the above general scenario may be typical for a grand minimum.
A very important feature of sunspot activity during the Maunder minimum was its strong north-south asymmetry, as sunspots were only observed in the southern solar hemisphere during the end of the Maunder minimum (Ribes and Nesme-Ribes, 1993;Sokoloff and Nesme-Ribes, 1994). This observational fact has led to intensive theoretical efforts to explain a significant asymmetry of the sun's surface magnetic field in the framework of the dynamo concept (see the review by Sokoloff, 2004, and references therein). a) According to Stuiver and Quay (1980); Stuiver and Braziunas (1989). b) According to Eddy (1977a,b). c) According to Goslar (2003). d) According to Usoskin et al. (2007). e) Exact duration is uncertain.  Table 1.

Grand minima on a multi-millennial timescale
The presence of grand minima in solar activity on the long-term scale has been mentioned numerously (e.g., Eddy, 1977a;Usoskin et al., 2003c;, using the radioisotope 14 C data in tree rings. For example, Eddy (1977b) identified major excursions in the detrended 14 C record as grand minima and maxima of solar activity and presented a list of six grand minima and five grand maxima for the last 5000 years (see Table 1). Stuiver and Braziunas (1989) and Stuiver et al. (1991) also studied grand minima as systematic excesses of the high-pass filtered 14 C data and suggested that the minima are generally of two distinct types: short minima of duration 50 -80 years (called Maunder-type) and longer minima collectively called Spörer-like minima. Using the same method of identifying grand minima as significant peaks in high-pass filtered ∆ 14 C series, Voss et al. (1996) provided a list of 29 such events for the past 8000 years. A similar analysis of bumps in the 14 C production rate was presented recently by Goslar (2003). However, such studies retained a qualitative element, since they are based on high-pass-filtered 14 C data and thus implicitly assume that 14 C variability can be divided into short-term solar variations and long-term changes attributed solely to the slowly-changing geomagnetic field. This method ignores any possible long-term changes in solar activity on timescales longer than 500 years (Voss et al., 1996). This approach, based on physics-based modelling (Section 3), allows for the quantitative reconstruction of the solar activity level in the past, and thus, for a more realistic definition of the periods of grand minima or maxima. A list of 27 grand minima, identified in the quantitative solar-activity reconstruction of the last 11,000 years, shown in Figure 17, is presented in Table 1 (after Usoskin et al., 2007). The cumulative duration of the grand minima is about 1900 years, indicating that the sun in its present evolutionary stage spends ∼ 1 / 6 (17%) of its time in a quiet state, corresponding to grand minima. Note that the definition of grand minima is quite robust.
The question of whether the occurrence of grand minima in solar activity is a regular or chaotic process is important for understanding the solar-dynamo machine. Even a simple deterministic numerical dynamo model can produce events comparable with grand minima (Brandenburg et al., 1989). Such models can also simulate a sequence of grand minima occurrences, which are irregular and seemingly chaotic (e.g., Jennings and Weiss, 1991;Tobias et al., 1995;Covas et al., 1998). The presence of long-term dynamics in the dynamo process is often explained in terms of the α-effect, which, being a result of the electromotive force averaged over turbulent vortices, can contain a fluctuating part (e.g., Hoyng, 1993;Ossendrijver et al., 1996) leading to irregularly occurring grand minima (e.g., Brandenburg and Spiegel, 2008). All these models predict that the occurrence of grand minima is a purely random "memoryless" Poisson-like process, with the probability of a grand minimum occurring being constant at any given time. This unambiguously leads to the exponential shape of the waiting-time distribution (waiting time is the time interval between subsequent events) for grand minima. Usoskin et al. (2007) performed a statistical analysis of grand-minima-occurrence time ( Table 1) and concluded that their occurrence is not a result of long-term cyclic variations, but is defined by stochastic/chaotic processes. Moreover, waiting-time distribution deviates significantly from the exponential law. This implies that the event occurrence is still random, but the probability is nonuniform in time and depends on the previous history. In the time series it is observed as a tendency of the events to cluster together with a relatively-short waiting time, while the clusters are separated by long event-free intervals (cf. Section 4.2). Such behavior can be interpreted in different ways, e.g., self-organized criticality or processes related to accumulation and release of energy. This poses a strong observational constraint on theoretical models aiming to explain the long-term evolution of solar activity (Section 4.5.1). However, because of the small number of events (27) this result is only indicative and needs further investigation.
A histogram of the duration of grand minima from Table 1 is shown in Figure 20. The mean duration is 70 year but the distribution is bimodal. The minima tend to be either of a short (30 -90 years) duration similar to the Maunder minimum, or rather long (> 100 years), similar to the Spörer minimum, in agreement with earlier conclusions (Stuiver and Braziunas, 1989). This suggests that grand minima correspond to a special state of the dynamo. Once falling into a grand minimum as a result of a stochastic/chaotic, but non-Poisson process, the dynamo is "trapped" in this state and its behavior is driven by deterministic intrinsic features.

The modern episode of active sun
We are presently living in a period of very high sun activity with a level of activity that is unprecedentedly high for the last few centuries covered by direct solar observation. The sunspot number was growing rapidly between 1900 and 1940, with more than a doubling average group sunspot number, and has remained at that high level until now (see Figure 1). Note that growth comes entirely from raising the cycle maximum amplitude, while sunspot activity always returns to a very low level around solar cycle minima. While the average group sunspot number for the period 1750 -1900 was 35±9 (39±6, if the Dalton minimum in 1797 -1828 is not counted), it stands high at the level of 75±3 since 1950. Therefore the modern active sun episode, which started in the 1940s, can be regarded as the modern grand maximum of solar activity, as opposed to a grand minimum (Wilson, 1988b).
Is such high solar activity typical or is it something extraordinary? While it is broadly agreed that the present active sun episode is a special phenomenon, the question of how (a)typical such upward bumps are from "normal" activity is a topic of hot debate.

Grand maxima on a multi-millennial timescale
The question of how often grand maxima occur and how strong they are, cannot be studied using the 400-year-long series of direct observations. An increase in solar activity around 1200 AD, also related to the Medieval temperature optimum, is sometimes qualitatively regarded as a grand maximum (Wilson, 1988b;de Meyer, 1998), but its magnitude is lower than the modern maximum (Usoskin et al., 2003c). Accordingly, it was not included in a list of grand maxima by Eddy (1977a,b).
A quantitative analysis is only possible using proxy data, especially cosmogenic isotope records. Using a physics-based analysis of solar-activity series reconstructed from 10 Be data from polar  Usoskin et al. (2003c stated that the modern maximum is unique in the last millennium. Then, using a similar analysis of the 14 C calibrated series,  found that the modern activity burst is not unique, but a very rare event, with the previous burst occurring about 8 millennia ago. An update (Usoskin et al., 2006a) of this result, using a more precise paleo-magnetic reconstruction by Korte and Constable (2005) since 5000 BC, suggests that an increase of solar activity comparable with the modern episode might have taken place around 2000 BC, i.e., around 4 millennia ago. The result by  has been disputed by Muscheler et al. (2005) who claimed that equally high (or even higher) solar-activity bursts occurred several times during the last millennium, circa 1200 AD, 1600 AD and at the end of the 19th century. We note that the latter claimed peak (ca. 1860) is not confirmed by direct solar or geomagnetic data. However, as argued by Solanki et al. (2005), the level of solar activity reconstructed by Muscheler et al. (2005) was overestimated because of an erroneous normalization to the data of ground-based ionization chambers (see also McCracken and Beer, 2007). This indicates that the definition of grand maxima is less robust than grand minima and is sensitive to other parameters such as geomagnetic field data or overall normalization. Keeping possible uncertainties in mind, let us consider a list of the largest grand maxima (the 50 year smoothed sunspot number stably exceeding 50), identified for the last 11,400 years using 14 C data, as shown in Table 2 (after Usoskin et al., 2007). A total of 19 grand maxima have been identified with a total duration of around 1030 years, suggesting that the sun spends around 10% of its time in an active state. A statistical analysis of grand-maxima-occurrence time suggests that they do not follow long-term cyclic variations, but like grand minima, are defined by stochastic/chaotic processes. The distribution of the waiting time between consecutive grand maxima is not as clear as that for grand minima, but also hints at a deviation from exponential law. The duration of grand maxima has a smooth distribution, which nearly exponentially decreases towards longer intervals. Most of the reconstructed grand maxima (about 75%) were not longer than 50 years, and only four grand minima (including the modern one) have been longer than 70 years. This suggest that the probability of the modern active-sun episode continuing is low 5 (cf. .

Related implications
Reconstructions of long-term solar activity have different implications in related areas of science. The results, discussed in this overview, can be used in such diverse research disciplines as theoretical astrophysics, solar-terrestrial studies, paleo-climatology, and even archeology and geology. We will not discuss all possible implications of long-term solar activity in great detail but only briefly mention them here.

Theoretical constrains
The basic principles of the occurrence of the 11-year Schwabe cycle are more-or-less understood in terms of the solar dynamo, which acts, in its classical form (e.g., Parker, 1955), as follows. Differential rotation Ω produces a toroidal magnetic field from a poloidal one, while the "α-effect" associated with the helicity of the velocity field produces a poloidal magnetic field from a toroidal one. This classical model results in a periodic process in the form of propagation of a toroidal field pattern in the latitudinal direction (the "butterfly diagram"). As evident from observation, the solar cycle is far from being a strictly periodic phenomenon, with essential variations in the cycle length and especially in the amplitude, varying dramatically between nearly spotless grand minima and very large values during grand maxima. The mere fact of such great variability, known from sunspot data, forced solar physicists to develop dynamo models further. Simple deterministic numerical dynamo models, developed on the basis of Parker's migratory dynamo, can simulate events, which are seemingly comparable with grand minima/maxima occurrence (e.g., Brandenburg et al., 1989). However, since variations in the solar-activity level, as deduced from cosmogenic isotopes, appear essentially nonperiodic and irregular, appropriate models have been developed to reproduce irregularly-occurring grand minima (e.g., Jennings and Weiss, 1991;Tobias et al., 1995;Covas et al., 1998). Models, including an ad hoc stochastic driver (Choudhuri, 1992;Schmitt et al., 1996;Ossendrijver, 2000;Weiss and Tobias, 2000;Mininni et al., 2001;Charbonneau, 2001;Charbonneau et al., 2004), are able to reproduce the great variability and intermittency found in the solar cycle (see the review by Charbonneau, 2005). A recent statistical result of grand minima occurrence  shows disagreement between observational data, depicting a degree of self-organization or "memory", and the above dynamo model, which predicts a pure Poisson occurrence rate for grand minima (see Section 4.3). This poses a new constraint on the dynamo theory, responsible for long-term solar-activity variations.
In general, the following additional constraints can be posed on dynamo models aiming to describe the long-term (during the past 11,000 years) evolution of solar magnetic activity.
• The sun spends about 3 / 4 of its time at moderate magnetic-activity levels, about 1 / 6 of its time in a grand minimum and about 1 / 5 − 1 / 10 in a grand maximum. Modern solar activity corresponds to a grand maximum.
• Occurrence of grand minima and maxima is not a result of long-term cyclic variations but is defined by stochastic/chaotic processes.
• Observed statistics of the occurrence of grand minima and maxima display deviation from a "memory-less" Poisson-like process, but tend to either cluster events together or produce long event-free periods. This can be interpreted in different ways, such as self-organized criticality (de Carvalho and Prado, 2000, e.g.,), a time-dependent Poisson process (Wheatland, 2003, e.g.,), or some memory in the driving process (Mega et al., 2003, e.g.,).
• Grand minima tend to be of two different types: short minima of Maunder type and long minima of Spörer type. This suggests that a grand minimum is a special state of the dynamo.
• Duration of grand maxima resemble a random Possion-like process, in contrast to grand minima.

Solar-terrestrial relations
The sun ultimately defines the climate on Earth supplying it with energy via radiation received by the terrestrial system, but the role of solar variability in climate variations is far from being clear. Solar variability can affect the Earth's environment and climate in different ways (see, e.g., a review by Haigh, 2007). Variability of total solar irradiance (TSI) measured during recent decades is known to be too small to explain observed climate variations (e.g., Foukal et al., 2006;Fröhlich, 2006). On the other hand, there are other ways solar variability may affect the climate, e.g., an unknown long-term trend in TSI (Solanki and Krivova, 2004;Wang et al., 2005) or a terrestrial amplifier of spectral irradiance variations (Shindell et al., 1999). Alternatively, an indirect mechanism also driven by solar activity, such as ionization of the atmosphere by CR (Usoskin and Kovaltsov, 2006) or the global terrestrial current system (Tinsley and Zhou, 2006) can modify atmospheric properties, in particular cloud cover (Ney, 1959;Svensmark, 1998). Even a small change in cloud cover modifies the transparency/absorption/reflectance of the atmosphere and affects the amount of absorbed solar radiation, even without changes in the solar irradiance. Since the CR flux at Earth is modulated not only by solar activity, but also by the slowly changing geomagnetic field, the two CR modulation mechanisms are independent and act on different timescales, thus giving one the opportunity to study the CR effect on Earth separately from solar irradiance (de Jager and Usoskin et al., 2005b). Accordingly, improved knowledge of the solar driver's variability may help in disentangling various effects in the very complicated system that is the terrestrial climate (e.g., de Jager, 2005;Versteegh, 2005). It is of particular importance to know the driving forces in the pre-industrial era, when all climate changes were natural. Knowledge of the natural variability can lead to an improved understanding of anthropogenic effects upon the Earth's climate.
Studies of the long-term solar-terrestrial terrestrial are mostly phenomenological, lacking a clear quantitative physical mechanism. Even phenomenological and empirical studies suffer from large uncertainties, related to the quantitative interpretation of proxy data, temporal and spatial resolution (Versteegh, 2005). Therefore, more precise knowledge of past solar activity, especially since it is accompanied by continuous efforts of the paleo-climatic community on improving climatic data sets, is crucial for improved understanding of the natural (including solar) variability of the terrestrial environment.

Other issues
The proxy method of solar-activity reconstruction, based on cosmogenic isotopes, was developed from the radiocarbon dating method, when it was recognized that the production rate of 14 C is not constant and may vary in time due to solar variability and geomagnetic field changes. Neglect of these effects can lead to inaccurate radiocarbon (or more generally, cosmogenic nuclide) dating, which is a key for, e.g., archeology and Quaternary geology. Thus, knowledge of past solar activity and geomagnetic changes allows for the improvement of the quality of calibration curves, such as the IntCal (Stuiver et al., 1998;Reimer et al., 2004) for radiocarbon, eventually leading to more precise dating.
Long-term variations in the geomagnetic field are often evaluated using cosmogenic isotope data. Knowledge of source variability due to solar modulation is important for better results.

Summary
In this section, solar activity on a longer scale is discussed, based on recent reconstructions.
According to these reconstructions, the sun has spent about 70% of its time during the Holocene, which is ongoing, in a normal state characterized by medium solar activity. About 15 -20% of the time the sun has experienced a grand minimum, while 10 -15% of the time has been taken up by periods of very high activity.
One of the main features of long-term solar activity is its irregular behavior, which cannot be described by a combination of quasi-periodic processes as it includes an essentially random component.
Grand minima, whose typical representative is the Maunder minimum of the late 17th century, are typical solar phenomena. A total of 27 grand minima have been identified in reconstructions of the Holocene period. Their occurrence suggests that they appear not periodically, but rather as the result of a chaotic process within clusters separated by 2000 -2500 years. Grand minima tend to be of two distinct types: short (Maunder-like) and longer (Spörer-like). The appearance of grand minima can be reproduced by modern stochastic-driven dynamo models to some extent, but some problems still remain to be resolved.
The modern level of solar activity (after the 1940s) is very high, corresponding to a grand maximum, which are typical but rare and irregularly-spaced events in solar behavior. The duration of grand maxima resembles a random Possion-like process, in contrast to grand minima.
These observational features of the long-term behavior of solar activity have important implications, especially for the development of theoretical solar-dynamo models and for solar-terrestrial studies.

Solar Energetic Particles in the Past
In addition to galactic cosmic rays, which are always present in the Earth's vicinity, sometimes sporadic solar energetic particle (SEP) events with a greatly enhanced flux of less energetic particles in the interplanetary medium also occur (e.g., Klecker et al., 2006). Strong SEP events mostly originate from CME-related shocks propagating in the solar corona and interplanetary medium, that lead to effective bulk acceleration of charged particles (e.g., Cane and Lario, 2006). Although these particles are significantly less energetic than GCRs, they can occasionally be accelerated to an energy reaching up to several GeV, which is enough to initiate the atmospheric cascade. Peak intensity of SEP flux can be very high, up to 10 4 particles (with energy > 30 MeV) per cm 2 per sec. In fact, the long-term average flux (or fluence) of SEP is mostly defined by rare major events, which occur one to two times per solar cycle, with only minor contributions from a large number of weak events Smart, 1990, 2002). As an example, energy spectra of GCR and SEP are shown in Figure 21 for the day of January 20, 2005, when an extreme SEP event took place. Such SEPs dominate the low-energy section of cosmic rays (below hundreds of MeV of a particle's kinetic energy), which is crucial for the radiation environment, and play an important role in solar-terrestrial relations. For many reasons it is important to know the variations of SEPs on long-term scales.
It is not straightforward to evaluate the average SEP flux even for the modern instrumental epoch of direct space-borne measurements (e.g., Mewaldt et al., 2007). For example, estimates for the average flux of SEPs with an energy above 30 MeV (called F 30 henceforth) for individual cycles may vary by an order of magnitude, from 8 cm −2 s −1 for cycle 21 up to 60 cm −2 s −1 for cycle 19 . Moreover, estimates of the SEP flux were quite uncertain during the earlier years of space-borne measurements because of two effects, which are hard to account for (e.g., Reeves et al., 1992;Tylka et al., 1997). One is related to the very high flux intensities of SEPs during the peak phase of events, when a detector can be saturated because of the dead-time effect (the maximum trigger rate of the detector is exceeded). The other is related to events with high energy solar particles, which can penetrate into the detector through the walls of the collimator or the detector, leading to an enhanced effective acceptance cone with respect to the "expected" one. Since the SEP fluence is defined by major events, these effects may lead to an underestimate of the average flux of SEPs. The modern generation of detectors are better suited for measuring high fluxes. The average F 30 flux for the last five solar cycles (1954 -2006) is estimated at about 35 cm −2 s −1 Shea et al., 2006).
The method of cosmogenic isotopes in terrestrial archives does not allow quantitative estimates of SEP events because of atmospheric and geomagnetic rigidity/energy cutoffs, which are too high for less energetic particles of solar/interplanetary origin. While the SEP effect is negligible in 14 C records, the 10 Be record may contain information on extreme SEP events in the past . In particular, an additional SEP-related enhancement of the 10 Be signal can appear around some solar activity maxima leading to an intermittent 5.5-year quasi-periodicity (McCracken et al., 2002). However, such peaks can be studied only a posteriori, i.e., for known SEP events (e.g., the Carrington flare), which can be associated with strong peaks in 10 Be, but the relation is not one-to-one. Therefore, SEP events cannot be unambiguously identified from the 10 Be data.
Since energy spectra of the two populations of cosmic rays are significantly different (very soft for SEP and much harder for GCR), they can be separated by their energy. Accordingly, for earlier times, when direct measurements were not possible, one needs a natural spectrometer in order to distinguish between cosmic particles of solar and galactic origin. Such a spectrometer should be able to naturally separate lower and higher energy components of the cosmic-ray spectrum and "record" them in different archives. Fortunately, there are such natural spectrometers, which allow us to study not only galactic cosmic rays but also solar energetic particle flux in the past.  (Mewaldt, 2006).

Lunar and meteoritic rocks
One spectrometer that is able to separate cosmic rays is lunar (or meteoritic) rocks. Figure 22 depicts an example of 14 C measured in a lunar sample (Jull et al., 1998). The dotted line shows the expected production of radiocarbon by GCR. The production increases with depth due to the development of a nucleonic cascade in the matter, initiated by energetic GCR particles, similar to the atmospheric cascade. Less energetic particles of solar origin produce the isotope only in upper layers of the rock, since their low energy does not allow them to initiate a cascade. On the other hand, thanks to their high flux in the lower energy range, the production of 14 C in the upper layers is much higher than that from GCR. Thus, by first measuring the isotope activity in deep layers one can evaluate the average GCR flux, and then the measured excess in the upper level yields an estimate for the SEP flux in both integral intensity and spectral shape. The result is based on model computations and therefore is slightly model dependent but makes it possible to give a robust estimate of the GCR and SEP in the past.
A disadvantage of this approach is that lunar samples are not stratified and do not allow for temporal separation. The measured isotope activity is a balance between production and decay and, therefore, represents the production (and the ensuing flux) integrated over the life-time of the isotope before the sample has been measured. However, using different isotopes with different life times, one can evaluate the cosmic-ray flux integrated over different timescales.
Estimates of the average SEP flux F 30 on different timescales, as obtained from various isotopes measured in lunar samples, are collected in Table 3. The average F 30 flux for the last five solar cycles (1954 -2006) is consistent with the average flux estimated in the past for longer timescales from 10 3 to 10 7 years (cf. Reedy, 2002).

Nitrates in polar ice
Another natural spectrometer is the Earth's atmosphere, where the penetration depth of cosmic rays depends on their energy. As an illustration, Figure 23 shows ionization due to cosmic rays in the polar atmosphere. One can see that ionization due to GCR is present in the entire atmosphere  Figure 22: Measured (dots) and calculated (curves) 14 C activity in a lunar sample 68815 (Jull et al., 1998). The big diamond implies contamination of a thin surface layer by 14 C implanted from solar wind. The dotted curve represents the expected production due to GCR, while the solid curve is the best fit SEP+GCR model production.  Figure 23: Altitude profile of the cosmic-ray-induced ionization (CRII) rate in the polar atmosphere, computed using the model by Usoskin and Kovaltsov (2006). Dashed and dotted curves depict the average CRII for solar maxima and minima, respectively. The solid curve denotes the average CRII due to solar energetic particles for the day of January 20, 2005.
with a minor solar-cycle variation. The ionization effect of a SEP event is minor in the troposphere even for the extreme event of January 20, 2005. On the other hand, the ionization rate is greatly enhanced (by orders of magnitude) in the polar stratosphere during strong SEP events (Bütikofer et al., 2008). Although such enhanced stratospheric ionization cannot be directly recorded in natural archives to be measured later, there is an indirect way to restore the history of extreme SEP events in the past. Cosmic-ray-induced ionization can lead to essential chemical changes in the polar stratosphere with enhanced production of "odd nitrogens" NO y (e.g., N, NO, NO 2 , NO 3 , 2N 2 O 5 , BrONO 2 , ClONO 2 , HO 2 NO 2 , and HNO 3 ) (see, e.g., Jackman et al., 1990Jackman et al., , 1993Vitt et al., 2000). Due to ionization of the ambient air and subsequent dissociation of O 2 and N 2 , energetic particles precipitating into the atmosphere facilitate formation of the odd nitrogen. The abundance of odd nitrogen dramatically responds to variations in the SEP flux on the background of a smooth 11-yr cycle due to GCR variations (Thomas et al., 2007). Odd nitrogen is long-lived during the polar night, and some species, such as nitric acid HNO 3 , can be effectively transported down and finally stored in polar snow/ice. By measuring such chemical species in a stratified polar-ice archive, one can obtain a record of the physical-chemical conditions of the polar atmosphere (Zeller et al., 1986;Zeller and Dreschhoff, 1995). Under some basic reasonable assumptions (stable atmospheric properties, no mixing in snow, etc.), the abundance of compounds, related to the stratospheric chemistry such as nitrate ions (NO − 3 ), provides a unique and long-term (potentially up to 10 5 years) record of the radiation environment and climate in polar atmosphere. The concentration of nitrates has been measured in polar ice from both the Southern (South Pole, e.g., Dreschhoff and Zeller, 1990) and Northern (Greenland, e.g., Zeller and Dreschhoff, 1995;Dreschhoff and Zeller, 1998) polar caps. The Greenland series of ultra-high resolution (better than 30 equispaced samples per year) provides a unique record in which (groups of) individual strong SEP events can be identified for as far back as 1562.
The nitrate abundance measured in polar ice depicts a clear seasonal cycle with a peak during the local summer because of increased snow sublimation under sunlight (Zeller et al., 1986). More- Figure 24: Cumulative probability of a large solar energetic particle event to occur (after Mc-Cracken et al., 2001b). The black histogram, extrapolated from the blue line corresponds to the directly observed SEP events (Reedy, 1996). Arrows (extrapolated from the brown line) depict an upper limit obtained from the analysis of lunar rock (see Table 3) assuming that the entire fluence has been produced within a few extreme events. Diamonds represent the result derived from the nitrate data for 1561-1950(McCracken et al., 2001b. over, the nitrate series contains clear time stamps of major volcano eruptions, which are apparent in the measured signal. Altogether it provides a solid basis for absolute dating of the samples with a time resolution within one year (McCracken et al., 2001b). Thus, the nitrate concentration in well-dated polar-ice cores provides a unique opportunity to evaluate the flux of SEP in the past, before instrumental observations. High time resolution allows the separating of (groups of) individual SEP events, such as the famous Carrington event in September 1859 (Shea et al., 2006;Thomas et al., 2007). It has been suggested that nitrate enhancements caused by strong SEP events can be reliably distinguished both from the relatively slow and shallow variations of galactic cosmic rays and from meteorologically-derived nitrate peaks.
The first studies based on the correlation between peaks of nitrate concentration and measures of solar activity were qualitative and aimed to search for periods of active/quiet sun (e.g., Zeller et al., 1986;Dreschhoff and Zeller, 1998) or even for a supernova event (Dreschhoff and Laird, 2006). Later, the method was explored more fully (McCracken et al., 2001a,b), which made it possible to identify large SEP events since 1560 and evaluate their fluence (see Table 1 in McCracken et al., 2001b). Only events with F 30 fluence exceeding 10 9 particles per cm 2 , which is estimated as a threshold for the nitrate signal, are identified in this way. This analysis leads to a statistical estimate of the occurrence frequency of large SEP events ( Figure 24). Thus-obtained statistics of SEP-event occurrence complements directly observed data and that reconstructed from lunar samples, indicating a kind of break in the distribution with the fluence above 10 10 protons/cm 2 . Note that estimates based on lunar-rock samples (arrows in Figure 24) provide only an upper bound for SEP fluence, since they are based on an extreme assumption that the net fluence was produced in a few extreme events. Therefore, extra-strong SEP events with the fluence exceeding 10 12 cm −2 are very improbable.
However, the method has a drawback as noted by McCracken et al. (2001b). The measured nitrate effect depends not only on the total fluence F 30 but also on the energy spectrum of SEP, namely, high energy particles produce more ionization and thus lead to more nitrate effect than lower-energy particles. Accordingly, the measured nitrate signal can be converted into the SEP fluence only by assuming a fixed energy spectrum of SEP by normalizing all data to, e.g., the event of August 1972(McCracken et al., 2001bShea et al., 2006). However, the actual energy spectrum of SEP remains a free parameter in such a reconstruction . Moreover, the nitrate record is only an indirect proxy of strong SEP events whose efficiency is not 100% (Palmer et al., 2001).

Summary
In this section, estimates of the averaged long-term flux of solar energetic particles (SEPs) are discussed.
Measurements of cosmogenic isotopes with different life times in lunar and meteoritic rocks allow one to make rough estimates of the SEP flux over different timescales. The directly spaceborne-measured SEP flux for past decades is broadly consistent with estimates on longer timescales -up to millions of years.
Measurements of nitrates in polar ice make it possible to reconstruct strong SEP events for nearly the past five centuries.
An analysis of various kinds of data suggests that the distribution of the intensity of SEP events has taken a break, and the occurrence of extra-strong events (with the F 30 fluence exceeding 10 12 ) is unlikely.

Conclusions
In this review the present knowledge of long-term solar activity on a multi-millennial timescale, as reconstructed using the indirect proxy method, is discussed.
Although the concept of solar activity is intuitively understandable as a deviation from the "quiet" sun concept, there is no clear definition for it, and different indices have been proposed to quantify different aspects of variable solar activity. One of the most common and practical indices is sunspot number, which forms the longest available series of direct scientific observations. While all other indices have a high correlation with sunspot numbers, dominated by the 11-year cycle, the relationship between them at other timescales (short and long-term trends) may vary to a great extent.
On longer timescales, quantitative information of past solar activity can only be obtained using the method based upon indirect proxy, i.e., quantitative parameters, which can be measured nowadays but represent the signatures, stored in natural archives, of the different effects of solar magnetic activity in the past. Such traceable signatures can be related to nuclear or chemical effects caused by cosmic rays in the Earth's atmosphere, lunar rocks or meteorites. The most common proxy of solar activity is formed by data from the cosmogenic radionuclides, 10 Be and 14 C, produced by cosmic rays in the Earth's atmosphere and stored in independently-dated stratified natural archives, such as tree rings or ice cores. Using a recently-developed physics-based model it is now possible to reconstruct the temporal behavior of solar activity in the past, over many millennia. The most robust results can be obtained for the Holocene epoch, which started more than 11,000 years ago, whose stable climate minimizes possible uncertainties in the reconstruction. An indirect verification of long-term solar-activity reconstructions supports their veracity and confirms that variations of cosmogenic nuclides on the long-term scale (centuries to millennia) during the Holocene make a solid basis for studies of solar variability in the past. However, such reconstructions may still contain systematic uncertainties related to unknown changes in the geomagnetic field or climate of the past, especially in the early part of the Holocene.
Measurements of nitrates in polar ice allow the reconstruction of strong solar energetic particle (SEP) events in the past, over the five past centuries. Together with independent measurements of the concentration of different cosmogenic isotopes in lunar and meteoritic rocks, it leads to estimates of the SEP flux on different timescales. Directly space-borne-measured SEP flux for recent decades is broadly consistent with estimates on longer timescales -up to millions of years, and the occurrence of extra-strong events is unlikely.
In general, the following main features are observed in the long-term evolution of solar magnetic activity.
• Solar activity is dominated by the 11-year Schwabe cycle on an interannual timescale. Some additional longer characteristic times can be found, including the Gleissberg secular cycle, de Vries/Suess cycle, and a quasi-cycle of 2000 -2400 years. However, all these longer cycles are intermittent and cannot be regarded as strict phase-locked periodicities.
• One of the main features of long-term solar activity is that it contains an essential chaotic/ stochastic component, which leads to irregular variations and makes solar-activity predictions impossible for a scale exceeding one solar cycle.
• The sun spends about 70% of its time at moderate magnetic activity levels, about 15 -20% of its time in a grand minimum and about 10 -15% in a grand maximum. Modern solar activity corresponds to a grand maximum.
• Grand minima are a typical but rare phenomena in solar behavior. Their occurrence appears not periodically, but rather as the result of a chaotic process within clusters separated by 2000 -2500 years. Grand minima tend to be of two distinct types: short (Maunder-like) and longer (Spörer-like).
• The modern level of solar activity (after the 1940s) is very high, corresponding to a grand maximum. Grand maxima are also rare and irregularly occurring events, though the exact rate of their occurrence is still a subject of debates.
These observational features of the long-term behavior of solar activity have important implications, especially for the development of theoretical solar-dynamo models and for solar-terrestrial studies.

Acknowledgements
I thank the editorial board of the Living Reviews in Solar Physics for the invitation to prepare this review and organizational aid. This work would be impossible without encouraging discussions, as well as direct and indirect help, from my colleagues Gennady Kovaltsov and Sami Solanki, whose invaluable support is acknowledged with my greatest gratitude. I have also greatly benefited from lively and open discussions and debates with Edouard Bard, Jürg Beer, Narendra Bhandari, late Giuliana Cini Castagnoli, Paul Charbonneau, Vincent Courtillot, Cornelis de Jager, Gisella Dreschhoff, Erwin Flückiger, Peter Foukal, Yves Gallet, Agnés Genevey, David Hathaway, Monika Korte, Natalia Krivova, Bernd Kromer, Devendra Lal, Michael Lockwood, Ken McCracken, David Moss, Kalevi Mursula, Raimund Muscheler, Hiroko Miyahara, Heikki Nevanlinna, Valery Ostryakov, Alexander Ruzmaikin, Dieter Schmitt, Manfred Schüssler, Margaret Shea, Don Smart, Ian Snowball, Dmitry Sokoloff, Willie Soon, Carla Taricco, Jóse Vaquero, and many others. I thank Evguenia Usoskina for help with editing this text. I am happy to acknowledge those colleagues who work in this field and make it living and vibrant.