Current status and possible extension of the global neutron monitor network

The global neutron monitor network has been successfully used over several decades to study cosmic ray variations and fluxes of energetic solar particles. Nowadays, it is used also for space weather purposes, e.g. alerts and assessment of the exposure to radiation. Here, we present the current status of the global neutron monitor network. We discuss the ability of the global neutron monitor network to study solar energetic particles, specifically during large ground level enhancements. We demonstrate as an example, the derived solar proton characteristics during ground level enhancements GLE $\#$5 and the resulting effective dose over the globe at a typical commercial jet flight altitude of 40 kft ($\approx$ 12 200m) above sea level. We present a plan for improvement of space weather services and applications of the global neutron monitor network, specifically for studies related to solar energetic particles, namely an extension of the existing network with several new monitors. We discuss the ability of the optimized global neutron monitor network to study various populations of solar energetic particles and to provide reliable space weather services.

ability. The energy loss of a neutron during elastic collision increases with decreasing the atomic mass, therefore the moderator is selected to contain a significant amount of low mass nuclei e.g. Hydrogen. The lead producer, surrounds the moderator, aiming production of more neutrons by inelastic interactions in a thick target. Therefore, the producer is built by high atomic mass material. The outermost layer of the NM represents a moderator, namely the reflector, which has a double purpose: first, it rejects the low energy neutrons result from interaction(s) of the very local surroundings from penetrating in the NM, secondly, it allows to keep the produced in the lead neutrons inside the monitor.
The introduction of a NM as a continuous recorder of CR intensity followed the design by Simpson et al. (1953). During the International Geophysical Year (IGY) 1957-1958 a 12 tube neutron monitor was constructed, but other configurations have been also used (Simpson, 1957;Shea and Smart, 2000b;Simpson, 2000). The IGY neutron monitor was used world-wide as a detector to study CR variations. Lately, in the mid-sixties, the design of the IGY NM was optimized resulting in increased counting rate (Hatton and Carmichael, 1964;Carmichael, 1968;Hatton, 1971). This second generation of NM design is known as NM64 or supermonitor (for details see Simpson, 2000;Stoker et al., 2000, and references therein). Recently, mini-NMs have been installed at several stations, exhibiting good performance, specifically at low cut-off rigidity and high-altitude locations (Poluianov et al., 2015).
The count rate of a NM provides reliable information about CR flux variations at the top of the Earth's atmosphere, both long-term (e.g. the 11-year sunspot cycle and the 22-year solar magnetic cycle), and short-term as Forbush decreases, diurnal CR variations and transient phenomena such as recently observed anisotropic cosmic ray enhancements (for details see Gil et al., 2018). NMs data are used to derive spectral and angular characteristics of GLEs and high-energy SEPs, specifically in the high-energy range and over the whole event timespan (e.g. Shea and Smart, 1982;Cramp et al., 1997;Bombardieri et al., 2006;Vashenyuk et al., 2006b;Mishev et al., , 2018b. The information retrieved from NMs is essential to assess important topics related to space weather, such as exposure to radiation of aircrew(s), henceforth exposure, and the influence of CRs on atmospheric chemistry (e.g. Bazilevskaya et al., 2008;Vainio et al., 2009;Usoskin et al., 2011;Mironova et al., 2015).
In order to offer a useful tool, specifically for space weather purposes, the global NM network shall provide coverage of the entire sky and real-time data access (e.g. Mavromichalaki et al., 2011). Here, we discuss the current status of the global NM network and present a plan for its extension, aiming an optimization of its performance as a space weather tool.

