© S. Musset et al., Published by EDP Sciences 2023

1 Introduction

The Sun accelerates electrons, protons, and ions to suprathermal energies in a variety of transient events, ranging from small flares to coronal mass ejections (CMEs). On rare occasions, the particle energies become relativistic, which means some GeV to a few tens of GeV for nucleons. These particles are of particular astrophysical interest because they challenge our understanding of the acceleration process.

Relativistic nucleons can be detected on the Earth through the cascades of secondary particles that they trigger in the atmosphere (reviews in Lopate, 2006; Belov et al., 2010). Neutron monitors on the Earth are the standard equipment to measure the cascades triggered by particles with energies above a few hundred of MeV, depending on the geographical location (Bütikofer, 2018). These relativistic nucleons, called cosmic rays, have either a galactic or a solar origin. Galactic cosmic rays are registered at all times by neutron monitors, with slow variations of their intensity over time induced by the varying magnetic fields of the Sun and the Earth. On the other hand, the contribution from relativistic solar protons is only seen during rare events called ground-level enhancements or ground-level events (GLEs). Seventy-three GLEs have been detected between 1942 and 2021.

Besides their astrophysical interest, cosmic rays are relevant to space weather, since the secondaries produced in the atmosphere are the source of enhanced radiation doses at flight altitudes. For this reason, the radiation doses received by aircrew are monitored. A real-time space weather service was established under the auspices of the International Civil Aviation Organization (ICAO) in 2019, which among other space weather disturbances covers the monitoring and real-time assessment of radiation doses.

The time profile of a GLE is usually simple, with a rise to levels that were found to range from just above the galactic cosmic ray level to several tens of that level. The count rate rises within a duration of a few minutes to several hours, and then generally decays on longer time scales. Detailed studies of GLE time profiles revealed a close correlation between the rise time and the decay time, with a two to three times longer decay than rise (Moraal et al., 2015; Strauss et al., 2017). But Moraal, McCracken and coworkers (McCracken et al., 2008, 2012; Moraal & McCracken, 2012) also found that in some cases, typically when the parent active region on the Sun was not too far from the root of a Parker spiral through the Earth, rather impulsive and sometimes very strong initial peaks occurred with the rise and initial decay on minute time-scales.

In this article, we re-visit the time profile of GLEs, in an attempt to derive general features that could be used in space weather applications. The relationship between decay time and rise time established by Moraal et al. (2015) and Strauss et al. (2017) is a hint that this should be possible. To this end we present, after a short description of neutron monitor data and their analysis in Section 2, the evaluation of a median GLE time profile (Sect. 3) and a statistical analysis of rise and decay times (Sect. 4). In Section 5, we search for statistical relationships between the GLE times and the parent eruptive solar activity. The results are discussed with respect to a physical interpretation and possible space weather applications in Section 6.

2 Observations and data analysis

The time histories of neutron monitor count rates during historical GLEs are publicly available at the GLE database at the University of Oulu¹ and the Neutron Monitor Database (NMDB) hosted by the University of Kiel². The GLE database has more events, but many historical GLEs were only observed with a cadence of 5 min, while GLEs of solar cycle 23 onwards are provided with 1-minute integration by NMDB. In the first step of our analysis, we selected GLEs using time profiles with 5-minute integration.

Neutron monitor count rates depend on the mass of air above the instrument, which determines the absorption of the cascade triggered by the primary energetic particle (see Hatton, 1971; Bütikofer, 2018, and references therein). The count rates, therefore, have to be corrected, and this is done for each station by multiplying the raw count rate with a factor , based on an empirically well-established linear correlation between the logarithm of the raw count rate and the pressure P measured at the same time. The pressure is the average pressure determined at each neutron monitor station. To compare the count rates of different neutron monitors, one has to correct them to a common level, such as sea level. The quantity λ, called the attenuation length, describes the absorption of the cascade. Since the interaction of the primary particle with the atmosphere is energy-dependent, the attenuation length depends on the energy spectrum of the incoming energetic particles. The attenuation length used at each neutron monitor is determined empirically from the correlation over long time intervals (weeks and more) between the logarithm of the raw count rate and the atmospheric pressure. It hence refers to the galactic cosmic ray spectrum. The empirical values in use at high-latitude stations range from 136 to 141 hPa in the southern hemisphere and from 133 to 139 hPa in the northern hemisphere³. It is well known that the energy spectra of GLEs are steeper so that the attenuation length is shorter. Following McCracken (1962) we correct the raw count rates of the individual neutron monitors to sea level, using the method of two attenuation lengths.

The galactic background is determined either as an average level before the GLE or by a linear or parabolic fit to the count rate time histories before and after the GLEs, along the procedure of Usoskin et al. (2020). The standard pressure correction of each neutron monitor is used in this evaluation. The fractional count rate increase at sea level due to the GLE then reads

(1)

N_r designates the total count rate before pressure correction, the index G designates the galactic component, i.e. the background. The raw background count rates are recalculated from the constant, linear or parabolic fit of the galactic background determined from the data corrected for each neutron monitor. hPa is the average pressure at sea level. The attenuation length λ is defined by

