On the Pulse Shape of Ground Level Enhancements

TL;DR

This study analyzes the pulse shape of cosmic ray ground level enhancements, revealing a linear relation between decay and rise times and emphasizing the importance of interplanetary transport effects in shaping GLE temporal profiles.

Contribution

It introduces a model linking GLE pulse shapes to interplanetary scattering, highlighting how transport effects can dominate observed profiles and challenge traditional event classification.

Findings

A linear relation $ au_d \\approx 3.5 \\tau_r$ is observed between decay and rise times.

Interplanetary scattering can fully determine GLE pulse shapes under certain conditions.

Transport effects can obscure initial acceleration signatures, affecting event classification.

Abstract

We study the temporal intensity profile, or pulse shape, of cosmic ray ground level enhancements (GLEs) by calculating the rise $(\tau_\mathrm{r})$ and decay $(\tau_\mathrm{d})$ times for a small subset of all available events. Although these quantities show very large inter-event variability, a linear dependence of $\tau_\mathrm{d} \approx 3.5 \tau_\mathrm{r}$ is found. We interpret these observational findings in terms of an interplanetary transport model, thereby including the effects of scattering (in pitch-angle) as these particles propagate from (near) the Sun to Earth. It is shown that such a model can account for the observed trends in the pulse shape, illustrating that interplanetary transport must be taken into account when studying GLE events, especially their temporal profiles. Furthermore, depending on the model parameters, the pulse shape of GLEs may be determined entirely by interplanetary scattering, obscuring all information regarding the initial acceleration process, and hence making a classification between impulsive and gradual events, as is traditionally done, superfluous.

Tables

Table 1. Table 1. : A summary of the GLE events selected for this study. The approximate solar longitude of the transient event and the maximum increase is taken from McCracken et al. ( 2012 ) . \ilabel \thearticle Tab:events˙table

GLE no.	Sol. long.111Ideal magnetic connection between the source of the SEPs and Earth would correspond to a source located at $\approx {60}^{\circ }$ W in terms of solar longitude.	Max. increase 222The maximum increase of the omni-directional counting rate.	Rise time	Decay time	NM stations 333ALTR – Alert, APTY – Apatity, BRBY – Barentzburg, CALG – Calgary, CAPS – Cape Schmidt, DPRV – Deep River, FSMT – Fort Smith, GSBY – Goose Bay, INVK – Inuvik, MCMD – McMurdo, MRNY – Mirny, MWSN – Mawson, NRLK – Norilsk, OTWA – Ottawa, PWNK – Peawanuk, TERA – Terra Adelie, THUL – Thule, TXBY – Tixie Bay
	(∘ W)	(%)	(min)	(min)	( GV)
30	40	19	25.2	75	ALRT, APTY, DPRV, GSBY, INVK,
					MCMD, MRNY, OULU, TXBY, SANAE
31	68	28	5.4	19.8	ALRT, APTY, DPRV, GSBY, INVK,
					MCMD, OULU, TXBY, SANAE
38	86	17	16.2	84.6	ALRT, APTY, DPRV, GSBY, INVK,
					MSWN, OULU, TERA, THUL, TXBY,
					SANAE, CALG
39	130	30	6	16.2	ALRT, APTY, DPRV, GSBY, INVK,
					MCMD, MWSN, OULU, TERA, THUL,
					TXBY, SANAE, CALG
42	105	246	78	158.4	APTY, CALG, CAPS, DPRV, INVK,
					MCMD, MWSN, MRNY, OTWA, OULU,
					TERA, THUL, TXBY, SANAE
43	-9	41	69	443.4	CALG, CAPS, DPRV, GSBY, INVK,
					MCMD, MWSN, OTWA, OULU, TERA,
					THUL, SANAE
45	57	91	75.6	241.8	APTY, CALG, CAPS, DPRV, GSBY,
					INVK, MCMD, MWSN, OTWA, OULU,
					TERA, THUL, SANAE
48	76	11	6	75	CALG, CAPS, DPRV, GSBY, INVK,
					MWSN, MCMD, OTWA, OULU, THUL,
					TXBY, SANAE
52	70	21	25.2	80.4	APTY, CALG, DPRV, GSBY, INVK,
					MCMD, OULU, TERA, THUL, TXBY,
					SANAE
