Issue 
A&A
Volume 647, March 2021



Article Number  A132  
Number of page(s)  16  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202040058  
Published online  22 March 2021 
New reconstruction of eventintegrated spectra (spectral fluences) for major solar energetic particle events^{⋆}
^{1}
University of Oulu, Oulu, Finland
email: koldobskiy@gmail.com
^{2}
National Research Nuclear University MEPhI, Moscow, Russia
^{3}
St. Petersburg State University, St Petersburg, Russia
^{4}
University of Turku, Turku, Finland
^{5}
Ioffe PhysicalTechnical Institute, 194021 St. Petersburg, Russia
Received:
3
December
2020
Accepted:
8
January
2021
Aims. Fluences of solar energetic particles (SEPs) are not easy to evaluate, especially for highenergy events (i.e. groundlevel enhancements, GLEs). Earlier estimates of eventintegrated SEP fluences for GLEs were based on partly outdated assumptions and data, and they required revisions. Here, we present the results of a full revision of the spectral fluences for most major SEP events (GLEs) for the period from 1956 to 2017 using updated lowenergy flux estimates along with greatly revisited highenergy flux data and applying the newly invented reconstruction method including an improved neutronmonitor yield function.
Methods. Low and highenergy parts of the SEP fluence were estimated using a revised spaceborne/ionospheric data and groundbased neutron monitors, respectively. The measured data were fitted by the modified Band function spectral shape. The bestfit parameters and their uncertainties were assessed using a direct Monte Carlo method.
Results. A full reconstruction of the eventintegrated spectral fluences was performed in the energy range above 30 MeV, parametrised and tabulated for easy use along with estimates of the 68% confidence intervals.
Conclusions. This forms a solid basis for more precise studies of the physics of solar eruptive events and the transport of energetic particles in the interplanetary medium, as well as the related applications.
Key words: Sun: particle emission / Sun: flares / Sun: activity / solarterrestrial relations
The reconstructed fluences in tabulated form and the corresponding bestfit parameters are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/647/A132
© ESO 2021
1. Introduction
In addition to Galactic cosmic rays (GCR) continuously bombarding the Earth with slightly variable flux, sporadic solar eruptive events, such as flares and/or coronal mass ejections (CMEs) may cause dramatic (by many orders of magnitude) enhancements of energetic particle fluxes near Earth, which are called solar energetic particle (SEP) events (Vainio et al. 2009; Desai & Giacalone 2016; Klein & Dalla 2017). SEP events occur quite frequently during maxima and early declining phases of solaractivity cycles and very seldom during the minimum phase (e.g., Bazilevskaya et al. 2014). Sometimes, eruptive events can be sufficiently energetic to accelerate particles to relatively high energy (> 400 MeV) so that they can penetrate the Earth’s magnetosphere and atmosphere, inducing atmospheric nucleonic cascades, which can be measured by groundbased detectors (e.g, Shea & Smart 2012). Such events are known as groundlevel enhancements (GLEs), which are consequently numbered from #1, which took place in February 1942, to the most recent #72 in September 2017. The first four GLEs were recorded by ionisation chambers (Forbush 1946) and cannot be quantitatively assessed, but those starting from GLE #5 (23Feb1956) were recorded by the network of groundbased neutron monitors (NMs) with a possibility of quantifying their intensities and spectral parameters (e.g., Mishev et al. 2018). All available information of NM data for GLEs #5−72 is collected in the International GLE Database^{1} (IGLED, see Usoskin et al. 2020a).
Studies of SEP events are important for different reasons. On one hand, solar eruptive events are wellobserved processes of energeticparticle acceleration (Vainio & Afanasiev 2018), which can be studied in detail using a multimessenger approach, complementing particle data with observations in different wavelengths (e.g., Plainaki et al. 2014; Cliver 2016; Kocharov et al. 2017). For this purpose, the peak flux intensity and detailed temporal variability of the particle flux are important as signatures of the acceleration process in the solar corona and the interplanetary medium (e.g., Desai & Giacalone 2016; Kong et al. 2017). Accordingly, numerous studies were focused on peak fluxes of SEPs and corresponding acceleration and transport processes (e.g., Kouloumvakos et al. 2015; Kocharov et al. 2017). On the other hand, enhanced fluxes of energetic particles affect the radiation environment near the Earth (e.g., Webber et al. 2007; Mishev et al. 2015), making not only the peak fluxes but also the fluence (eventintegrated flux) and its spectral shape of significant importance, especially for extreme events (e.g., Cliver et al. 2020). We emphasise that SEP fluences can not be used for the detailed study of SEP acceleration processes, because (i) the observations at 1 AU are also modified by transport, and (ii) different energies in the fluence spectrum can be dominated by different acceleration mechanisms or by the same mechanism operating under different conditions. It is evident from the proton timeintensity profiles alone that the fluence at MeV and 10 MeV energies is often dominated by acceleration at interplanetary shocks (e.g., Reames 1999). However, the question is more open at 100 MeV and GeV energies peaking much earlier, with possible contributions from flares and/or coronal shocks as the main candidates to account for the acceleration (see Cliver 2016, and references therein). Even if the same CMEdriven shock were responsible for the acceleration of 10 MeV and 1 GeV protons, the former would typically be accelerated mainly in the solar wind and the latter in the corona, and there is no reason to suggest that the spectral form of the fluence would reveal something common about the accelerator properties.
Composing an eventintegrated energy (or rigidity) spectrum (spectral fluence) of a SEP event is a difficult task. Direct measurements of the spectrum covering the highenergy (> 400 MeV) range can be done with modern spaceborne magnetic spectrometers PAMELA (Payload for Antimatter Matter Exploration and Lightnuclei Astrophysics – Adriani et al. 2014) operated from 2006 to 2016, and AMS02 (Alpha Magnetic Spectrometer – Aguilar et al. 2018), which has been in operation since 2011. SEP events analysed using PAMELA data were presented by Bruno et al. (2018), while data from AMS02 are not available yet (Bindi 2017). However, these instruments are located aboard loworbiting satellites, and thus spend most of their time inside the Earth’s magnetosphere and can detect lowenergy solar particles only intermittently (5−10 min per halforbit), leading to essential uncertainties in both SEP peak fluxes and fluences, especially for impulsive events. Additionally, uncertainties of the SEP fluxes were large during the earlier years (e.g., Reeves et al. 1992; Tylka et al. 1997) due to the saturation of the detectors by strong particle fluxes and the possible penetration of highenergy particles into the detector through the walls of the collimator, leading to an enhanced effective acceptance.
Thus, for most of the events, one has to combine data from different instruments, including lowenergy (< 400 MeV) spaceborne detectors located beyond the magnetosphere, and energyintegrating (above 400 MeV) groundbased NMs. The reconstruction of the spectrum is typically done by fitting a prescribed spectral shape to the data and finding the bestfit parameters of this shape along with their uncertainties. The first consistent reconstruction of combined spectral fluences for major GLE events was made by Tylka & Dietrich (2009) and updated by Raukunen et al. (2018). It was based on a combination of lowerenergy data from different spaceborne detectors and the highenergy tail based on groundbased NM datasets. The spectral fluence was estimated by fitting the Bandfunction spectral shape (double power law with an exponential junction – Band et al. 1993) to the integral rigidity spectral fluence. The spaceborne data were collected from different sources, while NM data were analysed by applying the NM yieldfunction by Clem & Dorman (2000). Spectral fluences were presented as tabulated parameters of the Bandfunction approximation for each GLE event. This dataset has been extensively used in numerous studies for solar and space physics (Cliver et al. 2020; Anastasiadis et al. 2019; Herbst et al. 2019), but it has become obsolete and requires an essential revision. First, the highenergy spaceborne data have been essentially revisited and corrected for known errors (Raukunen et al. 2020). Second, the data of the NM network for all GLE events have been revisited (Usoskin et al. 2020a): apparent errors were corrected, and the variable GCR background was taken into account. Furthermore, a most recent and directly verified using AMS02 data NM yield function (Mishev et al. 2013, 2020) and a new effectiveenergy analysis method were developed (Koldobskiy et al. 2018, 2019). Moreover, the Band function is not an optimum parametric shape for the GLE energy (or rigidity) spectrum, which requires a rolloff at the highest energies. Here, we present a complete revision of the reconstructions of the SEP spectral fluences for major SEP events (with GLE) using most uptodate knowledge of the SEP measurements on ground and in space, new models, and a modified spectral form including a rolloff at high energies.
2. Datasets
All sources of the energy/rigidity integral fluences used in this study are described below.
2.1. Spaceborne and ionospheric data
Data for the rigidities below 1 GV (energy < 430 MeV) were taken from spaceborne or ionospheric (in the presatellite era) data. For the period since 1989 (i.e. GLEs # 40−72), we used all publicly available data from the GOES (Geostationary Operational Environmental Satellite) energetic particle sensor (EPS) and high energy proton and alpha detector (HEPAD) datasets^{2} (Onsager et al. 1996; Sellers & Hanser 1996). The fluences at the lowenergy channels in this study, that is, > 30 MeV, > 50 MeV, > 60 MeV, and > 100 MeV, were calculated directly from the EPS dataset, but the higher energy HEPAD data were revised using a socalled bowtie method (Van Allen et al. 1974). The nominal HEPAD channels are wide in energy, and they have responses that vary significantly within the channels, and in some cases even outside the nominal channel range. In the bowtie analysis, calibrated channel responses were folded with an assumed spectral form (powerlaw), and by varying the spectral index within a realistic range for SEP events, optimal effectiveenergy and geometricfactor values were found for each channel. The method is explained in full detail in Raukunen et al. (2020). In addition, the cleaning of data spikes (single points with spuriously increased flux) and background subtraction was performed on both GOES datasets. We note that energy bounds of the HEPAD channels changed slightly after 1995 (GLE # 53) due to differences in the calibration procedures used for GOES6 (until 1995) and GOES8 onwards (after 1995). GLEs # 71 and 72 were measured by GOES13, with 10−100 MeV channels represented by two detectors facing east and west, so that the fluences were taken as the average of these two measurements. Resulting omnidirectional fluences (in units of 10^{5} cm^{−2}) are presented in Table 1. Data channels with rigidity above 1 GV were not used in the fitting procedure (Sect. 3) because of their lower reliability.
Eventintegrated omnidirectional integral fluences F(> E) (in units of 10^{5} cm^{−2}) obtained here for GLEs # 40−72.
We also used the SEP spectrum measured by the PAMELA experiment^{3} (Bruno et al. 2018) for the GLE #71 (17May2012). For years before 1989, we used fluences from several sources based on different spacecraft and experiments (King 1974; Reedy 1977; Goswami et al. 1988; Feynman & Gabriel 1990; Jun et al. 2007; Webber et al. 2007). The exact sources of the lowenergy data for each event are specified in the legends of the panels of the figure in Appendix A.
2.2. Groundbased data
Data for rigidity above 1 GV (energy > 430 MeV) were used from a recent reconstruction of SEP fluences based on the detrended data of the NM network from IGLED (Usoskin et al. 2020a). For some events and some NMs (mostly lowlatitude ones), only an upper limit of the enhancement can be set (see vertical blue bars in Fig. 1). A robust estimation of the SEP fluence using NM data requires the accurate identification of the GLE signal itself. Since GLEs # 6, 7, 9, 14, 15, 17, 34, 54, 57, and 68 were weak and lacking sufficiently clear signals, we did not analyse these events. Originally, the SEP event integration period was deduced from highenergy NM detrended data (Usoskin et al. 2020a) and was applied to the lowenergy data where we had this opportunity (i.e. starting from GLE #40). For earlier GLEs, we used all available data, which sometimes have bigger integration periods. It was taken into account as an additional uncertainty in Step 2 below.
Fig. 1.
Example of the fitting of the modified Band function (Eqs. (1) and (2)) to the data points (Sect. 3) for GLE #43, 19Oct1989. Blue crosses and vertical lines denote the highenergy NM data with uncertainties, and the upper limits, respectively (Usoskin et al. 2020a), while red circles correspond to lowenergy data from GOES6 (see Table 1). The grey line with shading depicts the bestfit MBF, along with its 68% confidence interval. 
3. Fitting procedure
For each event, the integral fluence points in both highenergy (Sect. 2.2) and lowenergy (Sect. 2.1) ranges were fitted with a prescribed spectral shape, as described below. Here, we used a modified Band function (MBF), which includes an exponential rolloff term, similar to the EllisonRamaty spectral shape (Ellison & Ramaty 1985):
where F(R) ≡ f(> R) is the omnidirectional fluence (in units of cm^{−2}) of particles with a rigidity greater than R, R is expressed in gigavolts, parameters γ_{1}, γ_{2}, R_{1}, R_{2}, and J_{2} are defined by fitting, and other parameters can be calculated as
This function is constructed in such a way that it and its first derivative are continuous, providing a smooth junction between the two parts. We used a fitting method based on a threestep Monte Carlo procedure further developed after Usoskin et al. (2020a), as illustrated in Fig. 1. The method includes L Monte Carlo iterations.
3.1. Step 1: Fitting the highenergy (NM) part
First, we fitted the highenergy part of the spectral shape (Eq. (2)) using the data from the NM network (blue crosses in Fig. 1) for each event with M spectral points with uncertainties, and m points, corresponding to the upper limits of the fluence (vertical blue lines in the figure). The uncertainties are both statistical (countrate statistic) and systematic, related to the ‘effectiverigidity’ (or bowtie) method ones. For the lth iteration, we simulated a set of randomised ‘exact’ fluence F^{*}(R^{*}) values, applying the same procedure as in Usoskin et al. (2020a), so that for each point k (1...M) the value of was randomly and uniformly taken inside R_{k} error bars; the value of was computed using Eq. (4) from Usoskin et al. (2020a), where the scaling factor K_{k} was randomly taken inside the error bars using the uniform distribution, and the GLE integral intensity X_{k} (in units of % h above the GCR background over the entire duration of the event, see Asvestari et al. 2017) was randomly taken inside the error bars, defined as σ_{Xk} = max[1; 0.1 ⋅ X_{k}]% h using the normal distribution. In order to avoid a bias towards more numerous polar NMs during the fitting procedure, we considered only one randomly selected point (NM) in each rigidity bin of 0.4 GV widths. For several weak events (# 18, 35, 53, 58, 63, 64), the bin width was reduced to 0.1 GV to keep the number of fitted points reasonable. We checked that reducing the bin size does not lead to any significant bias in the results. Then, this set of points F^{*}(R^{*}) was fitted by the spectral shape (Eq. (2)) applying a nonlinear leastsquares method (scipy.optimize.curve_fit function in Python), and the bestfit parameters J_{2}, γ_{2}, and R_{2} were found based on the minimisation of the logarithmic residual D:
where F_{fit} is the value computed using the fitfunction (Eq. (2)) for the rigidity R^{*}. We additionally checked that the obtained bestfit parameters are physically reasonable, that is, the obtained function does not have a positive second derivative anywhere in the highrigidity (> 1 GV) range since the differential SEP fluence is not expected to increase with R in this range. This is quantified as the condition that R_{2} > 0, γ_{2} > 0 or, for γ_{2} ≤ 0, we required that . We also required the formal fit to exceed none of the m upper limits. Fits that did not satisfy these conditions were discarded, and the corresponding Monte Carlo simulation was redone without counting it in the statistic. The obtained bestfit set of parameters (J_{2, l}, γ_{2, l,}, and R_{2, l}) was fixed and used in Step 2.
3.2. Step 2: Fitting the lowenergy part
The set of N lowenergy (R < 1 GV) fluences F^{*}(R^{*}) was obtained in a similar way to Step 1, assuming a fixed value of the uncertainty being 10% (20% for the preGOES/HEPAD era before 1989)^{4} of the tabulated value and applying the uniform distribution. Only points with R > 168 MV (corresponding to energy above 15 MeV for protons) were considered, since SEPs with the energy < 15 MeV may have different sources and different spectral parameters (e.g., Reames 1999; Cliver 2016).
Next, each of these F^{*}(R^{*}) values was divided by the extrapolated function obtained in Step 1 to form a set of values X(R^{*}):
The obtained set of points X(R^{*}) was fitted using the same nonlinear leastsquares method as in Step 1, with the following function:
which can be obtained by dividing the lowrigidity part of MBF (Eq. (1)) by its highrigidity part (Eq. (2)). We checked that the obtained fit parameters are mathematically reasonable, so the bestfit function does not have a positive derivative anywhere, since the integral fluence cannot increase with R. This condition is quantified as γ_{1}/R + 1/R_{1} > 0 for the rigidity range from 137 MV (equivalent to 10 MeV energy) to R_{b}.
Thus, for each lth iteration, a set of parameters J_{2, l}, γ_{2, l}, γ_{1, l}, R_{2, l}, R_{1, l} was calculated for an analysed event. Then, the formal value was computed as the merit function between the fitted curve 𝔉 and the data points F(R):
This set of parameters and the value of were recorded for the lth iteration, and then a new (l+1)st iteration started.
This procedure was repeated L = 5000 times, and the bestfit parameters corresponding to the minimum (among all L iterations) value of were saved and used in the next step of the procedure. We also checked points for outliers. If any data point contributed more than 100 to the total χ^{2} of the bestfit option (i.e. laying beyond 10σ from the bestfit curve), such a data point was discarded and the fit redone for that event. Only one outlier of this kind was found – the point corresponding to the fluence of protons with energy > 360 MeV for GLE #24 from Webber et al. (2007). The set of parameters corresponding to the was selected as the bestfit set and used in Step 3.
3.3. Step 3: Evaluation of the uncertainties of the parameters
Next, we performed an additional Monte Carlo study of the uncertainties of the obtained bestfit parameters. The value of each parameter was varied randomly (and independently of each other) so that the new value of a parameter P (any of the five parameters of the MBF) was taken as
where P_{0} is the bestfit value found in Steps 1 and 2, and r is a normally distributed pseudorandom number with a zero mean and the standard deviation of 0.5. All five parameters of the MBF were simultaneously and independently randomised in this way. If these parameters were not rejected by the physicalcriteria checks (described in Steps 1 and 2), the formal χ^{2} value (Eq. (7)) was calculated for the corresponding MBF and the data points. If the obtained χ^{2} value did not exceed , the corresponding set of parameters was recorded as being within a 68% confidence interval (c.i.) for a fiveparameter model (e.g., Chap. 15.6 of Press et al. 2007), otherwise it was discarded and the simulation was redone. If a value of was obtained during this step, it was assigned as a new , and the values of the bestfit parameters P_{0} reset, and Step 3 was restarted anew. This procedure was repeated 10 000 times for each analysed event, which involved (including the discarded iterations) 3 × 10^{8} iterations for GLE #5, shown in Fig. 2.
Fig. 2.
Distribution of the parameters and their pairwise correlations for the MBF fitting of the GLE #5 obtained for 10 000 iterations (see Step 3 in Sect. 3). The red dots in the bottom panels correspond to the minimum , while the horizontal red dashed line denotes the 68% confidence interval for the bestfit parameter values. Colour intensity corresponds to the point density. The diagonal panels depict histograms of the parameter values’ distribution. 
Figure 2 shows an example of such an analysis. Similar plots were constructed and analysed for all the studied events. The bottom row of panels depicts the relations between the χ^{2} value and the value of each of the five parameters. Each panel contains 10 000 points satisfying the condition of (see above). One can see that the relations have the inverted (generally asymmetric) bellshaped profile, clearly defining the uncertainties for each of the parameter. The corresponding ranges of the MBF fits are shown as grey shaded areas in Fig. 1 and Appendix A. However, it would be incorrect to provide 68% c.i. uncertainties for each parameter independently, since some of them are tightly interrelated (e.g., γ_{1} and R_{1}, and γ_{2}, R_{2}, and J_{2}), as can be seen in Fig. 2. We found that the two parts of the MBF (Eqs. (1) and (2)) can be considered independent (see Fig. 2 (b, c, d, e, g, h)), but within each part, the parameters’ values are highly correlated (a, f, i, j). As the key parameters, we consider γ_{1} and γ_{2}, which are not fully independent in the statistical sense (the Pearson’s correlation coefficient between them is 0.2), but the power of their common variability is only 4%, and thus their uncertainties can be assumed as roughly independent. Confidence interval (68%) uncertainties for γ_{1} and γ_{2} are given in the CDS table. Uncertainties of other parameters are related to these two key parameters and cannot be considered independent (Pearson’s correlation coefficients are statistically significant, ranging from 0.86 to 0.96). Accordingly, the uncertainties of other parameters can be calculated from those of γ_{1} and γ_{2}. The overall quality of the spectral fit can be evaluated via the range δF defined as the average ratio:
where F_{up}(R) and F_{low}(R) are the upper and lower bound values of the MBF for all combinations of the parameters within the 68% c.i. range ( – see Fig. 2) for each value of rigidity R; this is depicted as the upper and lower envelopes of the grey area in Fig. 1, and the averaging is performed in the (logarithmically sampled) rigidity range from the lowest rigidity R_{s} for the spaceborne data point (typically 230 MV), to the highest rigidity R_{n} of nonzero NM data points. The obtained bestfit MBF parameters (R_{2} values exceeding 100 GV are given as +∞), along with the 68% uncertainties of the key parameters γ_{1} and γ_{2}, as well as the values of δF for the analysed events are summarised in Table 2 and provided in machinereadable format at the CDS. A table of the numerical values of SEP fluences calculated using the obtained bestfit parameters is also available at the CDS. We note that the obtained bestfit spectra are not recommended to be extrapolated to energies below 30 MeV (i.e. beyond the energy range used for the fit) as this can lead to unphysical results in the lowenergy range.
4. Conclusions
In this work, eventintegrated fluences of solar energetic particles were reevaluated for most major SEP events (GLEs) using updated lowenergy flux estimates, greatly improved highenergy flux data, and the newly developed reconstruction methods. The earlier estimates (Tylka & Dietrich 2009; Raukunen et al. 2018) were essentially revisited here, providing an accurate parametrisation of the rigidity spectra, with the spectral parameters being tabulated in Table 2 and provided at the CDS. In particular, it was shown earlier (Usoskin et al. 2020a, Fig. 6) that the Band function spectral shape can lead to a significant overestimate of the highenergy tail of the spectrum, and a strong rolloff (assumed to be exponential here) is required. Accordingly, for the parametrisation of the spectral fluence shape, we propose a modified Band function (Eqs. (1)–(2)), which is a combination of the standard Band and a modified EllisonRamaty functions. The spectral fluences evaluated with the revisited datasets and a new method form a solid basis for more precise studies of the physics of solar eruptive events and the transport of energetic particles in the interstellar medium.
The new spectral fluences will be also useful in various applications of SEP influence on terrestrial effects such as cosmicrayinduced ionisation (e.g., Jackman et al. 2008; Usoskin et al. 2011; Duderstadt et al. 2016), radiation hazards (e.g., Feynman et al. 1993; Jiggens et al. 2014; Mishev et al. 2015; Raukunen et al. 2018), or cosmogenic isotope production (Webber et al. 2007; Kovaltsov et al. 2014; Mekhaldi et al. 2015). The latter is of particular importance for studies of the reference SEP events (Cliver et al. 2014; Usoskin et al. 2020b), and, accordingly, to a more precise assessment of historical extreme solar particle storms (Miyake et al. 2019).
Experimental data is available at the SSDC cosmic ray database: https://tools.ssdc.asi.it/CosmicRays/
This is an adhoc orderofmagnitude estimate based on the possibility of drawing a smooth curve throughout the experimental points within the error bars, or, on other words, keeping the bestfit merit function (Fig. 2) of the order of unity per degree of freedom.
Acknowledgments
Detrended data of GLE recorded by NMs were obtained from the International GLE database https://gle.oulu.fi. PIs and teams of all the groundbased neutron monitors and spaceborne experiments whose data were used here are gratefully acknowledged. This work was partially supported by the Academy of Finland (project No. 321882 ESPERA). Collection of the data from cosmicray/ionospheric experiments before 1989 and the creation of full dataset used in this work was supported by the Russian Science Foundation project no. 207210170. Development of the fluence fitting procedure was supported by the Russian Science Foundation project no. 206746016. The work in the University of Turku was performed in the framework of the Finnish Centre of Excellence in Research of Sustainable Space funded by the Academy of Finland (grant no. 312357). The authors benefited from discussions within the ISSI International Team work (HEROIC team) and ISWATCOSPAR S102 team.
References
 Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2014, Phys. Rep., 544, 323 [Google Scholar]
 Aguilar, M., Ali Cavasonza, L., Alpat, B., et al. 2018, Phys. Rev. Lett., 121, 051101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Anastasiadis, A., Lario, D., Papaioannou, A., Kouloumvakos, A., & Vourlidas, A. 2019, Philos. Trans. R. Soc. A, 377, 20180100 [Google Scholar]
 Asvestari, E., Gil, A., Kovaltsov, G. A., & Usoskin, I. G. 2017, J. Geophys. Res.: Space Phys., 122, 9790 [Google Scholar]
 Band, D., Matteson, J., Ford, L., et al. 1993, ApJ, 413, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Bazilevskaya, G. A., Cliver, E. W., Kovaltsov, G. A., et al. 2014, Space Sci. Rev., 186, 409 [NASA ADS] [CrossRef] [Google Scholar]
 Bindi, V. 2017, Adv. Space Res., 60, 753 [Google Scholar]
 Bruno, A., Bazilevskaya, G. A., Boezio, M., et al. 2018, ApJ, 862, 97 [CrossRef] [Google Scholar]
 Clem, J., & Dorman, L. 2000, Space Sci. Rev., 93, 335 [Google Scholar]
 Cliver, E. W. 2016, ApJ, 832, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Cliver, E. W., Tylka, A. J., Dietrich, W. F., & Ling, A. G. 2014, ApJ, 781, 32 [Google Scholar]
 Cliver, E. W., Hayakawa, H., Love, J. J., & Neidig, D. F. 2020, ApJ, 903, 41 [Google Scholar]
 Desai, M., & Giacalone, J. 2016, Liv. Rev. Sol. Phys., 13, 3 [Google Scholar]
 Duderstadt, K. A., Dibb, J. E., Schwadron, N. A., et al. 2016, J. Geophys. Res. Atmos., 121, 2994 [Google Scholar]
 Ellison, D. C., & Ramaty, R. 1985, ApJ, 298, 400 [Google Scholar]
 Feynman, J., & Gabriel, S. 1990, Sol. Phys., 127, 393 [NASA ADS] [CrossRef] [Google Scholar]
 Feynman, J., Spitale, G., Wang, J., & Gabriel, S. 1993, J. Geophys. Res., 98, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Forbush, S. E. 1946, Phys. Rev., 70, 771 [NASA ADS] [CrossRef] [Google Scholar]
 Goswami, J., McGuire, R., Reedy, R., Lal, D., & Jha, R. 1988, J. Geophys. Res., 93, 7195 [Google Scholar]
 Herbst, K., Grenfell, J. L., Sinnhuber, M., et al. 2019, A&A, 631, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jackman, C. H., Marsh, D. R., Vitt, F. M., et al. 2008, Atmos. Chem. Phys., 8, 765 [Google Scholar]
 Jiggens, P., ChavyMacdonald, M.A., Santin, G., et al. 2014, J. Space Weather Space Clim., 4, A20 [CrossRef] [EDP Sciences] [Google Scholar]
 Jun, I., Swimm, R. T., Ruzmaikin, A., et al. 2007, Adv. Space Res., 40, 304 [Google Scholar]
 King, J. 1974, J. Spacecr. Rocket., 11, 401 [Google Scholar]
 Klein, K.L., & Dalla, S. 2017, Space Sci. Rev., 212, 1107 [Google Scholar]
 Kocharov, L., Pohjolainen, S., Mishev, A., et al. 2017, ApJ, 839, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Koldobskiy, S. A., Kovaltsov, G. A., & Usoskin, I. G. 2018, Sol. Phys., 293, 110 [CrossRef] [Google Scholar]
 Koldobskiy, S. A., Kovaltsov, G. A., Mishev, A. L., & Usoskin, I. G. 2019, Sol. Phys., 294, 94 [CrossRef] [Google Scholar]
 Kong, X., Guo, F., Giacalone, J., Li, H., & Chen, Y. 2017, ApJ, 851, 38 [Google Scholar]
 Kouloumvakos, A., Nindos, A., Valtonen, E., et al. 2015, A&A, 580, A80 [CrossRef] [EDP Sciences] [Google Scholar]
 Kovaltsov, G. A., Usoskin, I. G., Cliver, E. W., Dietrich, W. F., & Tylka, A. J. 2014, Sol. Phys., 289, 4691 [Google Scholar]
 Mekhaldi, F., Muscheler, R., Adolphi, F., et al. 2015, Nat. Commun., 6, 8611 [NASA ADS] [CrossRef] [Google Scholar]
 Mishev, A., Usoskin, I., & Kovaltsov, G. 2013, J. Geophys. Res.: Space Phys., 118, 2783 [Google Scholar]
 Mishev, A., Adibpour, F., Usoskin, I., & Felsberger, E. 2015, Adv. Space Res., 55, 354 [Google Scholar]
 Mishev, A., Usoskin, I., Raukunen, O., et al. 2018, Sol. Phys., 293, 136 [Google Scholar]
 Mishev, A., Koldobskiy, S., Kovaltsov, G., Gil, A., & Usoskin, I. 2020, J. Geophys. Res.: Space Phys., 125, e2019JA027433 [Google Scholar]
 Miyake, F., Usoskin, I., & Poluianov, S. 2019, Extreme Solar Particle Storms: The Hostile Sun (Bristol, UK: IOP Publishing) [Google Scholar]
 Onsager, T., Grubb, R., Kunches, J., et al. 1996, GOES8 and Beyond (SPIE), 2812, 281 [Google Scholar]
 Plainaki, C., Mavromichalaki, H., Laurenza, M., et al. 2014, ApJ, 785, 160 [Google Scholar]
 Press, W., Teukolsky, S., Vetterling, W., & Flannery, B. 2007, Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd edn. (Cambridge: Cambridge University Press) [Google Scholar]
 Raukunen, O., Vainio, R., Tylka, A. J., et al. 2018, J. Space Weather Space Clim., 8, A04 [CrossRef] [EDP Sciences] [Google Scholar]
 Raukunen, O., Paassilta, M., Vainio, R., et al. 2020, J. Space Weather Space Clim., 10, 24 [EDP Sciences] [Google Scholar]
 Reames, D. V. 1999, Space Sci. Rev., 90, 413 [Google Scholar]
 Reedy, R. 1977, in Lunar and Planetary Science VIII, ed. R. Merril (Houston: Lunar and Planetary Institute), 825 [Google Scholar]
 Reeves, G., Cayton, T., Gary, S., & Belian, R. 1992, J. Geophys. Res., 97, 6219 [Google Scholar]
 Sellers, F. B., & Hanser, F. A. 1996, in GOES8 and Beyond, ed. E. R. Washwell (SPIE), 2812, 353 [Google Scholar]
 Shea, M. A., & Smart, D. F. 2012, Space Sci. Rev., 171, 161 [CrossRef] [Google Scholar]
 Tylka, A., & Dietrich, W. 2009, 31th International Cosmic Ray Conference (Lodź, Poland: Universal Academy Press) [Google Scholar]
 Tylka, A., Dietrich, W., & Boberg, P. 1997, IEEE Trans. Nucl. Sci., 44, 2140 [Google Scholar]
 Usoskin, I. G., Kovaltsov, G. A., Mironova, I. A., Tylka, A. J., & Dietrich, W. F. 2011, Atmos. Chem. Phys., 11, 1979 [Google Scholar]
 Usoskin, I., Koldobskiy, S., Kovaltsov, G. A., et al. 2020a, A&A, 640, A17 [EDP Sciences] [Google Scholar]
 Usoskin, I., Koldobskiy, S., Kovaltsov, G., et al. 2020b, J. Geophys. Res.: Space Phys., 125, e27921 [Google Scholar]
 Vainio, R., & Afanasiev, A. 2018, in Particle Acceleration Mechanisms, eds. O. Malandraki, & N. Crosby, Astrophys. Space Sci. Lib., 444, 45 [Google Scholar]
 Vainio, R., Desorgher, L., Heynderickx, D., et al. 2009, Space Sci. Rev., 147, 187 [Google Scholar]
 Van Allen, J. A., Baker, D. N., Randall, B. A., & Sentman, D. D. 1974, J. Geophys. Res., 79, 3559 [Google Scholar]
 Webber, W., Higbie, P., & McCracken, K. 2007, J. Geophys. Res., 112, A10106 [Google Scholar]
Appendix A: Reconstructed integral spectra
Eventintegrated fluences for the analysed GLE events are presented in the plots below. Each plot is similar to Fig. 1 of the main text in both style and notations. The GLE number and date are shown on the top of the plots. Symbols represent data from different sources, as specified in the legend: blue crosses and vertical lines correspond to data from neutron monitors and upper estimates, respectively, along with their 68% confidence intervals (Usoskin et al. 2020a); coloured dots correspond to spaceborne/ionospheric data with the source indicated in the legend (‘this work’ refers to Table 1); and the dark curve depicts the bestfit modified Band function (Eqs. (1) and (2); exact values of the parameters are available in Table 2), while the light grey shading represents the 68% confidence interval for the fits (see Step 3 of Sect. 3).
Fig. A.1.
Integral fluences of SEP reconstructed for GLEs considered in this work (the GLE number and date are given in the header of each panel). Notations are similar to Fig. 1 of the main text. 
Fig. A.1.
continued. 
Fig. A.1.
continued. 
Fig. A.1.
continued. 
Fig. A.1.
continued. 
Fig. A.1.
continued. 
Fig. A.1.
continued. 
Fig. A.1.
continued. 
All Tables
Eventintegrated omnidirectional integral fluences F(> E) (in units of 10^{5} cm^{−2}) obtained here for GLEs # 40−72.
All Figures
Fig. 1.
Example of the fitting of the modified Band function (Eqs. (1) and (2)) to the data points (Sect. 3) for GLE #43, 19Oct1989. Blue crosses and vertical lines denote the highenergy NM data with uncertainties, and the upper limits, respectively (Usoskin et al. 2020a), while red circles correspond to lowenergy data from GOES6 (see Table 1). The grey line with shading depicts the bestfit MBF, along with its 68% confidence interval. 

In the text 
Fig. 2.
Distribution of the parameters and their pairwise correlations for the MBF fitting of the GLE #5 obtained for 10 000 iterations (see Step 3 in Sect. 3). The red dots in the bottom panels correspond to the minimum , while the horizontal red dashed line denotes the 68% confidence interval for the bestfit parameter values. Colour intensity corresponds to the point density. The diagonal panels depict histograms of the parameter values’ distribution. 

In the text 
Fig. A.1.
Integral fluences of SEP reconstructed for GLEs considered in this work (the GLE number and date are given in the header of each panel). Notations are similar to Fig. 1 of the main text. 

In the text 
Fig. A.1.
continued. 

In the text 
Fig. A.1.
continued. 

In the text 
Fig. A.1.
continued. 

In the text 
Fig. A.1.
continued. 

In the text 
Fig. A.1.
continued. 

In the text 
Fig. A.1.
continued. 

In the text 
Fig. A.1.
continued. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.