(2)

McCracken (1962) gives the attenuation lengths λ_G = 135 hPa for the galactic contribution and λ_S = 98.1 hPa for the GLE part of the energetic particle spectrum. For a constant background, equation (1) corresponds to equation 6.5 of Bütikofer (2018) and to equation 11 of McCracken (1962). A constant solar attenuation length does not capture the variety of spectral slopes observed in the different GLEs. Duggal (1979) compiled values in the range (98–104) hPa for nine GLEs, and Ahluwalia and Xue (1993) evaluated 104 hPa for high-latitude neutron monitors in the 1989 Sep 29 GLE, with a linear increase with increasing vertical cutoff rigidity. In addition to variations between GLEs, the energy spectrum and therefore the attenuation length vary during a given GLE, which is not captured by simple models. We use McCracken’s value as a standard parameter that should give a more realistic result for the solar component than the use of the galactic attenuation length.

The time profiles of the count rates corrected to sea level are examined, and the neutron monitor with the strongest peak response is selected for further analysis, which is the one with the closest viewing direction to the interplanetary magnetic field. Only GLEs with peak fractional excess of at least 10% above the background level are considered. The selected neutron monitors and the GLE parameters are listed in Table 1. These neutron monitors are usually at high latitudes, where the low-energy cutoff of the energetic particle spectrum is determined by the atmospheric cutoff, rather than by the magnetic field of the Earth. Typical cutoff rigidities of the worldwide network of neutron monitors have been evaluated by Shea and Smart (2001). The atmospheric cutoff of 433 MeV (Mishev & Poluianov, 2021) corresponds to a magnetic rigidity of 1.0 GV for protons. In some GLEs, the strongest response was observed at sites with higher cutoff rigidity, such as Hobart (1.9 GV) on 1990 May 24 and Kerguelen (1.1 GV) on 1978 May 07 and 1982 Dec 07. The GLEs on 2003 Oct 28 and 29 are not listed in Table 1, because they occurred on top of a complex background that we could not determine satisfyingly.

3 A normalised GLE time profile

We first assess the common and particular aspects of GLEs by constructing a normalised time profile. The onset time t_st of each GLE is evaluated in two ways. A GLE threshold is set to be three standard deviations above the background. The first estimate of the start time is the midpoint between the latest instant where the count rate is below the GLE threshold and the following instant. Even if individual count rates may have exceeded the threshold before, they are not considered as part of the GLE, because they are isolated points. This estimate of the onset time is in column 3 of Table 1. The second estimate is based on a linear fit to the logarithm of the count rate in the early rise phase. This can in general only be computed for GLEs observed with a time resolution of 1 or 2 min. The start time is the intersection of the backward-extrapolated fit with the pre-event background. This estimate of the start time is listed in column 4, together with the statistical uncertainty of the fit. The peak time t_p (col. 5) is extracted from the 5-minute data or, when 1-minute data are used, from smoothed data over five points. This is to employ a consistent definition between the different datasets and to reduce the sensitivity of the determination to statistical count rate fluctuations. The time to maximum Δt_m = t_p – t_st is defined as the difference between peak time and start time (cols. 3 or 4 and 5). The start time derived from the extrapolation of the early rise of the curve (col. 4) is used when available. The time to maximum and the peak excess count rate at sea level are listed in columns 2 and 3 of Table 2, where the entries are ordered by increasing time to a maximum. A unified time axis of GLEs is built by referring the time of each GLE to its start time and normalising by the time to maximum: t_n = (t – t_st)/Δt_m. The count rate excess is normalised to its peak value in the 5-min data (hence: smoothed data if the 1-min cadence is used; this can lead to normalised 1-min count rates that exceed unity, e.g. on 2005 Jan 20). The normalised fractional count rate excess is plotted in Figure 1 as a function of the normalised time for all GLEs of Table 1. The top panel shows the individual GLEs, the bottom panel the median profile plus and minus the mean absolute deviation of the median. Since the steps of the normalised time are different for different GLEs, the normalised times of the individual GLEs were linearly interpolated onto a common time axis. The median GLE profile, which has a maximum near 0.92, was normalised to unity, and the mean absolute deviations were increased accordingly.

Fig. 1 Time histories of the normalised relative count rate excess of GLEs. The origin of the time axis is the start of the GLE, and the time is normalised by the time to maximum. Top: Individual GLEs, as noted in the right margin, plotted with different colours. The curve with a large excess in the decay phase is GLE 2000 Jul 14 (see text). Bottom: Median of the profiles in the top panel (plus signs, which merge into a continuous fat line before the rise and in the decay) plus and minus its mean absolute deviation (grey vertical bars). The median was in addition smoothed over 55 points (time step 0.55 units) to reduce noise, and normalised to unity. The red solid line is the best fit of a modified Weibull function to the median profile. The dashed lines show the fit plus and minus an analytically estimated statistical uncertainty (see text, Eqs. (3) and (4)).