59	7	31	22.8	96	APTY, CALG, GSBY, INVK, MCMD,
					MSWN, OULU, THUL, TXBY, SANAE
60	85	80	25.8	73.2	APTY, CALG, CAPS, FSMT, MCMD,
					NAIN, NRLK, OULU, PWNK, TERA,
					THUL, TXBY, SANAE
65	-2	12	22.8	40.2	APTY, BRBY, CAPS, FSMT, INVK,
					MCMD, NAIN, NRLK, OULU, PWNK,
					TERA, THUL, TXBY, SANAE
69	58	430	3	16.2	APTY, BRBY, CALG, CAPS, FSMT,
					INVK, MCMD, MWSN, NAIN, NRLK,
					OULU, TERA, THUL, TXBY, SANAE
70	24	27	33.6	70.2	APTY, BRBY, CAPS, FSMT, INVK,
					MCMD, MWSN, NAIN, NRLK, OULU,
					PWNK, TERA, THUL, TXBY, SANAE

Equations

$$\text{ilabel}\text{thearticle}Eq:TPE\frac{\partial f(z,\mu ,t)}{\partial t}=-\mu v\frac{\partial f}{\partial z}-\frac{1-{\mu }^2}{2L}v\frac{\partial f}{\partial \mu }+\frac{\partial }{\partial \mu }\left(D_{\mu \mu }(z,\mu )\frac{\partial f}{\partial \mu }\right),$$

$$L^{-1}(z)=\frac{1}{B(z)}\frac{\partial B(z)}{\partial z}$$

$${\mathbf{e}}_z\equiv \frac{\text{mathbfit}B}{B}=\cos \Psi {\mathbf{e}}_r-\sin \Psi {\mathbf{e}}_{\varphi }$$

$$\text{mathbfit}B=B_0{\left[\frac{r_0}{r}\right]}^2\left({\mathbf{e}}_r-\tan \Psi {\mathbf{e}}_{\varphi }\right),$$

$$\tan \Psi =-\frac{B_{\varphi }}{B_r}=\frac{\Omega r\sin \theta }{V},$$

$$\text{ilabel}\text{thearticle}Eq:de{f}_{z}z(r)={\int }_0^r\mathrm{d}z(r)={\int }_0^r\mathrm{d}r\sqrt{1+{\tan }^2\Psi }$$

$$\mathrm{d}\varphi =-\frac{\Omega }{V}\mathrm{d}r.$$

$$D_{\mu \mu }(z,\mu )=D(z)\left(1-{\mu }^2\right)\left\{{\left|\mu \right|}^{q-1}+H\right\}$$

$$\text{ilabel}\text{thearticle}Eq:de{f}_{l}ambda{\lambda }_{||}(z)=\frac{3v}{8}{\int }_{-1}^{+1}\frac{{\left(1-{\mu }^2\right)}^2}{D_{\mu \mu }(z,\mu )}\mathrm{d}\mu ,$$

$${\lambda }_{rr}={\lambda }_{||}{\cos }^2\Psi .$$

$$\text{ilabel}\text{thearticle}Eq:parametrized{\lambda }_{rr}={\lambda }_0{\left[\frac{r}{r_0}\right]}^{\alpha },$$

$$F(z,t)=\frac{1}{2}{\int }_{-1}^{+1}f(z,\mu ,t)\mathrm{d}\mu ,$$

$$\frac{1}{3}\frac{\partial f}{\partial t^{\prime }}+\mu v\frac{\partial f}{\partial z}$$

$$\frac{1}{3}\frac{\partial f}{\partial t^{\prime }}+\frac{1-{\mu }^2}{2L}v\frac{\partial f}{\partial \mu }$$

$$\frac{1}{3}\frac{\partial f}{\partial t^{\prime }}$$

$$A=3\frac{{\int }_{-1}^1\mu f\mathrm{d}\mu }{{\int }_{-1}^{1}f\mathrm{d}\mu }.$$

$$\frac{\partial f}{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{\kappa }_{rr}\frac{\partial f}{\partial r}\right).$$

$${\kappa }_{rr}={\kappa }_0{\left[\frac{r}{r_0}\right]}^{\alpha },$$

