Retrospective analysis of GLEs and estimates of radiation risks

28 February 2017 marked 75 years since the first confident registration of solar cosmic rays (SCRs), i.e., accelerated solar particles with energies from about 10 to ~10 10 eV. Modern state of the problems related to the studies of Ground Level Enhancements (GLEs) of relativistic SCRs is critically analyzed based on available direct and proxy data. We are also taking into account extremely large fluxes of non-relativistic solar energetic particles (SEPs). Both kinds of SCR events are of great astrophysical and geo-scientific (geophysical) interests. A number of the GLE properties (total statistics, occurrence rate, longitude distribution, ranking of GLEs, a number of specific GLEs – so-called ‘‘rogue’’ SEP events etc.) are discussed in some detail. We note also the problems of GLE identification (definition) by groundbased observations, the difficulties in the studies of weak (‘‘hidden’’, or sub-) GLEs etc. One of serious challenges to the problem of radiation hazard in space is a lack of a clear, unambiguous relation between the fluxes (fluences) of relativistic SCR and non-relativistic SEPs. Special attention is paid to the recent debate on the validity, origin and properties of the ‘‘ancient’’ events AD775, AD994, AD1859 (Carrington event) and BC3372. We demonstrate that, in spite of existing uncertainties in proton fluences above 30 MeV, all of them are fitted well by a unique distribution function, at least, with the present level of solar activity. Extremely large SEP events are shown to obey a probabilistic distribution on their fluences with a sharp break in the range of large fluences (or low probabilities). The studies of this kind may be extended for periods with different levels of solar activity in the past and/or in the future. Dose rates at aircraft altitudes are also demonstrated during some GLEs. Several examples of using the SCR data and GLE properties in radiation prediction schemes are considered.


Introduction
Solar energetic particles (SEPs) represent an important aspect of solar-terrestrial physics, as well as are among the three main components of space weather (geomagnetic field disturbances; enhanced energetic particle dosages; ionospheric electron density disturbances, e.g., Song, 2001). As well established long ago, during energetic solar phenomena, high-energy particles are generated in extended energy range -from about 1 MeV to 10 Ä 100 GeV. Historically, they were called Solar Cosmic Rays (SCRs) because of their close similarity to the true cosmic rays (CR) of galactic origin -Galactic Cosmic Rays (GCRs). The initial observations of GCRs relied upon measurements of secondary particles (muons) generated in the Earth's atmosphere. Events with relativistic SCRs are those where solar protons contain enough energy to initiate nuclear cascades in the atmosphere that can penetrate to the surface of the Earth.
According to modern ideas, SEPs are produced by the rapid release of magnetic energy during solar eruptions. Among them solar protons is one of significant populations of interplanetary particle fluxes that determine important properties of space weather. At the Earth's orbit these protons are registered as a solar proton event (SPE). The most homogeneous solar proton data base has been obtained from the CR records. From 1934, ionization chambers (IC), and, since 1951, neutron monitors (NM), recorded occasional increases in the CR intensity that were associated with the release of GeV protons from the Sun (e.g., Shea & Smart, 1992).
Continuous measurements of sea level ionizing radiation using ionization chambers began in the 1920's, but the validity of the observed intensity variations was doubtful because of Measurement, Specification and Forecasting of the Solar Energetic Particle Environment and GLEs atmospheric effects and instrument instability (see, e.g., Simpson, 1990). Later on, Compton et al. (1934) developed an ionization chamber (IC) of general purpose wherein the average cosmic ray background ionization was nulled out, so current variations above and below the ambient null were represented as time-intensity variations. Just this improved installation has played a crucial role in a discovery of SCRs at the beginning of 1940's.
Although there was evidence that observers in the 1920-1930's had recorded intensity increases, which were due to solar flares (e.g., Elliot, 1952;Chupp, 1996), the intensity increases of 28 February and 7 March 1942 (Lange & Forbush, 1942) and, especially, of 25 July 1946 (Forbush, 1946) registered in the USA and associated with solar flares, first drew attention to the importance of high-energy particles from the Sun. Similar observations and research works were going on in Europe, and the European researchers have reached similar conclusions (e.g., Elliot, 1952). In fact, due to these observations, two important astrophysical phenomena, namely, Ground Level Events of two kinds were discovered for the first time. According to Simpson (1990), the first kind is rapid intensity increase of SCRs (Ground Level Enhancement, or GLE), and the second one is rapid GCR intensity decrease (or Forbushdecrease, FD), both of them being caused by energetic solar phenomena.
Below we present an analytical (critical) review of the main features of GLEs that may be important for the Space Weather and Space Climate (SWSC) goals: 1) specific features of GLE registration; 2) general relation of GLEs to solar activity and reversals of global magnetic field (GMF) of the Sun; 3) occurrence (registration) rate of GLEs; 4) longitude distribution of the GLE sources; 5) specific ''rogue'' SEP events and GLEs; 6) small (''hidden'') events, or sub-GLEs; 7) new definition of GLE. The reliability of GLE analysis based on NM data is considered. The main characteristics of GLEs are discussed in the context of the radiation hazard in space and at aircraft altitudes. In this context, we give some examples of calculation of radiation doses at aircraft altitudes during the real GLEs. In particular, data on the event of 15 April 2001 are used along two actual flights, as well as ambient dose equivalent during GLE of 20 January 2005 computed by different groups for three assumed flights. The discrepancies in the dose estimates by different methods turned out to be depressingly large.
A number of methodological problems of the GLE studies are also critically considered, in particular, so-called ''GLE05 scenario'' in application to ''ancient'' large SEP events and GLEs, as well as possible using of SCR data in the prediction schemes. Taking into account recent proxy data on some large ''ancient'' SEP events, we suggest a new distribution function for the proton fluences above 30 MeV. We briefly discuss also a probability of super-flares at present Sun.