The figure shows that the time profiles of GLEs display a large amount of similarity. The normalised GLE time profiles define an overall median curve with a decay that is longer than the normalised time to maximum in agreement with a visual inspection. Events with a symmetric time profile, like 2005 Jan 20 and 1989 Oct 22, have shorter relative decay than average but do not show up as a distinct population in the plot. An event with particularly long normalised decay or particularly high excess count rate during the decay is GLE 2000 Jul 14. It appears as a singular case in our sample. This will be further discussed in Section 6.1. The decay profile is sensitive to the shape of the background. Our analytical representation by a constant, linear, or quadratic background is an approximation evaluated by eye. Therefore part of the variation of the decay time profiles may be induced by an inappropriate estimation of the background.

The smoothed median GLE profile intersects the value 0.5 at times 0.34 (rise) and 2.79 (decay). These values were derived from linear fits to the logarithm of the smoothed median profile using ten data points around the times where the mean rise and decay profiles go through 0.5. The ratio of the decay time (peak to 50%) over the rise time (50% to peak) is 2.7. This is smaller than the ratio of 3.5 reported by Strauss et al. (2017). The time profiles are plotted until 11 times the time to maximum. The relative statistical uncertainty of the median time profile increases about linearly from a value near 0 at the time of the maximum to 100% at 10 normalised time units later.

Ji et al. (2014) and Kahler and Ling (2017) explored the possibility to represent SEP profiles at energies above 10 MeV by analytical functions. They concluded that the best-suited one was a modified Weibull function (Eq. (2) of Kahler & Ling, 2017). When expressed in normalised time units t_n and normalised to its maximum, the modified Weibull function reads

(3)

The modification refers to the fact that α < 0. While the Weibull function is usually employed as a probability distribution, the modified version is used in the present context as a convenient description of a time profile with a rapid rise and slow decay. We found that this function describes well the decay phase of the median GLE profile, but is too steep in the rise phase. The bottom panel of Figure 1 shows the median GLE profile (plus signs) and the mean absolute deviation (vertical grey error bars) together with the modified Weibull function (red solid line) inferred from a non-linear least-squares fit. Only the decay phase (t_n > 1) was used for the fit. The best representation of the data is obtained for α = −0.880. The fit function is confounded with the fat symbols used to represent the median GLE profile.

During the decay of the GLE (t_n  > 1) the relative mean absolute deviation can be fit by a straight line

(4)

with c₀ = −0.031 ± 0.005 and c₁ = 0.0983 ± 0.0007. This fit can be used in an analytical estimate of the uncertainty in the decay phase. The two dashed red lines in the bottom panel of Figure 1 show the fit C_W(t_n) plus and minus the uncertainty ΔC_W = (c₀ + c₁t_n)C_W(t_n).

4 Correlation analysis

In a second approach, we analyse the correlation between rise and decay times for all GLEs of our sample. We approximate the GLE time profiles with an exponential rise and an exponential decay, using the following fitting function:

(5)

where A_r, A_d, B_r and B_d are constants. With the rise time τ_r, the decay time τ_d and the peak time t_p, seven variables have to be adjusted during the fit. The constants B_r and B_d represent constant backgrounds but allow that constant background to be different before and after the peak time, as a way to approximate the variation of background seen in a few events. Since the time profiles fitted are background subtracted, B_r should be very close to zero. Here and in the following the expression “rise time” refers to the time constant of the exponential fit, whereas “time to maximum” is the difference between the peak time and the start time of the GLE.

A few examples of fits to the profiles are shown in Figure 2. As can be seen, the fitting procedure is efficient for both strong and weak events. The fit does not represent precisely the time profile, however, it provides a good estimation of the rise and decay times for most of the events. With our goal to develop a tool that could be used in space weather forecasting, we prioritised this simple and generic fitting procedure over a more complex fit that would need to be adapted on a case-by-case basis to more accurately represent the time profiles.

Fig. 2 Time profiles for three GLEs, on 2001 Apr 15 (left), 1989 Nov 15 (middle) and 1982 Dec 7 (right). The ordinate displays the fractional count rate excess above the galactic cosmic ray level. The observed profile is in black; the two-exponential profile fitted to the observation is shown by the red continuous line, and the scaled median profile of GLEs established in Section 3 is shown by the dashed blue line.

The fit results were inspected by eye and in four instances, were found not satisfactory. These events are the GLEs on 1971 Sep 1, 1991 Jun 15, 1989 Oct 19, 1989 Oct 24. In these cases, the number of points observed with adequate time resolution was too small to constrain the pre-event background. As a consequence, both the pre-event background level and the rise time were overestimated. For these events, we made the additional assumption of the constant B_r being zero and performed the fit with this fixed parameter. Additionally, for the GLE 1978 May 7, the rise time was found to be of 8 min, and the data have a 5-minute cadence. While the fit looks reasonable, the rise time is therefore poorly constrained, and the uncertainty associated was too important. Therefore, we do not include this event in the analysis of the correlation between the rise time and decay time of the GLEs.