Performance and current status of the global neutron monitor network
Over the years, it was demonstrated that the global NM network is a powerful tool to study primary CR variations, transient phenomena, SEPs, and to provide data, which form an important input for space weather applications (e.g. Bütikofer, 2018b). In reality, the NM network as a whole, together with the geomagnetic field, represents a giant spectrometer, which allows one to observe the variations of the primary CRs, because NMs placed at various rigidity cut-offs are sensitive to different parts of CR spectrum. In addition, multi-vantage-point registration, specifically of SEPs, makes it possible to reveal the anisotropy of CRs in the vicinity of Earth, since the viewing cone of each station is a function on its location, particle rigidity, and angle of incidence of the arriving particle.
The global NM network presently consists of about 50 stations spread over the world, for details see Fig.1, where the NM stations with the corresponding rigidity cut-off are shown (Moraal et al., 2000;Mavromichalaki et al., 2011). Here the computation of the rigidity cut-off over the globe was performed with the MAGNETOCOSMICS code using the IGRF magnetospheric model corresponding to the epoch 2015 (Desorgher et al., 2005;Thébault et al., 2015).
The sensitivity of a NMs to primary CR is determined by the geomagnetic and atmospheric shielding. The rigidity cut-off is a function of the geomagnetic location of the monitor, while the thickness of the atmospheric layer above a given NM determines the atmospheric cut-off, since the primary CR must possess minimum energy (≈430 MeV nucleon −1 for the sea level) to induce an atmospheric cascade, whose secondary particles reach the ground (e.g. Grieder, 2001). The atmospheric cut-off plays an important role in polar NMs, specifically those at the sea level, since the geomagnetic rigidity cut-off is small in the polar regions. Several high-altitude polar NMs, e.g. SOPO/SOPB and DOMC/DOMB are more sensitive to primary CR, specifically SEPs, than mid-and high rigidity cut-off NMs. Therefore, the rigidity range of the global NM network is determined by the atmospheric cut-off at polar regions, which posses the lower rigidity cut-offs, accordingly by the highest geomagnetic cut-off at about 17 GV in the magnetic Equator.
Besides, polar NMs possess better angular resolution, which is important for the GLE analysis. With this in mind, a concept of the spaceship Earth, an optimized network consisting only of polar stations was proposed (Bieber and Evenson, 1995). However, one can see that the present NM network provides good coverage of arrival directions, and almost symmetric response (see Fig.1), but several gaps exist, as discussed below.

