Consistency of the average ﬂux of solar energetic particles over timescales of years to megayears

Aims. Solar energetic particles (SEPs) have been measured directly in space over the past decades. Rare extreme SEP events are studied based on terrestrial cosmogenic proxy data for the past ten millennia. Lunar rocks record the average SEP ﬂuxes on the megayear timescale. The question of whether the SEP ﬂuxes averaged over di ﬀ erent timescales are mutually consistent is still open. Here we analyze these di ﬀ erent datasets for mutual consistency. Methods. Using the data from directly measured SEPs over the past decades and reconstructions of extreme SEP events in the past, we built a distribution function of the occurrence of annual SEP ﬂuences for SEPs with energies above 30, 60, 100, and 200MeV. The distribution function was ﬁt with the Weibull and other types of distributions, and the long-term average SEP ﬂux was computed and compared with the megayear SEP ﬂux estimated from lunar data. Results. In contrast to the current paradigm, the direct space-era data are not representative of the long-term averaged SEP ﬂux because they are only 20–55% of it, while the major fraction was formed by rare extreme SEP events in the past. The combined statistics of direct and proxy data are fully consistent with megayear lunar data, implying that our knowledge of the whole range of the SEP ﬂuxes, from frequent weak to rare extreme events, is now consistent


Introduction
The majority of charged relativistic particles in the vicinity of Earth, known as cosmic rays, are galactic cosmic rays (GCRs), which are always present and have a moderate solar-modulation variability (e.g., Potgieter 2013; Aguilar et al. 2021). However, the lower-energy ( 100 MeV) part of the cosmic-ray flux is dominated by solar energetic particles (SEPs), which are accelerated in solar eruptive processes such as flares and coronal mass ejections (e.g., Desai & Giacalone 2016;Vainio et al. 2009). Although the energy of SEPs is lower than that of GCRs, their fluxes are much stronger and are highly variable in time. SEPs determine the near-Earth radiation environment, which is crucial for modern technology (e.g., Dyer et al. 2004), and they pose hazards for spacecraft, aircraft, and even habitability (Thomas et al. 2013;Yamashiki et al. 2019). Moreover, SEPs form the only measurable index of solar eruptive events in the past and thus can expand the range of our knowledge on solar physics from tens to millions of years (Reedy 1999;Usoskin et al. 2006;Cliver et al. 2022).
Three independent ways are currently used to study SEP fluxes on different timescales. One way (see Sect. 2.1) is based on direct measurements of SEPs in space (e.g., Shea & Smart 2000;Reedy 2012;Vainio et al. 2013) by various spacecraft since the 1960s, but these estimates earlier suffered from large uncertainties related to the difficult intercalibration of various space-borne detectors. Primary SEP detectors on board the Geostationary Operational Environmental Satellite (GOES) family have recently been recalibrated (Raukunen et al. 2022), which reduced the uncertainties of the SEP flux measurements since 1984.
Another way (Sect. 2.2) to study past SEP fluxes is by cosmogenic isotopes (CIs) 14 C, 10 Be, and 36 Cl produced by energetic particles in the Earth's atmosphere, which are stored in natural stratified archives such as tree trunks or polar ice and are measured in modern laboratories (e.g., Usoskin et al. 2006;Beer et al. 2012;Usoskin 2017). Normally, CI production is dominated by GCR fluxes, and the SEP contribution cannot be resolved (Usoskin et al. 2020;Mekhaldi et al. 2021). However, as discovered recently (e.g., Miyake et al. 2012Miyake et al. , 2013Usoskin et al. 2013;Paleari et al. 2022;Brehm et al. 2022), our Sun can occasionally produce extreme solar particle events (ESPEs) that lead to an enhanced (several orders of magnitude greater than for SEP events that were directly measured during the space era) CI production notably exceeding that by GCRs on the annual scale. The average SEP flux estimated from the direct data alone does not contain ESPEs because they did not occur during the space era.
The third way (Sect. 2.3) to assess the average SEP flux on a very long timescale is related to CIs measured in lunar samples brought to Earth by the Apollo missions in 1969-1972(e.g., Fink et al. 1998Jull et al. 1998;Reedy 1999;Nishiizumi et al. 2009;Poluianov et al. 2018). Since the Moon is not protected from energetic particles by the atmosphere and magnetosphere, 1.9 +0.9 −0.7 -Notes. Columns are: (1) flux notation; (2) average SEP flux for SC 22-24 (1987-2019) as taken from Table D1 of Raukunen et al. (2022); (3) the same but for SC 24 only (2009-2019); (4) F ∞ (see Eq. (2)) for the combined GOES and CI-based data (green curves in Fig. 1); (5) F * ∞ -the same as (3) but with GOES data, for SC 24 only; (6) average SEP flux reconstructed from lunar rocks for timescales of ≈10 6 years (Poluianov et al. 2018). Error bars are full-range for the lunar data and 1σ for other datasets. Error bars for Cols. 2 and 3 are <1% and not shown. The data are shown in Fig. 3.
CIs in lunar rocks can serve as a rough cosmic-ray spectrometer (Miyake et al. 2019;Poluianov et al. 2018) over the isotope lifetime, which may reach several megayears. The GCR contribution to lunar CI production can be estimated using the data from deep layers of lunar soil cores (Rancitelli et al. 1975;Nishiizumi et al. 1984;Poluianov et al. 2018), which enables the quantitative estimation of the lower-energy (30-80 MeV) average SEP flux using the data of the shallower lunar soil layers. This method cannot resolve individual SEP events (including ESPEs) nor the time variability, but it provides a robust estimate of the long-term average SEP flux. All these three methods assess the average SEP flux on different timescales and in independent ways. The current paradigm is that the megayear-averaged spectrum of SEPs is consistent with that over the past decades (e.g., Reedy 2012;Poluianov et al. 2018). In this Letter, we critically revise the existing paradigm and address the outstanding question of whether the SEP fluxes measured or derived by different methods on different timescales agree with each other and whether we can obtain a consistent view of the average flux of SEPs produced by the Sun in its current evolution state. This question is crucially important for studies in the fields of solar and stellar physics, terrestrial and planetary sciences, astrobiology, space weather, and modern technology.