Using the results of the fit for the GLEs in our sample, we looked into the possible correlation between rise and decay times. For this sample, Spearman’s rank correlation is 0.72, and Kendall’s rank correlation is 0.56, indicating a significant correlation (p-values are 2 × 10⁻⁴ and 3 × 10⁻⁴, respectively). We performed a linear regression using a total least squares fit method accounting for uncertainties on rise and decay times. The linear relationship is described by τ_d = Aτ_r + B, where A and B are constants. The regression gives A = 13.0 ± 1.8 and B = −15 ± 36 min. Figure 3 shows decay times versus rise times for the GLEs as reported in columns (4) and (5) in Table 2. The relation found from the linear regression is shown by the orange continuous line, and the uncertainty on parameters A and B is represented by the shaded area around this line. We performed the same fit on the median profile calculated in Section 3. The normalised rise time was found to be 0.38 ± 0.3 h, and the normalised decay time was 2.77 ± 0.9 h. The relation is shown by the blue dashed line.

Fig. 3 Rise and decay times of the GLEs, derived from the fit to the lightcurves (see text for details). The orange continuous line shows the result of the linear regression on this data, with the uncertainty on the fit parameter represented by the shaded area. The blue dashed line shows the result of the fit to the normalised median profile.

As was to be expected, there is a strong correlation between the time to maximum of the GLEs, i.e. the interval between the start and the peak, and the rise time derived from the exponential fit (see col. 4 in Table 2). However, the regression using the individual GLEs yields a longer decay time than the normalised median GLE profile.

5 GLE rise times and parameters of the parent activity

This section addresses the question of whether the GLE times are related to the parameters of the parent solar activity. We consider the soft X-ray (SXR) burst as a measure of overall energy release and the coronal mass ejection (CME) as an indicator of the expanding magnetic structures and possible shock waves. Solar SXR emission is continuously monitored by the Geostationary Operational Environmental Satellites (GOES) operated by NOAA. Lists of solar bursts provided by the National Centers for Environmental Information⁴ contain estimates of the start and peak times and the peak fluxes in the 0.1–0.8 nm wavelength range. For the GLEs studied here, the catalogue also lists the flare locations inferred from independent optical or EUV observations. Since the catalogued onset times often include early activity that is not unambiguously related to the eruptive event of interest, we used data witha 3-second resolution⁵ to refine the determination. Similar to the GLE onset, we determined the start of the SXR burst by the backward extrapolation to the pre-event level of a straight-line fit to the SXR time profile before the main peak. The first detection of the CME and the speed in the plane of the sky were taken from the SoHO/LASCO-based catalogue at the CDAW Data Center⁶ and the related halo CME catalogue⁷. Information on CMEs in 1989 was taken from the Solar Maximum Mission (SMM) CME catalogue at the High Altitude Observatory⁸. In two events (1989 Aug 16 and Nov 15) SMM observed a CME, but the velocity measurements in the catalogue are quoted as relating to prominence features, not to the front of the CMEs. These data are not used here.

The GLE parameters are compared with the parameters of solar eruptive activity in Table 2. Columns 6–8 are the parameters of the SXR burst. The start time inferred from the backward extrapolation is listed in column 6, with the value in the GOES catalogue within parentheses. Column 7 lists the time to maximum (i.e. difference between peak time and start time), again with the value from the GOES catalogue within parentheses. The peak flux is in column 8. Column 9 displays the heliographic latitude and longitude from the GOES catalogue. When CME observations of the event are available, they are reported in columns 10 and 11, with the first detection above the occulter (10) and the reported speed measured in the plane of the sky (11). The projected speed of the 2005 Jan 20 CME in the CDAW catalogue is 882 km s⁻². This appears very slow, given the projected height of the CME during the flare (Pohjolainen et al., 2007; Masson et al., 2009; Klein et al., 2014). The projected speed in the table is the one inferred by Gopalswamy et al. (2012). Simnett (2006) and Grechnev et al. (2008) argued for a lower speed, but this does not fundamentally change the conclusions we will draw.

Column 12 gives the three-dimensional speed V_3D inferred from the plane-of-sky speed V_pos and the CME cone model (Gopalswamy et al., 2010):

(6)

where θ is the angle between the CME axis and the plane of the sky, ω the half-width of the CME, which Gopalswamy et al. (2010) express by an empirical relationship with the propagation speed of the CME in the plane of the sky. If the CME is supposed to propagate radially outward from the flare site, the angle θ is determined by the longitude Φ and latitude Θ of the flare, and by the heliocentric latitude of the solar equator, B₀⁹:

(7)

The ordering with respect to the time to maximum of the GLE also provides some ordering of the SXR peak flux and the CME speed: the faster-rising GLEs tend to be accompanied by the faster CMEs and the weaker SXR bursts. This is also reflected by Spearman’s and Kendall’s rank correlation coefficients listed in the last two lines of Table 2, with the associated two-sided significance of their deviation from zero within parentheses.

