Application of the Verified Neutron Monitor Yield Function for an Extended Analysis of the GLE # 71 on 17 May 2012

Intense solar activity was observed in May 2012. A notable ground level enhancement (GLE) was registered on 17 May 2012 by several space‐borne instruments as well as on ground by neutron monitors (NMs). This event is known as GLE # 71. Here, we derived the spectral and angular characteristics, and apparent source position of the solar protons during the GLE # 71, employing verified newly computed NM yield function and sophisticated unfolding procedure. We considerably improved the previously derived information about the spectra and angular distribution, namely, the precision, time span, and time resolution of the derived characteristics, specifically during the event onset and late phase. A comparison with direct measurements, with the Payload for Antimatter Matter Exploration and Light‐nuclei Astrophysics (PAMELA) experiment, of the particle fluence was performed, and good agreement between NM and direct space‐borne data analysis was achieved. Subsequently, we computed the effective dose rates in the polar region at several altitudes during the event using the derived rigidity spectra of the solar protons as a reliable input for the corresponding radiation model. The contribution of the galactic cosmic rays and solar protons to the exposure is explicitly considered. We computed the integrated exposure during the event and discussed the exposure of crew members/passengers to radiation at several altitudes.


Introduction
A systematic study of high-energy particles of solar origin, known as solar energetic particles (SEPs), allows one to reveal their acceleration mechanism(s). Besides, it is important to study such events in order to understand their occurrence probability, important not only in solar physics but also for space weather (e.g., Pulkkinen, 2007, and references therein). Occasionally, fluxes of SEPs result from eruptive events on the Sun, such as solar flares and/or coronal mass ejections (CMEs), can enter the Earth's atmosphere (e.g., Cliver et al., 2004;Desai & Giacalone, 2016;Reames, 1999). Usually, fluxes of lower-energy SEPs (below several hundred MeV) are greater than that of corresponding galactic cosmic ray (GCR) background by several orders of magnitude. SEP events can last several hours, even in some cases tens of hours.
When the energy of SEP spectra exceeds about 1 GeV/nucleon, it results in the development of a complicated hadron-electromagnetic-muon cascade in the Earth's atmosphere due to interactions of the accelerated solar ions with the atmospheric constituent(s) (for details see Dorman, 2004;Grieder, 2011, and references therein). Normally, the low-energy part of the SEPs is absorbed in the atmosphere, dissipating its energy dominantly by ionization. In the same time, particles with energy of 1 GeV/nucleon and greater produce secondaries, which can reach the ground level, can be registered by the ground-based neutron monitors (NMs), resulting in the so-called ground level enhancements (GLEs) (e.g., Aschwanden, 2012;Poluianov et al., 2017;Shea & Smart, 1982;Stoker, 1995). In this case, an increase in the count rates of ground-based detectors, for example, NMs, can be observed (Hatton, 1971;Simpson, 2000).
Precise measurements of GLE particles can be performed with space-borne instruments. Over the years, space probes, such as Geo-stationary Operational Environmental Satellites (GOES) or Solar-Heliospheric Observatory (SoHO), have been used for continuous recording of SEP events (e.g., Domingo et al., 1995;Menzel & Purdom, 1994). Yet most of the space-borne instruments are focused on measurements mainly on the low-energy range of SEPs, that is, below some 100 MeV/nucleon. Therefore, they put apparent constraints for the precise study of GLEs. Two instruments have notably enhanced energy channels and can detect higher energies, for example, GOES/High Energy Proton and Alpha Detector (HEPAD) with energy range to about of 700 MeV/nucleon and SoHO/Electron Proton and Helium Instrument (EPHIN) up to 500 MeV/nucleon (for details see Kühl et al., 2017). Besides, Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics (PAMELA) which was in operation in the period June 2006 to January 2016 and Alpha Magnetic Spectrometer (AMS-02) in orbit since 2011, which were/are designed to measure directly more energetic particles, are able to register high-energy SEPs (see for details Adriani et al., 2016;Aguilar et al., 2015). However, the low orbiting probes are not well suited for SEP study, because of the limited sensitivity to low-energy cosmic rays inside the geomagnetosphere.
Hence, it is necessary to employ ground-based instruments; hereby GLEs can be studied using the worldwide NM network (Hatton, 1971;Mavromichalaki et al., 2011;Papaioannou et al., 2014;Simpson et al., 1953;Stoker et al., 2000). The sensitivity of an NM to the primary CR is determined by the geomagnetic and atmospheric shielding. The rigidity cutoff is a function of the location of the instrument, while the thickness of the atmospheric layer above a given instrument, namely, the altitude above sea level (a.s.l.), determines the atmospheric cutoff. The high-altitude polar NMs, for example, located at the South Pole standard and lead free NMs named SOPO/SOPB, respectively, and located on the Antarctic plateau at Dome C station named DOMC/DOMB, are most sensitive to SEPs (e.g., Poluianov et al., 2017). Naturally, the global NM network situated inside the magnetic field of the Earth represents a giant spectrometer. NMs located at various rigidity cutoffs are sensitive to different parts of CR spectrum and particle's arrival direction because the acceptance of each station at the border of the magnetosphere is a function on its location, particle rigidity, and the angle of incidence of the incoming SEPs. Besides, NMs located at different altitudes a.s.l. reveal different response versus primary SEPs. Yet the majority of those studies were fulfilled without the necessary systematic, for example, GLEs are usually analyzed using different model assumptions and methods resulting in bias and caveats of the derived spectra and angular distribution (e.g., Bütikofer & Flückiger, 2013Desai & Giacalone, 2016;Shea & Smart, 2012).
Here, we perform a revised analysis of 17 May 2017 GLE # 71, using newly computed and verified NM yield function for several altitudes and compare the derived fluence with direct measurements by the PAMELA space probe. A space weather application of the derived SEP characteristics is also demonstrated, namely, assessment of the exposure to radiation at commercial jet flight altitude(s).