Datasets
Solar energetic particles have traditionally been quantified via the integral particle flux F(> E) (in units of cm −2 s −1 ) above the given energy E. This flux definition differs from the standard definition, and should strictly speaking be called the omnidirectional intensity (see Chap. 1.6.3 in Grieder 2001), but it has been used in the field historically. Typical energies E for the measured integral flux of SEPs are 10, 30, 60, 100, and 200 MeV (Shea & Smart 2000;Raukunen et al. 2022), with the corresponding integral fluxes denoted henceforth as F 10 through F 200 , respectively.

Direct measurements
We used the annually averaged SEP integral fluxes based on the GOES dataset (Raukunen et al. 2022) in the energy range from 30-200 MeV for 1987-2019 (solar cycles (SCs) 22-24), Table 2. ESPE integral fluxes reconstructed from terrestrial CI data (Mekhaldi et al. 2015;Brehm et al. 2022;Paleari et al. 2022 Notes. The years of the events are given in the first column. Error bars mark the 68% confidence intervals. shown in Table 1. This dataset is an estimate of the regular SEP flux, including the occurrence of weak-to-strong SEP events, but excluding ESPEs. The period between the onset of the space era and the end of SC 23 in 2009 corresponds to the Modern Grand solar maximum (Usoskin et al. 2003;Solanki et al. 2004) of an unusually active Sun, with the cycle-averaged SEP flux in SC 22-23 about five times stronger than during the moderate SC 24 (Table 1). This is typical for the long-term solar activity (Usoskin 2017).

ESPEs in terrestrial CIs
Extreme solar particle events were discovered and are studied based on CIs that were recorded in terrestrial archives such as tree trunks or the polar ice (e.g., Usoskin et al. 2006Usoskin et al. , 2020Cliver et al. 2022); see Appendix A. We reconstructed integral ESPE fluxes using the recently developed effective-energy method (Kovaltsov et al. 2014;Koldobskiy et al. 2022), which estimates the SEP fluxes F 234 , F 236 , and F 60 based on the global production of 14 C (Brehm et al. 2022) and based on the production and deposition of 10 Be and 36 Cl in the polar caps (Mekhaldi et al. 2015), respectively. A short description of the method is presented in Appendix B. For the ESPE dated 7176 before the common era (BCE), the SEP flux values were taken from Paleari et al. (2022). The obtained ESPE flux values are summarized in Table 2.
Since the CI method detects only ESPEs but no regular SEP events (Usoskin et al. 2020;Mekhaldi et al. 2021), it cannot evaluate the total average SEP flux. Thus, this information is complementary to the regular SEP flux assessed from direct data (Sect. 2.1).

Megayear SEP flux estimates from lunar samples
In contrast to terrestrial proxies, lunar CI data lack a temporal resolution and thus cannot resolve individual ESPEs, but work as a simple particle spectrometer providing an integral cosmicray spectrum averaged over the CI lifetime (see Appendix A). The SEP flux averaged over long timescales in the energy range 30-80 MeV was recently estimated, as shown in Table 1, using data of cosmogenic 26 Al (lifetime 1.03 Myr) measured in lunar rocks (Poluianov et al. 2018), without an a priori assumption on the spectral form (see the brief description in Appendix C). The lunar-CI-based SEP fluxes should include both the regular (Sect. 2.1) and ESPE components (Sect. 2.2) because their L22, page 2 of 9 long-time averaging guarantees that rare ESPEs are covered, potentially even those with an occurrence rate <10 −4 year −1 .

Methods
Terrestrial CI data have an annual time resolution at best. To compare the statistics of proxy-based and direct datasets, we averaged the directly measured data to annual fluxes. Thus, possible recurrent SEP events produced by the same solar active regions are not distinguished and are seen as one annual event. We used GOES data for 33 years (Sect. 2.1) and CI data for eight ESPEs during 12 000 years of the Holocene (Sect. 2.2) to calculate the complementary cumulative distribution function (CCDF), that is, the probability that the average flux for a randomly chosen year exceeds the given value. For each annual flux data point, we estimated (following the method of Usoskin & Kovaltsov 2012) the median value of the probability W and its 68% confidence interval (c.i.) σ W . Since the number of ESPE flux reconstructions F 30 -F 100 is limited to three events (Table 2), their occurrence probability was assessed using F 200 data with a larger statistic. Figure 1 shows statistics of the annual SEP integral fluxes F 30 , F 60 , F 100 , and F 200 in the CCDF form, where data of regular SEP events and ESPEs are depicted by circles and stars, respectively. Conservative upper limits for GOES and terrestrial CI data are shown by filled symbols. They are defined under the assumption that a flux that is twice as high as the highest observed flux has confidently not appeared during the given period (33 years of the space era, or 12 000 years of the Holocene).
The CCDF of the SEP event occurrence is often described by the Weibull distribution (Weibull 1951;Gopalswamy 2018), where W is the probability of the flux F (e.g., F 200 ) for a randomly chosen year to exceed the given value F, k (dimensionless), and λ (in cm −2 s −1 ) are the distribution parameters. We fit the experimental data ( Fig. 1) with the Weibull distribution using the iterative Monte Carlo procedure. The procedure is described in Appendix D.
The best-fit parameters and their uncertainties are shown in Table 3. While the scaling parameter λ has large uncertainties and varies essentially with the SEP energy, the shape of the distribution (parameter k) is robustly defined and stable over all energies. This suggests that ESPEs are driven by physical processes that scale with the event strength (flux). The obtained parameters are interrelated, as discussed in Appendix D. The fit also included GOES-based data limited to the moderate cycle SC 24 as presented in the right-hand side (RHS) block, Cols. 4-5 of Table 3.
To validate the fit, we also tried fitting the obtained CCDF ( Fig. 1) with other distribution types. Neither purely exponential nor power-law distributions can formally fit the data. The Band function (Band et al. 1993) can provide a poor formal fit to the CCDF, but because it uses four free parameters compared to the two parameters in the Weibull distribution, it is rejected by the Akaike and Bayes information criteria (Chakrabarti & Ghosh L22, page 3 of 9

Results and discussion
The combined statistics of the direct and terrestrial-CI proxy data are nearly perfectly fit by the Weibull distribution (Fig. 1). Knowing (or estimating) the CCDF of the integral SEP flux, we can compute the full average flux (i.e., averaged over an infinite exposure time) as where ω(F) is the differential Weibull distribution. However, if the dataset is limited so that the highest observed flux is F , as is the case for the space era (Fig. 1), for example, the formal average flux F F (replacing ∞ with F in Eq. (2)) might be significantly underestimated compared to the true value because the contribution from rare strong-flux events is missed. The ratio F F / F ∞ as a function of F is shown in Fig. 2 for F 60 using the Weibull distribution given in the RHS block of Table 3. Direct space-era observations are limited to weakmoderate SEP fluxes (as marked in Fig. 2) and represent only ≈40% of the total flux F ∞ . For other energies, this fraction ranges from 20% to 55% (up to 71% for F 30 , caused by larger ESPE flux uncertainties in this energy range). However, the inclusion of the ESPEs data accounts for 95% of F ∞ . Thus, in contrast to the current paradigm, the directly measured SEP  Table 1. flux during the space era is not representative of the long-term flux because a great part of the ESPE contribution is missing. Next, we compare in Fig. 3 the F ∞ values estimated from the GOES and terrestrial CI data with the megayearaveraged SEP fluxes assessed independently from the lunar data (Poluianov et al. 2018) for the energy range 30-200 MeV. The average SEP fluxes that were directly measured during the space era are systematically lower than the megayear-averaged lunar data. This implies a significant discrepancy, particularly for the moderate SC 24, which represents the typical long-term solar activity (Usoskin et al. , 2016. However, the average fluxes F ∞ based on the combined space-era and CI-proxy dataset are fairly consistent with the lunar-based data as the shaded areas overlap (see also Table 1). Since the average SEP flux was high during SCs 22-23 because of the Modern Grand maximum, we also calculated F * ∞ , which is similar to F ∞ , but based on direct observations during SC 24 alone. The agreement between F * ∞ and the SEP fluxes from lunar data is very good, implying that the combined statistic of SEPs during a moderate SC and rare ESPEs is fully consistent with the megayear-averaged SEP spectrum.

Conclusions
In this Letter, we have analyzed average SEP fluxes in the energy range 30-200 MeV, estimated with different methods and datasets for different timescales: direct spacecraft measurements over the past few decades, ESPEs reconstructed from terrestrial CI-proxy data on the multi-millennial timescale, and lunar CI data over a megayear timescale. We conclude that direct data for the past decades are not representative of the longterm SEP flux. They contribute only 20-55% to it. Rare ESPEs, which were absent during the space era but were reconstructed with the terrestrial CI proxy over the Holocene, contribute 40-80% to the long-term SEP flux. This is a major contribution. The combined statistics of the direct and proxy-based SEP fluxes are fully consistent with the megayear-averaged SEP fluxes that were reconstructed independently from lunar samples.
Thus, we showed for the first time to our knowledge that two different components of the SEP flux, that is, the regular cycle variability and rare extreme events, contribute significantly to the long-term SEP flux, as shown by a good agreement with the lunar data. Neither of these datasets alone can represent the SEP flux variability, which is needed, for instance, to assess the radiation environment for planetary evolution and habitability and for the planning of long-lasting space missions.

Appendix A: Peculiarities of energetic particle measurements with terrestrial and lunar cosmogenic isotopes
Here, we provide an oversimplified illustration ( Figure A.1) of the formation of cosmogenic-isotope records in terrestrial and lunar natural archives. More details can be found elsewhere (Beer et al. 2012;Usoskin 2017). Fig. A.1a illustrates the formation of terrestrial isotopes exemplified in 14 C in tree rings. We first consider, in the lefthand side of the plot, two cosmic-ray particles impinging on the Earth's atmosphere at time t 1 : CR 1 with energy E 1 and CR 2 with energy E 2 < E 1 , as depicted by the blue dots. An energetic CR 1 particle can initiate a fully developed nucleonic cascade in the atmosphere and produce 14 C atoms near the ground. This cosmogenic 14 C becomes absorbed by living trees and is used to build a tree ring, corresponding to time t 1 (middle panel). Less energetic CR 2 particle can induce only a partial cascade and produce 14 C atoms in the middle stratosphere. These 14 C atoms are transported by the atmospheric circulation and can also eventually be absorbed by trees in the same annual ring. We now consider a similar pair of cosmic-ray particles, but a pair that impinges on Earth at time t 2 (as illustrated in the RHS of the plot), which is some years later. The 14 C atoms produced by CR 1 and CR 2 particles are stored in a way similar to that shown on the left-hand side, in a tree ring corresponding to the time t 2 . Thus, terrestrial cosmogenic 14 C provides a temporal resolution (different tree rings), but no spectrometry for CRs. If there were an enhanced flux of cosmic rays during a specific year, it would have been observed as a peak in the 14 C concentration in the corresponding tree ring. This allows us to study ESPEs and the times of their occurrence using terrestrial cosmogenic-isotope proxies. The scheme is similar for other isotopes, for example, 10 Be or 36 Cl, which are measured in stratified and independently dateable polar ice cores. The use of multiproxy data enables a rough estimate of the spectra because different isotopes have different energy thresholds for their production.
A sketch of the formation of the lunar cosmogenic isotope record is shown in Figure A.1b. We consider the same pair of cosmic-ray particles, that is, CR 1 and CR 2 , that impinges on the lunar surface at time t 1 (left-hand side). The more energetic CR 1 initiates a nucleonic cascade in the soil or rock, eventually producing an atom of 26 Al at some depth in situ. After production, the 26 Al atom remains exactly at the same depth at which it had been produced because there is no transport inside the rock. The less energetic CR 2 produces a 26 Al atom at a shallower depth, where it remains since then. Similarly, the CR 1 -CR 2 pair arriving at the moment t 2 (RHS) produces 26 Al atoms at the same depths, depending only on the cosmic-ray energy, regardless of the arrival time. The measured concentration of 26 Al is defined by the balance between production and decay rates. Depth in the soil or rock roughly represents the energy of the incoming particles, thus serving as an energy-integrating spectrometer. Thus, the lunar isotope provides an energy resolution, but no time resolution. Accordingly, we can robustly assess the long-term average spectrum of cosmic rays without resolving individual events.
L22, page 6 of 9 Primary cosmic-ray particles CR 1 and CR 2 with energies E 1 and E 2 < E 1 are denoted by blue dots, nuclear reactions by stars, secondaries by arrows, and produced cosmogenic isotopes by red dots. Two different moments of time t 1 and t 2 > t 1 are shown on the left-and right-hand sides, respectively.

Appendix B: Reconstructing the integral flux of ESPEs
To estimate the ESPE integral flux, we used the effective-energy method (also known as the bow-tie method) based on an assumption that there is an energy value, called effective energy E eff , that causes the CI response Q (i.e., the global production of 14 C or polar deposition of 10 Be and 36 Cl) for SEP events to be directly proportional to the integral flux of SEP particles with an energy that exceeds the effective energy, for a reasonably broad range of spectral shapes and parameters. The effective-energy method was validated, and the values of E eff and K eff were determined for different CIs by Koldobskiy et al. (2022). We used this approach, enriched with Monte Carlo modeling to assess the uncertainties as described below, and CI-proxy data Q from Mekhaldi et al. (2015) and Brehm et al. (2022) to calculate integral fluxes F 60 , F 234 , and F 236 , corresponding to 36 Cl, 14 C, and 10 Be, respectively, for the ESPE events listed in Table 2 in the main text. The values of Q need to be corrected for the changing geomagnetic field, which is typically represented by its virtual axial dipole moment, VADM (e.g., Nevalainen et al. 2013). To account for the realistic geomagnetic conditions corresponding to the exact time of the ESPE, we used the verified correction method (see Equation 5 in Koldobskiy et al. 2022) to reduce the value of Q to the standard present-day conditions, where M 0 = 7.75 · 10 22 A·m 2 is the present-day VADM, and the factor γ can be found in Table 2 of Koldobskiy et al. (2022). To account for notable uncertainties in the VADM reconstruction, we used the Monte Carlo modeling applied to two recent VADM series by Knudsen et al. (2008) and Panovska et al. (2018). For each iteration, we first randomly chose these VADM series. Then we randomly selected the exact value of VADM M(t), corresponding to the time t of the ESPE, within the reported uncertainties as where M * (t) and σ M (t) are the listed VADM values and their uncertainty for the time t, and r is a pseudo-random number with zero mean and unit dispersion. Then, the flux F was calculated using Equation B.1, where the value of K eff was randomly selected, similarly to Equation B.3, but for the values taken from Table 2 of Koldobskiy et al. (2022). The flux values calculated in this way were saved for each iteration. After 1000 iterations, we calculated the median ESPE integral flux values and their 68% confidence intervals. Because 14 C and 10 Be have similar effective energies (234 and 236 MeV, respectively), the corresponding fluxes were averaged. The necessary data for both 36 Cl and 10 Be exist for only three ESPEs (Table 2), which enables us to estimate fluxes in different energies. The other five ESPEs have only 14 C measurements, and only F 234 can therefore be estimated for them.
The values of F 200 , F 100 , and F 30 were inter-and extrapolated using the values of F 234 and F 60 . By assuming that the ESPE integral flux has a spectral shape similar to strong regular SEP events (Usoskin et al. 2020;Paleari et al. 2022), we based our inter-and extrapolations on the recent reconstructions of integral fluxes for the 58 strongest recently observed SEP events (Koldobskiy et al. 2021 (Nishiizumi et al. 2009) and the deep lunar-soil core Apollo-15 (R75, N84) (Rancitelli et al. 1975;Nishiizumi et al. 1984). Areas denoted SEP and GCR approximately show the amount of 26 Al produced by SEP and GCR, respectively.
The measured activity (concentration) of a cosmogenic isotope (CI) in a lunar rock or lunar soil has a distinct depth profile in which the isotope content continuously decreases over the depth. The example shown in Fig. C.1 is a compilation of the 26 Al measurements in lunar rock 64455 and the deep soil core Apollo-15 (Rancitelli et al. 1975;Nishiizumi et al. 1984Nishiizumi et al. , 2009). The double-slope shape is caused by the sum of CIs produced by GCRs and SEPs. Because GCRs are more energetic and thus produce CI in deeper layers (see Appendix A), CI measured in deep soil cores at depths > 10 g/cm 2 can be used to estimate the average GCR flux, which forms the background production of CI in shallower layers. The fraction of CI produced by SEP is calculated as the excess above the GCR-related CI background for each depth. This enables estimating the average energy spectrum of the SEP flux, as illustrated in Fig. C.1. The uppermost layers are corrected for surface erosion.
Full details of how exactly GCR and SEP average fluxes are reconstructed from the CI measurements can be found in Poluianov et al. (2018). The reconstructed integral SEP spectrum is shown in Fig. 3 of the main text, along with the full-range uncertainties. . (D.1) The best fit was found by a scan over the grid space of parameters. The values of k * and λ * and the corresponding minimum χ 2 were stored for each iteration.
After 10000 iterations, we analyzed the distribution of χ 2 values. The set of parameters giving the minimum value χ 2 min was considered as the best-fit solution with the 68% c.i. of the parameters defined as χ 2 ≤ χ 2 min + 2.30. In Figure D.1 we show an example of the distributions of the Weibull parameters versus χ 2 value and the parameter correlation plot for F 30 obtained for GOES SC 24 and terrestrial CIs. The best-fit solution is denoted with the red stars, and blue dots denote solutions within the 68% c.i. The corresponding values of the parameters are listed in Table 3 (first line, RHS block) in the main text. As shown in Fig. D.1 (RHS block), the obtained parameters k and λ appear to be strongly correlated.
The results are qualitatively similar for other energies and are listed in Table 3 in the main text. Distributions of Weibull parameters λ and k vs. χ 2 and parameter correlation plot for F 30 fit obtained using GOES SC 24 and terrestrial CI proxy data. Red stars correspond to the best-fit parameter values (χ 2 min ), and blue dots have χ 2 ≤ χ 2 min + 2.30, as shown with the dashed red line.