The four fastest-rising GLEs where CME measurements are available all have CMEs with 3D speeds above 2000 km s⁻¹, whereas the five GLEs with the slowest rise are all accompanied by CMEs that are slower than 2000 km s⁻¹. But the sample is small. For the time-to-maximum analysis, it can be completed by the 2003 Oct 28 GLE, which was discarded from the previous analysis because of the uncertain background, especially during the decay phase of the GLE. Assuming a constant pre-event background we find the time to maximum in the range of 32–37 min, which is longer than the median time to maximum of our sample of 18 min. The halo CME catalogue reports speeds of 2459 (projected) and 3128 km s⁻¹ (3D cone model), respectively, which would place this GLE in the group of fast CMEs, whereas the time to maximum places it in the second half of the sample.

The soft X-ray peak fluxes tend to be higher for GLEs that rise slowly, but there are again exceptions. The rise time of the soft X-ray burst is not ordered by the GLE time to maximum. Neither is the decay time (not listed in the Table).

6 Discussion

Between 1971 and 2022 the Sun produced 51 GLEs. 25 had a fractional peak count rate excess corrected to sea level above 10%. 23 events, where the galactic cosmic ray background could be reliably identified, were analysed in the present paper. Their time profiles are found to display a remarkable degree of similarity, with a correlation between the decay time and rise time of the individual GLEs, inferred from exponential fits. A median time profile can be determined using a time axis normalised by the time to maximum.

6.1 The median GLE time profile

The existence of a well-defined average or median GLE time profile is not necessarily expected, given that times to maximum range from 6 min to 3 h in the event sample of the present work. But the normalised excess count rate during the decay of GLE 2000 Jul 14 is clearly above the other events in Figure 1. This means that the count rate is enhanced by some means. The effect may be purely statistical, with the 2000 Jul 14 GLE showing an extreme limit, or it may be due to particular physical conditions, which are usually not realised.

We evaluate the difference between the observed time profile of the GLE and the median time profile in order to identify a possible additional particle source in the decay phase. Figure 4 shows the fractional increase of the GLE at the Thule neutron monitor (black curve), the median GLE decay time profile represented by the modified Weibull function (blue-shaded surface), scaled to the observed peak fractional increase, and the difference between the two time profiles (red). The Weibull function representation has been chosen here because it permits tracking the decay over a much longer duration than the median GLE profile in this GLE, which has a particularly short start-to-peak time. The difference profile is the time profile of the hypothetical additional particle source. Particles from this source would be first detected between 10:45 and 10:50 UT, about 16 min after the start of the event (Table 1). The peak is attained between 11:10 and 11:15 UT.

Fig. 4 Time histories of the fractional count rate excess of the GLE on 2000 Jul 14 (black), its decay profile represented by the modified Weibull function (blue curve and shaded surface), and the difference (red curve).

Bieber et al. (2002) and Bombardieri et al. (2006) analysed the anisotropy of the GLE. Bieber et al. (2002) found it to drop abruptly from high values before 10:35 UT to a constant low value between 11 UT and the end of the event (see their Fig. 5). The authors conclude that the particle transport with a pitch-angle scattering in a Parker-type heliospheric magnetic field cannot explain the rapid anisotropy decrease. They propose that the anisotropy decreases due to the arrival of sunward-travelling particles that were reflected back to Earth at a magnetic mirror point near 1.3 AU. The timing of the reflected particles gives a plausible physical interpretation to our finding of an unusually high fractional excess in the decay phase since the sudden decrease of the anisotropy and the rise of the excess above the median GLE profile coincide. Bieber et al. (2002) relate the magnetic mirror to local magnetic compression by a CME that had been observed near Earth about 20 h before the GLE start.

Backward reflection of particles can, of course, happen in other GLEs, since the majority occur under disturbed solar wind conditions (e.g., Masson et al., 2012). We compared the observed normalised time profiles of all GLEs with the median profile, scaled to the times of the individual GLEs. In several events, the profiles show evidence of two decay phases, a fast initial one and a slower following one. In some the initial decay is not far from the median profile, and the slower second decay appears as an excess with a rather well-defined onset. This is especially the case on 1982 Dec 07 (right panel of Fig. 2) and 1989 Oct 22. In these cases, Cramp et al. (1997a, 1997b) found that the anisotropies were better modelled when an interplanetary obstacle beyond the Earth scattered particles back. The arrival of the first backscattered particles could be timed on 1989 Oct 22. It coincides within 5 min with the time when the individual GLE profile started to exceed the median one. The 1990 May 21 GLE has a similar feature, but we found no report of model calculations of this event in the literature. Possible hints towards reflected protons were also reported in the 2005 Jan 20 GLE (Bombardieri et al., 2008; Bieber et al., 2013). Our interpretation of the deviation of the 2000 Jul 14 GLE from the median profile is that the contribution of backscattered protons was particularly large, while it may be hidden in the uncertainty of the decay profile in other events. The fluence represented by the difference profile in Figure 4 is found to be 74% of the fluence of the GLE during its decay phase. Since Asvestari et al. (2017) show that the integrated count rate of GLEs observed by neutron monitors at sea level is about proportional to the fluence of particles above 200 MeV and especially above 800 MeV (their Fig. 4), the 74% ratio between the fluences during the decay of the 2000 Jul 14 GLE can be related to the fraction of reflected particles to the total number of particles detected. This value is not far from the fraction of 85% of the particles that Bieber et al. (2002) estimated to be reflected back to the Earth in the 2000 Jul 14 GLE.