$$f\propto t^{3/(\alpha -2)}\exp \left[\frac{-r^2{\left(r_0/r\right)}^{\alpha }}{(2-\alpha {)}^2{\kappa }_{0}t}\right].\text{ilabel}\text{thearticle}Eq:genera{l}_{s}olution$$

$$t_{\max }=\frac{r^2{\left(r_0/r\right)}^{\alpha }}{3(2-\alpha ){\kappa }_0},$$

$$f_{\max }\propto {\left(\exp \left(1\right)\cdot t_{\max }\right)}^{3/(\alpha -2)},$$

$$\frac{f}{f_{\max }}={\left[\frac{t}{t_{\max }}\exp \left(\frac{t_{\max }}{t}-1\right)\right]}^{3/(\alpha -2)}.\text{ilabel}\text{thearticle}Eq:rewritten-solution$$

$$\frac{t_{1/2}}{t_{\max }}=2^{(2-\alpha )/3}\exp \left(1-\frac{t_{\max }}{t_{1/2}}\right)\text{ilabel}\text{thearticle}Eq:transcendenta{l}_{f}unction$$

$$x$$

$$y(x,\alpha )$$

$$\alpha =-2$$

$$\alpha =-1$$

$$\alpha =0$$

$$\alpha =1$$

Full text

{opening}

1 Introduction

Ground level enhancements (GLEs) are sudden increases in the cosmic ray (CR) intensity as observed at Earth’s surface, in recent times mostly through neutron monitors (NMs). Since the first observation of a GLE on 28 February 1942, such events have been observed seventy-one times (see the review by McCracken et al., 2012). A database of all GLE events is available at ftp://cr0.izmiran.rssi.ru/COSRAY!/FTP˙GLE/ (hosted by E. Eroshenko), https://gle.oulu.fi/, as well as http://usuarios.geofisica.unam.mx/GLE˙Data˙Base/ with details of the latter described by Moraal and Caballero-Lopez (2014). It is generally believed that transient solar eruptive episodes, such as coronal mass ejections (CMEs) and solar flares, are responsible for producing the enhanced CR flux detected during a GLE event, and hence these particles can be classified as solar energetic particles (SEPs). If the primary SEPs reach Earth with kinetic energies $\geq$ 500 MeV per unit charge, or rigidity $\geq$ 1 GV, they can produce atmospheric showers of secondary particles that can reach ground level and be registered by NMs. See also the work of, amongst others, Shea and Smart (1994), Miroshnichenko, De Koning, and Perez-Enriquez (2000), Oh et al. (2012), and Caballero-Lopez and Moraal (2016).

Despite the large amount of SEP observations by both ground and space based CR detectors over the years, the relationship between flares and CMEs, and their role in accelerating particles to relativistic energies during major solar events remains an exigent scientific challenge. Most studies, e.g. Reames (1999), classify SEPs into two distinct classes, impulsive and gradual events. Particles from these events are believed to be accelerated by different processes in different transient structures; impulsive SEPs are associated with magnetic reconnection in solar flares, while gradual SEPs are thought to be accelerated mainly by diffusive shock acceleration in the shock fronts of large CMEs. Having an accurate classification of SEP events may therefore give insight into where and how these particles were accelerated (e.g. Reames, 2013). The usual classification, based on the shape of the observed intensity time profile, the so-called pulse shape, is furthermore motivated by differences in the composition and charge state of the SEPs (e.g. Mewaldt et al., 2012), for example, impulsive events have higher ionization states. For ground-based observations of GLE, however, both the composition and change state information is lost and the GLEs are classified exclusively on their pulse shape. In Figure \iref\thearticlefig_compare_GLEs, as an example, we show the pulse shape of GLE 45 (blue curves) and GLE 69 (red curves) for all available NM stations with a cut-off rigidity below 1 GV. It is clear that GLE 45 was observed to be more gradual than the more impulsive GLE 69, but it is impossible to extract information solely from this temporal profile about the acceleration mechanisms responsible for each particle event without taking interplanetary transport into account. Once the SEPs are accelerated, presumably near the Sun, they still have to propagate (mainly) along the turbulent interplanetary magnetic field to reach Earth. Depending on the level and nature of the underlying turbulent fluctuations and how efficiently they can disrupt (i.e. scatter) the particle trajectories, the SEP temporal profile reaching Earth may be significantly different from that near the acceleration region.