Analysis of GLE # 71 Using NM Records
GLE # 71 occurred on 17 May 2012 as a result of eruptive processes in the active region NOAA 11476, which produced a CME and a moderately strong flare (class M5.1) at 01:25 UT. The active region was on the west side at N07 W88, so that the Earth was well magnetically connected to the eruption core. The event followed moderately strong M5.1 solar fare, which started at 01:25 UT and peaked at 01:47 UT. The solar flare was accompanied by a halo CME with estimated linear speed of 1,582 km⋅s −1 (e.g., Papaioannou et al., 2014). On the ground, the worldwide NM network recorded a weak enhancement of the count rates of several stations around 01:50 UT. Notably greater signals have been recorded by APTY, OULU, and SOPO/SOPB NMs (Figure 1), while the bulk of NMs registered only marginal count rate increases, which implied a large anisotropy of the arriving SEPs, specifically during the event onset. The SEP flux was observed above the background of GCRs for about 2 h. This event was extensively studied (e.g., Anastasiadis et al., 2019;Battarbee et al., 2018;Kocharov et al., 2018;Mishev, Kocharov, et al., 2014;Plainaki et al., 2014). A good summary of the observations is given in Papaioannou et al. (2014). However, despite the general agreement, several differences in the derived SEPs characteristics were reported, mostly due to the application of different assumptions, methods of data analysis, and models (for details see Bütikofer & Flückiger, 2015).
The spectral and angular characteristics of GLEs in the energy range ∼0.3-20 GeV/nucleon can be derived by modeling of the global NM network response and corresponding optimization of the model parameters over  (Table 1), the corresponding color lines, and numbers indicate the NM stations and asymptotic directions, which are plotted in the rigidity range ∼1-3 GV, while SOPO is plotted in the rigidity range ∼0.7-3 GV. The cross corresponds to the interplanetary magnetic field direction on the basis of ACE satellite measurements. The lines of equal pitch angles relative to the derived anisotropy axis are plotted for 15 • , 30 • , 45 • , 60 • , and 90 • for sunward directions (solid lines) and 165 • , 150 • , and 120 • for anti-Sun direction (dashed lines). the recorded count rate increases. Several methods have been developed over the years (e.g., Bombardieri et al., 2006;Cramp et al., 1997;Shea & Smart, 1982;. Here we employed method, whose detailed description and applications are given elsewhere and briefly summarized below (Mishev, Kocharov, et al., 2014;Mishev & Usoskin, 2016;. The NM response to primary CRs can be computed by integral of the product of the primary CR spectrum J(P, t) with the NM yield function S(P, h). Thus, the count rate of an NM at a given altitude h and time t is expressed as where P c is the local geomagnetic rigidity cutoff (Cooke et al., 1991;Smart et al., 2006), h is the atmospheric depth (or altitude), accordingly, S i (P, h) (m 2 sr) represents the NM yield function for primaries of particle type i (protons and/or -particles), and J i (P, t) (GV m 2 sr s) −1 is the rigidity spectrum of primary particle of type i at time t. The integration in Equation 1 is over the particle's rigidity, while the summation is over the types of the primary CRs.
Here, we modeled the NM response using a new NM yield function computed for several altitudes by Mishev et al. (2013Mishev et al. ( , 2020, which is fully consistent with the experimental latitude and altitude surveys and was recently validated by achieving good agreement between model results and space-borne with AMS-02 and ground-based NM measurements (for details see Gil et al., 2015;Koldobskiy, Bindi, et al., 2019;Lara et al., 2016;Mishev et al., 2020;Nuntiyakul et al., 2018;Usoskin et al., 2017). Therefore, the response of each NM was modeled with a yield function corresponding to the exact station's altitude a.s.l., which allowed us to reduce model uncertainties related to the application of the double-attenuation-lengths method, that is, normalization of high-altitude NM count rates to the sea level, as well as to employ different response(s) of NMs located at similar geographic positions, but at different altitude(s) (e.g., IRKT, IRK2, and IRK3).
Thus, through several consecutive steps: retrieving the necessary NM records, K p index data, and direction of the interplanetary magnetic field (if available); computation of asymptotic viewing directions and cutoff rigidities of all NMs used for the data analysis; making an initial guess of the optimization procedure; and performing optimization using modeled and recorded NM responses over a selected space of unknown parameters, we are able to analyze the GLE # 71 and to derive the spectral and angular characteristics of SEPs, as well as the apparent source position. A flow diagram summarizing all the steps, data inputs, and used models for GLE analysis and subsequent computation of the exposure to radiation, accordingly discussion related to initial guess influence on the optimization procedure, is given in .
Here, the magnetospheric computations were performed with the MAGNETOCOSMICS code, explicitly considering the measured K p index prior and during the event (Desorgher et al., 2005). A combination of the IGRF geomagnetic model as the internal field model (Thébault et al., 2015) and the Tsyganenko 89 model as the external field (Tsyganenko, 1989) was employed, which allows one to compute straightforwardly the rigidity cutoffs and asymptotic directions (see Figure 2) of NM stations used in our analysis (Kudela & Usoskin, 2004;Kudela et al., 2008;Nevalainen et al., 2013). The optimization was carried out by employing a method similar to Levenberg (1944) and Marquardt (1963), namely, by additional regularization (for details see Aleksandrov, 1971;Golub & Van Loan, 1980;Golub et al., 1999). Our method allows one to derive a reliable solution even in the case of ill-posed problem(s), usually resulted from complicated spectral and angular distributions (Aster et al., 2005;Mavrodiev et al., 2004;Mishev et al., 2005Tikhonov et al., 1995).
Using detrended records from 28 NMs (the details given in Usoskin et al., 2020) retrieved from the GLE database (for details see https://gle.oulu.fi; the full list of the NM stations with the standard abbreviations is given in Table 1), we derived the spectral and angular characteristics of SEPs during the GLE # 71 on 17 May 2012. In contrast to the previous studies, here we expanded considerably the time interval and the precision of the derived information (Balabin et al., 2013;Mishev, Kocharov, et al., 2014). The best fit is derived using a modified power law rigidity spectrum similarly to Cramp et al. (1997), Bombardieri et al. (2006), and Vashenyuk et al. (2008): where J || (P) is the particle flux with given rigidity P in GV arriving from the Sun along the axis of symmetry, which is defined by the geographic coordinate angles Ψ and Λ (latitude and longitude). In Equation 2, is the power law spectral exponent at rigidity P = 1 GV; accordingly, is the rate of the spectrum steepening. The angular distribution of the arriving SEPs is depicted by complicated pitch angle distribution (PAD) with shape similar to that considered by Cramp et al. (1997), namely, superposition of two Gaussians: where is the pitch angle, 1 and 2 are parameters corresponding to the width of the PAD, and B and ′ are parameters corresponding to the contribution of the second Gaussian, including direction nearly opposite to the derived axis of symmetry.
The merit function (normalized residual) determining the fit quality is given by where m is the number of NM stations and ΔN i N i is the relative NM count rate increase of the ith NM station. In general, steady convergence and robust solution are achieved when  ≤ 5% . Moreover, additional criteria have been used similarly to . Note that for weak events  can be about 8-14%.
Modeled and measured count rate increases of NMs with the maximum response for selected periods during the event are shown in Figure 3. The quality of the reconstruction is similar for the other NM stations used in the analysis. The quality of the fit is also depicted with the contour plot of the residual  for the best fit solutions against geographic latitude and longitude (Figure 4). We emphasize that an erroneous determination of the anisotropy axis direction normally leads to an incorrect assessment of the spectral and angular characteristics of SEPs (for details see Bütikofer & Flückiger, 2015). Here, one can see that the derived fit has reasonable quality, despite the assumed model of a complicated angular distribution of SEPs.
Accordingly, several derived spectra and PADs are presented in Figure 5; the details are given in Table 2. We can distinguish three phases of the event: initial (01:50-02:25 UT), when a relatively hard spectrum, a constant increase of SEP flux, and complicated PAD were derived; main phase (02:25-03:05 UT), when a steady softening of the SEP spectra and decrease of the SEP flux were derived; accordingly, a tendency of decrease of the steepening of the spectrum was observed; and late phase of the event (after 03:05 UT), when a pure power law SEP spectrum and wider a nearly isotropic in the very late phase PAD were revealed. The phases of the event can be quantitatively related to different acceleration stages of SEPs (e.g., Anastasiadis et al., 2019;Kocharov et al., 2017). An illustration of the time evolution of the derived spectra and PAD throughout the event is given in the supporting information.
We emphasize that using verified NM yield function and the corresponding sophisticated method for the data analysis, we considerably improved the derived information about the SEP spectra and PAD (time span,   Table 2. The black solid line denotes the GCR flux, which corresponds to the time period of the event occurrence. Time (UT) refers to the end of the corresponding 5-min interval over which the data are integrated. time resolution, and precision), specifically during the initial and late phases of the event, compared to the previous studies (e.g., Mishev, Kocharov, et al., 2014). Moreover, we reduced the derived SEP and PAD uncertainties and performed verification of the fit, not considered thoroughly before (Figure 4) (Mishev, Kocharov, et al., 2014).
Very good agreement of the derived fluence during GLE # 71 with the direct measurements made by PAMELA space probe was achieved (Bruno et al., 2018) (Figure 6); the details are given in . We emphasize that PAMELA detectors allowed one to distinguish the SEP signal on the GCR background only for rigidities less than about 2 GV, while information about higher-energies SEP can be retrieved only using the global NM network and/or AMS-02 data. Note that the direct measurements were fitted with Ellison-Ramaty spectral shape (e.g., Ellison & Ramaty, 1985), while for our model based on NM data, Equation 2 was employed, which exhibited smaller roll-off of the spectrum. Besides, good agreement with SOHO/EPHIN measurements fitted with pure power law, specifically in the energy range 300-700 MeV/nucleon, was also achieved (e.g., Kühl et al., 2015). Moreover, a very similar complicated angular distribution of the SEPs was revealed by direct PAMELA observations (the details are given in Adriani et al., 2015). Therefore, we claim that our method of analysis using a

10.1029/2020SW002626
verified NM yield function corresponding to the specific altitude of the station a.s.l. within reliable optimization is validated with good agreement by direct measurements.

Space Weather Effects During GLE # 71
During GLEs increased intensity of CRs of solar origin is observed (e.g., Miroshnichenko, 2018, and references therein). This can result in important space weather issues, specifically at flight altitudes (e.g., Shea & Smart, 2000, 2012, and references therein). During polar intercontinental flights, where the geomagnetic shielding is marginal, aircrews are exposed to enhanced radiation field, which can be significant during GLEs (Shea & Smart, 2000).
Several models aiming at assessment of the secondary CR radiation, namely, exposure to radiation, henceforth exposure, on aircrews have been developed in the last two decades (e.g., Copeland, 2017;Ferrari et al., 2001;Latocha et al., 2009). While a general agreement within about 10-20% between different models is achieved, specifically for the assessment of GCRs effects (for details see Bottollier-Depois et al., 2009;Yang & Sheu, 2020, and references therein), an important uncertainty, up to an order of magnitude, in the computation of the exposure during GLEs is reported (for details see Bütikofer & Flückiger, 2013. We note that during GLEs, the exposure is a superposition of the contributions of GCRs and SEPs. Hence, a nontolerable by the community discrepancy mostly due to the considered SEP spectra, that is, input for the models, is the main concern for assessment of this specific space weather effect (for details see Bütikofer & Flückiger, 2015;Tuohino et al., 2018, and references therein).
Here, we employed the derived spectra of GLE # 71 (Table 2) and an appropriate exposure model, which agrees well with reference and experimental data; the details are given elsewhere  to assess the effective dose at several altitudes during this event. In this study, we employed a numerical model for computation of the effective dose, that is based on precomputed yield functions, whose full description and applications are given elsewhere . In the model, the effective dose rate at a given atmospheric altitude (depth) h induced by a primary CR particle is computed by integral of the product of the CR particle spectrum with the corresponding yield function: where J i (T) is the differential energy spectrum of the primary CR arriving at the top of the atmosphere for the ith component of CRs (proton or -particle) and Y i is the corresponding effective dose yield function. The integration is over the kinetic energy T above T(P cut ), which is defined by the local cutoff rigidity P cut and over the sold angle Ω.
The effective dose yield function Y i is defined as where C j (T * ) is the fluence to effective dose conversion coefficient for secondary particle of type j (neutron, proton, , e − , e + , − , + , − , + ) with energy T * and F i, j (h, T, T * , , ) is the fluence of secondary particles of type j, produced by a primary CR particle of type i (proton or -particle) with a given primary energy T arriving at the top of the atmosphere from zenith angle and azimuth angle .
In the model, the conversion coefficients C j (T * ) are considered according to Petoussi-Henss et al. (2010). We emphasize that employment of different conversion coefficients C j (T * ) (e.g., ICRP, 1996) would lead to an increase of the assessed exposure of about 20%, which is considerably below the other model uncertainties (e.g., Copeland & Atwell, 2019;Yang & Sheu, 2020).
Using the derived spectra during GLE # 71 (Table 2 and Equation 5), we computed the effective dose rate at several altitudes, the results are presented in Figure 7, and the details are given in Table 3. The contribution of GCRs to the exposure at typical flight altitudes of about 30-50 kft a.s.l. varies between 5.7 and 14.4 μSv h −1 , while the superimposed exposure due to SEPs and GCRs is in the range 7.9-59.4 μSv h −1 . We stress that the exposure during GLEs can reach peak values for a relatively short period, mostly corresponding to the peak SEP flux. Usually, for assessment of the exposure, computations are performed over a period corresponding to the whole event or for a period corresponding to the flight duration in the polar region.
In Figure 8, we present the distribution of the effective dose over the globe at an altitude of 50 kft a.s.l., integrated over the whole GLE # 71. This computation corresponds to the most conservative scenario, integration over the whole time span of the event (about 2 h) at an altitude of 50 kft (≈15,200 m a.s.l.), which is representative for a polar flight considering the specifics of the polar atmosphere (for details see Mironova et al., 2015, and references therein;Mishev & Velinov, 2014). The computation of the rigidity cutoff over the globe was performed with the MAGNETOCOSMICS code (Desorgher et al., 2005). Despite the derived anisotropy, we conservatively assumed an isotropic angular distribution of the GLE particles similarly to Copeland et al. (2008) and . Naturally, the exposure is maximal in the polar region, where the received dose is slightly greater than that due to sole GCRs in quiet radiation conditions. An illustration of the time evolution of the exposure at an altitude of 35 kft a.s.l., accordingly the altitude dependence of the exposure during the peak phase of the event, is given in the supporting information.
The exposure was maximal in a polar region, rapidly diminishing in the low-latitude region. Moreover, the exposure dropped down significantly at lower altitudes. The integrated over the event exposure at an altitude of 50 kft a.s.l. was ≈50 μSv, which is less than that of the double due to GCRs. At lower altitudes, the integrated exposure was about 50% less, namely, ≈25 μSv at an altitude of 35 kft a.s.l. Here we would like Note. The contribution of the GCRs to the exposure is given separately in the first row. to point out that the pilots received annually around 3 mSv (e.g., Bennett et al., 2013). Therefore, relatively weak event such as GLE # 71 on 17 May 2012 does not represent an important space weather issue.

Summary and Discussion
In this study, we derived the rigidity spectrum and PAD of high-energy SEPs in the course of GLE # 71 on 17 May 2012. The analysis was performed employing a robust optimization method, using data from the global NM network, and a newly computed and verified NM yield function for the exact altitude of the stations. In contrast to the previous studies, we expanded considerably the time interval and the precision of the obtained information. Good agreement of the derived SEP fluence with direct measurements by the PAMELA space-borne detector was achieved. Moreover, reasonable agreement with other space-borne instruments is also observed. Thus, this specific method for data analysis of strong SEP events using NM records is verified at least for a single event and gives a good basis for further quantification of solar proton acceleration (for details see Kocharov et al., 2018, and references therein).
On the basis of the reconstructed spectra and employing a Monte Carlo-based model, we assessed the exposure for crew members/passengers at several typical commercial jet flight altitudes as well as at high-mountain ground level, specifically in the polar region. The altitude and time evolution throughout the event of the exposure are presented in great detail. Map of the global distribution of the event integrated effective dose at an altitude of 50 kft is also shown, which allowed one to assess the received exposure during the event to be about 50 μSv at an altitude of 50 kft and 25 μSv at an altitude of 35 kft a.s.l., respectively.
We emphasize that exposure during GLEs is studied usually a posteriori, since they occur sporadically and naturally differ from each other in spectra, duration, angular distribution of SEPs, and magnetospheric conditions, which govern the space weather impact. Therefore, exposure during GLEs is studied case by case. Yet the exposure at cruise flight altitudes during GLEs can be notably enhanced (e.g., Tuohino et al., 2018), one can see that in the case of GLE # 71 on 17 May 2012, the received exposure is not so significant. While the analysis of a GLE event and assessment of the corresponding space weather effects follow the described here standard procedure, in reality detailed studies of SEPs using jointly ground-based and space-borne instruments are performed on case-by-case basis, since they differ not only in characteristics but also in the amount and quality of the available experimental records (for discussion see Jiggens et al., 2019;, and references therein). The presented here analysis GLE # 71 and resulting exposure at flight altitude is comparable with other moderate GLEs.