The correlation between rise time and decay time and the existence of a well-defined median time profile are consistent with the correlation between a more restrictive definition of rise time and decay time by Moraal et al. (2015) and Strauss et al. (2017) and generalise their result. Quantitatively there is some difference, in that Strauss et al. (2017) inferred a ratio of 3.5 between the times of rise from the 50% level to the peak and decay from peak to 50%, whereas our median profile displays a ratio of 2.7. We find the same ratio of 2.7 when calculating the median or the average profile of the GLEs analysed by Strauss et al. (2017), except for the 2003 Oct 28 event. In the present study, the ratio of the decay time over rise time inferred from the correlation between the times of individual events is also higher than the ratio in the median GLE profile. This is seen by comparing the orange and blue lines in Figure 3. The least-squares correlation is pushed to longer decay times by the relatively larger role played by the events with slow decay, including the outlier on 2000 Jul 14.

6.2 How to understand the relationship between rise time and decay time?

In view of the well-established idea that the decay of GLEs is shaped by interplanetary propagation, the relationship with the rise time that we inferred through two different methods suggests that the rise is also largely determined by the interplanetary propagation, rather than the coronal acceleration, of the relativistic protons and nuclei. Moraal et al. (2015) and Strauss et al. (2017) substantiated this view by analytical evaluations and numerical simulations of diffusive transport. Their conclusion is that pitch angle diffusion in the interplanetary medium governs the time profiles in all events, and that the variation of scattering mean free paths creates substantial differences in rise times between individual GLEs.

If this interpretation is correct, no relationship will be expected between the rise time of GLEs and the parameters describing the particle acceleration. We examined the traditional indicators, namely soft X-ray emission and CME velocity. Indeed, no statistical relationship was found between the times to maximum of the GLE and the associated soft X-ray burst. The rise time of a soft X-ray burst, which reflects the heating of the plasma in the flare region, does of course not need to be physically related to the time scales of the particle acceleration processes. But there are hints that the faster-rising GLEs accompany faster CMEs as well as SXR bursts with weaker peak fluxes. The rank correlations appear significant, with probabilities of a chance coincidence of one to two percent. However, the 2003 Oct 28 GLE, which could not be used for the definition of the average GLE profile because of a strongly-varying galactic cosmic-ray background, is a counter-example because it combines a very fast CME with a rather slow GLE rise.

A correlation between the rise time or time to maximum of a GLE with the CME speed, if it were confirmed, could be interpreted in different ways. In the scenario of particle acceleration at the CME shock, the fast rise would be ascribed to the high CME speed, which presumably means a strong shock. But if the rapid rise is due to a particularly efficient particle accelerator, it is not easy to understand why under these conditions the decay of the GLE profile is fast, too. A different interpretation of a relationship between GLE rise time and CME speed could be a leaky-trap model of the GLE particles: if these particles were initially trapped in the erupting magnetic structures, for instance, the expanding flux rope of the CME, a rapid expansion would produce a fast increase of particle numbers in the magnetic trap because of the rapid increase of the mirror ratio in the flux rope, but also a fast energy loss of the confined particles through strong betatron and Fermi deceleration. Suggestive indications that GLE particles could leak out of such a trap when the expanding flux rope reconnects with neighbouring open magnetic structures were recently advocated for the small 2021 Oct 28 GLE (Klein et al., 2022), in line with scenarios studied with numerical simulations (Masson et al., 2013, 2019). Comparative studies with long-duration nuclear gamma-ray events as detected by Fermi/LAT (Share et al., 2018; Ajello et al., 2021) could help to shed light on this issue, because in the leaky-trap scenario the gamma-ray emission would come from protons dumped into the chromospheric footpoints of the confining flux rope (Mandzhavidze & Ramaty, 1992), and a similar relationship between the durations of the rise and decay might be expected.

The second significant correlation we found links more slowly-rising GLEs with stronger soft X-ray bursts. Since faster CMEs usually go along with stronger soft X-ray bursts (Salas-Matamoros & Klein, 2015), this correlation is distinctly different from the anti-correlation between CME speed and GLE rise time. Here again, the statistical relation includes outliers, for instance, the weakest and the strongest soft X-ray burst are in the third and fourth quartile of the sample ordered by the GLE time to the maximum, respectively. We mention this correlation for the record, but a physical interpretation remains to be found. Both statistical relationships with solar activity appear weak enough to argue that the interplanetary propagation probably contributes to shape even the rise of the GLE profile, but they do suggest a remaining link with the parent activity.

6.3 An application for space weather forecasting?