In this study we will show that the interplanetary transport conditions may alter the pulse shape in such a way as to remove any source information by the time it reaches Earth, possibly making the classification of GLEs into different classes superfluous.

2 Characterising the Pulse Shape

\ilabel\thearticle

Sec:charactering_the

The data base, which contains all available NMs observation, of all 71 GLEs since 1942, has been decsribed by McCracken et al. (2012). In this study, we select only 14 GLEs, based on the following 5 criteria: i) in order to remove any energy effects, we only use NMs with cut-off rigidities below 1 GV, as all these NMs are sensitive to the same particle energy range. As we will show later in this paper, particle propagation is very dependent on the transport coefficients, which are believed to be energy dependent, making a direct comparison of the pulse shape of different energy particles (i.e. NMs with different cut-off rigidities) impossible. ii) We select only GLEs that were observed simultaneously by more that 10 NM stations. It is well-known that NMs are sensitive to the arrival direction of particles due to deflection by the Earth’s magnetosphere (e.g. Shea and Smart, 1990). For an anisotropic event, such as a GLE where particles are predominantly streaming towards Earth, different stations at the same cut-off rigidity may respond differently to the same SEP pulse. Note, for example, the difference in peak intensity of GLE 69 (shown in Figure \iref\thearticlefig_compare_GLEs) as observed by different NMs due to this geographical effect of essentially viewing particles with different pitch-angles. By using more NM stations, we try, as far as possible, to remove this observational bias. In Figure \iref\thearticlefig_viewing_angle we show the calculated asymptotic viewing directions for a number of NM stations used in this study. Note the nearly uniform longitudinal coverage near the equatorial regions. iii) The selected GLEs have amplitudes more than 10% and iv) consists of both small and large events v) showing both gradual and impulsive characteristics. Table \iref\thearticleTab:events_table summarizes all the events that were selected and used in this study. We consider this selection as a good subset of all available 71 events. Lastly, before quantifying the pulse shape, we average the time profiles of the separate NMs to remove the observational pitch-angle bias and refer to this as the omni-directional time profile. This is also done in order to facilitate a later comparison to modelled solutions.

In Figure \iref\thearticlefig_define_times, we demonstrate how the pulse shape is characterized in this study: for each event, the time profile of the omni-directional intensity is calculated and the time of maximum intensity, $t_{\max}$, recorded. We then search for times when the intensity is half of its maximum value and note the corresponding times $t^{a}_{1/2}$ and $t^{b}_{1/2}$. The rise time is now defined as $\tau_{\mathrm{r}}:=t_{\max}-t^{a}_{1/2}$ and the decay time as $\tau_{\mathrm{d}}:=t^{b}_{1/2}-t_{\max}$. Here we deviate from Moraal, McCracken, and Caballero-Lopez (2016), and usual SEP studies such as Dresing et al. (2014), by not calculating the so-called onset time (defined as the time when the observed intensity breaches the background level) and defining the rise time in terms of this quantity. Observationally, the onset-times are biased towards larger events, or events with fast rising intensities (i.e. for impulsive events), that breach the background quicker (e.g. Xie et al., 2016), while, from a modelling point-of-view any chosen background level in a model is somewhat arbitrary (e.g. Wang and Qin, 2015). Our re-definition of the rise time removes this bias, while sufficiently quantifying the temporal slope of the increasing phase of an SEP event.