Extension of the global NM network
High-energy CRs are not deflected by the Earth's magnetic field. Therefore, NMs can record highenergy CRs, propagating almost along a straight line, determined by the latitude and longitude of the geographic position of the station. The situation is more complicated for low-energy particles, which are more strongly deflected. Thus, a NM is characterized by his asymptotic direction, i.e., the direction from which particles impinge on a given point in the atmosphere of the Earth arriving at the border of the magnetosphere. It depends on the location, particle incidence angle and rigidity (for details see Bütikofer, 2018a,b, and references therein). As a result, a NM is sensitive to a certain segment of the sky. While for the continuous recording of the isotropic GCR intensity, the asymptotic direction of a NM is not important, it is crucial for registration of GLEs, because SEPs reveal essential anisotropy, specifically during the event onset. Therefore, gaps in asymptotic directions of the global NM network can compromise the registration of GLEs, accordingly the corresponding analysis and alert services. The present situation of operational polar NMs allows one to derive a comprehensive picture of GLE characteristics and provide alert systems (see Figs.1, 2 and Table 1). However, a gap in the asymptotic directions of Arctic NMs is observed, precisely in the longitude range 130-250 • in the northern polar region. While the South polar NMs provide good coverage of the sky, those at North exhibit gaps (Fig.2). One can see that the majority of NMs are looking towards the Equator, i.e., NMs in the North hemisphere are looking southward, while those in Antarctica except DOMC, are looking northward. In addition, as was recently discussed, the high-altitude polar NMs such as DOMC and VSTK are more sensitive to SEPs . Therefore, there is a need for a NM, which is a counterpart of DOMC, i.e., high-altitude, low rigidity cut-off NM located in the North hemisphere close to the geomagnetic pole, as well as several stations to cover the gap and/or to improve the sensitivity, specifically in a low energy range.
For example, if a GLE with narrow angular distribution of the particle flux occurs (see the pitch angle distribution in Fig.3) with anisotropy axis located in the polar region of the northern hemisphere, e.g. at 150 • E, it would not be registered by the existing NMs, because the rapidly diminishing from the apparent arrival direction particle flux (see the contours of equal pitch angle which also depict the particle flux intensity in the upper panel of Fig.2 and the pitch angle distribution (PAD) of GLE #5 in Fig.3). According to the current definition, a GLE event is registered when there are simultaneous statistically significant enhancements of the count rates of at least two differently located NMs including at least one station near to sea level and a corresponding enhancement in the proton flux measured by space-borne instrument(s) (for details see Poluianov et al., 2017). Therefore, the global NM network could not see a possible event, similar to GLE #5, which poses major space weather thread (the strongest recorded GLE) occurring in the northern hemisphere (see the upper panel of Fig.2).
The existing gap can be filled, by an extension of the NM network with a NM at Severnaya Zemlya (SEVZ) (for details see the lower panel of Fig. 2 and Table 1) and by reopening of the presently non-operational, but previously existed NMs: Alert (ALRT) and Heis Island (HEIS). In addition, as a counterpart of DOMC, we propose a possible location of new NM on the Summit polar station in the Greenland plateau (Table 1), whose asymptotic direction is also given in Fig.2. Such an extended network of polar stations would provide almost global coverage in the maximal NM response rigidity range of 1-5 GV and nearly to symmetric response of NMs from both hemispheres. Here, the computations were performed with the PLANETOCOSMICS code employing the IGRF magnetospheric model corresponding to the epoch 2015 (Desorgher et al., 2005;Thébault et al., 2015).
The extension of the global NM network involved several steps: -Determination of the gaps in the current network and possible locations for new stations (we selected only places with an existing facility providing power supply and data transfer); -Computation of the asymptotic directions of the new stations; -Comparison of performance between current and extended network; -Estimation of the necessary funds and drafting the corresponding proposal;

Services and applications provided by the extended global NM network
Here, We present several abilities of the global NM network, related to space weather services and solar physics research.

Registration and analysis of GLEs
Registration of a GLE can provide an early alert for the onset of SEP event, which is specifically important for various space weather services (for details see Kuwabara et al., 2006a,b). Accordingly, alert systems, based on NM records have been developed Mavromichalaki et al., 2018;Dorman et al., 2019). Most of those alert systems are based on a good coverage of the arrival direction of GLE particles by the global NM network since a given number of stations shall exhibit a count rate increase. Therefore, an extended global NM network will provide a reliable basis for the corresponding alert service(s). Besides, the spectral and angular characteristics of strong SEP events, viz. GLEs in the energy range ∼ 0.3-20 GeV nucleon −1 , can be derived by modeling of the global NM network response. Methods for analysis of GLEs using NM data have been developed over the years, usually based on modeling of the global NM network response and optimization of a set of unknown model parameters n over the experimental data points corresponding to the number of NM stations (e.g. Shea and Smart, 1982;Cramp et al., 1997;Bombardieri et al., 2006;Vashenyuk et al., 2006b). In general, the relative count rate increase of a given NM during GLE can be expressed as: where N is the count rate due to GCR averaged over two hours before the event's onset (e.g. Usoskin et al., 2015), which can be also variable in case of a long event occurred during a Forbush decrease, ∆N(P cut ) is the count rate increase due to solar particles. J sep is the rigidity spectrum of i (proton or α-particle) component of SEPs, usually only protons are taken into account, accordingly J GCR i (P,t) is the rigidity spectrum of the i component (proton or α-particle, etc...) of GCR at given time t, G(α(P,t)) is the pitch angle distribution of SEPs, otherwise, for GCRs the angular distribution is assumed to be isotropic, accordingly, A(P) is a discrete function with A(P)=1 for allowed trajectories (proton with rigidity P can reach the station) and A(P)=0 for forbidden trajectories (proton with rigidity P cannot reach the station). Function A is derived during the asymptotic cone computations. P cut is the minimum rigidity cut-off of the station, accordingly, P max is the maximum rigidity of SEPs considered in the model, whilst for GCR P max = ∞. S k is the NM yield function for vertical and for oblique incidence SEPs (Clem, 1997). The contribution of oblique SEPs to NM response is particularly important for modeling strong and/or very anisotropic events, while for weak and/or moderately strong events it is possible to consider only vertical ones and using S k for an isotropic case, which considerably simplifies the computations (Mishev and Usoskin, 2016a). The background due to GCRs can be computed using a convenient model, e.g., the force-field model with the corresponding local interstellar spectrum, considering explicitly the modulation potential (Usoskin et al., 2005;Vos and Potgieter, 2015). The optimization can be performed over the set of model parameters n by minimizing the difference between the modeled and measured NM responses using a convenient method (Tikhonov et al., 1995;Mavrodiev et al., 2004;Aster et al., 2005;Mishev et al., 2005). The modeling of the global network NM response can be performed using the corresponding NM yield function, which establishes a connection between the primary CR flux at the top of the Earth's atmosphere and the count rate of the device. Since the secondary CRs, resulting from the primary CR induced cascade in the Earth's atmosphere, can reach the ground level and eventually be registered by a NM, the yield function incorporates the full complexity of the atmospheric cascade development including secondary particle propagation in the atmosphere and the efficiency of the detector itself to register the secondaries (e.g. Clem and Dorman, 2000, and references therein). The NM yield function can be determined by parameterization of experimental data, namely latitude survey(s) (e.g. Nagashima et al., 1989;Raubenheimer et al., 1981;Dorman et al., 2000) or can be assessed using Monte Carlo simulations of CR propagation in the atmosphere (e.g. Debrunner and Brunberg, 1968;Clem and Dorman, 2000). Recently, essential progress of Monte Carlo simulations of CR propagation in the atmosphere was achieved, which resulted in several newly computed NM functions Flükiger et al., 2008;Mishev et al., 2013;Mangeard et al., 2016). A recently computed NM yield function by Mishev et al. (2013Mishev et al. ( , 2020 is fully consistent with the experimental latitude surveys and was validated by achieving good agreement between model results and measurements, including space-borne data (Gil et al., 2015;Nuntiyakul et al., 2018;Koldobskiy et al., 2019b).
As an example, we present the derived spectra and PAD of GLE #5, which was the largest event ever observed by the global NM network. It occurred on 23 February 1956 and was registered by various ground-based detectors (ionization chambers, NMs and muon telescopes) and recently was reassessed (e.g. Vashenyuk et al., 2008;Usoskin et al., 2020). This event was very anisotropic. Significant asymmetry between the count rate increases recorded by several European NMs, namely Leeds (LEED), Stockholm (STHM) and Weissenau (WEIS) and American ones, namely Chicago (CHGO), Calgary (CALG) and Ottawa (OTWA) was observed. The stations in Europe revealed rapid and very large NM count rate increases, while those in North America were with considerably delayed maximum and smaller count rate enhancements, for details see gle.oulu.fi. The derived SEPs spectra and PAD are shown in Fig.3. The relativistic solar proton spectra were very hard, specifically during the event's onset initial phase, whilst a narrow PAD was revealed. The SEP spectra remained hard (with nearly exponential shape) during the whole event, in contrast to other GLEs (e.g. Miroshnichenko, 2018, and references therein).
The extended NM network allows to significantly improve the optimization procedure, namely it results in reduction of the residual D, which is defined as: where m is the number of NM stations, ∆N i N i is the relative NM count rate increase for the i NM station.
A robust optimization process and reliable solution are achieved when D ≤ 5 %, a criterion usually fulfilled for strong events, whilst for moderately strong and weak events D can be about 8-12 %. We emphasize that a solution can be obtained even in the case of D ∼ 20-30 %, but with considerably larger uncertainties. Usually, it is necessary to possess about 2(n-1) data (NM stations), n is the number of unknowns in the model, in order to be able to unfold the model parameters (e.g. Himmelblau, 1972;Dennis and Schnabel, 1996;Mavrodiev et al., 2004). Thus, it is sufficient to retrieve information from 15-20 NMs, specifically those in a polar region, whilst the mid-latitude stations provide the boundary conditions. However, this number of stations is reasonable in case of not complicated PAD and unidirectional SEP flux, such as GLE # 59 or GLE # 70 (for details see Mishev and Usoskin, 2016a;. In case of more complicated PADs and/or bi-directional SEP flux, e.g., GLE #69 or GLE #71 (for details see Mishev et al., , 2018b, the amount of required information considerably increases, leading to about 30-35 NM records necessary to perform a reliable analysis. Here, we examined the performance of the extended, actual and reduced NM network for an analysis of several GLEs, the details are given in Table 2. One can see that the extended NM network results in a notably smaller D compared to the actual number of NMs used for the analysis, whilst a reduction of the number of NMs leads to a considerable reduction of the ability of the global NM network to provide a reliable GLE analysis. The additional data used for the analysis with the extended NM network are based on forward modeling including realistic noise similarly to Mavrodiev et al. (2004) employing the derived spectra and PAD during the actual analysis. We note, that the extended analysis is performed with all polar stations from Table 1, which encompasses the extended network, who are added to the actual analysis (a partial overlapping exists for some events, since several NMs from Table 1 are used also for the actual analysis). For the analysis with the reduced NM network we removed about 5-10 NMs with moderate response.