Irrespective of the interpretation, the dependence of the decay time on the rise time of GLEs may be of use in space weather activities. The prediction of the time profile of a GLE can, for instance, enable a radiation service for civil aviation, as it presently exists under the auspices of the ICAO¹⁰, to estimate the duration of an alert for high radiation doses due to the particle cascade generated in the Earth’s atmosphere. An estimate of the duration of an alert will be a significant additional capability from an operational point of view. It could be implemented into models that use observations of GLEs by neutron monitors as their basic observational ingredient, such as the SIGLE_RT model in France, based on the former SIGLE model developed for post-event analysis (Lantos & Fuller, 2004), or the WASAVIES model in Japan (Kubo et al., 2015).

Several authors devised schemes for the description and forecasting of the evolution of SEP events. Ji et al. (2014) and Kahler and Ling (2017) explore a description of SEP time profiles at energies above 10 MeV by simple parameterised functions and find that a modified Weibull function allows for a satisfying description, through two parameters describing the shape and duration. Kahler and Ling apply fits to three energy channels in 14 large SEP events. They find satisfying fits (as seen from plots of the logarithm of the intensity) of both the rise phase and the early decay phase. Kubo et al. (2015) model SEP events using the focussed transport equation, which describes the propagation of particles along the interplanetary magnetic field considering magnetic focussing by its laminar part and pitch-angle scattering by its turbulent component. The particles are assumed to be injected near the Sun with an inverse Gaussian time profile. The application of numerical solutions to several SEP events shows a very good description of the observed time profiles near 165 and 433 MeV. They also show that the observation of the rise and maximum phase of the SEP event can be used to determine the model parameters and predict the decay phase.

The UMASEP scheme developed at the University of Malaga (Spain) is a sophisticated set of models for the real-time forecasting of SEP events and their time evolution. The original model predicts SEP events at energies above 10 MeV (Núñez, 2011). In this model the evaluation of the correlations, with some time shift consistent with the expected travel times, between the derivatives of the soft X-ray flux and of the particle intensity measured by a space-borne detector is used to estimate whether a magnetic connection exists between the spacecraft and the Sun. A SEP event is considered to be in progress when this correlation is sufficiently high and when a SXR burst above class C7 occurred. The SXR fluence and the connectivity parameter are then used to predict the forthcoming SEP flux. The model uses a number of empirical rules, based on historic SEP events, to determine thresholds and free parameters. This model was recently adapted to the prediction of relativistic solar particle events (energies above 500 MeV, Núñez et al., 2017) at GOES and of GLEs.

The correlation established in the present work between the rise time and the decay time of a GLE, as well as the median time profile and the simple analytical description of the decay and its statistical range by a modified Weibull function, allow one to predict the decay phase of a GLE once its rise phase has been observed. In the following, we examine the accuracy of possible predictions and describe the potential operational use of these predictions.

Using the median profile of the GLEs and the linear relation between rise and decay time that we established in this paper, we calculated the predicted decay times of the GLEs analyzed in this study, based on the knowledge of their rise time. A comparison between the observed and predicted decay times is shown in Figure 5. For the 22 GLEs in our sample (recall that the 1978 May 07 GLE was not considered in the analysis of Sect. 4), the prediction done with the median profile naturally leads to the decay time being overestimated for 11 events, and underestimated for 11 events. When the prediction is done using the linear relation between the rise and decay times, the decay time is overestimated in 13 events and underestimated in 9 events. We calculated the relative difference between the predicted and observed decay times by dividing the difference between observed and predicted values by the observed value of the decay time. The distribution of the relative differences is displayed in Figure 6 as a function of measured decay time. Several conclusions arise from this display of the relative difference between predictions and observations:

• the relative difference between predicted and observed decay times is on average higher for shorter decay times;

• for the longest measured decay times, both methods of prediction led to an under-estimate of the decay time;

• both methods display relative differences between predictions and observations of similar magnitude.

Fig. 5 Comparison between measured and predicted decay times of the GLEs. Left: predictions using the median profile of GLE; right: predictions using the linear relation between rise and decay times.

Fig. 6 Relative difference between the predicted and the measured decay time (in percentage), shown as a function of measured decay time. Predictions made using the median profile are shown in orange, and predictions made using the linear regression are shown in blue. Triangles directed up are used for predictions overestimating the decay time, and triangles directed down are used for predictions underestimating the decay time.

The main uncertainties in these predictions arise from 1) the background determination and 2) the deviation of the time profile from the exponential shape, in particular for the events that seem to present contributions from two different populations of energetic protons, as the GLE 2000 Jul 14 discussed in Section 6.1. As such profiles could be the result of a strong contribution to GLEs by particles that were reflected back to the Earth, our method of prediction would underestimate the duration of the decay phase for these events. The practicability of the method has to be tested under realistic operational conditions, and a few considerations are expressed below.