Figure \iref\thearticlefig_data_with_fit shows our calculations of the decay and rise times for the 14 selected GLE events as the data points. Although there is still a significant spread in these observed values, a linear trend is discernible. By fitting a linear curve (through the origin) to the data, we find $\tau_{\mathrm{d}}=(3.5\pm 0.8)\tau_{\mathrm{r}}$, where the error is simply taken to be a standard deviation. The regression curve, with its error, is shown on the figure as the red line and shaded region respectively. All the data-points are fairly well described by the linear fit, with the possible exception of that corresponding to GLE 48 (the circled data-point) which has an exceedingly long decay time. A closer inspection of GLE 48’s temporal profile shows a second peak (which might be a second event) some $\approx 4$ hours after the first. It is possible that this second increase influences our calculation of the decay time which was calculated only for the first event. Simply stated, our analysis therefore indicates that the decay time is roughly 3 to 4 times longer than the rise time for a broad range of time-scales containing both impulsive and gradual events. This value differs from that reported by Moraal, McCracken, and Caballero-Lopez (2016), i.e. $\tau_{\mathrm{d}}\approx 2\tau_{\mathrm{r}}$, although they used a different definition of $\tau_{\mathrm{r}}$, but similarly also found a linear relationship between $\tau_{\mathrm{r}}$ and $\tau_{\mathrm{d}}$. In the next sections, we interpret this seemingly universal dependence in terms of a particle transport model.

3 Comparison with Theory

3.1 The Transport Model

The propagation of a gyro-tropic distribution of SEPs can be described by the Roelof (1969) transport equation (TPE),

[TABLE]

which is valid for motion along a magnetic field directed along $\mathbf{e}_{z}$. In this expression, $\mu$ is the cosine of the particle’s pitch angle, so that the parallel speed becomes $v_{||}=v\mu$, and focussing due to the diverging Parker (1958) magnetic field is given in terms of the focussing length

[TABLE]

and pitch-angle diffusion is included via the pitch-angle diffusion coefficient $D_{\mu\mu}$. Although some analytical solutions of Equation (\iref\thearticleEq:TPE) do exist (see e.g.Effenberger and Litvinenko (2014)), these are notoriously inaccurate during the rising phase of a SEP event, and as such, we opt to solve the TPE numerically as outlined in Appendix \iref\thearticleSec:appe_A.

Usually, all transport quantities are specified, not in terms of $z$, but rather in terms of spherical coordinates, and as such, we need to convert between these two coordinate systems, $z\mapsto(r,\phi)$. The unit vector $\mathbf{e}_{z}$ is defined to point along the magnetic field, ${\mathbfit B}$, so that

[TABLE]

in terms of spherical coordinates with $\Psi$ being the magnetic field spiral angle, i.e. the angle between the magnetic field and the radial direction. For a Parker (1958) magnetic field,

[TABLE]

where $B_{0}$ is some reference values at $r_{0}$, we have

[TABLE]

which, at Earth in the equatorial plane $(\theta=\pi/2)$, gives $\Psi\approx 45^{\circ}$ when using observed values of the solar wind speed $V\approx 400$ km.s*-1* and the solar rotation rate $\Omega\approx 25$ days. As the magnitude of $z$ is now defined as the length along the spiral field, we may calculate

[TABLE]

where we have used

[TABLE]

For ease, Equation (\iref\thearticleEq:def_z) is integrated numerically in the present model. Due to the curvature of the magnetic field, we find $z(r=1\mathrm{AU})\approx 1.2\mathrm{AU}$ and $z(r=3\mathrm{AU})\approx 5.9\mathrm{AU}$.

Following Dröge et al. (2010), we parametrize the pitch-angle diffusion coefficient as

[TABLE]

where $q=5/3$ is the spectral index of the inertial range of the turbulence power spectrum and $H=0.05$ is a parameter to account for non-linear or dynamical processes which scatter particles though $\mu=0$ (see e.g. Dröge et al., 2010). Here $D(z)$ is a function which we will use to obtain the required magnitude of $D_{\mu\mu}$. In order to compare our results to analytical approximations, and to get a more intuitive feeling for the level of particle scattering, we specify the parallel mean free path $\lambda_{||}$ in the model rather than $D_{\mu\mu}$. These two quantities are related through (Hasselmann and Wibberenz, 1968)

[TABLE]

while $\lambda_{||}$ is also related to the effective radial mean free path $\lambda_{rr}$ (in the absence of perpendicular diffusion) through

[TABLE]

We parametrize $\lambda_{rr}$ as

[TABLE]