Space weather purposes -exposure during GLEs
The increased intensity of CRs during SEP events, leads to an important space weather issue, namely exposure at flight altitudes (e.g. Mewaldt, 2006;Pulkkinen, 2007;Shea and Smart, 2012, and references therein). During intercontinental flights over the sub-polar and polar regions, aircrews are exposed to non-negligible radiation field due to secondary particles, which can be significantly enhanced during major GLEs (Spurny et al., 1996(Spurny et al., , 2002Shea and Smart, 2000a). Assessments of the exposure during GLEs requires detailed information of SEP spectra as an input for a relevant model for computation of the exposure (e.g. Ferrari et al., 2001;Latocha et al., 2009;Copeland, 2017). The value of the merit function D obtained for the analysis of several GLEs (main phase of the event) as a function of the number of the used NM stations. Columns 1-2 correspond to the number and date of the GLE, while columns 3-5 correspond to D and number of the used stations (in the brackets) for extended NM network, actual NM network used for the analysis and the reduced NM network, respectively. N.A. depicts the case when the SEP spectra cannot be unfolded. The details for the analysis of the presented GLEs are given in Mishev and Usoskin, 2016a;Kocharov et al., 2017;Mishev et al., 2018b) as well as in this work. Here we present as an example the exposure to radiation at flight altitude during the strongest ever observed GLE. The computation was performed using a numerical model Mishev et al., 2018a). The effective dose rate at a given atmospheric depth h induced by a primary CR particle is computed by convolution of the exposure yield function with the corresponding primary CR particle spectrum: where J i (T ) is the differential energy spectrum of the primary CR arriving at the top of the atmosphere for i−th component (proton or α−particle) and Y i is the effective dose yield function for this type of particles. The integration is over the kinetic energy T above T cut (P c ), which is defined by the local cut-off rigidity P c for a nucleus of type i, T cut,i = Z i A i 2 P 2 c + E 2 0 − E 0 , where E 0 = 0.938 GeV/c 2 is the proton's rest mass.
Accordingly, the effective dose yield function Y i is: where C j (T * ) is the coefficient converting the fluence of secondary particles of type j (neutron, proton, γ, e − , e + , µ − , µ + , π − , π + ) with energy T * to the effective dose, F i, j (h, T, T * , θ , ϕ) is the fluence of secondary particles of type j, produced by a primary particle of type i (proton or α−particle) with given primary energy T arriving at the top of the atmosphere from zenith angle θ and azimuth angle ϕ. The conversion coefficients C j (T * ) are considered according to Petoussi-Henss et al. (2010). We note, that employment of different conversion coefficients C j (T * ) (e.g. ICRP, 1996), would lead to increase of the exposure assessment of about 20 %, which is considerably below the other model uncertainties (e.g. Copeland and Atwell, 2019;Yang and Sheu, 2020).
Using the derived rigidity spectra for GLE #5 (Fig.3) and Eq. (3), we computed the effective dose rate at a typical altitude for an intercontinental commercial jet flight of 40 kft (12 190m) a.s.l., altitude representative for a polar flight over a polar atmosphere (for details see Velinov, 2010, 2014;Mironova et al., 2015, and references therein), similarly to Mishev and Usoskin (2018); Copeland and Atwell (2019). Here, we would like to stress that the exposure during GLEs can usually reach peak values considerably greater than the GCR background, but for a relatively short period. Therefore, it is more relevant to integrate the exposure over a certain period, naturally related to the flight duration. However, during GLE #5, the derived SEP spectra remained hard even after the event initial and main phase of the event, i.e., for a relatively long period, which is comparable with a polar flight duration. The distribution of the effective dose over the globe at an altitude of 40 kft a.s.l., integrated over the first four hours after the event onset during GLE #5 is presented in Fig.4. One can see that the exposure is significant in a polar region, where the received dose is considerably greater than the suggested annual limit for occupational workers of about 6 mSv (e.g. EURATOM, 2014). The received dose for the population integrated over four hours in the polar region, which is a typical time span of flight in this region, is about an order of magnitude greater than the recommended of 1 mSv (e.g. EURATOM, 2014). The accumulated exposure is significant even at mid-and high-rigidity cut-off regions, because of the very hard SEP spectra.

