New Method of Assessment of the Integral Fluence of Solar Energetic (> 1 GV Rigidity) Particles from Neutron Monitor Data

A new method to reconstruct the high-rigidity part (≥ 1 GV) of the spectral fluence of solar energetic particles (SEP) for GLE events, based on the world-wide neutron monitor (NM) network data, is presented. The method is based on the effective rigidity Reff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{\mathrm{eff}}$\end{document} and scaling factor Keff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$K_{\mathrm{eff}}$\end{document}. In contrast to many other methods based on derivation of the best-fit parameters of a prescribed spectral shape, it provides a true non-parametric (viz. free of a priori assumptions on the exact spectrum) estimate of fluence. We reconstructed the SEP fluences for two recent GLE events, #69 (20 Jan. 2005) and #71 (17 May 2012), using four NM yield functions: (CD00 – Clem and Dorman in Space Sci. Rev.93, 335, 2000), (CM12 – Caballero-Lopez and Moraal in J. Geophys. Res.117, A12103, 2012), (Mi13 – Mishev, Usoskin, and Kovaltsov in J. Geophys. Res.118, 2783, 2013), and (Ma16 – Mangeard et al. in J. Geophys. Res.121, 7435, 2016b). The results were compared with full reconstructions and direct measurements by the PAMELA instrument. While reconstructions based on Mi13 and CM12 yield functions are consistent with the measurements, those based on CD00 and Ma16 ones underestimate the fluence by a factor of 2 – 3. It is also shown that the often used power-law approximation of the high-energy tail of SEP spectrum does not properly describe the GLE spectrum in the NM-energy range. Therefore, the earlier estimates of GLE integral fluences need to be revised.