where $\lambda_{0}$ is a reference value at $r_{0}=1$ AU. Note that the same parametrization is used in Appendix \iref\thearticleSe:appendix_2, keeping in mind that $\kappa_{rr}=v/3\lambda_{rr}$, in order to compare with the analytical approximations. By specifying $\lambda_{0}$ and $\alpha$, we can then follow Equations (\iref\thearticleEq:def_lambda) - (\iref\thearticleEq:parametrized) in reverse to estimate $D(z)$ which is the input to the model. The magnitude of $\lambda_{0}$ can already give us an indication as to the behaviour of SEP intensities: at Earth, $L\approx 1$ AU, and the length of the field line connecting Earth to the Sun is $z\approx 1.2$ AU. Therefore, if $\lambda_{0}\ll L,z\approx 1$ AU, we will have very efficient particle scattering and the resulting distribution will be nearly isotropic. If not, i.e. $\lambda_{0}\approx 1$ AU, or larger, the resulting distributions will be highly anisotropic.

After solving Equation (\iref\thearticleEq:TPE) to obtain the distribution function $f(z,\mu,t)$ we can calculate the omni-directional intensity by simply averaging over pitch-angle:

[TABLE]

and compare the modelled time profile of $F(z=1.2\mathrm{AU},t)$ to that derived from GLE observations. See also the recent modelling, using a similar approach, by Bütikofer et al. (2016).

3.2 General Modelling Results

To be compatible with the GLE observations, we solve Equation (\iref\thearticleEq:TPE), for 2 GV protons, by specifying isotropic impulsive injection at the inner boundary of the model, $f(z=0.05\mathrm{AU},\mu,t)=\delta(t)$. Examples of the resulting time profiles of $F$, at Earth, are shown in the left panel of Figure \iref\thearticlefig_model_characteristics for a fixed value of $\alpha=0$ and different values of $\lambda_{0}$. As expected, the level of particle scattering (quantified through the magnitude of $\lambda_{0}$) has a large influence on the temporal profiles. For large values of $\lambda_{0}$ (weak scattering), we note a very fast rise in intensities, followed by a sharp decrease. For smaller values of $\lambda_{0}$ (stronger scattering), the intensities rise more gradually. We are therefore able to control the temporal profile of the distribution at Earth by adjusting the amount of scattering that the particles experience. Hence, the injection profile (here, simply a delta-function) is not the only controlling factor of the pulse shape of GLEs; interplanetary transport conditions can also have a large influence. In the right panel of Figure \iref\thearticlefig_model_characteristics, we keep $\alpha=0$ and $\lambda_{0}=0.05$ AU constant and show the distribution function at different pitch-angles. Of course, the particles that stream along the magnetic field, i.e. these with $\mu=1$, reach Earth first, and their profile shows a very quick rise (i.e. a short rise time). Particles mostly gyrating perpendicular to the field ($\mu=0$), take longer to reach Earth but also have a longer rise time. This is an additional motivation for us to use only the omni-directional intensity to quantify the temporal profile (for both the model and the observations) of the GLE events, thereby removing the dependence of the rise time on pitch-angle, or, in the observational sense, the dependence of the rise time on the geographical position of the NMs.

Following our standard definition, we calculate the rise time from the modelled solutions, and present these in Figure \iref\thearticlefig_rise_vs_lambda. As alluded to in Figure \iref\thearticlefig_model_characteristics, larger values of $\lambda_{0}$ lead to shorter rise times. What is however very interesting, is the linear dependence of $\tau_{\mathrm{r}}$ for small values of $\lambda_{0}$; when $\lambda_{0}$ is less than $\approx 0.1$ AU, $\tau_{\mathrm{r}}$ decreases linearly with time, with the slope of the linearity dependent on the value of $\alpha$. These small values of $\lambda_{0}$ must result in fairly isotropic distributions and, as we will show in the next section, this linear dependence is characteristic of the isotropic scenario where pitch-angle scattering is sufficiently strong. As shown in the next section, a similar linear dependence is observed for $\tau_{\mathrm{d}}$ in this limit.

3.3 Comparison between Theory and Observations