Registration of solar neutrons
The global NM network provides a good opportunity to study solar neutrons (e.g., Usoskin et al., 1997;Dorman, 2010;Artamonov et al., 2016). During solar eruptions, accelerated high-energy ions can interact with matter in the solar atmosphere, resulting in in-situ production of different types of secondary particles, e.g. γ-rays and neutrons (for details see Hurford et al., 2003;Dorman, 2010, and references therein). Of specific interest are neutrons, the so-called solar neutrons (e.g., Lingenfelter et al., 1965, and references therein). Since the solar neutrons are neutral, they propagate straight to the Earth, therefore bringing direct information of the acceleration site. If the energy of solar neutrons is greater than about 100 MeV, they can induce a nucleonic cascade in the Earth's atmosphere and can be registered by NMs. The sensitivity of a NM to solar neutrons is greater when the atmospheric depth in the solar direction is smaller, because the atmosphere attenuates the flux of secondary nucleons in the cascade. An optimal location is high-altitude, close to the Equator (Usoskin et al., 1997). In order to improve this capability, it is recommended to extend the current network with at least two high-altitude NMs, namely one located at the Canary Islands and the other in New Zealand, and to re-open the Haleakala (HLEA) NM, details given are in Table 1, (see also Artamonov et al., 2016). We note that the Canary Island NM is under construction (Private communication).

Conclusions
We discussed the current status and application of the global neutron monitor network to study solar energetic particles, specifically for space weather purposes, namely alerts, assessment of SEP characteristics and the corresponding computation of the exposure to radiation at flight altitudes.
As an example, we presented the ability of the global NM network data to be used for derivation of the spectra and angular distribution of SEPs during the strongest GLE event of the observational era: GLE #5 and the related in the course of the event effective dose over the globe. In order to improve those capabilities, we propose to reopen four previously operational NMs, namely ALRT, HEIS, HLEA and VSTK (see Table 1, stations below the dashed line) and to build four new stations: CANI, NZLD, SEVZ, SUMT (see Table 1, stations below the dashed-dashed line). Hence, covering several existing gaps and improving its sensitivity specifically in the low energy range, the global NM network will be a useful tool to study various populations of solar particles and will be a useful instrument for space weather services.
Besides, in order to keep operational those capabilities of the global NM network, we would like to stress that even a partial reduction of the number of existing NMs would considerably influence the usage of the global NM network as a convenient tool for space weather services. Since at present the existence and continuous functioning of several NM stations are under question, the support of the network from governments, foundations(s) and space flight operators is crucially needed.