Historical outline
The original ionization chambers (IC) and counter telescopes are now classified as muon detectors, in particular, standard muon telescope (MT). These detectors respond to primary high-energy (relativistic) protons (>4 GeV) interacting at the top of the atmosphere. It is worth to mention that an IC of the Compton's construction (Compton et al., 1934) was prone to uncontrolled drifts in the current thus making it very difficult to provide a long-term variability of GCRs (e.g., McCracken & Beer, 2007). Later on, this problem was taken into account, and in the late 1940s, in a much more sophisticated Yakutsk design (Shafer & Shafer, 1984), this defect was significantly reduced. By the Yakutsk IC, in particular, one of the first large GLEs (namely, GLE04 on 19 November 1949) has been registered (e.g., Krasilnikov et al., 1955).
In the 1950's, development of the CR neutron monitor for secondary neutrons and protons (nucleonic component) lowered the detection threshold to >450 MeV primary protons (e.g., Simpson, 1957). A number of standard neutron monitors (NM of IGY type) was deployed for the International Geophysical Year (IGY, 1957(IGY, -1958, and due to those instruments, in particular, the outstanding SCR event of 23 February 1956 has been recorded. Soon after the NM design was improved by Carmichael (1968) with the development of the so-called ''super'' neutron monitor (now the conventional name for a ''super-NM'' is NM64).
The present-day worldwide network for continuous CR registration includes~50 stations equipped mainly with NM64 super-monitors ( Fig. 1), the data of which form the international NM database (NMDB, http://www.nmdb.eu). It should be stated that NMDB is not a complete dataset of the worldwide NM network (not all stations provide data to NMDB). The complete dataset of all NM data is collected at World Data Center on Cosmic Rays (WDC-CR) (http://cidas. isee.nagoya-u.ac.jp/WDCCR/). A number of ground MTs of different design make it possible to register SCRs arriving at large angles to the vertical. Several underground MTs are also used to register extreme hard relativistic events, such as the GLEs of 23 February 1956 and 29 September 1989. Groundbased observations are satisfactorily completed with the balloon measurements in the stratosphere (e.g., Stozhkov et al., 2009), as well as with the network of solar neutron telescopes (SNTs) (e.g., Flückiger et al., 1998), which register the arrival of secondary neutrons generated by primary accelerated ions in collisions with nuclei in the lower layers of the solar atmosphere.
In total, seventy GLEs have been registered from February 1942 to December 2006 (solar cycles 17-23). The list of 70 GLEs may be found, e.g., in (Miroshnichenko & Pérez-Peraza, 2008;Miroshnichenko, 2014). Similar list (up to GLE of 10 September 2017) is given at the website of IZMIRAN (http://cosrays.izmiran.ru). The official International GLE database, stored in Oulu (http://gle.oulu.fi), contains information on the 72 GLEs observed so far (see details below). For the convenience of researchers, at the turn of the 1990's all GLEs received official numbers, starting from 28 February 1942 (GLE01). The last event in cycle 23 of solar activity (SA) was observed on 13 December 2006 (GLE70). In current cycle 24 (started in January 2009), proton activity of the Sun was manifested with a delay: the first GLE in this cycle occurred only on 17 May 2012 (GLE71). The second, rather weak, but officially recognized GLE event of cycle 24 occurred on 10 September 2017 (GLE72).
At present, many researchers believe that most of SEPs are generated by shock waves from Coronal Mass Ejections (CMEs) (e.g., Reames, 1999Reames, , 2017. As to relativistic solar protons (RSP), current thinking is that GLEs are produced mainly due to fast particle acceleration caused by magnetic reconnection processes in solar flares (e.g., Somov, 2013a,b). At the same time, we do not exclude that some fraction of SEPs may be accelerated up to relativistic energies by CME-driven coronal/interplanetary shocks (e.g., Li et al., 2013Li et al., , 2016. In any case, from the point of view of the SWSC, we have serious grounds to consider the GLE particles as an extension of general SEP spectrum into relativistic energy range.

GLE main properties and databases
The events with the presence of relativistic solar protons (i.e., GLEs) are very diverse in their characteristics, starting from their magnitudes and ending with the peculiarities of the energy spectrum (all aspects of the definition of a GLE are discussed below in separate Sect. 10). The first GLEs (before 1956) were registered at sparse stations equipped only with ICs and/or MTs, which were mainly intended for measuring one hard (muon) component. Since the effective registration energy of NMs is lower than that of MTs and ICs, the latter detectors are less sensitive to SCRs.

Features of GLE registration
It is interesting to note that the last event in cycle 23 (GLE70 on 13 December 2006) was registered not only at the worldwide NM network but also with non-standard ground detectors, specifically, with the URAGAN muon hodoscope (Timashkov et al., 2007). Moreover, this GLE was also registered with the IceTop extensive air shower (EAS) detector, which is a component of the IceCube neutrino telescope in Antarctica (Abbasi et al., 2008). All GLEs registered from 1942 to 2012 are listed, e.g., in , see also the website of IZMIRAN (http://cosrays.izmiran.ru) and the International GLE database (http://gle.oulu.fi).
A special technique is used to derive the GLE characteristics at the Earth's orbit (e.g., Shea & Smart, 1982;Vashenyuk et al., 2011;Miroshnichenko et al., 2013; see also Mishev et al., 2017). This technique takes into account the anisotropy of SCRs fluxes that approach the Earth, steep energy spectrum of SCRs, and the high NM sensitivity. At the same time, some weak GLEs (~1 Ä 10% by 5-minute data of NMs) were registered only at high-latitude or polar stations that are most suitable for registration of weak events. For example, on 29 October 2015 two neutron monitors on the Antarctic plateau (Dome C) at high altitudes (3233 m elevation) have observed an increase in the count rate caused by SCRs  whereas no other polar neutron monitor station at or near sea level have registered a change in the count rate.
If the energy of primary protons is E < 100 MeV (R < 0.44 GV), neutron monitors practically do not respond to these protons due to the atmospheric absorption of the secondary neutrons (this is so-called ''atmospheric cutoff'' at the rigidity R a ), the maximum of the NM response to GCR being within 1-5 GV. It means that all high-latitude (polar) NM stations start to record secondary neutrons efficiently from the same rigidity of the primary protons about 1 GV (E p % 433 MeV), irrespective of the NM nominal (calculated) ''geomagnetic cutoff rigidity'', R c . As it fortunately happened, a rigidity % 1.0 GV is approximately the midway between the low-rigidity satellite interval (e.g., GOES) and relativistic rigidity range (Smart & Shea 1996), and it turned out to be a convenient reference point as a characteristic rigidity cutoff at the polar NM stations. It is pertinent to note that the NM response to SCR is determined by the ''effective rigidity'' of their spectra (e.g., Miroshnichenko, 2001, Chapter 9), and the maximum response of a polar NM to SCRs is lower than that to GCRs (see below).
According to the recent studies (e.g., Evenson, 2011;Evenson et al., 2011), below a geomagnetic cutoff of about 0.6 GV (E p % 200 MeV), atmospheric absorption determines the cutoff. At the same time, the high altitude NM stations at the South Pole (Antarctic Plateau, SOPO/SOPB at Amundsen-Scott station, 2835 m elevation, R c = 0.11 GV) have a number of important advantages: a uniform energy response, an excellent directional sensitivity, and a focusing of obliquely incident primaries. A neutron monitor at the South Pole is in operation since 1964 http://www.nmdb.eu). After the start of operation of the two mini-NMs at Dome C/B (Concordia Research Station, Central Antarctica) in 2015, the worldwide NM network became more sensitive to the registration of SEP events, because high elevation polar NMs possess lower compared to the sea level atmospheric cutoff, R a . As a result, the high altitude polar NM stations are able to detect SEP events, which would not be registered by the other (nearly sea level) NMs (threshold energy of about 300 MeV instead of 433-450 MeV). Some features of locations of geophysical Although the South Pole and Dome C stations are located on the Antarctic Plateau, they have quite different asymptotic acceptance cones (Poluianov et al., 2015). Dome C has the most poleward acceptance cone among all the stations: particles with rigidities between 1 and 20 GV are coming from latitudes between À75°and À90°. Other Antarctic stations, including South Pole, receive particles from a wider range of latitudes between À20°and À85°. Since less energetic cosmic rays are more numerous, most of the particles detected at those stations are coming from the equatorial plane, while Dome C has its cones almost totally polar-oriented even for the lowenergy part of spectrum.
As to the two mentioned mini-NMs, they were installed at the point with rather specific conditions (Dome C, Central Antarctica, 75°06 0 S, 123°23 0 E, 3233 m a.s.l.). The site has unique properties ideal for CR measurements, especially for detection of SEPs: very low ''nominal'' cutoff rigidity <0.01 GV, high elevation and, as already noted above, poleward asymptotic acceptance cones pointing to geographic latitudes >75°S. The instruments consist of a standard neutron monitor and a ''bare'' (lead-free) neutron monitor. As known, a neutron monitor without lead has enhanced sensitivity to low-energy nuclei and, consequently, to weaker cascades and softer primary particles (e.g., Oh et al., 2012), but at the cost of reduced total instrument efficiency . The instrument operation started in mid-January 2015. Several interesting events, including SPE of 29 October 2015 have been registered (Poluianov et al., 2015;Mishev et al., 2017). The data sets are available at the websites https://gle.oulu.fi and http://www.nmdb.eu.
The website cosmicrays.oulu.fi/ is the official international database of NM count rates during GLEs caused by solar energetic particles. After being created and maintained by L.C. Gentile, M.A. (Peggy) Shea, D.F. Smart, and M.C. Duldig, this database has been moved in 2014 to Finland (Sodankylä Geophysical Observatory, University of Oulu, http:// cosmicrays.oulu.fi/, data manager I.G. Usoskin). The European Community's Seventh Framework Programme (the FP7 project) which has developed the neutron monitor database (NMDB) started in 2008. The NM database becomes operative in August 2008 (website http://www.nmdb.eu, data manager C.T. Steigies, Kiel, Germany). All the GLEs, starting from number GLE05 are archived at the International GLE database http://gle.oulu.fi .

Statistics of GLEs
Many researchers have for decades attempted to streamline the statistics of the GLEs. One of the first lists of such events was drawn up in the work by Shea & Smart (1993). Later, the list of events was gradually replenished with various new data, along with their statistical and physical analyses (e.g., Miroshnichenko & Pérez-Peraza, 2008;Andriopoulou et al., 2011;Belov et al., 2010;Miroshnichenko et al., 2013;Asvestari et al., 2017;; see also the web-sites http://cosmicrays. oulu.fi/GLE.html; http://www.nmdb.eu; http://cosrays.izmiran.ru). For obvious reasons, early events in the history of the SCR study are of particular interest. Table 1 shows the numbers of registered GLEs per cycle, together with available data on the peculiarities of solar activity (SA) and reversals of global magnetic field (GMF) of the Sun for the entire period of GLE observations in solar cycles 17-24. Data on the SA level in the past were taken from Venkatesan (1958), Duggal (1979), McCracken (2007, and for cycle 24 -from the site https://solarscience.msfc.nasa.gov/ SunspotCycle.shtml (D. Hathaway's predictions). The last column of Table 1 contains necessary information on the detectors used to record SCRs, as well as on some features of solar activity in different cycles. For example, starting from the cycle 18, it was established a presence of the second maximum in the sunspot numbers. This effect was called ''Gnevyshev Gap'' after the astronomer who initiated investigation of the double-peak L.I. Miroshnichenko: J. Space Weather Space Clim. 2018, 8, A52 structure of the solar activity cycle (see Gnevyshev, 1977 and references therein).
Based on the available data on 72 GLEs (Table 1), we can assume that in the early years of observations some weak GLEs were not registered on technical and methodological reasons. If the average occurrence rate of the GLEs is g~1.0 yr À1 , the number of omitted events in the cycles 17 and 18 could be considerable (Miroshnichenko et al., 2012), up to 8 events per cycle (see Table 1). On the other hand, in the cycles 20-23 a number of GLEs per cycle (from 12 to 16 events) is practically independent of the cycle magnitude (by sunspot numbers), starting from the cycle 20, when a technique for the GLE registration was considerably improved, namely, the worldwide network has been equipped by the super-monitors on the NM64 neutron counters. From the beginning of the cycle 20 (October 1964) up to December 2008 (44 years) in total 56 GLEs have been registered, i.e., their occurrence rate was about 1.27 per year. This empirical fact seems to be one of the most important for the understanding of SCR production.
The current cycle 24, however, severely violated this outlined pattern. Indeed, a prolonged minimum of cycle 23 ended in December 2008, and cycle 24 (started in January 2009) proceeds very slowly (flabbily): sunspot formation, solar flare and proton activities are generally at a rather low level. Thus, only two rather weak GLEs (GLE71 on 17 May 2012 and GLE72 on 10 September 2017) were registered during almost ten years of the cycle (up to October 2018). Of course, it cannot be excluded that in the descending phase of the cycle, the Sun will produce 1-2 strong GLEs as it happened during the descending phase of the cycle 23. But such a possibility seems highly unlikely. On the other hand, due to new installation at Dome C (Antarctic) several sub-GLEs have been registered since 2015 (for details see Sect. 8).

Diversity of GLE properties
Due to diversity of the GLE properties, their hierarchy (e.g., by magnitude) strongly depends on the key classification parameters -the intervals of data averaging, range of energies (magnetic rigidities) under consideration, spectrum shape, SCR intensity and/or fluence, in particular, on so-called ''GLE strength'' , etc. An important effect on the magnitude of a GLE has also the source location due to the effect of interplanetary magnetic field (IMF) connectivity. For example, a considerable flux of relativistic particles does not mean that a powerful flux of protons will also be observed in the non-relativistic range. Similarly, an abundant flux of non-relativistic particles cannot unambiguously indicate that the flux of relativistic SCRs is large (e.g., Mewaldt et al., 2007;Miroshnichenko, 2014). As a result of this spectral peculiarity (a sharp spectrum steepening with increasing of particle energy), it is very difficult to model SCR acceleration processes (e.g., Priest & Forbes, 2000;Somov, 2013a,b;Miroshnichenko & Pérez-Peraza, 2008;Miroshnichenko, 2014 and construct models for predicting radiation hazard (e.g., Miroshnichenko, 2003a).
As a characteristic example, in Figure 3 we present the spectra for several observed SPEs, including a number of well-known GLEs (Mewaldt et al., 2005a(Mewaldt et al., ,b, 2007Wang, 2009). The left panel shows so-called ''time-of-maximum'' differential spectra for the three events measured in different ranges of energy at/near the Earth's orbit. For comparison, we show the data of three spacecraft for the event of 5 December 2006. These last measurements have been carried out at significant distance from the Earth (in particular, onboard spacecraft ACE in the Lagrange point L1). The middle panel shows integral intensity spectra for several GLEs in the energy range between 10 MeV-10 GeV (Wang, 2009) along with expected upper limit spectrum (ULS) by Miroshnichenko (1996) (dashed curve). The right panel demonstrates integral energy spectra of so-called fluences (Mewaldt et al., 2005a(Mewaldt et al., ,b, 2007 for several outstanding SPEs-GLEs (the fluence is the proton flux integrated during the whole event under consideration). This characteristic of SPE is important, first of all, for the studies of geophysical effects of SCRs (e.g., for the calculations of the production rate of cosmogenic isotopes in the terrestrial atmosphere) as well as for the estimates of radiation hazard near the Earth's orbit.
One can clearly see significant differences between the events -both in the form (slope) of their spectra (e.g., left panel in Fig. 3) and in the magnitude of the fluence in each event (right panel). The spectra have a clearly pronounced variable slope, and with increasing energy the spectrum becomes steeper, i.e., the exponent c depends on the energy. An explanation for this behavior has not yet been found. It can only be stated that starting from the break-off energy the spectrum can be described by a power-law function with exponential cutoff (Ellison & Ramaty, 1985). According to Mewaldt et al. (2007), Tylka & Dietrich (2009 and many other researchers, such spectra are better described by a specific double power-law function proposed by Band et al. (1993).
In the paper by Atwell et al. (2010) the authors review and discuss two methods for describing the proton energy spectra of several major GLEs and compare the ''old method'' of exponential extrapolation in proton rigidity with the ''new method'' based on the Band function fit, which combines both GOES satellite data with ground-based NM data to completely describe the entire energy spectrum from 10 MeV to several GeV. In particular, they use the Band function fitted spectra for four major SPEs-GLEs (February 1956, August 1972, October 1989. The Band-function describes the integral rigidity spectrum by a double power law in rigidity with a smooth exponential junction in-between. This approximation describes the integral spectrum by a double power law in rigidity with a smooth roll over in-between. A tremendous work has been performed by Tylka & Dietrich (2009 to make a Band-function fit to the most GLE events (59 from 72), using both NM data for high-energy tail and in-situ space-borne measurements for the lower part of the spectrum.
Concerning event fluence distributions, it should be noted that the right panel in Figure 3 covers a set of data from about 50 years (several solar cycles), thereby tacitly assuming that these data are representative of a much longer timescale which is by no means guaranteed (see, e.g., Miroshnichenko & Nymmik, 2014). Therefore it is appropriate to emphasize that solar cycles differ widely, and this may prove to be important, in particular, for evaluating fluences from ''ancient'' SPEs (Sects. 6 and 7).

GLE ranking
As it was noted, the historical beginning of SCR observations was the GLE of 28 February 1942. By now (October 2018), 72 GLEs have been confidently recorded by the worldwide network of CR stations. Last of them (GLE72) has been registered on 10 September 2017 (https://gle.oulu.fi) with rather small amplitude (a few percent by 5-minute data) at several NM stations. These events characterize only one, relativistic part of the entire energy spectrum of SCR (kinetic energy about E ! 433 MeV/nucleon, or magnetic rigidity R ! 1 GV). It means that by ground based CR detectors only the relativistic range of the energy spectrum of the SCRs can be studied.
The spectra for the events in the solar cycle 23 were summarized in (Wang, 2009). After an outstanding relativistic event on 29 September 1989 (GLE42), the list of Wang was completed with a very large event on 20 January 2005 (GLE69), which was mainly registered with NMs. Standard MTs apparently did not register any increase, although some non-standard muon detectors (e.g., D'Andrea & Poirier, 2005, project GRAND) registered statistically significant effects. Moreover, the GLE69 caused us (see Miroshnichenko et al., 2013) and independently another research group (Bieber et al., 2013) to revise the entire ''hierarchy'' of the events based on the maximal increase at relativistic energies (Smart & Shea, 1991). To achieve the reliable results of GLE ranking, the SCR data must also be evaluated for anisotropies, and identical time intervals should be used for comparisons. How be it, the events may be ranked now as shown in Table 2 (N.O. means ''no observations''). As one can see, special attention in Table 2 is paid to the interval of the data averaging. Table 2 shows the largest GLEs with important information on the interval of data averaging taken mainly from Dorman (1957) and the website of IZMIRAN (http://cosrays.izmiran.ru). In the earlier GLEs the time resolutions of the measurements were, e.g., 60 min, 30 min, and 10 min. Only later with the advent of computers it was possible to store the counts with higher time resolutions, e.g., 1-minute data. Especially, in the event of 20 January 2005 several polar NMs (Terre Adélie, Antarctica, South Pole and others) and special MT station (project GRAND, D'Andrea & Poirier, 2005) were able to register the data with 1-minute resolution.
It should be noted that the main purpose of Table 2 was just to demonstrate (identify) the dependence of the GLE ranking, first of all, from the intervals of data averaging. Therefore, there was no reason to reduce data in Table 2 to the same time resolution, just the opposite -it was necessary to give them at actual averaging. For example, if the SOPO-B NM (with no lead producer) detected a > 4800% increase by 5-minute data between 06:50-06:55 UT on 20 January 2005, then a 15-minute average for the interval between 06:45-07:00 UT will be about 2700% that is apparently less than maximum increase on 23 February 1956 by 15-minute data of the NM Leeds. Time-of-maximum differential spectra for several SPEs observed at the Earth's orbit (14 July 2000, 28 October 2003, 12 May 2006 and at some distance from the Earth (5 December 2006) (adapted from Mewaldt et al., 2007). Middle panel: Integral intensity spectra for several GLEs in the energy range between 10 MeV-10 GeV (from Wang, 2009) along with expected upper limit spectrum (ULS) by Miroshnichenko (1996) (dashed curve). Right panel: Spectra of proton fluences obtained for several large SPEs of 1956-2005 (adapted from Mewaldt et al., 2005aMewaldt et al., ,b, 2007. On the other hand, over the same 6-minute span on 20 January 2005, the NM64 rate at the sea level station at Terre Adélie (1-minute records) increased by about 4200% (i.e., a factor of 43) over the pre-event GCR background, while the rate at the sea level station at McMurdo (Antarctica) increased by a factor of about 2600%. The increase at the high-altitude NM64 (2820 m) station at South Pole was even greater, a factor of about 5500% (1-minute data). However, this distinction is owing to the unique location of South Pole NM64, as mentioned above, at both high latitude and high altitude. Corrected to sea level the increase at South Pole NM64 would have been ''only'' about 2300% (Bieber et al., 2013).
In the last column on the right of Table 2 we present the results of recalculation by Duggal (1979) for five ''early'' GLEs to possible amplitudes of count increase in the neutron component, under assumption that high-latitude (polar) NMs were in operation at that times. Event GLE69 occupies the second line in the list; GLE04 (19 November 1949) is at the third line; two first GLEs are at 5th and 6th lines, respectively. At the same time, the well-known GLE42 with a very hard spectrum occupies only the modest seventh line. However, if we start from the amplitude of increase at the NMs, an outstanding GLE05 still remains the event of rank 1. It is timely to note that, proceeding mainly from Duggal's estimates, Bieber et al. (2013) obtained an analogous distribution of GLEs by their rank. It should be specified that above discussion on the GLE ranking uncertainties is related only to the peak intensities. The ''GLE strength'' proposed by Asvestari et al. (2017) is free of these uncertainties (see also Sect. 6.1).
Here it would be appropriate to emphasize that the powerful GLEs are, as a rule, strong solar proton events (SPEs) that cause significant geophysical effects, in particular, production of large amount of nitrates (''odd nitrogen NO y '') and cosmogenic isotopes (such as 3 H, 7 Be, 10 Be, 14 C, 36 Cl). NO y is a generic term given to the various nitrate compounds generated in the atmosphere as a result of the ionization, dissociations, dissociative ionizations and excitations caused by the solar proton interactions with the air molecules. Cosmogenic isotopes are produced mainly by splitting the target nuclei of atmospheric oxygen 16 O, nitrogen 14 N, argon 40 Ar, and some others, due to proton interactions in the atmosphere.
As known, nitrates and cosmogenic isotopes can be accumulated and stored for a long time in some natural archives (polar ice, tree rings, deposits on the bottom of the seas and oceans, etc.). This makes it possible to use those proxy data for studying ''ancient'' SEP events. Note, however, that now the dominating current thinking is that the nitrate data (''nitrate hypothesis'', ''nitrate method'') cannot be used as a quantitative proxy for SEP (GLE) events. Some aspects of this controversial (disputable) issue are critically examined in more detail in Section 6.

Occurrence (registration) rate
Another interesting aspect, which characterizes the Sun as a star, was revealed as a result of a wavelet analysis of the GLE registration rate g, depending on the SA level (sunspot number SS) and solar cycle epoch (Pérez-Peraza et al., 2009;  1. For technical reasons, when registering the early GLEs, the averaging time intervals (given after the name of CR station) were different for different detectors and may be changed from 1 h (60 min) to 1 min. 2. Maximum amplitudes of the first 4 GLEs have been reduced to a 15-minute averaging time interval (Dorman, 1957). 3. Ionization chamber (IC) events have been normalized to correspond to the increase that a high-latitude neutron monitor (NM) would have observed at sea level (Duggal, 1979). 4. Neutron component records for the GLE04 in Manchester were not obtained with the help of a ''classical'' NM, but with the help of the neutron detector that consisted of boron-trifluoride filled counters embedded in a ''pile'' built of high purity graphite (Adams, 1950). 5. The effective vertical cutoff rigidities of the neutron monitors (in GV units) are indicated for the CR stations with maximal observed amplitude of the corresponding GLEs. Miroshnichenko et al., 2012). As known, SA is notable for nonlinear and chaotic behaviour, as well as for short-term and long-term variations on time-scales in the range of hours to several hundred (and even thousand) years, such as the well-known sunspot cycle and its longer-period modulations. Therefore, it seems reasonable to use the corresponding wavelet analysis methods in order to compare the observed variations in activity indicators with their mathematical models. Such a comparison will make it possible to better describe the behaviour of the toroidal and poloidal magnetic fields of the solar dynamo. Wavelet analysis is a powerful tool for revealing the predominant oscillation mode and for studying the oscillation time evolution by transforming nonlinear time series into a ''time-frequency'' space (e.g., Torrence & Compo, 1998). On the other hand, wavelet analysis is a complicated and modernized version of harmonic analysis, when the latter is performed together with the data compression of different data sets (e.g., solar and geophysical ones) in a wide region of applications. The simplest example of data compression is a removal of spaces (empty spaces) in the data sets. The wavelet analysis technique is often applied, in particular, to time series in order to rapidly discover short-term phenomena. Note that an ordinary time series makes it possible to adequately localize a signal in time but does not include any information on the oscillation frequency. On the other hand, the Fourier transform gives a high frequency resolution, but does not make it possible to localize a signal in time. In this connection, a wavelet is an ''optimal'' combination of the time localization and frequency characteristics of the studied oscillations.
Using the dates of the 70 GLEs and the wavelet Morlet scheme, we constructed the PWM time series (Pulse Width Modulation) for parameter g, which includes the statistically significant oscillation with a period of~11 years (Fig. 4). In this case, the oscillations of g to a certain degree are coherent with the time series of the parameters of the photosphere (sunspot number SS) and corona (coronal index CI). Since the main sources of large GLEs are solar flares, then it should apparently be expected that similar coherence will be present if one uses also so-called Solar Flare Index (SFI), which, along with the coronal index CI, serves as one of indirect (proxy) characteristics of the SA level (see, e.g., Gupta et al., 2007). In spite of the limitations of the GLE statistics and restrictions of the wavelet analysis method, these results can be interesting to understand the periodic phenomena in the solar dynamo, solar atmosphere, interplanetary medium, and galactic cosmic rays, as well as for long-term forecast of GLEs.
The tendency of GLEs to group mainly on the ascending and descending branches of solar cycles (Fig. 4, top panel) is suggested to be caused by the specific features in the spatialtemporal behaviour of the GSMF. As is known, this field reverses its sign precisely near SA maximums. In this context, we mention the results by Nagashima et al. (1991). These authors used MT and NM data for the 43 GLEs from 1942 to 1990 in order to analyze the above GLE tendency. They indicated that the solar proton eruptions (flares) that cause GLEs are basically absent (''forbidden'') during the cycle transient phase, when the GSMF reverses its sign. The absence of GLEs at a SA maximum is explained by a decrease in the particle acceleration efficiency during the GSMF reconfiguration rather L.I. Miroshnichenko: J. Space Weather Space Clim. 2018, 8, A52 than by the suppression of SCR escape processes due to strong magnetic fields.
Since certain periodicities found by Miroshnichenko et al. (2012) are coherent for parameters g, SS, and CI (and, probably, for SFI), therefore, during this stage of the study, the conclusion can be made that oscillations are synchronized in different layers of the solar atmosphere -from the photosphere to the corona. This can indicate that the SCR (GLE) generation involves extended areas in the solar atmosphere rather than being a local (isolated) process.
Based on these results, interesting prospects of longterm GLE prediction have been recently formulated by Miroshnichenko et al. (2013): a) the GCR oscillation distribution can be used to study a solar cycle, to predict SPEs, and in other practical applications in the space weather problem; b) oscillations in the SS (sunspots), CI (coronal index), and GCRs obey the hierarchy principle (e.g., GCR intensity is modulated by solar activity SS etc.); c) meanwhile, oscillations in the GLE occurrence rate are apparently of absolutely different nature, with the only statistically significant oscillation with a period of~11 years. Only this period can be confirmed by two means in the wavelet analysis -1) by the oscillation spectrum of parameter g (PWM series for the registration rate, see above) and 2) by the coherence evolution analysis for the series of g, SS, CI, and SFI (solar flare index). All other oscillations (with shorter periods) are close to or below the red-noise level in the wavelet method. Further discussion of the GLE forecasting issues is given in Section 8.

Longitude effect of GLEs
As it was noted above, total GLE statistics (72 events), accumulated for more than 75 years of ground-based observations, makes it possible to study some problems related to spatial-temporal variations of SA and properties of the GMF of the Sun. For example, it is of great interest to determine the frequency distribution GLEs as function of the heliolongitude of the GLE sources (flares). It was established long ago that interplanetary magnetic field (IMF) is a guiding factor in the process of SCR flux formation near the Earth's orbit (e.g., Dorman & Miroshnichenko, 1968;Toptygin, 1985). According to a current paradigm, the IMF consists of two components -quasi-regular and irregular ones -which affect the SCR particles in different ways. The relative role of each of them is determined by the particle energies. In general, the particle movement along the IMF lines of force is described by the model of guiding center. If the particle scattering on inhomogeneities (irregularities) of the IMF becomes essential, then it is necessary to use a more complex mathematical apparatus (e.g., Toptygin, 1985).
Although relativistic solar particles in their movement to the Earth, as a rule, do not undergo significant scattering (sometimes their transport path in the IMF can be compared with 1.0 AU), the probability of reaching the Earth obviously strongly depends on the angle of curvature of the Parker spiral of the IMF. This leads to a rather strong dependence of GLE registration rate g on the heliolongitude of the assumed source ( Fig. 5): the most of the sources are associated with a heliolongitude interval of~30°W Ä 90°W. It is striking, however, that in 12 cases out of 72 the SCR particles came to the Earth even from the behind-the-limb sources. The sources of two recognized GLEs of the cycle 24 (GLE71 and GLE72) have had the helio-coordinates N11W86 and S09W92, respectively.

GLEs and specific ''rogue'' SEP events
From the point of view of radiation risks, the studies of GLEs are of special interest during so-called ''uncontrolled'', or ''rogue'' SEP events. As applied to some unusual SCR (GLE) events, this name was first proposed by Kallenrode & Cliver (2001a,b). As noted by these authors, about once in a solar cycle, a SEP event occurs whose fluence in non-relativistic range near the Earth's orbit dominates that for the entire cycle (e.g., . Rogue SEP events also have been observed in the inner heliosphere (with Helios 1 on 4 November 1980 at 0.5 AU) and with Ulysses in March 1991 at 2.5 AU. Kallenrode & Cliver (2001a,b) focused on the solar and interplanetary conditions that gave rise to unusual SEP events. It was suggested to identify them by some leading criteria, such as long-lasting high intensities of SEPs between two converging shocks. Since the converging shocks appear to be the leading criterion, they called these particle phenomena ''rogue events'', in analogy to rogue ocean waves during which, due to the superposition of two wave fields, wave heights of more than 30 m can be acquired. Well-known events of this kind are, in particular, the SPEs registered on 14-17 July 1959 (GLE07), 4 August 1972 (GLE24), 19 October 1989 (GLE43), and 14 July 2000 (GLE59) (Kallenrode & Cliver, 2001a,b). Similar events apparently also took place on 12 November 1960 (GLE10)  (Chirkov & Filippov,1977;, 1978 and 12 October 1981 (GLE36) (Kuzmin et al., 1983).
The events in August 1972 and July 1959 are discussed in some details by Smart & Shea (1989) and Miroshnichenko (2014). In particular, Smart & Shea (1989) believe that the SEP fluxes observed on 4 August 1972 by Pioneer 9 spacecraft and surface detectors are the result of a large ejection and reacceleration of solar particles into a region of space where the converging interplanetary structures existed. Moreover, if satellite measurements had been available in July 1959, then the fluence from this activity episode would have been equal to or greater than the fluence observed in August 1972. In fact, the total fluence of the >10 MeV protons, U(!10 MeV) = 1.53 · 10 10 cm À2 , deduced for the July 1959 solar activity episode (e.g., McDonald, 1963;Webber, 1963;Feynman et al., 1988Feynman et al., , 1990, exceeded the measured fluence U(!10 MeV) = 1.1 · 10 10 cm À2 , observed for the August 1972 solar activity episode. At the same time, the relationship between the corresponding fluences for the !30 MeV protons turned out to be inverse: U(!30 MeV) = 3.31 · 10 9 cm À2 in July 1959 and U(!30 MeV) = 5.0 · 10 9 cm À2 in August 1972 (Smart & Shea, 1989;. Note, however, that the accuracy of determining the fluences is low and can be questioned.
In the paper by Levy et al. (1976) the GLE of 4 August 1972 is considered as an example of adiabatic Fermi acceleration of energetic particles between converging interplanetary shock waves. On the other hand, one of the key parameters in the studies of ''rogue events'' can be the shape (slope) of the spectrum of the accelerated particles and its temporal behaviour. Although on 4 August 1972 the increase was recorded on the Earth's surface (e.g., in Apatity its amplitude reached~15% in the 15-min NM data), it cannot be considered as an ordinary GLE. Using the NM data from several CR stations (Tixie Bay, Deep River, Uppsala, London, Khabarovsk),  have constructed the rigidity spectra for a number of relevant GLEs. It was shown that the evolution of the spectrum at R > 1.0 GV during the ''normal'' GLE, e.g., of 7 August 1972 (GLE25), is considerably different from that obtained for the GLEs of 4 August 1972, 17 July 1959, and 12 November 1960. The spectra of the ''rogue events'' are shifted into higher energy (rigidity) range and become harder in time during the events. For example, if the spectrum is presented in power-law form~R Àc , then the value of c have changed from 7.0 to 3.0 in the interval 02:00-08:00 UT on 17 July 1959. For the event of 4 August 1972 the spectrum hardening was observed after intensity maximum about 13:00 UT. In fact, the value of c between 14:00-18:00 UT has changed from 6.0 to 2.1. A hardening of the spectrum on 12 November 1960 took place in a more complicated way due to the presence of two maxima in this event (see Fig. 2 in . Taking into account disturbed interplanetary conditions during the August 1972 episode of solar activity,  suggested that the GLE24 was very likely caused by protons accelerated between two convergent interplanetary shock waves. The same appears to be true also for the event of 17 July 1959. Simple model calculations by  confirmed that a first order Fermi mechanism (see also Levy et al., 1976) operating in shock waves is the most promising one for the interpretation of the sources of ''rogue GLE events''. The above results on the spectrum hardening with time for the GLEs of 17 July 1959, 12 November 1960, and 4 August 1972 are observational evidence of the CR acceleration up to relativistic energies 0.5-3.0 GeV in the interplanetary medium. The best recent reviews in the field of shock acceleration are given by Priest & Forbes (2000, Chapter 13) and Somov (2013b), Part II, Chapter 9).
The unusual character of the events of 17 July 1959 and 12 November 1960 was noted long ago (e.g., Steljes et al., 1961;Carmichael, 1962). After the gigantic GLE05 of 23 February 1956 (see Table 2), rather weak GLE06 of 31 August 1956(McCracken, 1959a, and the blank years 1957 and 1958, the GLE07 of 17 July 1959 was abnormal in several aspects (Carmichael, 1962). One of the intriguing results by Carmichael (1962) was that no appreciable flux of accelerated particles of magnetic rigidity greater than 1.2 GV apparently reached the Earth from the Sun on 17 July 1959. This feature is also confirmed by the data on very soft rigidity spectra of the GLE24 of 4 August 1972 (e.g., Lockwood et al., 1974). The event of 12 November 1960 was noteworthy in that two plasma fronts, rather than one only, were in transit towards the Earth at the time of the flare (Carmichael, 1962). It seems that from a physical point of view the problem of unusual (possibly -''rogue'') GLEs deserves more careful study already at the present level of understanding.

Radiation at aircraft altitudes
As a result of GCR interactions with the atmosphere, there are produced a huge amount of secondary particles which together with the primary incident particles give rise to ionizing radiation exposure through the atmosphere. Intensity of secondary particles decreases with depth from the altitude of supersonic aircraft down to sea level (Fig. 6). In addition, solar energetic particles can contribute to the aircraft crew exposure through occasional SPEs. During such events an increase fluence of particles at aviation altitudes may be observed. In this context, of special interest is the validation of models for dose assessment of SEP events, using data from neutron ground level monitors (GLEs), in-flight measurement results obtained during a SPE and proton data measured by instruments onboard the satellites.

Early models and estimates
As known, cosmic radiation is included by the International Commission on Radiological Protection (ICRP) in radiation protection recommendations. The radiation exposures due to SEP events were estimated by many research groups (for a review see, e.g., Beck, 2007) with different models. Figure 7 shows a world dose map calculated for the radiation conditions (Lantos, 2006) during the GLE69 on 20 January 2005 based on the SiGLE model (Lantos & Fuller, 2004). Calculations have been made for the flight altitude of 40 kft (about 12,192 km).
Several years ago, Shea & Smart (2012) considered in detail the space weather aspects of practical importance in relation to an increased GLE occurrence rate in cycle 23. L.I. Miroshnichenko: J. Space Weather Space Clim. 2018, 8, A52 They presented the calculated radiation doses for aircrafts on the polar routes for each GLE event in the previous cycle. The space weather effects during large solar events in October and November 2003 are of special interest. The authors emphasize that it is important to use NM data (see also Sect. 9) in order to predict SPEs in the region of the most radiation hazardous SCR energies (several tens and hundreds of MeV). Such a prediction makes it possible to inform aircraft and spacecraft crews about an impending radiation hazard in advance.
As recently demonstrated, the expected computed effective dose rate at flight altitude during some major GLEs is highly dependent on assumed SEP spectra (Bütikofer & Flückiger, 2011, 2013. In addition, such characteristics of the GLEs, as SCR anisotropy and main flux arrival direction derived from ground-based NM measurements may considerably vary in time and space. As these characteristics are used as input parameters for the assessment of the radiation dosage, e.g., at flight altitudes, During the GLE69 the largest count rate increase was observed by the South Pole NM with more than 200% in the 5-minute data. This event is ranked as the second largest GLE after GLE05 during the last 50 years with gigantic count rate increases of several thousand percent at the south polar NM stations (Table 2). In addition, GLE69 exhibited strong north-south (N/S) anisotropy during the early phase of the event and showed a complex intensity-time profile in the data of the NMs of the worldwide network. Figure 8 demonstrates the computed equivalent rates during the main phase of GLE60 along two actual flights . In addition, the 5-minute data of the NM station Nain, which is located close to the flight routes during the GLE maximum phase, are plotted. The computations were based on the GLE parameters deduced by different research groups (Apatity, Athens, Australia, and Kiel). The Apatity group gives an accuracy of the SCR spectrum of ±25-30%. The comparison of the GLE analysis procedures and the achieved results does not allow to conclusively identify the causes for the differences in the deduced GLE characteristics.
As for the event of 20 January 2005, the same authors  have made similar calculations of the Fig. 6. Approximate values of ionizing radiation exposure at different altitudes for the polar regions due to interaction of galactic cosmic radiation with the Earth's atmosphere (Beck, 2007).  (Lantos, 2006). ambient dose equivalents (GCRs + SCRs) from the GLE onset until 08:00 UT based on the GLE parameters of the different groups for three assumed flights. The results are presented in Table 3. Although, e.g., the Apatity, Australia, and Bern groups use the same yield function, the derived GLE parameters differ significantly. More detailed information on the GLE analysis procedures is needed to find out the reasons for the large discrepancies in the results of GLE analysis and to improve the reliability of GLE analysis.
In a more recent paper, Bütikofer & Flückiger (2015) discussed in detail the various GLE analysis procedures and try to identify possible reasons for the differences in the results obtained. As they noted, the reasons for the spread of the GLE parameters deduced from NM data can be manifold and are at present still unclear. They include differences in specific properties of the various analysis procedures (e.g., NM response functions, different ways in taking into account the dynamics of the Earth's magnetospheric field), different characterizations of the solar particle flux near Earth as well as the specific selection of NM stations used for the analysis. The authors quantitatively investigated this problem for a time interval during the maximum phase of the GLE70 (13 December 2006). They presented and discussed the changes in the resulting GLE parameters when using different NM response functions, different model representations of the Earth's magnetospheric field as well as different assumptions for the solar particle spectrum and pitch angle distribution near Earth. Several NM yield functions have been published after 2013, when  made their analysis. The new yield functions (Mishev et al., 2013;Mangeard et al., 2016) may alter the earlier GLE analyses, but this assumption can only be checked by new calculations.

Advanced models
A new numerical model of estimating and monitoring the exposure of personnel due to secondary cosmic radiation onboard aircraft, in accordance with radiation safety standards as well as European and national regulations, has been developed quite recently by . The model aims to calculate the effective dose at flight altitude (about 12,000 km) due to secondary cosmic radiation of galactic and solar origin (GCRs and SCRs). In addition, the model allows the estimation of ambient dose equivalent at typical commercial airline altitudes in order to provide a comparison with reference data. The model is based on a straightforward full Monte Carlo simulation of the cosmic ray induced atmospheric cascade. The cascade simulation is performed with the PLANETOCOSMICS code (Desorgher, 2005(Desorgher, , 2007Desorgher et al., 2005). The flux of secondary particles, namely neutrons, protons, gammas, electrons, positrons, muons and charged pions was calculated.
A subsequent conversion of the particle fluence into the effective dose or ambient dose equivalent was performed as well as a comparison with reference data. An application of the model has been demonstrated, using a computation of the effective dose rate at flight altitude during the ground level enhancements of 20 January 2005 (GLE69), 13 December 2006 (GLE70) and 17 May 2012 (GLE71). The first two GLEs mentioned have been analyzed with taking into account the presence of two relativistic components in the SCR flux (prompt and delayed ones), and, correspondingly, two forms of derived spectra, as they were reconstructed by the Apatity research group (see, e.g., Vashenyuk et al., 2011).
It was shown that during the initial phase of GLE69 the effective dose rate at the flight altitude (39,000 ft or 12,000 m above sea level) was %150 lSv/h at North polar region and %1000 lSv/h at South polar region . During the late phase of this event the computed effective dose rate was still significant %200 lSv/h. Therefore, this event could significantly increase the potential biological risk of aircrew members and passengers. The computed effective dose rate    and references therein). The computed effective dose rate at flight altitude during the initial phase of GLE71 on 17 May 2012 (21.1 lSv/h) was greater than that for the late phase of GLE70. This value (21.1 lSv/h) is roughly three times greater than the average effective dose rate due to GCR for the same period of GLE71 (initial phase of the event, up to 02:00 UT).
The authors  used the new yield functions (Mishev et al., 2013). As noted above, the expected computed effective dose rate at flight altitude during some major GLEs is highly dependent on assumed SEP spectra. In addition, other model assumptions (Mishev & Velinov, 2010) lead to about 15% difference in secondary particle characteristics and, accordingly, in the effective dose rate. In this respect, the obtained results of computations are obviously in a satisfactory agreement with some previously reported early estimates (Bütikofer et al., 2008;Matthiä et al., 2009a,b).
In one of the latest works on this issue, the authors (Hubert & Aubry, 2017) are focused, first of all, on neutron fluxes as function of altitude. Secondly, atmospheric ionizing radiation impacts on biological doses during quiet period (GCRs) and extreme solar events (SCRs) are considered. On the basis of the modern models for GCRs and SCRs they calculated the ambient dose equivalents for different flight conditions. In particular, a GLE05 model was applied to London M New York flight dose calculations. Specific case of Antarctica is discussed because it combines both the high altitude and the very low magnetic field. Obtained results show that dose values vary drastically, on the one hand -with the route path (latitude, longitude, altitude), on the other hand -with the phasing of the solar event.
Based on these first results, it is important to quantify the impact of an extreme solar event. Thus, cartographies of ground ambient dose equivalent rate (in lSv/h) in Antarctica continent during GLE69 are presented in (Hubert & Aubry, 2017) for different time phases of the event: (a) onset T 0 , (b) T 0 + 15 min, and (c) T 0 + 30 min. Ambient dose equivalent rates increase drastically during the GLE event, particularly at the stations close to the specific topographic features of the Central Antarctica (http://www.gdargaud.net/Antarctica/ DomeCFAQ.html), such as Dome A, Dome B (Ridge B) and Dome C (see Fig. 2). This implies that the occupants of Concordia (France, Italy) and Vostok (Russia) stations suffered directly the GLE69 impact. However, this GLE is rather short, involving a total dose of the order of 50-100 lSv, that is low compared to the annual total dose.
Then, these results demonstrate that atmospheric ionizing radiation induced by CRs can be a problematic in Antarctica environment from the point of view of classical recommendations established for public. In the specific case of GLE69 event, calculations show limited impact on additional dose, in contrast to avionic altitudes. However, considering the limited number of monitored events and influences of the GLE anisotropy on dose assessment, this issue deserves special attention. In the authors' opinion (Hubert & Aubry, 2017), this highlights the importance of monitoring atmospheric ionizing radiation, more particularly extreme solar events, then to develop semiempirical and particle transport method for reliable calculation of dose levels.

Extreme SEP events and GLEs
In this section, we consider a number of actual problems related to extreme SEP events. First of all, it is necessary to present a modern concept of ''extreme SEP event'', especially considering the wide variety (diversity) of the SEP spectra and also proceeding from the requirements of practical cosmonautics. Then, we describe the observational situation with known ''ancient'' SEP events. Since it is naturally implied that possible ancient SEP events should be largest ones (highly likely -GLEs), in Section 7 we consider an applicability of so-called ''GLE05 scenario'' to SCR studies. With taking into account all data on the past extreme SEP events, we suggest a new distribution function of the occurrence probability of largest SEP events as a function of integral fluence U of protons with the energy of E ! 30 MeV, U(!30 MeV), over the entire available range of fluences. Finally, we discuss briefly a probability of super-flares at present Sun.

Concept of extreme SEP event
There are two different approaches (e.g., Panasyuk et al., 2018) to a definition of ''extreme event'' for solar cosmic rays, extreme SEP event, or extreme SPE. The first of them suggests that SPE may be considered as an extreme particle event when the worldwide network of ground-based CR stations has detected a significant flux and/or fluence of relativistic solar  . Earlier dose rate estimates are marked with asterisks (Bütikofer et al., 2008;Matthiä et al., 2009a,b protons in a given GLE. In this case, the question immediately arises of the magnitude of the flux (fluence) of relativistic SCRs. It is common to characterize the strength of a GLE as the peak intensity in percentage of the increase above the GCR background (see Table 2). However, for some purposes it is more useful to study the integral intensity of the events. The eventintegrated intensity I (in units %*h) is defined by Asvestari et al. (2017) as the integral of the excess above the GCR background over the entire duration of the event. It corresponds to the total fluence of SEPs with energy sufficient to cause an atmospheric cascade (several hundred MeV). It is necessary to emphasize that, while the peak intensity is important for the problem related to particle acceleration at/near the Sun and transport in the interplanetary medium, the total fluence is more relevant for the terrestrial effect (e.g., Usoskin et al., 2011;Mironova et al., 2015). Quite recently, Asvestari et al. (2017) suggested to use a ratio of two fluences, U(!30 MeV)/ U(!200 MeV), as a new characteristic of the GLEs.
An alternative definition of extreme SPE is based on the detection of a large enhancement of SEP flux at the Earth's orbit in the range of non-relativistic energies, more exactly -in the range of~10-100 MeV (for protons). Note that the SEP energy spectrum covers several orders of magnitude (e.g., Miroshnichenko, 2014) -from !10 6 to~10 11 eV, so relativistic particles are only a part of total SEP spectrum. Therefore, we use below a term ''SEP events'' as the most common one to describe the particles of solar origin. In the context of our motivations and goals of the present paper, one can suggest to select the extreme SEP events by the integral fluence of U(!30 MeV). This empirical choice is dictated mainly by the practice of space research, in particular, by the needs of practical cosmonautics.
Up to recent years, the so-called ''Carrington event'' (CE) of 1-2 September 1859 (AD1859) with the estimated integral fluence of U(!30 MeV) = 1.88 · 10 10 cm À2 (McCracken et al., 2001a,b) was suggested to be the most suitable for this purpose (e.g., Townsend et al., 2003Townsend et al., , 2006Smart et al., 2006;Miroshnichenko & Nymmik, 2014). In the paper by Smart et al. (2006) even a possible intensity-time profile of the >30 MeV solar protons associated with the Carrington flare of 1 September 1859 has been constructed. The absence of a statistically significant increase in the observed concentration of the cosmogenic nuclide 10 Be for 1859, in their opinion, indicates that the SCR produced in the CE had a soft spectrum, similar to other interplanetary shock-dominated events (e.g., GLEs in July 1959, November 1960, August 1972 and some others). If so, such an event has been considered as the ''worst-case'' reference for the radiation hazards in space. Indeed, the two nearest candidates for a role of ''worst-case'' are the events on 15 November 1960 (GLE11) and 4 August 1972 (GLE24). They were characterized by lesser values of U(!30 MeV), of about 9 · 10 9 cm À2 and 5 · 10 9 cm À2 , respectively .
The event AD1859 happened due to a powerful, so-called ''white'' solar flare (Carrington,1860). The white flares seem to have been associated with a very large SPE as originally suggested by McCracken (1959b).The importance of this flare could be characterized by its position on the solar disk relative to the Sun-Earth line (12°W of the central meridian), a large geomagnetic storm which was recorded 17 h later and an unusual display of the aurora which was reported as far south as the tropics (e.g., Botley, 1957;Vallance, 1992). Note that this flare, similarly to the outstanding flare of 23 February 1956 (solar cycle 19) happened at the ascending stage of solar cycle 10, two years before its maximum.
Of course, we have no direct evidence of an associated SEP event for the AD1859. But a certain correlation between U(!30 MeV) fluence and flare size based on modern data allowed to Cliver & Dietrich (2013) to obtain a best guess U(!30 MeV) value of~1.1 · 10 10 cm À2 (with the ± 1r uncertainty spanning a range from~10 9 -10 11 cm À2 ) for a composite (multi-flare plus shock) 1859 event. This value is approximately twice that of estimates/measurements -ranging from (5 Ä 7) · 10 9 cm À2 -for the largest SEP episodes (July 1959, November 1960, August 1972 in the modern era ).
These publications were preceded by numerous attempts to link the impulsive nitrate events observed in polar ice with large SPEs (e.g., Dreschhoff & Zeller, 1990;Zeller & Dreschhoff, 1995;Shea et al., 1999Shea et al., , 2006. The authors have examined impulsive nitrate ''spikes'' above the relatively large terrestrial nitrate background found in polar ice. Researchers consistently noted that observed impulsive nitrate depositions are delayed by only a few weeks relative to the occurrence of a large fluence SPE. Another group (Kepko et al., 2008(Kepko et al., , 2009 suggested that the first four large GLEs, starting from GLE01, are also manifested in the form of impulsive cosmogenic nitrate ''spikes'' above the relatively large terrestrial nitrate background found in a Summit, Greenland ice core drilled in the summer of 2004. However, a few years after, those publications gave rise to the ''Great Debates'' about the reliability of nitrate data for identification of ancient SPEs, including the event AD1859. As it turned out, these debates are not only of methodological interest, but are of a fundamental nature. According to the current thinking, dominating at present, the nitrate method is not suitable at all for studying past SEP events (e.g., Wolff et al., 2012Wolff et al., , 2016Duderstadt et al., 2014Duderstadt et al., , 2016aMekhaldi et al., 2015Mekhaldi et al., , 2017Sukhodolov et al., 2017;Katsova et al., 2018). In their opinion, this is partly due to equivocal detection of nitrate spikes in single ice cores and possible alternative sources, such as biomass burning plumes, especially in the case of AD1859 event (Wolff et al., 2012). This approach is strongly disputed by other researchers (e.g., Smart et al., 2014Smart et al., , 2016McCracken & Beer, 2015;Melott et al., 2016;. On the other hand, several years ago quite certain signatures of cosmic-ray increase in the years 774-775 (thereinafter the event AD775) were found by the 14 C annual data from tree rings analysis which was carried out in Japan, Europe, Russia and America (Miyake et al., 2012;Usoskin et al., 2013;Jull et al., 2014). Integral fluences U(!30 MeV) were independently estimated as 4.5 · 10 10 cm À2 , 8.0 · 10 10 cm À2  and (2.1 Ä 3.0) · 10 10 cm À2 (Miroshnichenko & Nymmik, 2014). The event AD775 has been also measured by polar ice-core annual proxies of 10 Be and by 36 Cl data with a cerca 5-year resolution (e.g., Sigl et al., 2015;Mekhaldi et al., 2015;Sukhodolov et al., 2017). It is timely to mention that ice cores from where the radionuclide data come are mainly from Greenland (project GRIP for the 36 Cl measurements, and projects NGRIP, Tunu, and NEEM -for the 10 Be measurements). Some 10 Be data come from the project WAIS Divide (Antarctica). In addition, some researchers, again by the 14 C data, have found two other large SEP events in the past, namely AD994 (Miyake et al., 2013) and even more ancient one BC3372 (Wang et al., 2017). Both the AD775 and AD994 events were also measured with 10 Be and 36 Cl data. Integral fluence U(!30 MeV) for the event AD994 was about 1.3 · 10 10 cm À2 (Mekhaldi et al., 2015(Mekhaldi et al., , 2017. The intensity of the event BC3372 is estimated (Wang et al., 2017) to be about 0.6 times as large as the AD775 event, i.e., about (1.5 Ä 1.7) · 10 10 cm À2 . It is easily seen that all those estimates are comparable to the estimate by McCracken et al. (2001a) for the event AD1859.
From the very beginning, the most researchers believed that the event AD775 is due to solar activity and represents an extreme large SPE or a short sequence of several SPEs (e.g., Thomas et al., 2013). On the other hand, some doubts appeared about the reliability of above fluence estimates and even about the nature of the event AD775 (e.g., Pavlov et al., 2013;Cliver et al., 2014;Liu et al., 2014), and the discussion on this matter continued for some time. The hypothesis by Liu et al. (2014) of a cometary impact was quickly turned down (Melott, 2014;Usoskin & Kovaltsov, 2015). A multi-proxy analysis (Mekhaldi et al., 2015;Sukhodolov et al., 2017) seems to put the end to this doubts: it was shown that the AD775 was a SEP event (or a short consequence of events) with a hard spectrum, close to that of GLE05 or even harder. The energy spectrum of the event of AD775 was assessed by Mekhaldi et al. (2015) as very hard, with the U(!30 MeV) being~(2.5 Ä 2.8) · 10 10 cm À2 (Mekhaldi et al., 2015(Mekhaldi et al., , 2017. This value seems to be currently the best possible estimate of the fluence U(!30 MeV) for the AD775. It may be seen that this value practically coincides with our earlier estimate of 2014 (see above) obtained from independent considerations.
Be that as it may, it is obviously that now the event AD775 can pretend for the role of the ''worst case''. On the other hand, based on modeling results and measurements in different ice cores, Sukhodolov et al. (2017) have concluded that, even with a fine temporal resolution and excluding post-depositional processes, a nitrate signal of the 775AD event is not observable in polar ice core records. Quite recently, Mekhaldi et al. (2017) presented new continuous high-resolution measurements of nitrate and of the biomass burning species ammonium and black carbon, from several Antarctic and Greenland ice cores. They have thoroughly investigated periods covering the two largest known SEP events -AD775 and AD994, as well as the Carrington event and the hard SEP event of February 1956, and found no coincident nitrate spikes associated with any of these benchmark events.

Limitations of nitrate method
As noted above, many researchers believe that, at present, the nitrate data from polar ice should be considered completely unsuitable for the identification of ancient SPEs and for quantitative estimates of their properties and should be dismissed. Moreover, as discussed by several authors (e.g., Beer et al., 2012;Usoskin & Kovaltsov, 2012), the very high SEP fluence proposed by McCracken et al. (2001a,b) for the outstanding space weather event AD1859 contradicts the available cosmogenic isotope data. In particular, Usoskin & Kovaltsov (2012) have analyzed annual radiocarbon data from different sources and annual 10 Be series from Berggren et al. (2009). They found that those data cannot resolve the Carrington peak. Therefore, Usoskin & Kovaltsov (2012) conclude that the cosmogenic data do not support the hypothesis of a very strong SPE related to the Carrington flare. One of the objections against this conclusion was mentioned above: it might happen due to a composite (multi-flare plus shock) nature of the 1859 event. Another argument is that, in analogy to the GLE05 (Smart et al., 2014), the fluence U(!30 MeV) for the AD1859 was not strong enough to produce large nitrate signal.
Two events from the list by Usoskin & Kovaltsov (2012), ca. 780 AD and 1460 AD, appear in different series that makes them strong candidates for extreme SPEs. From here these authors received a new strict observational constraint on the occurrence probability of extreme SPEs. Practical limits can be set as U(!30 MeV) % 1, 2 Ä 3, and 5 · 10 10 cm À2 (i.e. 10, 20 Ä 30, and 50 times greater than GLE05) for the occurrence probabilities of 10 À2 , 10 À3 , and 10 À4 yr À1 , respectively. Given these uncertainties, the results by Usoskin & Kovaltsov (2012) should be considered with a precision of up to an order of magnitude.
In our opinion, a brief discussion of potential possibilities and limitations of the nitrate method is of independent interest. In general, in the framework of the review article, it seems appropriate to describe a number of the methodological and physical issues associated with the SEP studies by the nitrate method and to give, at least, the outlines of the history of relevant discussions. For example, as shown in McCracken et al. (2001a,b), there were five events in the interval 1561-1950 that equal or exceed so-called ''streaming limit'' for the >30 MeV fluences. Of these SPEs there are good magnetic and solar records for the events of 1859 (CE), 1895 and 1896. The ''streaming limit'' effect (Reames & Ng, 1998Reames, 1999Reames, , 2013Reames, , 2017 seems to be caused by a possible resonance interaction between the particle flux and the plasma waves at the interplanetary shock. Although this effect is discussed since 1998 and seems do really work, up to now it still remains an elegant hypothesis, and no more. It is timely to emphasize that the streaming limit does not depend at all upon the nature of the particle source itself, but only on the intensity and spectra of the accelerated particles (Reames & Ng, 2010). As estimated by McCracken et al. (2001a,b), the streaming limit is effective at the fluences U(!30 MeV) = (6 Ä 8) · 10 9 cm À2 . At the same time, they recognized that the CE and several other SPEs displayed the fluences greater than their estimates.
All those values of the fluences U(!30 MeV) have been obtained by measurements of nitrate anomalies (excess) in the thin layers of the Greenland and Antarctic ices (McCracken et al., 2001a,b). The layers have a characteristic short timescale (<6 weeks) and are highly correlated with periods of major solar-terrestrial disturbance, the probability of chance correlation being less than 10 À9 . The authors did not estimate the errors of the fluence values, but indicate a distinguishable (minimal, background) fluence for the nitrate event of 31 August 1956 at the level of U(!30 MeV) = 5 · 10 8 cm À2 . As known, this event coincided with the small, but quite discernible GLE06 (McCracken, 1959a), with the fluence U(!30 MeV) = 2.5 · 10 7 cm À2 . With this likely threshold for the detection of SPEs in ice cores and with the fluence of 1.88 · 10 10 cm À2 for the CE, we get a ratio 5 · 10 8 cm À2 /1.88 · 10 10 cm À2 = 2.63 · 10 À2 . It means that an estimated error of the fluence determination for the Carrington event is below 3%. By the way, the time of occurrence of the Carrington nitrate event is estimated with confidence to be 1859.75 ± 0.2 (McCracken et al., 2001a).
Those data underwent the analysis and interpretation in a number of studies (e.g., Townsend et al., 2003Townsend et al., , 2006Smart et al., 2006;Kepko et al., 2008Kepko et al., , 2009Cliver & Dietrich, 2013;Miroshnichenko & Nymmik, 2014). In practice immediately, however, serious doubts appeared (e.g., Wolff et al., 2008Wolff et al., , 2012Schrijver et al., 2012) on the validity of the ''nitrate signal'' itself for the Carrington event and, respectively, on the reliability of all fluence estimates by McCracken et al. (2001a,b). Especially, Wolff et al. (2012) have analyzed 14 ice cores from Greenland and Antarctic focusing on the period around the CE and found that there is no signal in most of them, and only a few datasets from South/Central Greenland depict a nitrate peak which is, however, accompanied also by peaks in ammonium, formate, black carbon and vanillic acid of biomass burning plumes, primarily from North America. Their final conclusion was: ''The Carrington event not observed in most ice core nitrate records''.
In turn, Schrijver et al. (2012) stated: ''We have had to conclude that nitrate concentrations in polar ice deposits cannot, at present, be used to extend the direct observational records of SEP events to a longer time base without at least significantly more study''. In addition, they suggest to evaluate ''the probabilities of large-energy solar events by combining solar flare observations with an ensemble of stellar flare observations''. This rather pessimistic paper, in essence, moves the discussion into another, very ''unreliable'' field of solar super-flares (see Sect. 7.3).
Meanwhile, independent nitrate data have been obtained earlier with a high-resolution continuous flow analysis (CFA) of a Summit, Greenland ice core drilled in the summer of 2004 (Kepko et al., 2009). The authors additionally included data from Windless Bight, Antarctica providing an interhemispheric comparison. As noted above, during the period of 1940-1950 four GLEs were recorded by cosmic ray ionization chambers. All four of these were time associated with significant impulsive nitrate enhancements in both the Greenland and Antarctic ice core data within 1-2 months after the GLE (Kepko et al., 2009). This interhemispheric correlation is strong evidence for a global response. These results and observational materials presented earlier , in our opinion, still deserve attention as experimental support of the arguments by McCracken et al. (2001a,b) that there is a short time delay between the generation of massive amounts of nitrates throughout the polar atmosphere and their deposition in polar ice.
At the same time, there are serious doubts in validity of CFA method (Wolff et al., 2016) because of possible distortions of the nitrate signal (''smoothing'' effects). In addition to biomass burning plumes, Wolff et al. (2016) make a number of other caveats. The last ones are strongly determined by the conditions of chemicals' transportation and deposition. For example, if the chemicals are dry deposited, they may give a very thin layers between snowfalls. If they are deposited by wet deposition, then they may give a signal that initially has the width of a snowfall (typically millimeters to centimeters). Also, the month of snowfall may be important (seasonal effect). After snowfall, snow may be redistributed through wind, leading to mixing of different snowfall layers. Nitrate that is deposited as nitric acid is known to be mobile in the snow pack, so that even sharp peaks become smoothed through vapor redistribution. Finally, an influx of marine air containing high concentrations of sea salt will scavenge acidic nitrate from the atmosphere, leading to deposition of a nitrate spike even in the absence of any primary input of nitrate.
Undoubtedly, above doubts and disputes have serious reasons, and it seems to be not accidental that a principal debate on methodological topics is still continuing (e.g., Smart et al., 2014Smart et al., , 2016Duderstadt et al., 2014Duderstadt et al., , 2016aMelott et al., 2016;Wolff et al., 2016;Mekhaldi et al., 2017;Sukhodolov et al., 2017;Katsova et al., 2018;. In particular, Smart et al. (2014) did not accept the ''closure'' of the event AD 1859 by Wolff et al. (2012) and brought new arguments in favor of their nitrate method. Smart et al. (2014) demonstrated that the data used by Wolff et al. (2012) cannot detect impulsive nitrate events because of resolution limitations. The authors (Smart et al., 2014) suggested to reexamine the data from the top of the Greenland ice sheet at key intervals over the last two millennia, with attention to fine resolution, and replicate sampling of multiple species.
In an effort to understand the different interpretations between McCracken et al. (2001a) and Wolff et al. (2012), the authors (Smart et al., 2014) concluded that a correlation test between impulsive nitrate spikes in polar ice (Dreschhoff & Zeller, 994, 1998) and SPEs could be conducted by examining ice core nitrate deposition during the only years when both well-documented extremely large relativistic SCR events with hard spectrum occurred and ultrafine-resolution nitrate measurements in consolidated firn exist, namely the 14 year interval 1937-1951, as well as for extremely large SPEs between 1956 and 1989. They use these time periods to ascertain if the ice core data available for analysis contain evidence of short-duration impulsive nitrate enhancements that could be time associated with these large relativistic SCR events. They believe that this will allow further insight into polar depositional processes on a sub-seasonal scale, including atmospheric sources, transport mechanisms to the ice sheet, post-depositional interactions, and a potential SPE association. For example, the nitrate signal of GLE05 can be seen at submonthly resolution (see Fig. 1 from Smart et al., 2014). This signal is quite distinct, but not very impressive (rather small) because of a fluence of U(!30 MeV) in this case was not large -about 10 9 cm À2 only . It is noteworthy that independent estimates of the fluences for the GLE05 event also indicate a low value of U(!30 MeV) in this event (see Fig. 3, right panel, Mewaldt et al., 2007).
A critique of the Wolff et al. (2012) analysis by Smart et al. (2014) was replied by Wolff et al. (2016). These authors agree with Smart et al. (2014) that very sharp nitrate peaks may be detected in a system with very high resolution. However, in the opinion by Wolff et al. (2016), there is no evidence that these arise from SEPs. Liquid conductivity data, in fact, confirm that many sharp peaks are not acidic, making it likely that they are a by-product of influx of sea salt, biomass burning, or dust aerosol, rather than a consequence of a significant influx of nitrate. So, there seems no a priori way to distinguish SEP signals from other causes of nitrate enhancement. Wolff et al. (2016) again emphasize that the evidence that was used previously, regarding the coincidence in date of SEP events and nitrate spikes, is no longer valid. In any case, they have shown clearly that such spikes are not reliably deposited in adjacent cores. ''This makes the use of such sharp features to log the power and frequency of SEP events impractical''. At the same time, Wolff et al. (2016) do not doubt that the largest SEPs will cause an enhancement of nitrate in the middle atmosphere and in extreme cases -in the troposphere. However, they cannot see a plausible route for identifying and using nitrate deposition in snow to diagnose past SEPs. In their opinion, ''the use of 10 Be seems a much more promising path (e.g., Usoskin et al., 2013), although at much lower time resolution''.
Advocating the alternate use of cosmogenic isotopes, Wolff et al. (2016), however, omit recent work by McCracken & Beer (2015). These authors have identified all of the same events in the 1940-1956 time interval, that have been previously identified as impulsive nitrate events in Figures 1 and 7 of Smart et al. (2014). Recall that the 10 Be and 14 C records are annual; the 10 Be increases appear 12-15 months after the high-energy SPEs. As noted by Smart et al. (2016), including fine-resolution nitrate signatures together with the 10 Be and 14 C records would enhance the overall database of historic high-energy SPEs and might refine the dating and enable estimates of the spectral shape and fluence of past events. Wolff et al. (2016) mention that they await more detailed modeling studies, which are in progress. Such a study of the 23 February 1956 event has recently been conducted . The authors have computed the amount of nitrate available from this SPE using full spectrum modeling involving high-energy showers. The results indicate that nitrate production in the stratosphere is roughly 2-10 times greater than estimated by past analytic models. Summing the ionization from the surface to approximately 45 km at Summit gives the nitrate deposition found in the integrated winter 1956 GISP2-H peak. Therefore, this event can account for the measured nitrate with timely (~2 months) transport downward to the surface via mechanisms involving the polar vortex, denitrification, and polar stratospheric cloud. The authors come to the conclusion that the nitrate spike in GISP2-H data is mostly acidic and likely due to the event of 23 February 1956.
Recently, Duderstadt et al. (2016a) used the Whole Atmosphere Community Climate Model (WACCM) to calculate how large an event would have to be to produce enough ''odd nitrogen'' throughout the atmosphere to be discernible as nitrate peaks at the Earth's surface. It has been shown that SEP events cannot lead to a claimed very sharp nitrate spike (duration of about a fortnight) since nitrate is mostly produced in the stratosphere, so that the expected increase, if any, should have a duration of several months, up to a year. Finally, Sukhodolov et al. (2017), using state-of-the-art CR cascade and chemistryclimate models have successfully reproduced the observed variability of cosmogenic isotope 10 Be around AD775, in four ice cores from Greenland and Antarctica, thereby validating the models in the assessment of this event. It has been shown that there is no nitrate peak in five analyzed nitrate series from Greenland and Antarctica for the event of AD775. Note that two of series were at a high resolution, about 20 samples per year. Therefore, Duderstadt et al. (2016a) conclude that ''nitrate data cannot provide a quantitative index of SEP events, and the results based on this dataset should be dismissed''.
One characteristic example of discrepancies in the approaches to the study of ancient SPEs and differences in the results can be demonstrated by several recent works Sinnhuber, 2016;Duderstadt et al., 2016a,b). As noted by Sinnhuber (2016), both groups of authors Duderstadt et al., 2016a) make similar assumptions about the probable atmospheric ionization by large SPEs with hard spectra, and about the subsequent production of nitric oxides in the lower and middle atmosphere. Nevertheless, both come to opposing conclusions. Sinnhuber (2016) does not assess which of these two approaches is more realistic. In his opinion, they both have different problems attached to them, ''but do not seem, on the face of it, implausible''. He suggests to leave it to the authors, or to future research, to conclude which provides the more realistic assessment.
In their turn, Smart et al. (2016) noted that it is still incorrect to compare calculated percentage nitrate increases (Duderstadt et al., 2016a) relative to all NO y species in the atmospheric reservoir, because it consists largely of organic nitrates (Jones et al., 2011), which have a long residence time compared to HNO 3 (Ridley et al., 2000). If one uses a percentage enhancement, this dilutes the expected effect. Instead, calculations of the absolute amount of additional, SPE-produced, inorganic, nitric acid available for prompt deposition are consistent with the observed nitrate increase . Nevertheless, main conclusion from Duderstadt et al. (2016b) is again: ''Existing and previous studies that utilize nitrate peaks in the ice core record to identify individual SPEs are flawed''.
This, however, does not necessarily exclude a possibility to use nitrate records, for example, as indicators of past solar activity (Dreschhoff & Zeller, 1998), e.g., in the period of the Maunder minimum (Dreschhoff et al., 1997) and as proxy for the studies of centennial GCR variability (Traversi et al., 2016). In particular, Traversi et al. (2016) presented the first direct comparison of cosmogenic 10 Be and chemical species in the period of 38-45.5 kyr BP spanning the Laschamp geomagnetic excursion from the EPICA-Dome C ice core. A principal component analysis (PCA) allowed to group different components as a function of the main sources, transport and deposition processes affecting the atmospheric aerosol at Dome C. Moreover, a wavelet analysis highlighted the high coherence L.I. Miroshnichenko: J. Space Weather Space Clim. 2018, 8, A52 and in-phase relationship between 10 Be and nitrate at this time. The evident preferential association of 10 Be with nitrate rather than with other chemical species was ascribed to the presence of a distinct source, here labeled as ''cosmogenic''. The evidence that the nitrate record from the EPICA-Dome C ice core is able to capture the Laschamp event hints toward the possibility of using this marker for studying GCR variations and open up new perspectives in paleoclimatic studies.
It is appropriate to note here a great diversity of sources for the production of nitrates and minor constituents in the Earth's atmosphere (e.g., Krivolutsky & Repnev, 2012;Ogurtsov & Oinonen, 2014), multidisciplinary nature of the problem itself (e.g., Sigl et al., 2015) and existence of a number of unresolved and/or poorly studied issues/questions of SCR interaction with the terrestrial atmosphere (e.g., Miroshnichenko, 2008;Mironova et al., 2015;Duderstadt et al., 2016a).
Finally, we would like to cite several important notes by Beer et al. (2012, p. 260): ''. . .we note that one scientists' proxy may be another investigators' noise''. Indeed, as well-known, ''. . .the cosmogenic production rate is modulated by solar activity and on longer time scales by the geomagnetic field intensity. However, this production signal is not directly recorded in an archive. Before being archived, the cosmogenic radionuclides are subject to various transport processes which modify the original production signal (Chapter 13). A solar physicist may call these modifications ''noise'' because they disturb the solar proxy that is of interest. On the contrary, for an atmospheric scientist, it may be exactly these modifications that are of interest, because they may assist in the improvement of atmospheric mixing models and lead to a better understanding of stratosphere-troposphere exchange. This is a clear example of the old truth that ''somebody's noise is somebody else's signal''.
In dismissing nitrates, while acknowledging their production by SPEs in the middle atmosphere and troposphere, Wolff et al.'s (2016) main conclusion looks unproven, and their call for excluding nitrates from 10 C and 10 Be studies is premature. Be that as it may, the authors  reiterate and stand by their main point: ''Low time resolution analysis of polar ice cores cannot detect impulsive nitrate events.''

SEP events and space weather
Note that in all cases of the largest ancient SEP events we deal with several significant difficulties that can affect the results of analysis and their interpretation. First of all, the experimental information obtained from different proxies is rather heterogeneous; second, all estimates are highly dependent on the data quality (the averaging period, magnitude of statistical error, etc.); finally, the estimates of SEP characteristics can depend on the accepted model (scenario) of the event. On the other hand, the modern technological implications of such events may be extreme. As a result, in spite of some doubts about the reliability of above fluence estimates for the events AD775, AD1859 and others, this make us to use them (among the all relevant data) to construct new distribution function of SEP events by the values of U(!30 MeV) over the entire available range of fluences (Sect. 7.2).
It is important to note that the largest SPEs are produced, as a rule, by multiple sequenced solar flares and can be accompanied by a number of prominent geophysical disturbances (auroras, geomagnetic and ionospheric storms, etc.). Storm-time magnetic field depression changes the paths of the SEPs and extends the region of energetic particle penetration inside the magnetosphere. As a result, strong geomagnetic storms (GMS) can significantly enhance the effects of SPE. In particular, according to Tsurutani et al. (2003), the 1-2 September 1859 magnetic storm was the most intense geomagnetic disturbance in the recorded history. One presently used measure of such geomagnetic disturbances is the D st index, which measures the globally averaged change in the magnetic field near Earth's equator. Tsurutani et al. (2003) showed that the 1 September 1859 flare most likely had an associated extremely strong coronal mass ejection (CME) which led to the magnetic storm in the Earth's magnetosphere with D st~À 1760 nT. This is consistent with the response of DH = 1600 ± 10 nT, registered at Colaba magnetic observatory (near Bombay-Mumbai, India). However, more careful analysis by Siscoe et al. (2006) has suggested that an equilibrium solution of the equation describing the solar wind forcing of the magnetosphere would lead to D st~À 850 nT for the Carrington event.
It is generally accepted now that the September 1859 solarterrestrial disturbance was, in fact, the first recognized space weather event (e.g., Cliver & Svalgaard, 2004). For example, the auroras associated with the 1859 magnetic storm occurred globally and have been reported by many researchers (e.g., Tsurutani et al., 2003). Moreover, the associated magnetic disturbances produced numerous disruptions of telegraph transmissions, which attracted much public attention and were widely reported (e.g., Boteler, 2006). The first magnetic disturbance started in the evening of August 28 and telegraph operations were disrupted in North America and Europe until the next morning. The second disturbance started with a sudden storm commencement at 04:40 UT on September 2 and a major disturbance followed immediately. Telegraph operation in Europe was severely affected by the initial magnetic disturbance. Instead, the North American telegraphs were affected by the second phase of the disturbance. Note that telegraph was introduced in the 1840s and underwent rapid expansion, so by 1859 there were lines across North America, Europe, and parts of Australia and Asia. Both the August 28 and September 2 disturbances had a widespread impact on telegraph operation.
It was also found that both the flare energy as well as the speed of associated CME were extremely high but not unique during the event on 1-2 September 1859. Other events with more intense properties have been detected; thus magnetic storm of this or even larger intensity may occur again. Because the data for the high-energy tails of solar flares and geomagnetic storms are extremely sparse, the tail distributions and therefore the probabilities of occurrence cannot be assigned with any reasonable accuracy. Finally, Tsurutani et al. (2003) concluded that a serious complication is a lack of knowledge of the saturation mechanisms of flares and magnetic storms.
Extreme geomagnetic storms can significantly extend the region of SEP penetration into the high-latitude magnetosphere (e.g., Panasyuk et al., 2018). For example, changes in the radiation environment at low extraterrestrial orbits (or LEO -Low Earth Orbit) during a solar proton event depend on the power of the solar event itself and the efficiency of the penetration of SEPs inside the Earth's magnetic field. The stronger a magnetic storm, the larger the area affected by penetrating SEPs will be observed. As known, sometimes such events (GLEs) can be detected on-ground by the network of neutron monitors. Such observations can be used to verify the magnetic field models (e.g., Desorgher et al., 2009).
On 23 July 2012, solar active region 1520 (~141°W heliographic longitude) gave rise to a powerful CME with an initial speed that was determined to be 2500 ± 500 km/s (Baker et al., 2013). The eruption was directed away from Earth toward 125°W longitude. STEREO-A sensors detected the CME arrival only about 19 h later and made in situ measurements of the solar wind and interplanetary magnetic field. The authors (Baker et al., 2013) address the question of what would have happened if this powerful interplanetary event had been Earthward directed. Using a well-proven GMS forecast model, they found that the 23-24 July event would certainly have produced a geomagnetic storm that was comparable to the largest events of the twentieth century (D st~À 500 nT). Using plausible assumptions about seasonal and time-ofday orientation of the Earth's magnetic dipole, the most extreme modeled value of storm-time disturbance would have been D st = À1182 nT. This is considerably larger than estimates by Siscoe et al. (2006) for the well-known Carrington storm of 1859 (see above). This finding has far reaching implications because it demonstrates that extreme space weather conditions such as those during March of 1989 or September of 1859 can happen even during a modest solar activity cycle such as the one presently underway. How be it, the event of 1859, undoubtedly, joins the annals of modern powerful geomagnetic storms such as 4 August 1972 that has had severe impacts on power grids, satellites, and communication systems.

Fluence distribution functions
In connection with the increasing of the duration and range of space missions, interest in the radiation safety of space crews and electronics of space vehicles has strongly increased in recent years. In this context, an occurrence rate of large SEP events in the past is of special importance from the point of view of their probability in the future. In this section, we summarize available data, estimates and models on this issue.

Largest SEP events and GLE05 scenario
About 11 years after the fundamental studies by McCracken et al. (2001a,b), the authors (Usoskin & Kovaltsov, 2012) have published a new list of 23 candidate events for large (powerful) ancient SPEs with the fluences of U(!30 MeV) ! 10 10 cm À2 , that could have occurred in the past at the time scale of 100-1000 years. The values of U(!30 MeV) have been estimated by proxy data on cosmogenic isotopes 14 C and 10 Be. It was assumed that the event-candidate developed according to the so-called ''scenario GLE05'', i.e., to the scenario of the outstanding event of 23 February 1956 which is still considered to be the largest in the flux and the fluence of particles in the region of relativistic energies (Table 2). Note, however, that the ''GLE05 scenario'' is vulnerable to criticism -first of all, because the spectrum of non-relativistic protons on 23 February 1956 was fairly flat (e.g., Bailey, 1959;Miroshnichenko, 1996Miroshnichenko, , 2001, and the fluence U(!30 MeV), as mentioned above, reached only~10 9 cm À2 (see also Fig. 3, right panel, Mewaldt et al., 2007). In turn, certain evidence was also obtained that the flux of relativistic solar protons during GLE05 consisted of two distinct components: prompt and delayed ones (Vashenyuk et al., 2008). At the same time, all events from the list  deserve the most thorough further study.
In our opinion, this discussion cannot be completed without taking into account, at least, one important point, namely, a great diversity (variety) of the slopes of observed SEP-GLE event spectra (e.g., Miroshnichenko, 1996). Quite recently, Asvestari et al. (2017) have shown that all the strong GLE events (the strength I > 100% · h) are characterized by a hard spectrum with the U(!30 MeV)/U(!200 MeV) ratio being 10-50. Meanwhile, if we consider all list of 71 GLEs (their Table 1), the U(!30 MeV)/U(!200 MeV) ratio, in fact, changes in rather wide range -from about 10 (GLE08 of 4 May 1960) to about 500 (GLE24 of 4 August 1972). On the other hand, the most of modern calculations of radionuclide concentrations are based mainly on the two SPE scenarios -GLE05 and GLE24 (e.g., . Note that the hard SEP-SCR spectrum assumes a flattening in the range of non-relativistic energies, as it happened, e.g., in the GLE05 (Bailey, 1959). The flattening inevitably leads to lower values of the proton fluence U(!30 MeV). In turn, this will affect some major geophysical effects (production of nitrogen oxides and cosmogenic isotopes, depletion ozone layer, etc.). On the other hand, some processes in the atmosphere, except for the ionization, have a resonance character, for example, a generation of the stable isotope 7 Be (Webber et al., 2007). These authors performed detailed calculations of the production rate of the 3 H, 7 Be, 10 Be, and 36 Cl cosmogenic isotopes, using the FLUKA (Monte Carlo code) software and taking into account recent data on the interaction crosssections for vertically incident protons with energies varying from 10 MeV to 10 GeV. This made it possible to study the isotope production due to SCRs and GCRs in the range of low energies, where the production rate is a very sensitive energy function. The results of new calculations show interesting patterns for SCR production at energies below~1 GeV.
As known, for the studies of ancient SEP events the proxy data on the cosmogenic isotopes 10 Be and 36 Cl (along with 14 C) are of the greatest interest. Based on the data for the SPE series of October-November 2003, Webber et al. (2007) have shown that the 10 Be production rate reached a maximum during these events when the proton energy was~100 MeV. At the same time, the differential yields for 7 Be and 36 Cl (Figs. 1 and 2 from Webber et al., 2007) show an excess around 30 MeV (0.24 GV). The 7 Be and 36 Cl isotopes are produced by splitting the target nuclei of atmospheric nitrogen 14 N and argon 40 Ar, respectively, due to proton interactions in the atmosphere. More precisely speaking, the 7 Be and 36 Cl isotopes were shown to produce more intensively at the energy of~25 MeV. In this sense we can speak on some kind of ''resonance effect''. If the SCR spectrum is steeper than in October-November 2003, the maximum of production rate will shift toward lower energies. If the spectrum is less steep (as was registered, e.g., on 20 January 2005, Mewaldt et al., 2005a,b), the production peak will shift to higher energies. L.I. Miroshnichenko: J. Space Weather Space Clim. 2018, 8, A52 Solar cosmic ray spectra have been reevaluated by Webber et al. (2007) for all of the major events that occurred since 1956. Finally, as it was shown, in terms of yearly production of 10 Be, only the GLE05 makes a major contribution, namely, about 40% of the total production (GCRs + SCRs) above 65°o f geomagnetic latitude in 1956. For 36 Cl the SCR contribution is generally larger because of a resonance near 25 MeV for production from 40 Ar. Note that the production above 65°latitude is largest in 1960 where the SCR production was~110% of the yearly production by GCRs. For complete latitudinal mixing the effects of SCRs are almost completely washed out for 10 Be.
In general, the different energy sensitivities of the cosmogenic isotope production along with the SCR energy spectrum itself suggest a hierarchy of responses to SCR events. For the events detected supposedly by their nitrate production (e.g., McCracken et al., 2001a,b; see Sects. 6.2 and 6.3), the peak response appears to be at~10 MeV, for SCR detection by 7 Be or 36 Cl production in the atmosphere the peak response is 15-25 MeV and for detection by 10 Be production the peak response is~100-200 MeV. ''Thus SCR events detected by one method may not be seen using another and vice versa'' (Webber et al., 2007).
In recent years, many researchers are strongly fascinating by the results of Tylka & Dietrich (2009). The authors have developed a new technique for analyzing data from the world-wide NM network. The method was used to derive absolutely normalized event-integrated proton spectra for 53 of the 66 GLEs recorded since 1956. Unfortunately, these authors (as far as we know) did not publish important methodological details of the Band function representation (Band et al., 1993) of solar proton spectra in GLEs. Also, they distinctly recognized: ''Compared to other techniques for analyzing GLEs, our method has many shortcomings. Most importantly, we do not take careful account of anisotropies that are generally observed at the beginnings of events; instead, we rely on the world-wide network and wider time-binning to average out these effects, conceptually analogous to what happens with a spinning satellite. We also neglect other factors, such as rapid evolution of the spectrum in the initial stages of the event and careful numerical integration over each station's rigidity dependent asymptotic viewing cone''.
In spite of those restrictions, quite recently, Raukunen et al. (2018) have presented two new statistical models of high energy solar proton fluences based on GLE observations during solar cycles 19-24. As the basis of their modeling, the authors utilized a four parameter double power-law function (known as the Band function) fits to integral GLE fluence spectra in rigidity. For each GLE, a least-squares fit to a power-law in rigidity is performed on the corrected fluences, and the value of spectral index c that gives the best fit is selected. As an example of the corrected NM fluences, the integral spectrum of GLE71 (17 May 2012) for the proton fluences is shown in Figure 9. The NM observations are shown in orange, GOES observations in green (MEPAD) and red (HEPAD), and the Band-fit spectrum in blue.
In the first model, the integral and differential fluences for protons with energies between 10 MeV and 1 GeV are calculated using the fits, and the distributions of the fluences at certain energies are modeled with an exponentially cut-off power law function. In the second model, Raukunen et al. (2018) used a more advanced methodology: by investigating the distributions and relationships of the spectral fit parameters the authors found that the parameters can be modeled as two independent and two dependent variables. Therefore, instead of modeling the fluences separately at different energies, one can model the shape of the fluence spectrum. The authors presented examples of modeling results and showed that the two methodologies agree well except for short mission duration (1 year) at low confidence level. Raukunen et al. (2018) also showed that there is a reasonable agreement between their models and three wellknown solar proton models -JPL91 (Feynman et al., 1993), ESP (Xapsos et al., 2000) and SEPEM (Crosby et al., 2015), -despite the differences in both the modeling methodologies and the data used to construct the models.

New distribution function
In the light of the above debate (Sect. 6.3) and existing uncertainties, we do not see serious obstacles to use all available data and estimates on confirmed ancient events AD775 and AD994, as well as on ''problematic'' event AD1859, align with data on some other events from the two lists compiled by McCracken et al. (2001a,b) and Usoskin & Kovaltsov (2012). In spite of ''hypothetical'' nature of some fluence estimates (e.g., for the AD1859 event), or model dependency of some of them (e.g., for the event AD775) we propose to apply all of them tentatively within a trial empirical model for the construction of modern distribution function for the SPEs on their fluences of U(!30 MeV). It would allow to understand the applicability of those data, at least, at the test level. We believe that such an approach would be also helpful to resolve old and new doubts discussed in Section 6.3. Moreover, Fig. 9. Event-integrated proton fluence spectrum for GLE71 in rigidity. NM observations are shown in orange, GOES observations -in green (MEPAD) and red (HEPAD), and the Band-fit spectrum in blue (Raukunen et al., 2018). L.I. Miroshnichenko: J. Space Weather Space Clim. 2018, 8, A52 we incline to consider the studies by Kepko et al. (2008Kepko et al. ( , 2009 and Smart et al. (2006Smart et al. ( , 2014Smart et al. ( , 2016, as well as the observation materials and ideas expressed therein, in favour of the methodology and results obtained by McCracken et al. (2001a,b), and as a base for further studies of relevant problems.
Therefore, all of the data described in Section 6.2 have been recently used by Miroshnichenko & Nymmik (2014), Panasyuk et al. (2018,  and in preset work to construct new distribution function of the occurrence probability of largest SEP events on the U(!30 MeV) values over the entire available range of fluences (Fig. 10). Basic features of our approach and main sources of fluence data have been described in detail earlier (e.g., Nymmik, 1999aNymmik, ,b, 2011Miroshnichenko & Nymmik, 2014). These authors undertook their own analysis of the fluence data for different solar cycles. They also compared available distribution functions for the integral fluences of U(!30 MeV), namely, JPL91 model by Feynman et al. (1993) and Emission of Solar Protons (ESP) model by Xapsos et al. (2000) with the model of SINP MSU (Nymmik, 1999a,b;Miroshnichenko & Nymmik, 2014). The results of comparison are shown at the left panel of Figure 10. Red curve corresponds to the modified SINP model.
To approximate the fluence data we applied a function obtained by Nymmik (1999b): W(!U) = (U/10 6 ) Àc · exp(ÀU/U 0 ), where c is a power-law index and U 0 -characteristic exponential constant. As it was shown by Lu et al. (1993), such a form of distribution function seems to be universal for the description of many manifestations of solar flares (peak fluxes and/or energy fluences in X-ray and radio wave bursts, in proton and electron emissions, etc.). Solid red line in Figure 10 marks this approximation which turned out to become reasonable with the following parameters: c = 0.32; U 0 = 7·10 9 cm À2 (Miroshnichenko & Nymmik, 2014). Those parameters are, in practice, identical to that obtained earlier (e.g., Nymmik, 1999b) by analyzing the events measured onboard the satellites only, namely, the value of U 0 = 6 · 10 9 cm À2 .
The right panel of Figure 10 covers all recent changes in the U(!30 MeV) estimates for the confirmed and/or hypothetical large SEP events in the past. In particular, if we were try to use the maximal U(!30 MeV) estimates as 4.5 · 10 10 cm À2 , 8.0 · 10 10 cm À2 , then the fluence data for the event AD775 would be on the same red curve (blue column), but the column would be slightly shifted lower and to the right in comparison with the proposed position of hypothetical fluence from the event AD1859. By the way, empirical estimate of U(!30 MeV)~1.1 · 10 10 cm À2 by Cliver & Dietrich (2013), within the error bars, apparently rests well on the red curve (Fig. 10, right panel). Full and open triangles at the right panel demonstrate our estimates of U(!30 MeV) based on the data by Kiraly & Wolfendale (1999) for the fluences of protons at energies ! 60 MeV, with the extrapolation into the past for 1 Myr and 100 Myr, respectively. It seems obvious that the estimates by Kiraly & Wolfendale (1999) are greatly overestimated. Difference in the energies of protons (between 30 and 60 MeV) makes this discrepancy even much more visible.
On the other hand, recently  presented a combined cumulative occurrence probability distribution of large SEP events based on three timescales: directly measured SEP fluences for the past 60 years; estimates based on the terrestrial cosmogenic radionuclides 10 Be and 14 C for the multi-millennial (Holocene) timescale; and cosmogenic radionuclides measured in lunar rocks on a timescale of up to 1 Myr. These three timescales yield a consistent  Feynman et al., 1993); Emission of Solar Protons (ESP) (model by Xapsos et al., 2000); SINP MSU -modified model (Nymmik, 1999a,b;Miroshnichenko & Nymmik, 2014). Right panel: Space era data (red dots) + Greenland ice data (blue rhombus) (Miroshnichenko & Nymmik, 2014), including Carrington event. Event AD775 -blue column, upper and lower estimates, 2.0 · 10 10 Ä 3.0 · 10 10 cm À2 (adapted in this work from Miroshnichenko & Nymmik, 2014, and from assessment by Mekhaldi et al., 2015). Full and open triangles demonstrate the extrapolation of integral fluences of U(!30 MeV) into the past for 1 Myr and 100 Myr, respectively (estimated by Miroshnichenko & Nymmik, 2014, from the data by Kiraly & Wolfendale (1999) for the protons ! 60 MeV). An orange asterisk corresponds to the upper estimate of U(!30 MeV) from .
L.I. Miroshnichenko: J. Space Weather Space Clim. 2018, 8, A52 distribution. The data suggest a strong roll-over of the occurrence probability, so that the SEP events with a fluence of U(!30 MeV) > 10 11 protons cm À2 yr À1 are not expected on a Myr timescale. Applying to this estimate our formulae for W(!U) (see above), we get the probability of such event about 1.5 · 10 À8 (orange asterisk). As one can see, this value is located, in practice, on the red curve (Fig. 10, right panel).
As follows from Figure 10, the modified model of SINP MSU seems to provide the best agreement with the data on ancient extreme SEP events. In this context, one important question arises: can our Sun produce a superflare with the fluence of U(!30 MeV) considerably larger, or even much larger, than the AD775 event? With our current understanding of the problem, there is no direct answer to this question.

Probability of super-flares at present Sun
Using Kepler mission data, Maehara et al. (2012) discovered 365 superflares (flares with energy ! 10 33 erg) on 148 solar-type stars. Among them, they found 14 superflares (Fig. 11) on Sun-like stars (slowly rotating G-type main sequence stars, which have rotational periods longer than 10 days and surface temperature of 5600 K < T eff < 6000 K). From this, they estimated that the occurrence rate of superflares with the energy of 10 34 erg is once in 800 yr, and that of 10 35 erg superflares is once in 5000 yr on Sun-like stars. Simple calculations by Shibata et al. (2013), in combination with the results of analysis by Maehara et al. (2012), may indirectly suggest that the Sun is able produce super flares of 10 34 erg once in 800 years (recall that the energy of the observed largest solar flares is usually about 100 times less). Recent publications (see below) seem do not support those preliminary estimates.
A detailed discussion of stellar flares is out of the scopes of the present review. Nevertheless, in the context with above considerations, in this section we, at least, briefly fix some relevant problems and ideas concerning the probability of solar superflares with present level of solar activity. Exhaustive statistics of the Kepler mission on superflares at other Sun-like stars with solar rotation rates (Maehara et al., 2012;Shibayama et al., 2013;Nogami et al., 2014;Wichmann et al., 2014, and many others) have greatly reinvigorated the discussion of whether or not such events could in principle also occur at the Sun.
Observations of superflares on Sun-like stars are challenging for solar physics in general and for the solar dynamo -in particular. On one hand, the energy of the solar magnetic fields is clearly insufficient to produce a superflare, if one assumes a reasonable efficiency of the flare production (e.g., Aulanier et al., 2013;Schmieder, 2017). On the other hand, cosmogenic isotope studies have identified, at least, two historical events (AD775 and AD994) that were shown to be 30-50 times stronger by the fluxes of SEPs than the most energetic SEP events observed instrumentally (e.g., Usoskin et al., 2013;Mekhaldi et al., 2015). It is suggested that these events may form the upper limit of the intensity of SEP events during the last eleven millennia, while the strongest stellar superflares are substantially stronger. In addition to the events AD775 and AD994, a new one ca. 3372 BC was found recently (Wang et al., 2017) from ancient buried tree rings in China.
The largest solar flares observed over the past few decades have reached energies of a few times 10 32 erg, possibly Fig. 11. Comparison between the occurrence rate of superflares on G-type stars and those of solar flares . The solid-line histogram shows the frequency distribution of superflares on Sun-like stars. The error bars in the histogram represent the square root of the event number in each bin. This distribution can be fitted by a power-law function with an index of %1.5 ± 0.3 (thick solid line). The dashed line, dotted line, and dot-dashed line indicate the power-law distribution of solar flares observed in EUV (Aschwanden et al., 2000), soft X-rays (Shimizu, 1995), and hard X-rays (Crosby et al., 1993), respectively. The red line corresponds to a power-law function with an index of 1.8 (from theoretical estimate by Aschwanden, 2012). Largest solar flare marked by blue oval.
L.I. Miroshnichenko: J. Space Weather Space Clim. 2018, 8, A52 up to 10 33 erg. Flares at active Sun-like stars reach up to about 10 36 erg. In the absence of direct observations of solar flares within this range, complementary methods of investigation (e.g., dynamo models, radionuclide data etc.) are needed to assess the probability of solar flares beyond those in the observational record. As known, recent data from the Kepler mission has revealed the occurrence of superflares at Sun-like stars which exceed by far any observed solar flares in release of energy. Meanwhile, the distribution of superflares over the studied star population was found very inhomogeneous (e.g., Candelaresi et al., 2014). This fact may indicate that they may have some specific dynamo features not applicable to the Sun (e.g., Kitchatinov & Olemskoy, 2016;Katsova et al., 2018). In addition, the predicted flare energy significantly exceeds a theoretical limit for the Sun (e.g., Aulanier et al., 2013;Livshits et al., 2015).
In particular, Aulanier et al. (2013) evaluated from observations that 30% of the area of sunspot groups are typically involved in flares. This is related to the strong fragmentation of these groups, which naturally results from sub-photospheric convection. When the model is scaled to 30% of the area of the largest sunspot group ever reported, with its peak magnetic field being set to the strongest value ever measured in a sunspot, it produces a flare with a maximum energy of~6·10 33 erg. Livshits et al. (2015) have presented their estimates based on the observations and the methods of the non-linear force-free field extrapolation. They obtained the upper estimate of the energy of flares, which are able to occur in a given large active region. For the Sun, this value does not exceed~3·10 32 erg. Because the average longitudinal magnetic fields of G stars (at an age of 1-2 Gyr) are ten-times stronger than the maximum magnetic field of the Sun as a star, the energy of solar-type stellar flares can only slightly exceed 10 34 erg. Note also that dynamo activity is known to decrease with the star's effective temperature (Kitchatinov & Olemskoy, 2016), which then leads to less frequent and less energetic flares.
According to Schmieder (2017), with our Sun as it is today, it seems impossible to get large enough sunspots in order to produce super-flares with energy >10 34 erg. The author considers different sunspot groups under the assumption that the magnetic field in the bipole (sunspot) may be used as boundary condition of the MHD simulation. This approach provides the flare energies of 10 35 -10 36 erg, as it is recorded for other stars by the Kepler satellite. Finally, it may be concluded that in order to produce stronger flares, the Sun-like stars should have a much stronger dynamo than the Sun and a rotation rate exceeding several days. In this context, the prediction of having extreme solar storms in 800 years would be very speculative.
Thus, the question on application of the Kepler stellar data to assessment of solar superflares remains open. It is interesting to note, however, that superflares at Sun-like stars, solar flares, microflares, and nanoflares (see Fig. 11) are roughly on the same power-law line with an index of 1.8 (solid red line) for a wide energy range from 10 24 erg to 10 35 erg.

Absolute fluxes of relativistic protons in GLEs
From the point of view of maximum possibilities of solar accelerator(s), the GLEs are of specific interest (Miroshnichenko, 2003b). At present, estimates of maximum absolute intensities (fluxes) of the relativistic solar protons, I max (>1 GV), in the so-called proton flux units pfu (1 pfu = 1 proton/(cm 2 s sr), are available, in practice, for all GLEs registered for the entire period of ground-based observations. The intensity values for 34 GLEs registered in the solar cycles 17 through 22 have been obtained earlier (Miroshnichenko & Pérez-Peraza, 2008). Here we complement them by our new estimates based on the data from several lists and/or Catalogues of SPEs (for details see Miroshnichenko & Yanke, 2016). Some estimates were made by the data on two relativistic components of SCRs, prompt (PC) and delayed (DC) ones, as they were defined long ago ; the results for the GLEs of 1956-2006 have been summarized in Vashenyuk et al., (2011) andMiroshnichenko et al., (2013). At last, in recent publications (Li et al., , 2015 the estimates have been obtained for one small GLE of cycle 24 registered on 17 May 2012 (GLE71). The event of 6 January 2014 was also taken into account as a ''candidate'' in the GLE72. So, finally we analyze data of 71 events (Fig. 12). It should be recalled that the estimates of I max (>1 GV) have been obtained with an accuracy factor of about 2.0 only, i.e., the absolute error of I max values is about 100%. As it turned out later, the event of 6 January 2014 was rather a sub-GLE, whereas the officially recognized GLE72 is now the event of 10 September 2017.
One interesting aspect here is that most of the values of I max (>1 GV) in Figure 12 are higher than the integral intensity of galactic cosmic rays (GCR) at R ! 1 GV, that is, I GCR (>1 GV) ! 1 pfu, as was estimated and plotted earlier (Miroshnichenko, 2003b). On the other hand, Figure 12 exhibits a rather flat integral size distribution of GLEs in the lowintensity range. This fact emphasizes once more the importance of studying the GLEs of low intensity (''hidden GLEs'') for better understanding of the SEP spectrum formation, especially in the range of relativistic energies (e.g., Li et al., 2015).

Weak and hidden GLEs
GLE events with relatively small amplitudes of increase on the Earth's surface are of certain interest. As two characteristic examples we note the SPEs of 19 and 24 September 1977 (GLEs No. 28 and 29, respectively). They were, apparently, the first SCR increases recorded on the Earth's surface in the solar cycle 21. According to the hourly data of the NM South Pole, the maximum amplitude of increase did not exceed 3.2% and 11.8% on September 19 and 24, respectively. Due to a small amplitude of effects and complicated helio-geophysical conditions in September 1977, the analysis of these events presents certain difficulties. But at the same time, as it was shown (Kepicova et al., 1982), certain quantitative information on the SCR spectrum may be obtained on the basis of the method of integral multiplicities, m(R, h), of generation of secondary CR particles (this term means the number of secondary particles, e.g., neutrons, generated by a primary particle of rigidity R and registered at the depth h, e.g., at the sea level). A similar analysis was carried out by Miroshnichenko (1990) for the three other small events: 21 August 1979, 10 April 1981, and 10 May 1981.
Note also a search for weak GLE54 (2 November 1992) by the data of Mexican NM at very high cutoff rigidity R c = 8.27 GV (Vargas-Cárdenas & Valdés-Galicia, 2012). They performed a search for GLEs on the full 5-minute database of the Mexico City NM using wavelet filters and two different statistical tests. A detailed analysis of the time series of 2 November 1992 allowed to find a previously unreported increment matching the onset time of the impulsive phase of GLE54, thus providing evidence of an effective detection of high energy SCRs. Solar cycle 23 was remarkable with a record number of GLEs (in total 16) in general, including several rather weak events, namely, GLE56, 57 and 58 (2 May, 6 May and 28 August 1998, respectively), as well as GLE66 (29 October 2003) and GLE68 (17 January 2005). Of special interest are several effectively registered and (probably) ''hidden'' (or unidentified) GLEs of the current cycle 24 (Li et al., 2015). Apparent paucity of GLEs in this cycle even inspired some authors (Atwell et al., 2015) to suggest a new terms for the events of this kind -''sub-GLE''. Moreover, in the light of the well-known discussions on the possibility of particle acceleration up to energies of ! 1.0 GeV by shocks in the solar corona (e.g., Zank et al., 2000;Berezhko & Taneev, 2003, 2013Li et al., 2013;Thakur et al., 2014) the studies of weak and sub-GLEs become especially interesting. Table 5 shows a tentative list of some SPEs of solar cycle 24 as was presented by Li et al. (2015), with important additions from other papers. In particular, we include the results by Makhmutov et al. ( , 2015 on international observations with the CARPET cosmic ray detector, as well as observations on several sub-polar and high-latitude neutron monitors (Belov et al., 2015a,b). Also, we demonstrate one of the recent attempts to predict GLEs in cycle 24 (Pérez-Peraza & Juárez-Zuñiga, 2015). At last, we consider the backside CME of 23 July 2012 with rather short Sun-Earth shock transit time (18.5 h) and associated SPE with a maximum intensity of >10 MeV protons about 5000 pfu (Gopalswamy et al., 2016). In some cases we used information (e.g., about CMEs) from the NOAA website (NOAA, 2018).
The data of Table 5 are not complete, rather controversial and some of them are questionable. For example, certain doubts cause the results of forecasting (Pérez-Peraza & Juárez-Zuñiga, 2015) for the events 13 and 15. According to these authors (their Table 3), the event 13 formally falls into the interval of prognosis between 1 February and 30 June 2015, but no real SPE (or GLE) was observed at that time (http://legacy-www.swpc.noaa. gov/weekly/pdf/prf2075.pdf), as far as we know. As to the event 15, they give the interval of prognosis between 11 August and 2 December 2015. On the other hand, quite recently Mishev et al. (2017) announced that two NMs (South Pole and Dome C) have registered an increase in the count rate caused by SCR. Nevertheless, despite of existing uncertainties, we believe that this information will be helpful for further thorough studies of weak and hidden GLEs (or sub-GLEs). It is especially true if we take into account a theoretical possibility that some portion of SEPs can be accelerated up to relativistic energies by CME-driven shocks. Such a possibility has been assumed, e.g., for GLE71 (17 May 2012), as it was considered by Li et al. (2013). The same assumption was applied to the interpretation of the sub-GLE occurred on 6 January 2014 (Thakur et al., 2014;Li et al., 2016). The event 16 deserves special attention: it occurred at the descending stage of the cycle 24 (close to its minimum), and from the point of view of space weather all disturbed period of 4-10 September 2017 is of certain interest.

SCR and GLEs in prediction schemes
In the middle of the 1980s, among numerous techniques and schemes of helio-geophysical predictions, many researchers started to consider the idea of using ground-based CR observations in order to make short-term forecasts of different phenomena. Specifically, they put forward several interesting proposals to use relativistic solar protons (R ! 1 GV) as SPE predictors in the non-relativistic range. Thus, Dorman et al. (1990) for the first time considered the possibility of diagnosing interplanetary medium and predicting the SPE onset based on the solution of the inverse problem of SCR propagation. Based on observations   (Kurt et al., 2018) up to the SCR (GLE) maximum in the Earth's orbit, it was assumed, first, to restore the SCR emission function and, then, to predict the SPE development for several hours ahead. Although the methodical aspects of such an approach were sufficiently justified, it remained unclear how the proposed scheme could be verified based on observational data. One of the difficulties was that a large flux of relativistic protons was not always accompanied by the same increase in the flux in the non-relativistic region (see, e.g., Fig. 2). Subsequently, many researchers considered the role of relativistic SCRs in prognostic schemes (Belov & Eroshenko, 1996;Dorman & Zukerman, 2003;Mavromichalaki et al., 2009;Vashenyuk et al., 2011;Veselovsky & Yakovchuk, 2011;Pérez-Peraza et al., 2011). Belov & Eroshenko (1996) developed an empirical approach to the determination of the solar proton spectrum near the Earth in the 10 MeV-10 GeV energy range, directly using observational data without any preliminary assumptions regarding the possible spectral shape. Their method also made it possible to reconstruct the intensity time profile for protons with any energy. Vashenyuk et al. (2011) tried to predict the form of the maximum flux spectrum for non-relativistic protons during SPEs using data on the spectrum of delayed component (DC) in the corresponding GLE. In this case, the spectrum in the 500 MeV range was considered as a natural smooth continuation of the DC spectrum. Using the GLE47 event (21 May 1990) as an example, these researchers indicated that the DC spectrum is in good agreement with its extrapolation into the region of low energies ( 430 MeV), where direct measurements on the Meteor spacecraft and stratospheric balloon measurements were performed. The authors proposed a reasonably limited calculation model, where data of~20 NMs are used. This model makes it possible to obtain real-time SCR spectra with accuracy sufficient for routine prediction and to automatically solve a number of space weather problems.
However, the next studies indicated that the scheme described above is still piloting (tentative) and cannot be used to solve the problem in detail. In particular, Veselovsky & Yakovchuk (2011) verified the application of the method in (Mavromichalaki et al., 2009). The method, based on NMDB, was developed in order to give early warning that solar protons with E p~1 0-100 MeV approached the Earth. A post-event analysis and comparison with the observations from 2001 to 2006 indicated (Veselovsky & Yakovchuk, 2011) that more than 50% of SPEs were omitted in the case of using such a prediction method. To increase the reliability of this method, it is necessary to use additional data on the state of solar and heliospheric activity.
A new tool for short-term predicting the occurrence of GLEs has been suggested quite recently by Núñez et al. (2017). This real-time tool, called HESPERIA UMASEP-500, is based on the detection of the magnetic connection, along which protons arrive in the near-Earth environment. This is done by estimating the lag-correlation between the time derivative of 1-minute soft X-ray flux (SXR) and the time derivative of 1-minute near-Earth proton fluxes observed by the GOES satellites. One of important advantages of the suggested scheme is that, unlike current GLE warning systems, this tool can predict GLE events before the SCR detection by any NM station. The prediction for the interval 1986-2016 is presented for two consecutive periods, because of their notable difference in the tool performance. For the 2000-2016 period, this prediction tool gave a probability of detection of 53.8% (7 of 13 GLE events), a false alarm ratio of 30.0%, and average warning times of 8 min with respect to the first NM station's alert and 15 min to the GLE Alert Plus's warning (for more details see original papers by Núñez, 2015;Núñez et al., 2017). The authors have tested the model by replacing the GOES proton data with SOHO/EPHIN proton data, and the results turned out to be similar in terms of probability of detection, false alarm ratio and average warning times for the same period.
An interesting new approach to the problem of GLE prediction has been suggested by Pérez-Peraza et al. (2011) who used the specific features of variations in the GLE registration rate (see Fig. 4). The nearest goal of these authors was to develop a method for predicting such events during cycle 24. A tentative prediction indicated that the next event (GLE71) would be registered between 12 December 2011 and 2 February 2012. This event factually occurred on 17 May 2012; it was small and was only observed at high latitudes. The maximum enhancement (~23% according to 5-min NM data) was registered at the South Pole station.
Further development of this idea was undertaken by Pérez-Peraza & Juárez-Zuñiga (2015) who used the results of wavelet analysis of the GLE registration rate (see Sect. 3.5 and Fig. 4) and applied so-called ''fuzzy logic tools''. Their predictions, however, were somewhat disappointing. For example, the authors predicted that the next GLE (after GLE71) would occur in the first half of 2015 (see Table 5 above), but nothing of this kind has happened at that time. In fact, so far (October 2018) the worldwide network of NMs has confidently registered only two rather small GLEs -on 17 May 2012 (GLE71) and 10 September 2017 (GLE72).
In recent years, to develop predictive schemes, researchers have increasingly turned to fluxes (fluences) of protons with energies ! 100 MeV (Veselovsky & Yakovchuk, 2011;Belov et al., 2015b) and ! 200 MeV . Fluence ratios U(!30 MeV)/U(!200 MeV) are also used (see, e.g., Asvestari et al., 2017). Note that such analysis is carried out ''inside'' the SCR spectrum that has already been formed near the Earth for this event. Therefore, its results, in our opinion, cannot be used reliably for serious prognostic recommendations. At the same time, such analysis is important to understand the nature of the ''kinks'' in the SCR spectrum (see, e.g., Miroshnichenko & Nymmik, 2014, and references therein).
In the context of above consideration, it seems to be timely to note some prospects of the problem of radiation hazard in the light of available data on the total weakening of the solar activity in the solar cycle 24 and in the next coming decades. Now the Sun's arrhythmia, in its different parameters and displays, is widely discussed at the level of observational data, empirical estimates and theoretical expectations (e.g., Popova et al., 2017;Panasyuk et al., 2018), but this topic would go far beyond the scope of our review.

New definition of GLE
From the 1970s up to present, the conventional definition of GLE events requires a detection of SEP by at least two L.I. Miroshnichenko: J. Space Weather Space Clim. 2018, 8, A52 differently located neutron monitors. At present, as a GLE definition is usually taken a statistically significant increase of the count rate of at least two differently located NMs, which is accompanied with an increase in data from space-borne instruments. In 2014 the authors (Li et al., 2015) suggested the term ''hidden GLE'' in order to indicate the problem of search for the weak GLEs that are yet to be discovered (found) in the data of the world network of NMs at the boundary of detection limit, with given interval of data averaging. This is especially interesting and important for the SEP events of the weak current cycle 24.
In July 2015, another research group (Atwell et al., 2015), instead of traditional GLE definition, introduced the division of GLEs into actual GLEs, sub-GLEs and sub-sub-GLEs. According to (Atwell et al., 2015), GLEs are extremely energetic SPEs having proton energies often extending into several GeV range and producing secondary particles in the atmosphere, mostly protons and neutrons, observed with ground NMs. Sub-GLE events are less energetic, extending into the several hundred MeV range, but without producing detectable levels of secondary atmospheric particles. Sub-sub-GLEs are even less energetic, with an observable increase in protons at energies greater than 30 MeV, but no observable proton flux above 300 MeV. Our term ''hidden GLE'' (Li et al., 2015) is close to the term ''sub-GLE''. With our understanding, the event on 6 January 2014 (Thakur et al., 2014;Li et al., 2016) was a typical sub-GLE. All this diversity in the terms, in our opinion, requires to return to the problem of GLE definition, taking into account the high efficiency (accuracy) of CR registration by NMs and new experimental capabilities for SCR observations (e.g., Poluianov et al., 2015) on high-altitude polar stations.
As noted in Section 3.1, some places at the Earth's surface are exceptionally well suitable for ground-based detection of SEPs, -high-elevation polar regions with negligible geomagnetic and reduced atmospheric energy/rigidity cutoffs. Cosmic-ray detectors located at high altitude in the polar regions possess higher sensitivity to low-energy cosmic rays than those located near sea level and/or at non-polar locations. There are two regions with these properties in the world: the top of the Greenland ice sheet (3205 m a.s.l), and the Antarctic plateau (average elevation of about 3000 m a.s.l). At present, there is no high-altitude cosmic-ray station in Greenland, but there are two NM stations in such locations on the Antarctic plateau (see Fig. 2 above): SOPO/SOPB (at Amundsen-Scott station, 2835 m elevation), and DOMC/DOMB (at Franco-Italian research station Concordia, 3233 m elevation). The neutron monitor station, called ''Dome C'', started continuous operation in early 2015. Since then, in particularly, a relatively weak SEP event on 29 October 2015 that was not detected by sea-level NM stations was registered by both SOPO/SOPB and DOMC/DOMB, and it was accordingly classified as a GLE. This would lead to a distortion of the homogeneity of the historic GLE list and the corresponding statistics.
Indeed, another small GLE happened on 10 September 2017, and it was officially classified as GLE72 (http://gle.oulu.fi). To address this issue, Poluianov et al. (2017) propose to modify the GLE definition so that it maintains the homogeneity, namely, a GLE event is registered when there are near-time coincident and statistically significant enhancements of the count rates of at least two differently located neutron monitors, including at least one neutron monitor near sea level and a corresponding enhancement in the proton flux measured by a space-borne instrument(s). Relatively weak SEP events registered only by high-altitude polar neutron monitors, but with no response from CR stations at sea level can be classified as sub-GLEs. New definition apparently does not affect the present list of GLEs, because all 72 ''official'' events (see, e.g., the site http://gle. oulu.fi) have been registered by at least one of near sea level NMs.
As far as follow from the paper (Poluianov et al., 2017), the authors suggest to define the GLEs and sub-GLEs as ''statistically significant enhancements of the count rates'' based on 5-minute NM data during the maximum phase of the event, as it was applied, e.g., by Mishev et al. (2013) in their analysis of the GLE71 (17 May 2012). This point should be emphasized in order to avoid possible misunderstandings with time interval of data averaging (see Sect. 3.3). Another important note concerns ''a corresponding enhancement in the proton flux measured by a space-borne instrument(s)''. This formulation should be added by distinct indication, at least, on the threshold energy and corresponding proton flux, as well as on suggested (possible) type of spacecraft (LEO, GOES, ACE etc.).
The GLE definition is crucial, e.g., for the studies based on GLE occurrence rate (see, e.g., Fig. 4). As a useful addition to this subject, we would propose to search for ''hidden GLEs'' (or ''sub-GLEs'') in the existing data of the global network of NMs in past years, especially, in current solar cycle 24. This would clarify some important issues in the physics of generation and transportation of SCR, as well as in their interaction with the Earth's atmosphere.

Conclusions and suggestions
As shown above, solar energetic particles (SEPs) represent an important aspect of solar-terrestrial physics, as well as are among the three main components of space weather. The results of our critical analysis allow to make conclusions of great physical, methodological and/or applied interest.
1. As it was mentioned above, there is no unambiguous relation between the fluxes (fluences) of relativistic SCR (!500 MeV) and non-relativistic SEPs (below 100 MeV). This fact is one of the serious challenges to the problem of estimation and prediction of radiation hazard in space. 2. Some GLEs (apparently closely related to so-called ''rogue'' SEP events) may have an interplanetary (not solar) origin, and their diagnostics and prediction deserve special studies. 3. On the other hand, recently, based on the techniques of five different research groups (Apatity, Athens, Bern, Kiel, Hobart), radiation doses were calculated for three real air routes from two GLEs (15 April 2001 and20 January 2005). The discrepancies in the dose estimates by different methods turned out to be depressingly large. This definitely indicates the existence of significant methodological problems in the analysis and use of NM data. The estimates diverge mainly because of inaccurate determination and consideration of the shape of the SCR intensity spectrum in the Earth's orbit at various times of the proton event. Some new progress in the dose calculations was achieved in the latest works by . 4. Special attention should be paid to the well-known dilemma ''Flare-CME'' from the point of view of GLE production by these two (alternative?) sources. At present, there are compelling evidences that both flares and CMEs are products of magnetic reconnection processes, usually occurring at the Sun in active magnetic field regions of great complexity. 5. Although the situation with weak GLEs remains uncertain, their study seems to be very promising, first of all, for understanding the mechanisms of SCR acceleration. In this context, the SCR community needs a new, more modern definition of the GLE event (see above special Sect. 10). This task seems to be quite feasible by the NMDB community efforts. 6. The results of the recent studies are expected to yield a basis for the reduction in the spread of the GLE parameters deduced from NM data. Recently, the future of the Neutron Monitor Data Base (NMDB) has been discussed at a special meeting ''10 years of NMDB'' (Athens, 20-23 March 2017, http://www.nmdb.eu). In particular, there were considered the software for reliable GLE alerts and SPE prediction with inclusion of real-time NM measurements; the definition of GLE and a new class of ''sub-GLE'' events (see also Sect. 10); some new suggestions on registration, identification, interpretation, physical mechanisms and possible sources of GLEs. 7. As to the methodological and physical aspects of the problem of ancient large SEP events (GLEs), the current state of research on the matter is difficult to recognize as satisfactory. It is enough to recall that: 1) all available data on the ancient SPEs are proxy, or indirect, ''by definition''; 2) both of the two methods (nitrate and multi-proxy radionuclide analysis) are burdened by a number of serious restrictions, and 3) fluence estimates are strongly depending on the event ''scenario'' that is being applied. 8. A promising idea (and difficult research task) is to quantify the relative contribution of SEPs into the total nitrate signal for a given SPE (GLE), e.g., for the welldocumented event of 20 January 2005 (e.g., Vashenyuk et al., 2011;Krivolutsky & Repnev, 2012). Within such an approach, it would be reasonable to reexamine thoroughly such SPE features as the N/S asymmetry of SEP invasion, their energy spectrum, the duration of the event, and to assess again the rates and absolute amounts of nitrate production, e.g., with the WACCM model. This will make it possible to turn difficulties and challenges into incentives for the development of our ideas about ancient proton events.
In 76 years of ground-based observations (from February 1942 to present -October 2018) a total 72 GLEs were confidently registered. GLE events are being intensively studied in many countries. For example, in 2012, in a special issue of Space Science Reviews (v. 171) a collection of 7 articles on various aspects of GLE research was published. Unfortunately, a number of poor-studied, forgotten and/or somewhat neglected questions and/or problems, as well as controversial issues of SCR-GLE studies remained outside of this review. In particular,  presented an overview of the observed properties of the GLEs and the two associated phenomena, viz., flares and CMEs, both being potential sources of particle acceleration. In turn, Kahler et al. (2012) examined statistically the associated flare, active region (AR), and CME characteristics of about 40 GLEs observed since 1976 to determine how the GLE e/p and Fe/O ratios, each measured in two energy ranges, depend on those characteristics. This interest in the problem undoubtedly reflects its fundamental nature and great practical importance.