We show the modelled relationship between $\tau_{\mathrm{r}}$ and $\tau_{\mathrm{d}}$ in Figure \iref\thearticlefig_data_with_model, as the symbols, for different choices of $\alpha$. As already discussed, by decreasing $\lambda_{0}$ in the model, both $\tau_{\mathrm{r}}$ and $\tau_{\mathrm{d}}$ increase, and in the limit of $\lambda_{0}\ll 1$ AU, these increases are linear. The same behaviour is illustrated in Figure \iref\thearticlefig_data_with_model where, for $\tau_{\mathrm{r}}>10$ min, a linear relationship is evident between $\tau_{\mathrm{r}}$ and $\tau_{\mathrm{d}}$ with the slope dependent on the value of $\alpha$. This linear dependence is indicative of isotropy, and this idea can be tested by comparing the modelled solution to the analytical approximations discussed in Appendix \iref\thearticleSe:appendix_2, which are derived for an isotropic particle distribution, and shown on Figure \iref\thearticlefig_data_with_model as the different dashed lines. Indeed, at larger values of $\tau_{\mathrm{r}}$ both the model and the approximations follows the same linear trend. Our results thus indicate that for GLEs with $\tau_{\mathrm{r}}>10$ min, we may approximate the particle distributions to be isotropic and use the approximations derived in Appendix \iref\thearticleSe:appendix_2. This linear trend of course breaks down below $\tau_{\mathrm{r}}<10$ min, as these events are highly anisotropic with $\lambda_{0}\approx 1$ AU. No clear dependence is noted for these events.

Also shown on Figure \iref\thearticlefig_data_with_model as the red band, is the linear dependence of $\tau_{\mathrm{r}}$ and $\tau_{\mathrm{d}}$ as derived form GLE observations in Section \iref\thearticleSec:charactering_the. We can now explain this linear dependence in terms of particle transport in the isotropic limit. Moreover, a comparison between modelled and observed results indicates the value $\alpha\approx 0$ ($\lambda_{rr}$ independent of radial distance) is consistent with the observations; the possible implications thereof are discussed in the next section.

4 Summary and Conclusions

By calculating the rise and decay times for a small subset of all GLE events, we find a linear relationship between these two quantities, conveniently summarized as $\tau_{\mathrm{d}}\approx 3.5\tau_{\mathrm{r}}$, i.e. the decay phase lasts, on average, 3 to 4 times longer than the rising phase of the event. This relationship seems to hold for a large range of $\tau_{\mathrm{d}}$ and $\tau_{\mathrm{r}}$ values (almost two orders of magnitude), suggesting that GLEs do not fall into two distinct classes of either impulsive or gradual events, but follow a continuous distribution of impulsive-like or gradual-like events. It is difficult to imagine how different acceleration mechanisms can lead to such a “universal” linear relationship, and hence we interpret this result as being indicative of interplanetary transport. We show that interplanetary scattering can significantly affect the temporal profile of GLEs as observed at Earth, in a sense obscuring the initial acceleration profile. In the limit of very effective particle scattering, i.e. when the resulting distribution is nearly isotropic, both the numerical and analytical solutions reproduce the observed linear trend. We thus conclude that interplanetary transport may have an extremely large effect on the observed pulse shape of GLEs and should not be completely ignored as is mostly done.

There are very large inter-event variations in the values of $\tau_{\mathrm{d}}$ and $\tau_{\mathrm{r}}$, although the values follow the same general trend. These large variations are most likely related to changes in the level of interplanetary scattering between events (see e.g. Dröge (2000)) and can be accounted for in our model by varying, for example, $\lambda_{0}$. A comparison between the observations and especially the analytical approximations, suggests, however, that a constant radial dependence of $\lambda_{rr}\propto r^{0}$ (i.e. $\lambda_{rr}$ being independent of radial position) best reproduces the observations for all events. Perhaps rather fortuitously, the same assumption is frequently used in many transport models (see the discussion by Dröge (2000)), giving us some confidence in our modelling approach and conclusions.

In future, we plan to extend our analyses to include more events, as our 14 selected events only represent $\approx 20\%$ of the total recorded GLEs. However, based on the findings of Moraal, McCracken, and Caballero-Lopez (2015), who included more events, the qualitative conclusions presented here are not expected to change much. Of particular interest will also be to perform an equivalent study for SEP events that are observed in-situ by spacecraft and that can be classified more rigorously into impulsive (flare accelerated) and gradual (CME accelerated) events.

Appendix A Notes on the Numerical Model

\ilabel\thearticle

Sec:appe_A

Following Strauss and Fichtner (2015), we integrate Equation (\iref\thearticleEq:TPE) numerically by first transforming it into a set of three, 1D convection and diffusion equations, using a first-order operator splitting method, to obtain

