Effective Energy of Cosmogenic Isotope (10Be, 14C and 36Cl) Production by Solar Energetic Particles and Galactic Cosmic Rays

Cosmogenic isotopes 14C, 10Be and 36Cl measured in datable natural archives provide the only known quantitative proxy for cosmic‐ray (CR) and solar‐activity variability before the era of direct measurements. Studies of relations between the measured isotope concentrations and CR variability require complicated modeling including the isotope production and transport in the terrestrial system. Here we propose a rough “effective energy” method to make quick estimates of the CR variability directly from the cosmogenic data using an approximate linear scaling between the measured isotope concentrations and the energy‐integrated flux of CR above the effective energy. The method is based on the thoroughly computed effective yield function presented here. A simple way to account for the variable geomagnetic field is also provided. The method was developed for both solar energetic particles (SEPs) and galactic cosmic ray (GCR) variability and is shown to provide a robust result within 20% and 1% accuracy, respectively, without an assumption of the specific spectral shape. Applications of the effective‐energy method to the known extreme SEP events and the secular GCR variability are discussed. The new method provides a simple and quick tool to assess the CR variability in the past. On the other hand, it does not supersede the full detailed modeling required for precise results.

concentrations correspond not to the local cosmic-ray-induced cascades but to a weighted flux of cosmic rays distributed over the globe. Generally, the production and transport of cosmogenic isotopes in the terrestrial atmosphere can be considered as an energy-integrating distributed cosmic-ray detector, which does not directly measure the energy spectrum but rather the integral response. The response is defined by the yield function of the detector (see details in Section 2). The ideal integral detector should have a step-like yield function characterized by a threshold energy E th so that the detector's efficiency is zero and unity for energies below and above E th , respectively. Such a detector could directly measure the integral flux of cosmic rays with energy above the threshold. The realistic yield function of cosmogenic isotope production is not step-like but grows with the particle's energy (Poluianov et al., 2016). However, in some cases, the so-called effective energy E eff exists that reduces the real detector to the ideal one with E eff being the threshold energy, with no or little dependence on the exact spectral shape.
The concept of effective energy was introduced to study GCR variability (Alanko et al., 2003; and SEPs  using data from the NM network as well as for cosmogenic isotope production in lunar rocks (Poluianov et al., 2018). The effective-energy concept was applied by Kovaltsov et al. (2014) to the production of cosmogenic isotopes 14 C and 10 Be in terrestrial archives by SEP events. In this paper, we further develop the concept of the effective energy for cosmogenic isotopes using updated yield functions (Poluianov et al., 2016) and revisited energy spectra of SEP/GLE events (Koldobskiy et al., 2021). We also consider GCR directly measured by AMS-02 (Aguilar et al., 2018a).

Modeling
The concept of effective energy E eff or rigidity is based on an assumption (and its validation) that the energy-integrated flux (in the number of particles per [cm 2 sr s]) of cosmic rays with energy above it J (>E eff ) is roughly proportional to the measured quantity Q: Since the intensity of SEP/GLE events is usually characterized by the energy-integrated omnidirectional fluence F (>E) (in units of particles per cm 2 ), this equation for SEP/GLE events is For cosmogenic isotopes, Q can be the depositional flux (in the framework of a specific transport model) or the production rate and can be represented as (see the full formalism in Section 3 of : where the summation is over different species of cosmic-ray particles, J i (E) is the differential energy spectrum of particles of the ith type with the atomic mass A i , Y eff,i (E, M) is the isotope's effective yield function for nucleons, which depends on the geomagnetic field characterized by the virtual axial dipole moment (VADM - Korte & Constable, 2005) M and transport of the isotope in the atmosphere. The differential energy spectrum is related to the omnidirectional integrated flux as All the computations here were made for the reference present-day VADM of M 0 = 7.75⋅10 22 A m 2 (Thébault et al., 2015). Corrections for other VADM values for both SEP and GCR are described in Section 3.3.
Here we computed the effective yield functions for the three isotopes applying the formalism described by . We used specific yield functions provided in Poluianov et al. (2016). Global production was considered for 14 C, while transport/deposition of 10 Be and 36 Cl in polar ice was modeled using parameterization by Heikkilä et al. (2009) and Heikkilä et al. (2013). The computed effective yield functions for primary cosmic ray protons are shown in Figure 1 and tabulated in the Supporting Information S1.

Differential Response Function
The product of the energy spectrum and the effective yield function, viz. the integrand of Equation 3, is called the differential response function D = J⋅Y eff and quantifies the sensitivity of the isotope production/deposition to the energy of primary cosmic-ray particles. As an illustration, Figure 2b shows the differential response functions for the GLE event #5 of 23 February 1956 which was the strongest directly observed event with the hardest known spectrum. One can see that 14 C and 10 Be have very similar responses with the highest sensitivity to SEPs of a few hundred MeV; 36 Cl is sensitive to lower-energy SEPs (several tens of MeV); while NM is sensitive to SEPs with the energy of around one GeV. Basing on the significant difference between the peaks for 36 Cl and 10 Be, a rough assessment method of the extreme SPE spectral characteristics was proposed by Webber et al. (2007) and applied for recent events (Mekhaldi et al., 2015;O'Hare et al., 2019). Figure 2c depicts the cumulative relative contributions to the production/deposition Q. Crossings of the cumulative curves with the horizontal line denote the median energies of the isotope production. The median energy varies significantly (not shown) between the analyzed GLE events, especially for 36 Cl, and thus cannot be used as effective energy of cosmogenic isotope production.

Effective Energy
In this Section, we determine the effective energy of cosmogenic isotope production by SEPs as defined by Equation 2.
Here we considered all GLE events with reliably defined spectra, viz. 58 events with the reconstructed proton energy spectra (see Table 2 in Koldobskiy et al., 2021). Contribution of helium and heavier species into the cosmogenic-isotope production is small (less than 5%) for SEPs (Mekhaldi et al., 2015(Mekhaldi et al., , 2021. The time resolution of the cosmogenic data is typically annual, and a sophisticated analysis can lead to, at best, seasonal accuracy of the isotope production (e.g., Büntgen et al., 2018;Sukhodolov et al., 2017;Uusitalo et al., 2018). Accordingly, individual SEP events cannot be distinguished in the cosmogenic-isotope data from a series of consecutive events originating from the same solar active region, as for example, took place in 1989. Therefore, we have combined SEP events into annual SEP fluences by summing up fluences of individual events over a calendar year. Thirty Cl. The values are tabulated in the Supporting Information S1. The curves were scaled, as indicated in the legend, for better visibility. For comparison, the yield function (in sr m 2 ) of a polar sea-level 6NM64  is also shown.

of 13
yearly SEP fluences were made as shown in Table 1. We used these annual fluences as an input for the production model (Equation 3).
For all annual-fluence spectra, we computed, using Equation 3, the corresponding productions Q of each isotope as shown in Table 1. One can see that the annual SEP-related production varies by up to three orders of magnitude between different years. Next, we calculated the scaling factor as a function of energy following Equation 2, viz. ( ) = (> )∕ . Figures 3a-3c show the relation between  and E for different cosmogenic isotopes. The relations for individual annual fluences are shown by gray curves. One can see that, despite the huge difference in the annual isotope production rates, all the gray curves cross at nearly the same location, making a "bow-tie" plot (e.g., Raukunen et al., 2020), implying the existence of the effective energy so that Equation 2 is valid irrespectively of the exact spectrum of SEPs. The effective energy E eff is defined as the energy where the relative dispersion of individual curves in Figure 3 is minimal. The relative dispersion was quantified as follows. Vertical slices of the  -vs-E plots ( Figure 3) were taken for different values of E. For each such slice, the distribution of ( ) values was constructed, and the mean  > and the dispersion  were computed. Their ratio, wiz. the relative dispersion of the scaling factor ∕ <  > was used as the merit function so that we defined the effective energy E eff as the value of energy E which minimizes the relative dispersion as shown in Figure 4.
The best-fit scaling factor eff is defined as the mean value over all individual curves, corresponding to E eff , viz. ( eff) > , and its 68% uncertainties are taken as ( eff) . The uncertainties of the E eff are defined as the 68% range of the energies at the best-fit eff for individual curves. The values and 68% confidence intervals of the effective energies and scaling factors are shown by the red crosses in Figures 3a-3c and given in the central block of Table 2 for the three isotopes studied here.
One can see that the E eff is defined quite accurately, within 11%. The effective energy appears to be 230-240 MeV for 14 C and 10 Be, in agreement with a previous study by Kovaltsov et al. (2014), but lower (≈60 MeV) for 36 Cl. The scaling factors are different by up to four orders of magnitude between different isotopes, from 44 for 14 C to ≈ 5 ⋅ 10 5 for 36 Cl but eff is very stable for each individual isotope and is defined with the accuracy of better than 25%.
Quite often, the ratio of F ( 30 MeV)/F ( 200 MeV) is considered as an estimate of a measure of the hardness of the SEP spectrum for extreme events (e.g., Asvestari, Willamo, et al., 2017;Cliver et al., 2020;Mekhaldi et al., 2021), where F ( 200 MeV) is roughly related to the effective energy of 14 C production (Kovaltsov et al., 2014), while F ( 30 MeV) is a typical quantity for SEP fluence based on direct in-situ measurements during the space era (e.g., Shea & Smart, 1990). However, as shown here, the effective energy of 36 Cl production/deposition is about 60 MeV, and F ( 30 MeV) cannot be reliably defined from cosmogenic-isotope data. Accordingly, the F ( 60 MeV)/F ( 240 MeV) ratio ought to be used instead.

Dependence on the Geomagnetic Field
Figure 3 corresponds to the current epoch with the geomagnetic field dipole moment of M = 7.75⋅10 22 A m 2 . However, extreme SEP events detectable in the cosmogenic isotopes occur during different times over millennia (e.g., Miyake et al., 2020), when the geomagnetic shielding can be different. Accordingly, we modeled the isotope production for different values of M ranging from (6-12)⋅10 22 A m 2 as corresponding to the VADM variability during the last 10 millennia of the Holocene (e.g., Korte et al., 2011;Usoskin et al., 2016). The effective energy E eff appears amazingly independent on the geomagnetic field in a wide range of the VADM values: E eff changes within 2% in the full range of the M-values that is much smaller than the 68% confidence intervals for E eff of 8-11% (Table 2). Thus, the effective energy is defined unambiguously and very accurately for SEPs for a broad range of the VADM values. The scaling factor changes more significantly, within ±20%, viz. larger than the error bars, following the corresponding changes in the production. The changes are systematic and can be accurately (within 0.2% in the range of M of [6-12]⋅10 22 A m 2 ) approximated (see Figure 5a) as:  Therefore, from the measured production/deposition of cosmogenic isotopes in relation to extreme SEP events or candidates, one can directly and straightforwardly assess (within ±20% for the known VADM) several spectral points of the SEP energy spectrum, without any assumptions of the exact spectral shape.

An Example of Application to the Extreme Solar Events in the Past
As an example of the application of the effective-energy concept for a simple assessment of the SEP events recorded in cosmogenic proxy data, we considered here extreme solar events discovered recently over the last millennia (see s review Miyake et al., 2020). Presently, three confirmed extreme SEP events are known: 774/5 CE originally discovered by Miyake et al. (2012) in Δ 14 C measured in a Japanese cedar tree and later confirmed in other datasets (Büntgen et al., 2018;Sukhodolov et al., 2017;Usoskin et al., 2013;Uusitalo et al., 2018); 660 BCE discovered by Park et al. (2017) and confirmed later (Büntgen et al., 2018;O'Hare et al., 2019); and 993 CE discovered by Miyake et al. (2013) and confirmed later (Mekhaldi et al., 2015;Sakurai et al., 2020). In addition, three candidates to extreme SEP events have been found recently in 14 C datasets and are pending confirmation with other isotopes: 1052 CE and 1279 CE were discovered by Brehm et al. (2021), and 5410 BCE by Miyake et al. (2021).
For all these events we have calculated the fluence F 234 ≡ F (>234 MeV) based on 14 C data following the methodology developed here, as shown in Table 3. First, from the published values of the 14 C production rates Q and geomagnetic dipole strength VADM M for each event, we computed the production rate 0 reduced to the standard geomagnetic conditions using Equation 5. Next, using Equation 2, we calculated the omnidirectional fluence F 234 , where the values of E eff and eff were obtained from Table 2. Finally, the obtained fluence F 234 was compared with that for the GLE#5 (1.38 ± 0.37)⋅10 8 cm −2 (Koldobskiy et al., 2021) as the R 1956 value. The events in Table 3 are sorted by their (corrected for the VADM value) strength. The strongest SEP event of 774 CE was approximately 70 times stronger than the GLE#5 making it an extreme and maximum known event. Other events take the R 1956 value from 15 to 70 (within uncertainties) reflecting the occurrence probability distribution of the extreme SEP events. It is interesting that three event candidates have the R 1956 value of 20-30 which approached the theoretical sensitivity limit of the 14 C method of R 1956 = 15 .
It should be noted that the R 1956 values in Table 3 are approximate as based only on 14 C data, while other isotopes may lead to slightly different values. A more detailed study will form a subject to a forthcoming work. Note. The computations were performed for the reference VADM M 0 = 7.75⋅10 22 A⋅m 2 and correspond to Figure 3. The shown uncertainties correspond to the 68% confidence interval.

Table 2 The Effective Energy E eff , the Scaling Factor eff , and the Parameter γ of the Relation for the Virtual Axial Dipole Moment (VADM, Equation 5), Calculated Here for the Three Cosmogenic Isotopes, for Solar Energetic Particle (SEP) (Central Block) and for Galactic Cosmic Ray (GCR) (Right-Hand-Side Block)
Figure 3. "Bow-tie" plot of the relation between  and E (Equation 1) for cosmogenic isotopes 14 C, 10 Be and 36 Cl. Lefthand-side (a-c) panels correspond to solar energetic particles (SEPs), where gray curves correspond to different SEP annual fluences (Table 1). Right-hand-side (d-f) panels correspond to galactic cosmic rays (GCRs), where gray curves correspond to GCR nucleonic spectra for different Bartels rotations. The geomagnetic dipole moment (VADM) M is set to 7.75⋅10 22 A m 2 . The red dots depict the best-fit values of eff and E eff (see Table 2), while red crosses denote the 68% confidence levels.

Effective Energy for Galactic Cosmic Rays
In this Section, we computed the effective energy for the production of the three cosmogenic isotopes by GCR. As the spectra of GCR, we considered direct measurements of protons, helium and heavier species performed by the AMS-02 experiment (Aguilar et al., 2017(Aguilar et al., , 2018a(Aguilar et al., , 2018b during the period of 2011-2017 with a 27-day (Bartels rotations, BRs) cadence. The data set covers 79 BRs from 15-May-2011 (start of BR # 2426) through 09-May-2017 (end of BR # 2506) corresponding to the period from the mid-ascending to the mid-descending phases of solar cycle #24 and characterized by moderate modulation of GCR. By considering directly measured spectra of GCRs, we avoid uncertainties related to the modeled GCR spectra Herbst et al., 2010). In contrast to SEPs consisting mostly of protons, GCR include also helium and heavier species whose contribution to the atmospheric effects can be significant (Koldobskiy, Bindi, et al., 2019). Since heavier species are modulated and shielded by magnetic fields differently with respect to protons of the same energy, they cannot be represented by a simple scaling of protons. Accordingly, we consider the GCR spectrum of not only protons but also heavier species reduced to that of nucleons with the energy expressed per nucleon.  The event year. b The global 14 C production Q (10 8 cm −2 ) corresponding to the measurements. c Reference to the 14 C data (Büntgen et

Differential Response Function
Similarly to SEPs, we computed the differential response function for GCR nucleons as shown in Figure 2 (panels d-f). Panel d depicts the energy spectrum of all nucleons of GCR defined as where J i (E) is the energy spectrum of ith specie of GCR with atomic mass A i , and E represents kinetic energy per nucleon. Panel e shows the differential response function. One can see that it peaks at much higher energies than that for SEPs reflecting the different spectra of GCR and SEPs. The shape of the response functions for both 10 Be and 36 Cl produced by spallation reactions are nearly identical to each other, while 14 C produced by (n,p) reactions is slightly more sensitive to CR in the energy range of 3-20 GeV. The response functions for all the studied isotopes peak at around 1.5 GeV energy, while the peak of the polar sea-level NM response corresponds to the energy of about 5 GeV. Figure 2f depicts the cumulative relative contribution, which defines the median energies, that are 4.5, 4.5, 5.9 and 11 GeV for the 10 Be, 36 Cl, 14 C and polar NM, respectively. The median energy for the polar NM is consistent with earlier estimates (Ahluwalia & Fikani, 2007). This is nearly an order of magnitude higher than those for SEPs.

Effective Energy
The "bow-tie" plots for  = n(> )∕ for GCR are shown in Figures 3d-3f, where each gray curve corresponds to one BR of AMS-02 data. One can see that the curves lie more compact than for SEP events (Figures 3a-3c) reflecting the fact that variability of GCRs is smaller and more regular than that for SEPs. The effective scaling eff and energy E eff are defined in the same way as for SEPs and listed in the right-hand-side block of Table 2.
Since the standard way of quantifying the GCR flux is the intensity in units of (cm 2 sr s GeV) −1 , the scaling eff has the dimension of sr −1 in contrast to SEP fluences where it is dimensionless. The effective energy was found very close to each other for 10 Be and 36 Cl, being about 2 GeV/nucleon. The effective energy for 14 C appeared slightly higher, viz. about 2.5 GeV/nucleon. This is related to the fact that 14 C is produced globally, while 10 Be and 36 Cl are more weighted, because of the atmospheric transport/deposition pattern, toward higher latitudes with lower geomagnetic cutoffs.

Dependence on the Geomagnetic Field
The dependence of the isotope's production rate by GCR on the geomagnetic field strength, quantified via the VADM M is shown in Figure 5b for 10 Be produced by GCR. It can be seen that the computed dependence is nearly perfectly described by formula 5 with the index γ = 0.364 (see Table 2). Similar dependencies for other isotopes (not shown) are equally well described by formula 5 with indices listed in Table 2. The dependence is stronger (γ = 0.546) for the globally produced 14 C than for mid/high-latitude dominated production and transport of 10 Be and 36 Cl.

An Example of GCR Intensity Variability Over the Last Millennium
As an example of the application of the effective energy approach to GCR, we show in Figure 6 estimate of the integral intensity J ( 2.48 GeV/nuc) of GCR over the last millennium, since 970 CE. The estimate is based on the recent reconstruction of the annual 14 C global production rate by Brehm et al. (2021) for the last millennium. Conversion to the standard geomagnetic conditions M 0 = 7.75⋅10 22 A m 2 was made by a Monte-Carlo approach (similar to Usoskin et al., 2021) by randomly selecting values of M for each year from four alternative archeomagnetic reconstructions (Hellio & Gillet, 2018;Nilsson et al., 2014;Pavón-Carrasco et al., 2014;Usoskin et al., 2016) and applying the correction using formula 5. Then the corrected production rate was converted, using Equation 1 with values E eff and eff from  (Usoskin, 2017), including grand solar minima ca 1050 (Oort minimum), 1300 (Wolf), 1400 (Spörer), 1700 (Maunder) and 1800 (Dalton). Because of the anthropogenic Suess effect (extensive use of fossil fuel which dilutes 14 C in the atmosphere - Suess, 1955) and bomb-effect (nuclear tests produce a large amount of 14 C), radiocarbon can be hardly used as a cosmic-ray index after the onset of the industrial revolution.
For comparison, we also show, as the red curve in Figure 6, the J ( 2.48 GeV/nuc) flux directly computed since 1951 from the world NM-network data using the methodology of Usoskin et al. (2017) and Koldobskiy, Bindi, et al. (2019). It can be seen that the level of GCR is significantly lower in the second half of the twentieth century than that for the millennium before. It is comparable only to a short period ca. 1350 CE between Wolf and Spörer grand minima. This systematically reduced flux of GCR corresponds to the Modern grand maximum of solar activity (Usoskin, 2017). On the other hand, the GCR flux was higher during the last two weaker solar cycles being comparable to that of the late nineteenth century, in agreement with the solar activity estimates (e.g., Hathaway, 2015). This indirectly confirms the validity of our approach and the overall agreement between indirect 14 C-based method to estimate cosmic-ray variability and direct estimates for the modern era.

Conclusions
Here we developed and proposed a simple and quick "effective energy" method to estimate the flux of cosmic-ray particles directly from concentrations of cosmogenic isotopes 14 C, 10 Be and 36 Cl measured in natural archives. Even though the method is simple, it is based on full detailed simulations of the production and transport of the cosmogenic isotopes in the terrestrial environment, including the thoroughly computed effective yield function tabulated in the Supporting Information S1.
The method is based on an assumption that there is such effective energy E eff that the relation between the flux of CR particles with energy above it F (> eff ) is directly proportional to the measured concentration of the isotope (Equation 3), irrespectively of the exact energy spectrum of CR particles. We first verified this assumption separately for SEPs and GCRs (Figure 3) using the data of recorded GLE events for the last 60 years for SEPs and the space-borne measurements by AMS-02 detector from 2011 to 2017 for GCR. We calculated the values of the effective energy E eff and the scaling factor eff separately for each isotope and each type of CRs, as presented in Table 2. The accuracy of the method is within 20% for SEP and within 1% for GCR. To account for the variable geomagnetic field quantified through the VADM concept, we provided a simple way (Equation 5) to reduce the measured data to the standard conditions corresponding to the modern epoch with VADM M 0 = 7.75⋅10 22 A m 2 .
As an illustration, the method was applied to the data of 14 C for six known extreme SEP events or candidates that allowed us to rank these events by strength in comparison to the strongest directly observed SEP event of 23-Feb-1956 (GLE#5), as summarized in Table 3. We also provided, as an example, the temporal evolution of the integral nucleonic flux of GCR J ( 2.48 GeV/nuc) over the last millennium ( Figure 6) as based on the recent data of the Figure 6. Time profile of the galactic cosmic ray integral nucleonic flux J ( 2.48 GeV/nuc), estimated using the 14 C annual production rate for the last millennium Usoskin et al., 2021). The 68% confidence level uncertainties are shown by the gray shaded area. The red line depicts the annual values of F ( 2.48 GeV/nuc) calculated from the NM network data since 1951.
10.1029/2021JA029919 11 of 13 annual 14 C global production rates. The estimated flux agrees well with that calculated from the direct NM data for the last decades, providing an indirect verification of the method.
Thus, we have developed an approximate method for a simple and quick estimate of the CR integral fluxes directly from the measured cosmogenic isotope measurements (production rate or depositional flux) for SEPs and GCRs. This makes a useful tool for rough assessments of the CR variability and spectra. However, the proposed method does not supersede the full modeling of the isotope production and transport that is still required for precise data analysis and interpretation.

Data Availability Statement
Spectra of solar energetic particles can be computed from the information presented by Koldobskiy et al. (2021), viz. parameters (Table 2) therein applied to Equations 1-3 therein. Spectra of GCR are available from AMS-02 data (Aguilar et al., 2018a), collected in the Cosmic Ray Database https://tools.ssdc.asi.it/CosmicRays (Di Felice et al., 2017).