The peak intensity during GLEs is reached in less than one hour in 17 events in our sample of 23 events, and often in less than 20 min. An exponential fit of the rising phase enables a first prediction of the decay time even before the peak. As soon as the peak is observed, the measured time to maximum and the measured peak count rate can be used to scale the median GLE profile to the observation, and to obtain a second prediction of the decay phase.The relations derived here used the time profiles of the monitor that registered the maximum increase of flux. Therefore, an operational service would 1) start the alert when the flux increases above a pre-set threshold in a given monitor; 2) as soon as the flux increases in time, identify the monitor in which the count rate increased the most; 3) measure the rise time with one or both methods presented in this paper; 4) make an estimation of the decay time and of the related duration of the event (possibly provide a range for the duration) and provide an estimation of the duration of the alert ahead of time. Typical alert durations would range between 2 h and 10 h after the peak of the event, as seen in Figure 5. This estimation of the duration of the alert can then be adjusted as the event proceeds through the decay phase until the end of the alert is decided based on the observations (count rate returning to a level below the pre-set threshold). The procedure will need to scan a range of neutron monitors to identify the one with the strongest peak. The real-time database NMDB (https://www.nmdb.eu) offers the necessary data.

7 Conclusion

In this study of 23 strong events that happened between 1971 and 2012, we show that relativistic solar particle events (GLEs) share similarities in their time profiles, with an evident correlation between rise and decay times. A median time profile was derived using normalised time profiles of individual GLEs. Departures from this median time profile have been found in a few events. Excesses above the median were especially found in GLEs where independent modelling analyses reported a strong contribution of back-scattered particles from ICMEs beyond 1 AU. We additionally showed that both the rise and decay times of GLEs can be described with an exponential fit of the time profiles. The rise times show some correlation with the peak soft X-ray flux and some anticorrelation with the CME speed, which suggests that the parent solar activity plays a role in shaping the GLE time profiles, in addition to interplanetary transport processes. Our findings demonstrate that a prediction of the decay time can be produced early during a GLE. This information could be used by space weather operation services that issue alerts related to high radiation doses for civil aviation to provide an estimate of their duration. The correlation between rise and decay times also provides additional information regarding the physical processes controlling the GLE time profiles and can give useful input to future studies and simulations of GLEs.

Acknowledgments

The authors acknowledge the providers of the neutron monitor data used in this work, namely the PIs of the worldwide neutron monitor network, the NMDB database hosted by the University of Kiel and the GLE database at the University of Oulu. The French Polar Institute IPEV is acknowledged for its operation of the neutron monitors at Kerguelen Island and Terre Adélie, and the personnel at these stations for their dedicated work under difficult conditions. Sophie Musset is supported by the European Space Agency Research Fellowship programme. This study received financial support from the French space agency CNES. The authors are grateful to many colleagues, especially R. Buetikofer, S. Dalla, P. Démoulin, E. Flueckiger, A. Papaioannou, and R.D. Strauss, as well as participants at the NMDB 2022 workshop, for helpful discussions, and to the International Space Science Institute (ISSI) Bern for a number of motivating team meetings. We thank the referees for their careful reading and constructive comments, and J. Watermann for his help with technical issues of the manuscript. The editor thanks Sergey Koldobskiy and an anonymous reviewer for their assistance in evaluating this paper.

References

All Tables

All Figures

Fig. 1 Time histories of the normalised relative count rate excess of GLEs. The origin of the time axis is the start of the GLE, and the time is normalised by the time to maximum. Top: Individual GLEs, as noted in the right margin, plotted with different colours. The curve with a large excess in the decay phase is GLE 2000 Jul 14 (see text). Bottom: Median of the profiles in the top panel (plus signs, which merge into a continuous fat line before the rise and in the decay) plus and minus its mean absolute deviation (grey vertical bars). The median was in addition smoothed over 55 points (time step 0.55 units) to reduce noise, and normalised to unity. The red solid line is the best fit of a modified Weibull function to the median profile. The dashed lines show the fit plus and minus an analytically estimated statistical uncertainty (see text, Eqs. (3) and (4)).

In the text

Fig. 2 Time profiles for three GLEs, on 2001 Apr 15 (left), 1989 Nov 15 (middle) and 1982 Dec 7 (right). The ordinate displays the fractional count rate excess above the galactic cosmic ray level. The observed profile is in black; the two-exponential profile fitted to the observation is shown by the red continuous line, and the scaled median profile of GLEs established in Section 3 is shown by the dashed blue line.

In the text

	Fig. 3 Rise and decay times of the GLEs, derived from the fit to the lightcurves (see text for details). The orange continuous line shows the result of the linear regression on this data, with the uncertainty on the fit parameter represented by the shaded area. The blue dashed line shows the result of the fit to the normalised median profile.
In the text

	Fig. 4 Time histories of the fractional count rate excess of the GLE on 2000 Jul 14 (black), its decay profile represented by the modified Weibull function (blue curve and shaded surface), and the difference (red curve).
In the text

	Fig. 5 Comparison between measured and predicted decay times of the GLEs. Left: predictions using the median profile of GLE; right: predictions using the linear relation between rise and decay times.
In the text

Fig. 6 Relative difference between the predicted and the measured decay time (in percentage), shown as a function of measured decay time. Predictions made using the median profile are shown in orange, and predictions made using the linear regression are shown in blue. Triangles directed up are used for predictions overestimating the decay time, and triangles directed down are used for predictions underestimating the decay time.

In the text