Introduction
While galactic cosmic rays (GCRs) always bombard the Earth's atmosphere, with the intensity being somewhat modulated by solar activity in the course of the 11-year solar cycle (see, e.g. a review by Potgieter, 2013), sometimes sporadic fluxes of solar energetic particles (SEPs) can impinge on the Earth's atmosphere, as caused by solar eruptive events like solar flares or coronal mass ejections (e.g., Vainio et al., 2009;Desai and Giacalone, 2016). During SEP events, fluxes of lower-energy particles (below several hundred MeV) can get enhanced over the GCR background by many orders of magnitude during several to tens of hours. It is important to study such events for different reasons, from purely academic, viz. studying solar eruptive events and probing the inner heliosphere, to very practical ones, since these fluxes pose serious radiation hazards for space-based technologies and even to high-latitude commercial jet flights (Gopalswamy, 2018;Shea and Smart, 2012).
Variability of SEPs is continuously monitored by space missions, such as GOES (Geostationary Operational Environmental Satellites), SoHO (Solar-Heliospheric Observatory), etc. over the last several decades. However, due to natural limitations, most space missions are able to measure mainly the low-energy range of particles, ≤ 100 MeV. A few instruments can detect higher energies, GOES/HEPAD (High Energy Proton and Alpha Detector) can extend the energy range to 700 MeV (P10 channel) and to the integrated flux above 700 MeV (P11 channel). In addition, the SoHO/EPHIN (Electron Proton and Helium Instrument) can measure SEPs up to 500 MeV (Kühl et al., 2017). Two missions are/were able to measure more energetic particles in space, PAMELA (Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics, Adriani et al., 2014) was in operation June 2006 through January 2016, while AMS-02 (Alpha Magnetic Spectrometer, Aguilar et al., 2017) is in operation since 2011. Despite their excellent performance and sensitivity, both these missions are not well suited for SEP monitoring because of their low orbits, whose major fraction is located inside the geomagnetic field and is thus protected from low-energy cosmic particles. Therefore, the only type of detectors able to continuously monitor the energy range above several hundred MeV is a ground-based neutron monitor (NM; see Simpson, 2000). On one hand, NM is an energy-integrating detector unable to directly measure the particle energy/rigidity spectrum. On the other hand, there is the world-wide network of NMs, located in different places with different geomagnetic rigidity cutoffs, which makes it possible to roughly assess the spectrum of energetic particles during SEP events. The key here is the knowledge of the yield function of a NM that quantifies the response of a NM to a monoenergetic unit flux of primary energetic particles on the top of the atmosphere (e.g., Clem and Dorman, 2000). Usually, the spectrum of SEPs is reconstructed parametrically, so that the best-fit parameters of a prescribed SEP spectral shape are defined by fitting the modeled responses of several NMs to the measured ones (e.g., Cramp et al., 1997;Vashenyuk, Balabin, and Stoker, 2007;Mishev, Kocharov, and Usoskin, 2014), explicitly considering also the SEP pitch-angle anisotropy, which can be large in the initial impulsive phase of the event. This method, while allowing for estimate of the time-variable spectral and angular distributions of SEPs during the events, is very laborious and not always stable, and may lead to large uncertainties (Bütikofer and Flückiger, 2015), mostly due to differences in NM yield functions. Of course, the NM-based estimates can be made only for hard-spectrum SEP events, which can initiate atmospheric cascades and be detected by ground-based NMs. This class of events is called Ground-Level Enhancements or GLEs (Poluianov et al., 2017). At present there are known 72 such events (a list can be found at https://gle.oulu.fi).
For practical applications, it is often sufficient to know not the peak flux and its temporal/angular distributions, but the integral fluence (flux integrated over the entire event). Determination of the event fluence is more robust and is usually done under an assumption of the isotropic distribution of SEP particles near Earth. A detailed method for that was proposed by Tylka and Dietrich (2009), who fitted a power law in rigidity tail of the Bandfunction spectral shape to the measured NM responses for most of the GLEs (see Raukunen et al., 2018, -called R18 henceforth). However, this method is parametric, viz. based on an explicit assumption of the power-law spectral shape. In addition, it uses an outdated yield function of Clem and Dorman (2000).
Here we propose a further development of the method by Tylka and Dietrich (2009), by introducing the effective rigidity of a NM, which enables one to make a non-parametric (viz. free of explicit assumptions of the spectral shape) reconstructions of the GLE integral fluence, based on the data from the NM network.

General Approach
A solid method of an assessment of the integral fluence of SEPs for a GLE event, using the data from the world-wide NM network was developed by Tylka and Dietrich (2009) and updated in R18. According to the method, the integral omnidirectional fluence of SEP is assumed as a power-law dependence of the proton rigidity (in the NM energy range of above several hundred MeV): (1) Then, the SEP fluence in the Earth's vicinity can be defined via the measured response of a given ith NM, characterized by the effective vertical geomagnetic cutoff rigidity, P c i , to the SEP event (GLE) as where the first term is the measured (asterisks denote the measured values henceforth) relative increase of the ith NM count rate due to SEPs, integrated over the entire duration of the event, with respect to the background count rate caused by GCR; the second term is the theoretically expected response of an ideal NM to GCR (see Equation 5); and the last term is the scaling (corresponding to the 'correction' term in Tylka and Dietrich, 2009) factor, which is defined as where γ is the spectral index and Y is the yield function of the ideal NM at the atmospheric depth (height) h i , and integration is over the particle's rigidity R. It is important that the value of the scaling factor depends on the SEPs' spectral shape, viz. the spectral index γ . Using data from several NMs with different geomagnetic cutoffs P c i , the best-fit values of F 0 and γ were found by means of the χ 2 method, along with their uncertainties. By applying this method, integral SEP fluences were estimated for most of the GLE events (R18). The method described above contains two important simplifications: first, it ignores the angular distribution of SEPs, and second, it is based on a prescribed spectral shape (power law in rigidity). While the former one is reasonable as the integral fluence is largely defined for the main isotropic phase of the event, the latter assumption makes the method parametric (viz., not the spectrum per se is estimated but parameters of a prescribed shape, without validation whether this shape is applicable) and may lead to a significant systematic uncertainty.
Here we propose a non-parametric method of SEP fluence assessment that is free of strict assumptions of the SEPs' spectral shape. The proposed method leads to a true nonparametric spectrum reconstruction rather than to estimation of some prescribed spectral parameters. We have modified the method of Tylka and Dietrich (2009) by introducing the concept of the effective rigidity R eff and the scaling factor K eff so that the SEP eventintegrated fluence is defined as (cf. Equation 2) The values of R eff i and K eff i are constant characteristics of each (ith) NM. They are calculated once and do not depend on the SEPs' spectral shape. Accordingly, this method allows one to directly relate the NM response to the integral SEP fluence for a GLE event, without any a priori assumptions on its spectral shape. In other words, instead of fitting a model to different observational points, one can plot each NM response as a single point (with error bars) on an F -R diagram. Details are discussed later in Section 2.3. The theoretically expected response of an ideal NM (characterized by the geomagnetic cutoff rigidity P c and the height h of its location) to GLE or GCR can be calculated as where Y j (R, h) is the yield function of the NM (located at height h) for primary cosmicray particles of type j (protons, helium, heavier species), and J j is the differential intensity of primary particles of type j at the Earth's orbit but outside the magnetosphere and atmosphere. The spectrum can be either GCR or SEP. In the case of SEP only protons are considered, so that there is no summation over the different types of primary particles. We used the force-field parameterization (Caballero-Lopez and Moraal, 2004;Usoskin et al., 2005) of the GCR spectra near Earth based on an updated local interstellar spectrum for protons (Vos and Potgieter, 2015) as constrained by recent measurements of Voyager (Stone et al., 2013) and PAMELA (Adriani et al., 2013) in-situ data. The methodology is described elsewhere (Koldobskiy, Kovaltsov, and Usoskin, 2018a;Koldobskiy et al., 2019). The values of the modulation potential φ in the force-field parameterization were calculated for each GLE pre-increase interval from polar NM data, using the methodology described in Usoskin et al. (2017).

NM Yield Function
The NM is a ground-based detector, where secondary nucleonic particles produced in the atmospheric cascade are detected instead of the primary energetic cosmic rays. Thus, the process of the atmospheric cascade needs to be properly modeled in order to study the cosmic-ray variability. Successful efforts in modeling the cosmic-ray induced atmospheric cascade were made over 60 years (e.g., Debrunner and Brunberg, 1968), but only during the last decades it became possible, thanks to the improving computer performance and development of appropriate full-target Monte-Carlo packages, to conduct detailed simulations. This led to the concept of the yield function (YF) of a NM defined as the response of the detector (in terms of counts) to the unit flux of primary cosmic-ray particles outside the Earth's atmosphere and magnetosphere (e.g., Clem and Dorman, 2000). At present, four yield functions for the standard NM64 detector (Simpson, 2000) are used in the research community: • CD00 (Clem and Dorman, 2000) YF was computed numerically as a first detailed Monte-Carlo simulation, using the FLUKA (FLUktuierende KAscade or Fluctuating Cascade - Fassò et al., 2001;Ballarini et al., 2006) package, of the cosmic-ray induced atmospheric cascade in the atmosphere; • CM12 (Caballero-Lopez and Moraal, 2012) YF was empirically constructed based on latitudinal surveys of a NM, and thus defined only for the rigidities below 15 GV, it was extended to higher energies/rigidities theoretically; • Mi13 (Mishev, Usoskin, and Kovaltsov, 2013) YF was computed using the PLANETO-COSMICS GEANT-4 simulation tool (Desorgher et al., 2005, considering, for the first time, the finite lateral size of the atmospheric cascade and the NM's electronic dead time; • Ma16 (Mangeard et al., 2016b,a) YF was also computed using the FLUKA package (release 2011; see Böhlen et al., 2014).
We note that the method of R18 was based solely on the CD00 YF, which may be outdated and does not agree with the NM latitudinal surveys (Caballero-Lopez and Moraal, 2012; Mishev, Usoskin, and Kovaltsov, 2013). Figure 1 shows these YFs and their ratios to the CD00 one. One can see that while they generally agree, within ±10%, with each other for high rigidities > 20 GV, there is a significant discrepancy between them in the low-rigidity range < 10 GV. Ma16 and CD00 YFs appear to be quite close for R < 2 GV, whereas CM12 and Mi13 YFs are significantly, by a factor of two -three, lower. This has been noted by Koldobskiy et al. (2019), who, using the direct in-situ cosmic-ray measurements by PAMELA and AMS-02 over a solar cycle, showed that Mi13 YF is the most consistent with the data, while CD00 and Ma16 tend to overestimate the sensitivity of a NM to low-energy cosmic rays. CM12 YF lies close to Mi13 but gives slightly larger uncertainties. That result was based on an analysis of GCR data related to a higher energy/rigidity range. Here we re-analyze the response of a NM to SEPs, which are less energetic and thus make it possible to study the low-energy part of the NM YF in greater details. We note that YF is usually defined for the intensity of primary energetic particles, J (see Equation 5), while SEPs are typically presented via the omnidirectional flux/fluence F . For the isotropic case, the two quantities are related as F = 4π · J (Grieder, 2001, Chapter 1.6).

NM Effective Rigidity for GLE
The NM is an energy-integrating detector that cannot directly measure the differential energy spectrum of cosmic rays. However, it can record the integral spectrum. The ideal integral particle detector would have a step-like YF (viz., zero below the threshold energy E th , and constant above it). The response of such an ideal detector is directly proportional to the integral flux of primary particles with energy above this threshold energy E th . While the YF of NM is not ideal, it is close to that (very sharp, nearly step-like rise, especially for non-polar NMs followed by a gradual increase roughly proportional to the energy; see Figure 1a). This makes it possible to define the effective energy/rigidity of a NM for the given typical spectrum of primary particles, GCR or GLE. The effective energy/rigidity, E eff /R eff , is such that the integral flux of primary particles above this threshold is nearly proportional to the count rate of the detector, viz. an analog of E th .
The concept of the effective energy/rigidity for NM has been used in application to study GCR variability (e.g., Alanko et al., 2003;. The same concept was applied also to such integral 'detectors' as production of cosmogenic isotopes in the Earth's atmosphere (Kovaltsov et al., 2014; and lunar rocks (Poluianov, Kovaltsov, and Usoskin, 2018). Such an effective energy/rigidity was recently introduced for the sea-level polar NMs for detection of GLEs (Koldobskiy, Kovaltsov, and Usoskin, 2018b).

Figure 1 Panel A:
Yield functions of NM64 sea-level neutron monitors for primary protons, used in this study (denoted in the legend as Mi13 (Mishev, Usoskin, and Kovaltsov, 2013), Ma16 (Mangeard et al., 2016b), CM12 (Caballero-Lopez and Moraal, 2012) and CD00 (Clem and Dorman, 2000)), normalized to the unity at 100 GV rigidity. Panel B: Ratios of the yield functions shown in panel A to that of CD00.
This approach makes it possible to assess the SEP energy-integrated spectrum directly, without making any a priori assumption of the spectral shape. Let us assume that there is an effective rigidity R eff such that the integral fluence of SEPs above this rigidity F (> R eff ), is proportional to the NM response to GLE N GLE : where K eff is a scaling factor for a given NM, which is ideally a constant irrespectively of the strength and energy/rigidity spectrum of the analyzed event (Koldobskiy, Kovaltsov, and Usoskin, 2018b). The expected response of a NM to GLE is calculated using Equation 5.
Here we assume that SEPs causing the GLE consist of protons only, otherwise a term accounting for heavier species needs to be considered. This assumption is valid for large SEP events, where helium composes typically about 0.01 of protons (Torsti et al., 2002), but not for GCR, where helium and heavier species should be explicitly considered. We consider the modified power law in rigidity as the spectral shape of SEP fluence (e.g. Cramp et al., 1997;Mishev et al., 2018): where R is rigidity in GV, γ is the spectral index, and δγ in GV −1 is the rate of the spectrum steepening. In the forthcoming analysis we varied the value of γ in a range from five to nine, which corresponds to typical GLE spectra. The rate of the spectrum steepening δγ was considered in the range from 0 (no steepening of the spectrum, viz. the power law) to 1 GV −1 (strong steepening). The analyzed range of γ and δγ corresponds to a wide range spectra of real SEP events. For each value of γ and δγ from these ranges we calculated the value of K for different values of R, where N GLE is defined using Equation 5. This forms a ribbon in the K-vs-R diagram (see Figure 2). We considered the vertical full-range width of the ribbon, for a given R as K.
Next we found such a value of R, called the effective rigidity R eff , which minimizes the value of δK(R) ≡ K/ K , where K is a mean value of K for the given value of R. The value of K is called the effective scaling factor K eff . Full-scale uncertainties for both values R eff and K eff were evaluated as illustrated by the red error bars in Figure 2. For the standard near-sea-level polar 1NM64 we found, for the Mi13 YF, that R eff = 1.43 +0.05 −0.11 GV and K eff = 5.44 +1.07 −1.45 (cm 2 count) −1 . The value of K eff is defined for 1NM64 throughout the paper.
In a similar way we have defined effective rigidities and scaling factors for different geomagnetic rigidity cutoffs P c and altitudes h, as shown in Figures 3 and 4 and summarized in Tables 1 and 2, respectively, for Mi13 YF. Since the Mi13 yield function was computed only for the sea-level NMs, we used the altitudinal dependence according to Flückiger et al. (2008) applied to the Mi13 yield function. One can see that the effective rigidity R eff is very close to the geomagnetic rigidity cutoff P c for low-and mid-latitude locations (P c > 3 GV) but saturates at 1.3 -1.5 GV (depending on the atmospheric depth) for high-latitude sites. The value of the K eff varies with the geomagnetic cutoff depicting a shoulder at high-latitude Figure 3 The effective rigidity R eff of a NM, using the Mi13 YF, as a function of the cutoff rigidity P c for different atmospheric depths as denoted in the legend. Error bars represent the full-range uncertainties. Values are tabulated in Table 1.

Figure 4
The same as Figure 3 but for the scaling factor K eff . Values are tabulated in Table 2. locations and a nearly exponential decrease with P c for low latitudes and mid-latitudes. This relation is shaped by two different processes, viz. the atmospheric cutoff (particles must possess sufficient energy of a several hundred MeV to initiate an atmospheric cascade reaching the ground) and the geomagnetic cutoff (particles must possess sufficient rigidity to be able to enter the atmosphere). While the geomagnetic cutoff dominates at low latitudes and mid-latitudes, the atmospheric cutoff becomes crucial at high latitudes.
It is important that the effective rigidity and scaling factor are defined robustly for a wide range of the geomagnetic rigidity cutoffs, and they are independent of the exact SEP spectrum, in a reasonable range of parameters. Thus, for each GLE and each NM, one can estimate, using Equation 6, the integral fluence F (> R eff ) of SEPs, and a set of such NMs with different values of R eff makes it possible to perform a non-parametric reconstruction (i.e., one without an explicit assumption on the spectral shape) of the event's integral spectrum.

Test of the Effective-Rigidity Method
In this section we test the R eff method for two well-studied events: GLE #69 and #71. The relative measured response of ith NM to a GLE, integrated over the entire event is defined Table 1 Values of the effective rigidity R eff (in GV) with the full-range uncertainties as a function of cutoff rigidities P c (rows) for different atmospheric depths (columns) as given in the top line. Computations were done for the Mi13 YF. See Figure 2. P c (GV) 700 g cm −2 800 g cm −2 900 g cm −2 1000 g cm where X i is the event-integrated relative intensity of GLE in units of percent-hours , where percents are with respect to the pre-increase NM count rate due to GCR. The value of X i was computed (see column 5 in Table 3) as time integration of the relative (%) increase in time profiles obtained from the international GLE database (https://gle.oulu.fi). These values were defined here and thus may differ from those used in R18. Geomagnetic cutoff rigidities P c were taken from the NMDB database. The theoretically expected response of a NM to GLE N GLE , which enters Equation 6, is whose statistical uncertainties have two independent sources: the accuracy of the determination of X i considered as 1% hr and the statistical uncertainty of the GCR count rate: Another source of uncertainties of the final reconstruction is also the uncertainties in the definition of R eff and K eff (Section 2.3). All the uncertainties were accounted for in the following analysis.

GLE #69, 20 Jan. 2005
As the first test of the method we considered GLE #69, which occurred on 20 January 2005 and was the strongest event over the last cycles and the second strongest ever directly observed. It was an impulsive event with the high peak (≈ 3350 and 4800% for the SOPB and SOPO NMs, respectively) and was characterized by a very strong anisotropy of its initial phase (e.g., Plainaki et al., 2007;Matthiä et al., 2009). Since the effects of anisotropy and temporal evolution of fluxes cannot be considered by the effective-rigidity approach, we can only provide a rough estimate, while the exact spectral reconstruction requires a laborious and complicated method with particle tracing (e.g., Plainaki et al., 2007;Mishev and Usoskin, 2016). For the comparison we use four NM yield functions as mentioned before. The corresponding values of the R eff and K eff were computed for all the YFs, but shown only for Mi13 in Table 3.
Thus, for each NM we can estimate one value of the integral fluence F (> R eff ), that yields a spectral estimate for a set of NMs with different cutoff rigidities (ignoring the anisotropy effects). This was done for the GLE #69 as shown in Figure 5 for the four different YFs, along with the integral fluence (using the Band-function parameterization) estimated for the same event by R18. One can see that the results obtained here by using the effective-rigidity concept (viz., no explicit fitting of the spectral shape) generally agree with the Band-function fitting R18 using NM network and spacecraft data and also ignoring the possible anisotropy, but there are also significant differences. In order to study it in more detail we show in Figure 6 the ratio of the fluences reconstructed here using different YFs from that provided by R18 for GLE #69. Two main features can be observed.
First, the overall shape of the spectrum, reconstructed here, is robust for different YFs and disagrees with the Band-function shape (note that in this energy/rigidity range the Band function is close to a simple power law with the index γ 2 , following the notations of R18): while all spectra merge at the high-rigidity tail of 6 -7 GV and low-rigidity head (< 1 GV), they have an essential excess in the mid-rigidity range of several GV. We note that the wide spread of low-rigidity points is caused by the anisotropic impulsive phase observed by polar NMs. This suggests that the single power-law tail of the Band-function shape is not well representative of such spectrum. Despite the discrepancy in the absolute calibration of the reconstructed spectrum, all YF-based calculations agree in that the spectrum rolls off Table 3 Results of the GLE #69 analysis: NM name (four-letter abbreviations are given according to the list at http://gle.oulu.fi); Geomagnetic cutoff rigidity P c ; Atmospheric depth; Effective rigidity R eff ; Scaling factor K eff (for the standard 1NM64 counter); GLE integral increase X; Estimated SEP integral fluence F SEP (> R eff ). The values correspond to the Mi13 yield function.   11). The integral fluence, reconstructed by Raukunen et al. (2018) in the Band-function form and using the CD00 YF, is depicted with the blue line.

Figure 6
Ratio of the integral fluences reconstructed for the GLE #69 (20 Jan. 2005) here, using different YFs from that provided by Raukunen et al. (2018) using the Band-function approximation.
at high rigidities, above 5 GV, with respect to a purely power-law high-energy tail of the Band function. In order to exclude a possibility that the roll-off is an artifact related to the 'training' of the method on the modified power-law spectral shape (Equation 7), we have tested the method also for the pure power law (δγ = 0) and found that it does not have any significant effect on the values of R eff and K eff . Accordingly, the roll-off of the spectrum is realistic and should be taken into account either using the modified power law (Equation 7) or the Ellison-Ramaty spectra shape (Ellison and Ramaty, 1985), which is a power law over energy with an exponential cutoff.

Figure 7
Reconstructed fluence for GLE #69. Reconstruction based on the effective-rigidity method using the Mi13 YF (red dots, identical to Figure 5a) is confronted with the full reconstruction (blue dashed line; see text) and that from R18 (green solid line, identical to that shown in Figure 5. Second, the mid-rigidity excess ranges for different yield functions between 1.5 -2 for CD00 and Ma16 YFs to 2.5 -4 for CM12 and Mi13 YFs. Overall, the fluence reconstructed here with the use of the CD00 and Ma16 YF is comparable to that of R18 which was also based on the CD00 YF. In the case of Mi13 and CM12, the reconstructed fluence above a few GV rigidity is a factor of 2 -3 higher than that proposed by R18. Considering that CM12 and Mi13 YFs better represent the low-energy part (Koldobskiy et al., 2019), particularly CM12 one which was empirically constructed from the latitudinal survey data, these estimates look more realistic. Figure 7 focuses on the spectral fluence reconstruction for GLE #69. The reconstruction (red dots) is based on the effective-rigidity method introduced here, using the Mi13 YF. The full reconstruction, based on the full inversion of the NM responses, which accounts for the anisotropy of the impulsive phase and temporal evolution of the SEP energy spectrum throughout the events (e.g., Mishev et al., 2018), using the same Mi13 YF, is shown as the blue dashed line. One can observe that, while the shapes of the spectra are nearly identical within the error bars, the full reconstruction agrees within 10% with respect to the method presented here. Since this event was strongly anisotropic during the impulsive phase and both methods use the same Mi13 YF, we consider this as the uncertainty related to the anisotropy of the SEP flux, since the method assumes the SEP flux to be isotropic. On the other hand, the R18 fluence is significantly lower than the full reconstruction, which is caused largely by the use of the CD00 YF (Bütikofer and Flückiger, 2015).

GLE #71: 17 May 2012
We have also applied the effective-rigidity method of spectrum reconstruction to GLE #71 (17 May 2012). While the event was moderate in strength, much weaker (peak 15 -17 % in South Pole, Apatity and Oulu NM count rates) than that #69, it was also highly anisotropic and relatively hard in the impulsive phase (Mishev et al., 2015). Peak fluxes and temporal evolution of SEPs for that event have been studied in detail (e.g. Mishev, Kocharov, and Usoskin, 2014;Plainaki et al., 2014), however, here we focus on the event-integrated fluence. Most importantly, the event was observed by the PAMELA experiment, which directly measured the SEP integral fluence (Bruno et al., 2018). In Figure 8 we confront the spectra reconstructed here, using different YFs, with that measured directly in space. One can see that the spectra, reconstructed using the R eff method, are fully consistent with the direct data (gray area) in the spectral shape, but those based on CD00 and Ma16 YFs tend to underestimate the low-energy range of the spectrum. On the other hand the Band-function approx-  (Bruno et al., 2018). imation (R18) systematically underestimates, by the factor 2 -3, the fluence in the rigidity range of about 2 GV (energy about 1 GeV). The full NM-based reconstruction based on Mi13 YF (dashed red line) demonstrates the excellent agreement with the spectrum reconstructed here, and a fair agreement, within the uncertainties, with the directly observed one.

Summary
Here we presented a new fast method for assessment of the high-rigidity part (above 1 GV) of the spectral fluence of SEPs for GLE events based on the data from the world-wide NM network. The method is based on the effective rigidity R eff and scaling factor K eff , calculated for each NM, where the SEP integral fluence is directly related to the NM response to the event. This method is non-parametric so that it provides an estimate of the spectrum without any explicit assumption of the spectral shape. This is an important difference to most of the other methods which fit a prescribed spectral shape into the data by finding the bestfit parameters for that shape, often without a proper validation. This method is simple and fast avoiding laborious computations needed for the full reconstruction, but it neglects the possible anisotropy of the impulsive phase of the event.
We tested the method for two recent GLE events, #69 (20 Jan. 2005), which was the second strongest observed one and had a very anisotropic SEP distribution during the impulsive phase of the event, and #71 (17 May 2012), which was a moderate, but also strongly anisotropic event, whose spectrum was directly measured by PAMELA instrument. For both events we found a good qualitative agreement between the spectral shape of the reconstructed here, on one hand, and that for fully reconstructed spectra, on the other hand. By confronting the spectra assessed here with those obtained via a full NM-based reconstruction, we estimated the uncertainty related to the anisotropy as being within 10%.
Next, we compared our present reconstructions with the directly measured spectrum for GLE #71. Only the reconstructions based on Mi13 and CM12 YFs appear quantitatively consistent with the measurements, while the results based on CD00 and Ma16 lead to an underestimate of the spectrum by a factor 2 -3 in the lower-rigidity range. This is most probably related to the fact that the latter two YFs tend to overestimate the NM response to lower-energy particles, as was found by Koldobskiy et al. (2019) from an analysis of GCR spectra measured directly by AMS-02 and PAMELA experiments. In particular, the CD00 YF is known to fail reproducing the observed NM latitudinal surveys (Clem and Dorman, 2000;Caballero-Lopez and Moraal, 2012;Mishev, Usoskin, and Kovaltsov, 2013) implying that it overestimates the sensitivity of a NM to the lower-energy particles. Since the model by R18 is based on the CD00 YF, it may underestimate the SEP fluence by a factor of two to three. It is also shown that the power-law approximation of the high-energy tail, used in particular in the Band-function parameterization, does not properly describe the form of the GLE spectrum in the NM-energy range. Therefore, the earlier estimates of GLE integral fluences need to be revised. A proper reconstruction of the SEP integral fluence for the known GLE events is planned for a forthcoming work.