[TABLE]

where $\mathrm{d}t^{\prime}=\mathrm{d}t/3$. The diffusion equation (Equation (\iref\thearticleEq:diff_equation)) is solved by an explicit Euler method, while the convection/advection equations (Equations (\iref\thearticleEq:adve_1) and (\iref\thearticleEq:adve_2)) are solved by a flux limiter corrected upwind scheme. For all equations, flux conserving boundary conditions are applied (see the discussion by Strauss and Fichtner (2015)). To validate this modelling approach, we compare our results to those of Dröge et al. (2010) in this section. These simulations, using exactly the same transport parameters as given by Dröge et al. (2010), are for 4 MeV protons with the results shown as a function of time at Earth’s position. In Figure \iref\thearticlefiga1 our model results are labelled as the finite-difference (FD) solutions, while the Dröge et al. (2010) results, computed by making use of stochastic differential equations (SDEs) are labelled accordingly. The left panel compares the omni-directional intensity and the right the first-order anisotropy, defined as

[TABLE]

The excellent agreement between the two independent models validates our modelling approach while giving us confidence is our calculated temporal profiles.

Appendix B (Semi-) Analytical Approximations

\ilabel\thearticle

Se:appendix_2

If the particle distribution can be approximated to be nearly-isotropic, we can follow Moraal, McCracken, and Caballero-Lopez (2016) and obtain $f$ by solving the spherical symmetric diffusion equation

[TABLE]

Assuming that the effective radial diffusion coefficient, $\kappa_{rr}$, can be parametrized as

[TABLE]

the solution of $f$ for impulsive injection at $r=0$ and $t=0$ is (Duggal, 1979)

[TABLE]

The time of peak, or maximum, intensity is evaluated as

[TABLE]

with a corresponding peak intensity of

[TABLE]

so that Equation (\iref\thearticleEq:general_solution) may be rewritten, in terms of these quantities, as

[TABLE]

This analytical approximation is shown, for a variety of different parameters, and discussed in more detail in Moraal, McCracken, and Caballero-Lopez (2015, 2016). Here, however, we are only interested in times where the distribution obtains a half of its maximum value, i.e. we are interested in finding $t^{a}_{1/2}$ and $t^{b}_{1/2}$ as defined in Figure \iref\thearticlefig_define_times. By setting $2f=f_{\max}$ in Equation (\iref\thearticleEq:rewritten-solution), we obtain the required transcendental equation

[TABLE]

that must be solved numerically to obtain the two values of $t^{a}_{1/2}<t_{\max}$ and $t^{b}_{1/2}>t_{\max}$. Note also that this expression is independent of the magnitude of the diffusion coefficient and only on its radial dependence through $\alpha$. Figure \iref\thearticlefiga2 illustrates how to find the roots (values of $x$ where $y=0$) of Equation (\iref\thearticleEq:transcendental_function) for different values of $\alpha$ by setting

[TABLE]

For example, $\alpha=0$ gives $t^{a}_{1/2}\approx 0.44t_{\max}$ and $t^{b}_{1/2}\approx 3.13t_{\max}$, which, substituted into our definitions of the rise and decay times (see again Section \iref\thearticleSec:charactering_the) leads to $\tau_{\mathrm{d}}\approx 3.8\tau_{\mathrm{{r}}}$. Our results are summarized as

[TABLE]

with the special case of $\tau_{\mathrm{{d}}}=\tau_{\mathrm{{r}}}$ for $\alpha=2$ describing a free-streaming scenario. Note that, for all values of $\alpha$, we obtain a linear relationship between $\tau_{\mathrm{{d}}}$ and $\tau_{\mathrm{{r}}}$. These resulting estimates, referred to as the isotropic solutions, are included in Figure \iref\thearticlefig_data_with_model as the dashed lines.

Acknowledgements

This work is based on the research supported in part by the National Research Foundation (NRF) of South Africa (grant no. 106049). Opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the NRF. OO acknowledges the support of the post-doctoral programme of the North-West University in South Africa. This paper is dedicated to the memory of the late Harm Moraal.

Disclosure of Potential Conflicts of Interest

The authors declare that they have no conflicts of interest.