We address two $\gamma$-ray events: 2012 January 23 and 2012 May 17 and approach the problem by, first, simulating the proton acceleration at the shock and, second, simulating their transport back to the Sun.

## 9.1 Introduction

The novel $\gamma$-ray observations by the Large Area Telescope (LAT) on the Fermi $\gamma$-ray Space Telescope spacecraft (Atwood et al. 2009), taken in a systematic way at unprecedented high energies, have presented a puzzle to the solar energetic particle (SEP) research community. More than two dozen >100 MeV $\gamma$-ray events were observed between 2008 and 2016 (see Chap. 8), many of which have properties that challenge the traditional idea that high-energy (>300 MeV) protons needed for the production of the $\gamma$-rays, via the pion-decay process deep in the chromosphere, are accelerated in solar flares (Ackermann et al. 2014). Specifically, the Fermi/LAT observations indicate that particles are precipitating to the solar atmosphere for up to a day after the impulsive phase of the flare, which is difficult to reconcile with a model of impulsive acceleration followed by particle trapping in the coronal magnetic field. On the other hand, coronal mass ejection (CME) driven shock waves can emit protons with energies above 300 MeV for several hours after the onset of the associated solar eruption as observed at 1 AU (e.g., the 2012 May 17 SEP event, see Chap. 8). Therefore, as an alternative view on the genesis of the long-duration $\gamma$-ray events, shock acceleration needs to be considered.

One of the challenges of the shock-acceleration hypothesis is that the SEP events observed in connection with the Fermi/LAT $\gamma$-ray events are not always very large, nor do they extend to very high energies when observed at 1 AU (see Chap. 8). Therefore, one of the key aspects to understand about these events is the spatial distribution of the accelerated particles at the CME-driven shock wave as well as the relation between interplanetary and interacting protons. Several factors contribute to this relation: (1) In-situ observations are local, i.e., performed in a particular interplanetary flux tube, whereas the observed high-energy $\gamma$-rays are produced over an extended emission region involving contributions from different field lines. (2) The energy spectrum of the particles accelerated at the shock is modified by transport effects when the particles propagate both downstream and upstream to reach the Sun and 1 AU. (3) Particles can modify their own transport conditions upstream of the shock due to amplification of Alfvén waves, so the fluxes observed at 1 AU can be partially decoupled from the fluxes at the shock. (4) Compressive and stochastic acceleration in the downstream region close to the CME can modify the spectrum of particles propagating toward the Sun.

We tackle the problem by conducting simulations of acceleration of protons in the shock and of their transport to 1 AU and back to the Sun for two long-duration $\gamma$-ray events: 2012 January 23 and 2012 May 17, of which the latter is associated with a ground level enhancement (GLE) of SEPs. In what follows, we outline in Sect. 9.2 the modelling techniques applied, then we present simulation results in Sect. 9.3 and their discussion as well as conclusions in Sect. 9.4.
## 9.2 Model Description

In this section, we outline the Shock-and-Particle (SaP) and Coronal Shock Acceleration (CSA) simulation models used to infer the proton spectrum at the shock and describe in detail the DownStream Propagation (DSP) model developed in the frame of High Energy Solar Particle Events foRecastIng and Analysis (HESPERIA) project and used to simulate proton transport from downstream of the shock to the Sun's surface.

### 9.2.1 Shock and Particle Model

The Shock-and-Particle (SaP) simulation model allows one to determine the injection rate of energetic particles from a propagating shock, using simulations of the shock propagation and particle transport, combined with the fitting of the simulation results to the observations. The particle injection rate $Q$ at the shock is defined as $Q = df(z_s, E, t)/dt$, where $f$ is the particle phase-space distribution function and $z_s(t)$ is the shock position along the magnetic field line, which is, of course, a function of time. In SaP, the Parker spiral magnetic field is assumed. Instead of $Q$, one can consider the injection rate integrated over the cross-section $A$ of the magnetic tube, $G = QA(z_s)$. See Chap. 4, Pomoell et al. (2015) and Aran et al. (2007) for further details.

### 9.2.2 Coronal Shock Acceleration Model

The Coronal Shock Acceleration (CSA) model is a Monte Carlo simulation model dealing with acceleration of ions in a coronal shock, taking into account ion-induced generation of Alfvén waves in the solar wind upstream of the shock. CSA simulates evolution of particles and Alfvén waves on a single radial magnetic field line. The shock is treated as a magnetohydrodynamic (MHD) discontinuity, the gas and magnetic compression ratios of which are computed through Rankine-Hugoniot relations, using the shock speed along the field line $V_s$, the shock-normal angle $\theta_{Bn}$ (the angle between the magnetic field vector **B** and the shock normal **n**) and the ambient solar wind parameters (the plasma density $n$, the magnetic field **B** and the temperature **T**) at the shock position. All these parameters vary as the shock propagates outward from the Sun, which is implemented in CSA by using analytic functions of time or heliocentric distance to describe such variations. The analytic functions are determined by a number of free parameters that have to be provided as input in a simulation (see Afanasiev et al. 2017 for details). For instance, the
variation of the magnetic field strength B with heliocentric distance r is implemented in CSA as

$$B(r) = B_0 \left(\frac{r_{\oplus}}{r}\right)^2 \left[1 + b_{rf} \left(\frac{R_{\odot}}{r}\right)^6\right], \quad (9.1)$$

where $R_{\odot}$ is the solar radius, $r_{\oplus} = 1$ AU, and $B_0$ and $b_{rf}$ are free parameters. The parameter $b_{rf}$ accounts for a super-radial expansion of the associated radial magnetic flux tube close to the solar surface. In HESPERIA, the free parameters were determined by fitting the analytic functions used to the data obtained by using techniques of semi-empirical modelling of the shock (see Rouillard et al. 2016).

The treatment of wave-particle interactions in CSA is based on quasi-linear theory, assuming only outward-propagating Alfvén waves, i.e., waves propagating away from the Sun. Particles experience elastic pitch-angle scattering in the wave frame, which is governed by the quasi-linear pitch-angle diffusion coefficient:

$$D_{\mu\mu} = \frac{\pi}{4} \Omega \frac{f_{\text{res}} P(f_{\text{res}})}{B^2} (1 - \mu^2), \quad (9.2)$$

where $P(f)$ is the wave power spectrum, $f_{\text{res}} = \Omega V/(2\pi v)$ is the resonant wave frequency, $v$ and $\mu$ are the particle speed and pitch-angle cosine as measured in the wave-rest frame, $V$ is the Alfvén wave propagation speed as measured in the solar-fixed frame and $\Omega$ is the ion cyclotron frequency. The expression given for the resonance frequency $f_{\text{res}}$ represents a simplified (pitch-angle-independent) resonance condition of particle pitch-angle scattering (see, e.g. Afanasiev et al. 2015). If considered in other reference frames, in particular in the plasma frame, pitch-angle scattering leads to the energy exchange between particles and waves. Pitch-angle scattering of particles off outward-propagating waves, as viewed in the plasma frame, can give rise to growth of waves. After a series of approximations (Vainio 2003), the wave growth is given as

$$\Gamma(f_{\text{res}}) = \frac{\pi}{2} \Omega \frac{p_{\text{res}} S(r, p_{\text{res}}, t)}{n V_{\text{A}}}, \quad (9.3)$$

where $S(r, p, t) = 2\pi \int_{-1}^{+1} p^2 \mu v f(r, p, t) d\mu$ is the particle streaming, $p_{\text{res}}$ is the resonant particle momentum corresponding to the resonance frequency $f_{\text{res}}$, $n$ is the proton density and $V_{\text{A}}$ is the Alfvén speed.

Particles are followed in the guiding-center approximation in the shock's upstream. The seed particle population is modelled by a kappa distribution in speed with an exponential cutoff energy $E_0$. The spatial distribution of seed particles is parameterised as $n_{\text{seed}} = \epsilon n$, where $n_{\text{seed}}$ is the density of seed particles and $\epsilon$ is the injection parameter. Particle injection at the shock as well as particle-shock interactions are modelled by testing particles hitting the shock from the upstream for reflection/transmission from/through the shock front and for their transport back
to the upstream, if transmitted. Alfvén waves are followed in the Wentzel-Kramers-Brillouin (WKB) approximation, with an additional diffusion term in frequency accounting for wave energy cascading. For a comprehensive description of CSA see Battarbee (2013).

### 9.2.3 DownStream Propagation Model

The DownStream Propagation (DSP) model is a Monte Carlo simulation model that has been devised to simulate proton transport from the shock's downstream back to the Sun. It is based on Parker's equation, which assumes quasi-isotropic particle distributions and, hence, diffusive transport:

$$ \frac{\partial f}{\partial t} + \mathbf{u} \cdot \nabla f - \frac{1}{3} p (\nabla \cdot \mathbf{u}) \frac{\partial f}{\partial p} = \nabla \cdot (\kappa \cdot \nabla f), \quad (9.4) $$

where $f = (d^6N/(d^3x d^3p))$ is the isotropic part of the particle distribution function, $\mathbf{u}$ is the velocity of the background fluid (solar wind plasma) and $\kappa$ is the diffusion tensor. The second term on the left-hand side describes advection of particles with the solar wind, the third term describes adiabatic cooling due to the solar wind expansion and the right-hand side term describes spatial diffusion of particles. Hence, DSP describes the propagation of particles in the test-particle approximation.

Assuming that particles are confined within a magnetic flux tube, which gives $\kappa = \kappa \mathbf{b} \mathbf{b}$ with $\mathbf{b}$ being a unit vector along the magnetic field, and that $\mathbf{u} \parallel \mathbf{b}$, Parker's equation can be reduced to

$$ \frac{\partial f}{\partial t} + u \frac{\partial f}{\partial z} - \frac{1}{3} p \frac{1}{A} \frac{\partial (Au)}{\partial z} \frac{\partial f}{\partial p} = \frac{1}{A} \frac{\partial}{\partial z} \left( A \kappa \frac{\partial f}{\partial z} \right), \quad (9.5) $$

where $z$ is the (curvilinear) coordinate measured from the solar surface along the field and $A(z)$ is the cross-sectional area of the flux tube. Changing to $F = 4\pi p^2 Af$, one can obtain

$$ \frac{\partial F}{\partial t} + \frac{\partial}{\partial z} \left[ \left( u + \frac{1}{A} \frac{\partial (A\kappa)}{\partial z} \right) F \right] - \frac{\partial}{\partial p} \left( \frac{1}{3} p \frac{1}{A} \frac{\partial (Au)}{\partial z} F \right) = \frac{1}{2} \frac{\partial^2}{\partial z^2} (2\kappa F) . \quad (9.6) $$

Equation (9.6) is equivalent to the following set of stochastic differential equations (SDEs):

$$ dz = \left( u + \frac{\kappa}{L} + \frac{\partial \kappa}{\partial z} \right) dt + \sqrt{2\kappa} dW_t, \quad (9.7) $$

$$ dp = -\frac{1}{3} p \left( \frac{u}{L} + \frac{\partial u}{\partial z} \right) dt, \quad (9.8) $$
<figure><img src="image_73.png" /><figcaption>**Fig. 9.1** Spatial simulation domain in the DSP model</figcaption></figure>

where $L(z) = A/(dA/dz)$ is the focusing length of the magnetic flux tube, $dW_t$ is a stochastic differential normally distributed with zero mean and variance $dt$. Note that the term $(\kappa/L)dt$ in Eq. (9.7) represents the net effect of focusing to a quasi-isotropic particle population (Kocharov 1996). These SDEs are solved for each Monte Carlo particle in the simulation, using the standard explicit Euler-Maruyama method.

The simulation is performed in an expanding 1-D spatial simulation box in the shock's downstream, i.e., the box is confined by the solar surface from one side and by the moving shock front from the other side, and, thus, expands with time along the magnetic field line (Fig. 9.1). For the plasma speed $u$ and the spatial diffusion coefficient $\kappa$ in the simulation box, the following linear models are adopted:

$$u(z, t) = \frac{u_0 z}{V_s t}, \qquad (9.9)$$

$$\kappa(z, p, t) = \frac{1}{3} \lambda_0 \frac{p}{m_p} \frac{V_s t - z}{V_s t}, \qquad (9.10)$$

where $u_0$ is the plasma speed immediately behind the shock front and $V_s$ is the shock speed along the field, both measured in the solar-fixed frame, $m_p$ is the proton mass and $\lambda_0$ is the particle mean free path at the solar surface. One can see from Eqs. (9.9) and (9.10) that the plasma speed decreases linearly from $u = u_0$ at the shock ($z = V_s t$) to $u = 0$ at the Sun ($z = 0$), and the mean free path recovers from $\lambda = 0$ at the shock to $\lambda = \lambda_0$ at the Sun. Hence the transport of particles from the upstream through the shock is purely advective. On the other hand, the linear dependence of the spatial diffusion coefficient $\kappa$ on $V_s t - z$ secures that no particles escape from the downstream beyond the shock.

The focusing length $L$ is specified in a form similar to the one that can be derived from Eq. (9.1):

$$L(z) = \frac{1}{2}(R_{\odot} + z) \left[ \left( 1 + \frac{z}{R_{\odot}} \right)^6 + b_f \right] \left[ \left( 1 + \frac{z}{R_{\odot}} \right)^6 + 4b_f \right]^{-1}, \qquad (9.11)$$

where $b_f$ is a free parameter.
The initial size of the simulation box is specified by the initial position of the shock $z_{s0}$. This gives the initial value for time $t$ as $t_0 = z_{s0}/V_s$. The shock speed $V_s$ is assumed constant and $u_0 = V_s(r_c - 1)/r_c$ with $r_c$ being the gas compression ratio taken constant as well.

To be consistent with the chosen form of $\kappa$, the amount of particles being injected into the simulation from the shock is determined by the net particle flux to the far downstream. This flux is given by

$$ \mathcal{F} = \int d^3 p \mu v f(p, \mu) = \int d^3 p' u_2 f'(p') , \qquad (9.12) $$

where $\mu$ is the pitch-angle cosine, $v$ is the particle speed, $u_2$ is the downstream solar wind speed; the unprimed symbols designate quantities as measured in the shock frame and the primed symbols in the downstream plasma frame. The second expression in Eq. (9.12) is derived under the assumption that the particle pitch-angle distribution as measured in the downstream plasma frame is isotropic. Then, it can be derived to the first order in $u_2/c$ ($c$ is the speed of light in vacuum), taking into account that $\lambda \to 0$ at the shock, that

$$ \frac{d\mathcal{F}}{dp} = 4\pi u_2 p^2 f_0(p), \qquad (9.13) $$

where $f_0(p)$ is the isotropic part of the particle distribution function at the shock (on the upstream side). This equation is used to relate the amount of particles injected to the downstream with the particle intensity at the shock, $j_s = p^2 f_0(p)$, which is the output of CSA simulations. Specifically, the number of particles per unit momentum injected downstream from a source (at the shock) of size $A_s(z)$ in time $dt$ is determined by

$$ \frac{dN}{dp} = 4\pi \frac{V_s}{r_c} A_s(z) j_s(t, E) dt . \qquad (9.14) $$

In case the particle injection rate $G(t, E)$ is given (output of SaP), the particle intensity at the shock is computed as

$$ j_s(t, E) = \frac{p^2 C \Delta r_g G(t - t_0, E)}{v r^2 \cos \psi(r)}, \qquad (9.15) $$

where $v$ is the particle speed, $C = 3.524$ is a constant determined by the numerical scheme implemented in SaP, $\Delta r_g$ is the spatial grid size in SaP simulations, $r$ is the heliocentric distance and it is taken into account that in SaP the injection rate $G$ is determined under the assumption of a Parker-spiral upstream magnetic field, i.e., $\cos \psi(r) = [1 + (\Omega_{\odot}r/u_{sw})^2]^{-1/2}$ with constant solar wind speed $u_{sw}$ ($\Omega_{\odot}$ is the angular speed of the solar rotation).
The proton injection into the downstream is implemented in the following way. We introduce an injection time step $\Delta t$ and deposit a certain number of Monte Carlo particles into the simulation box at each $\Delta t$. A Monte Carlo particle is a representative of a group of physical particles and is characterised by its weight $w$. To obtain weights of injected Monte Carlo particles at a given time $t_0$, we compute values of the spectrum $dN/dp$ at the momentum grid using Eq. (9.14) and interpolate those spectral values by power laws between the grid nodes. Then, the weight of a Monte Carlo particle with momentum $p$ between the grid points $p_j$ and $p_{j+1}$ is given by

$$w = \frac{1}{N_{\text{MC}}} \int_{p_j}^{p_{j+1}} \frac{dN}{dp} dp, \qquad (9.16)$$

where $dN/dp \propto p^{-q_j}$ and $N_{\text{MC}}$ is the number of Monte Carlo particles assigned to the interval $[p_j, p_{j+1}]$. Note that the particle momentum $p$ is randomly chosen from a $p^{-q_j}$ distribution in this interval and $q_j = -\ln(S_{j+1}/S_j)/\ln(p_{j+1}/p_j)$ with $S_{j(j+1)} = dN/dp|_{j(j+1)}$.

The Monte Carlo particles injected at a given time $t_0$ are placed into a small spatial 1-D volume $\Delta z = u_2\Delta t$ behind the front, to mimic their advection by the bulk plasma during $\Delta t$. Then, each of these particles is propagated by solving the SDEs until it hits the Sun, i.e., precipitates, or up to $t = t_{\text{max}}$. The particle transport time step $\delta t$ must fulfil the condition $\delta t \ll L^2/\kappa$. Furthermore, the linear dependence of $\kappa$ on distance from the shock dictates that Monte Carlo particles that might appear in the upstream side during a simulation have to be returned to the downstream side. This requires a reflective boundary condition for particles to be applied at the shock.

During the simulation, Monte Carlo particles that hit the Sun are collected and, based on their momenta, precipitation times and weights, the flux and time-integrated energy spectrum of absorbed protons are calculated. The latter is compared with the proton spectra derived from $\gamma$-ray observations.

## 9.3 Results

The simulations based on the models presented above were conducted for two long-duration >100 MeV $\gamma$-ray events occurred on 2012 January 23 and 2012 May 17. These $\gamma$-ray events are associated with substantial SEP events as observed at 1 AU (the 2012 May 17 event is also a GLE). The event characteristics are described in Chap. 8 (see also Rouillard et al. 2016). In what follows we first present modelling results for the 2012 May 17 event due to its association with a GLE and then results for the 2012 January 23 event.
## 9.3.1 2012 May 17 Event

### 9.3.1.1 Modelling of the SEP Event

We begin with results of the SaP modelling of the accelerated proton population at the shock, which utilises the observations of the shock and SEP event at 1 AU. Figure 9.2 presents results of simulations of the ambient solar wind and shock propagation to the Earth. The solar wind was simulated starting from the onset time of the SEP event. Table 9.1 provides values of the input solar wind parameters at $1.03 R_⊙$. The input parameters of the shock-driving disturbance are its density $ρ_{cme} = 0.3 × 10^{-13}$ kg m$^{-3}$, speed $v_{cme} = 1000$ km s$^{-1}$, and the shock-front shape parameter $Δφ = 0.5a_{cme}$ determined by the angular extent $a_{cme}$ of the disturbance (see Pomoell et al. 2015 for details). One can see that the simulations reproduce quite well the average characteristics of the solar wind at 1 AU prior to the shock arrival (note, however, that the temperature is somewhat underestimated) as well as the shock arrival time. At the same time, the observed jumps in the plasma characteristics are reproduced well with some underestimation of the magnetic field increase.

In order to derive the proton injection rate from the shock into the magnetic flux tube connecting the observer with the shock front, we fitted the proton

<figure><img src="image_74.png" /><figcaption>**Fig. 9.2** Simulated characteristics of the solar wind (blue lines) superposed on the observed characteristics (shown in orange and red) before and during the shock arrival at 1 AU for the 2012 May 17 SEP event. Shown are: the plasma density, the solar wind radial velocity, the magnetic field magnitude and the solar wind proton temperature. The temperature was measured by the Wind spacecraft and the other parameters by the ACE spacecraft. Time counts from the time of the interplanetary shock passage by the ACE spacecraft (marked by a dotted vertical line)</figcaption></figure>
<table><caption>**Table 9.1** Input solar wind parameters used in the SaP simulations of the 2012 January 23 and 2012 May 17 SEP events</caption><thead><tr><th>Event</th><th>ρ<sub>0</sub> (kg/m<sup>3</sup>)<sup>a</sup></th><th>T<sub>0</sub> (K)</th><th>S<sub>0</sub> (W/m<sup>3</sup>)</th><th>L (R<sub>☉</sub>)</th><th>γ</th><th>B<sub>r0</sub> (T)</th></tr></thead><tbody><tr><td>Jan 2012</td><td>1.169 × 10<sup>-13</sup></td><td>1.22 × 10<sup>6</sup></td><td>0.335 × 10<sup>-7</sup></td><td>0.735</td><td>1.55</td><td>-2.15 × 10<sup>-4</sup></td></tr><tr><td>May 2012</td><td>1.169 × 10<sup>-13</sup></td><td>1.18 × 10<sup>6</sup></td><td>0.35 × 10<sup>-7</sup></td><td>0.70</td><td>1.55</td><td>1.26 × 10<sup>-4</sup></td></tr></tbody></table>

<sup>a</sup>The parameters provided are the plasma density ρ<sub>0</sub>, temperature T<sub>0</sub> and the radial component of the magnetic field B<sub>r0</sub> at the heliocentric distance r<sub>0</sub> = 1.03 R<sub>☉</sub>; S<sub>0</sub> and L are the coronal heating function parameters (Pomoell et al. 2015); and γ is the adiabatic index.

<figure><img src="image_75.png" /><figcaption>**Fig. 9.3** Examples of the intensity-time profiles from the SEPEM reference data and the GOES/HEPAD detector (shown in red) with superposed fits (light blue lines). The short green horizontal lines indicate the pre-event background intensities that have been subtracted from the measured intensities. The two solid blue vertical lines indicate the times of the shock arrival at the ACE and GOES spacecraft positions. The dotted orange vertical line indicates the onset of the associated X-ray flare and the dotted purple line indicates the time of a change in the interplanetary magnetic field direction</figcaption></figure>

intensity-time profiles provided by the ACE/EPAM instrument (in the energy range 0.59–4.8 MeV), SEPEM reference data (Jiggens et al. 2012) (6–166.3 MeV) and GOES/HEPAD detector (330–700 MeV). The fitting was performed using SaP particle transport simulations, following the method described in Pomoell et al. (2015). Figure 9.3 shows examples of the fitted intensity-time profiles for several high-energy channels. Note that when performing the fitting, we focused on the first several hours of the event as the corresponding γ-ray event lasted for about 2 h.

Figure 9.4 shows the proton injection rates G(t) resulting from the fitting of the intensity-time profiles at selected energies. The obtained injection rates were used
<figure><img src="image_76.png" /><figcaption>**Fig. 9.4** Temporal evolution of the proton injection rate *G* at high energies. *Black lines* represent polynomial fits to the simulated *G* functions</figcaption></figure>

as input for the DSP simulations. The simulations reveal that the injection rate at high energies drops by more than one order of magnitude several hours after the X-ray flare onset, which agrees with rather short duration of the associated γ-ray event (~2 h).

9.3.1.2 Simulations of Proton Acceleration at the Shock

To carry out the CSA simulations of proton acceleration at the shock in this
event, we utilised the ambient plasma and shock parameters derived from the semi-
empirical modelling of the shock (Rouillard et al. 2016). Those parameters are the
plasma density *n*, the magnetic field strength *B*, the shock speed *V*<sub>s</sub> and the shock-
normal angle θ<sub>Bn</sub>, which were determined along individual magnetic field lines. The
parameters were fitted by the analytic functions of time/distance implemented in
CSA (Afanasiev et al. 2017). Figure 9.5 shows an example of the data obtained for
a single field line, superposed by the corresponding fits. In total, the data for over
100 field lines were available and fitted. Among the field lines possessing good-
quality fits, nine were chosen for simulations with CSA.

The detailed analysis of the simulation results presented in Afanasiev et al. (2017) reveals that the parameter mainly controlling the acceleration efficiency is the Alfvénic Mach number of the shock. Figure 9.6 shows examples of the evolution of the Alfvénic Mach number of the propagating shock for three different magnetic field lines and simulated proton energy spectra at the shock, corresponding to these field lines, obtained at *t* = 1000 s. The simulated spectra (for one of the simulated field lines) were then used as the other input for the DSP simulations.
<figure><img src="image_77.png" /><figcaption>**Fig. 9.5** Example data set obtained with the semi-empirical shock modelling approach in the 2012 May 17 event, plotted together with the corresponding fits. The *upper left panel* shows the shock-normal cosine cos θ<sub>Bn</sub> (“plus” symbols) and the inverse Alfvénic Mach number M<sub>A</sub><sup>-1</sup> (asterisks) vs. time; the *upper right panel* shows the shock speed V<sub>s</sub> along the magnetic field line vs. time. Note that time is counted from the moment when M<sub>A</sub> > 1.5. The *bottom left panel* shows the magnetic field magnitude B vs. radial shock position and the *bottom right panel* shows plasma density n vs. radial shock position</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>**Fig. 9.6** Alfvénic Mach number of the shock versus simulation time for three different magnetic field lines (**left panel**) and corresponding proton energy spectra at the shock at *t* = 1000 s, resulting from CSA simulations (**right panel**). Note the correspondence between the Alfvénic Mach number and the spectral cutoff energy</figcaption></figure>
### 9.3.1.3 Modelling of the Proton Transport Back to the Sun

The proton transport back to the Sun was simulated with the DSP model, assuming a radial magnetic field in the shock’s downstream. The particle source size at the shock was modelled as $A_s(z) = (R_\odot + z)^2\Omega_0$. The parameter $\Omega_0$, which can be interpreted as the global angular size of the shock, was taken to be 1 steradian. In fact, the realistic effective source size should be smaller because of the substantial difference in the particle acceleration efficiency along different field lines, as revealed by the CSA simulations for this event. We took this into account by considering an additional parameter, so-called filling factor, $a_{fill}$ that is the relative fraction of field lines at which high-energy, $\gamma$-ray-productive, protons can be produced. Based on the ascertained dependence between the particle acceleration efficiency and the Alfvénic Mach number magnitude, and available data on >100 individual field lines, we estimated that $a_{fill} = 0.1$. To compare the proton spectra resulting from the DSP simulations with the proton spectrum obtained from the Fermi/LAT observation, we multiply the simulated spectra by $a_{fill}$. The parameters of the observationally-derived spectrum (the total number of >500 MeV protons and the power-law spectral index of >300 MeV protons) were kindly provided to us by G. Share (see Chap. 8 for references). The other DSP model parameters were taken to be $z_{s0} = 1.6 R_\odot (1.15 R_\odot)$, in the case of CSA(SaP) input, $V_s = 1510$ km s$^{-1}$, $r_c = 3.6$ and $u_{sw} = 387$ km s$^{-1}$.

The DSP simulations were conducted for a set of values for the downstream transport parameters $\lambda_0$ and $b_f$. The resulting (integrated over the duration of the $\gamma$-ray emission; i.e., 2 h) energy spectra of protons hitting the Sun are presented in Fig. 9.7 along with the shock-injected spectrum and the observationally-derived

<figure><img src="image_79.png" /><figcaption>**Fig. 9.7** Time-integrated energy spectra of protons precipitated at the Sun, resulting from DSP simulations using the results of SaP (**left panel**) and CSA simulations (**right panel**) for the 2012 May 17 event. The spectra are shown by **blue and green lines** with the corresponding DSP model parameters indicated. The integration time is 2 h, which is approximately the duration of the >100 MeV $\gamma$-ray event. Also shown are the time-integrated spectrum of protons injected by the shock to the downstream (**black dashed line**) and the spectrum of interacting protons derived from the Fermi/LAT observation (**red line**)</figcaption></figure>
spectrum of interacting >300 MeV protons. It can be calculated that the number of >500 MeV protons injected from the shock in the CSA+DSP simulation exceeds by more than $10^4$ the corresponding number of protons derived from the observation. On the other hand, the number of absorbed high-energy protons is sensitive to the transport parameters. In particular, it can be easily reduced by increasing the parameter $b_f$ controlling the focusing length, which enhances proton mirroring from the flux tube base and adiabatic cooling, in accord with Eqs. (9.7) and (9.8). Note also that the DSP model completely neglects the possibility for particles in the downstream side of the shock to escape to the upstream side. This process, if taken into account, should decrease the number of absorbed particles as well.

In contrast, the spectrum of shock-injected protons obtained from the SaP+DSP simulation and the absorbed spectrum corresponding to $b_f = 0$ (radial flux tube) are in good correspondence with the observed spectrum. However, it should be noted that this correspondence holds only for rather idealistic conditions of the DSP model (no particle escape to the upstream) and one can expect a lack of high-energy protons in the simulations, if a more realistic downstream transport model is considered. The possible reasons of this result are discussed in Sect. 9.4.

## 9.3.2 2012 January 23 Event

### 9.3.2.1 Modelling of the SEP Event

The results of simulations of the ambient solar wind and shock propagation to 1 AU are shown in Fig. 9.8. The input parameters describing the initial shock-driven disturbance in this case are $\rho_{cme} = 0.3 \times 10^{-13}$ kg m$^{-3}$, $v_{cme} = 1650$ km s$^{-1}$, and $\Delta\phi = 0.25a_{cme}$ (see also Table 9.1 for the initial solar wind parameters). Like for the 2012 May 17 event, the modelling reproduces well the average characteristics of the solar wind at 1 AU prior to the shock arrival and the shock arrival time. The observed jump in the magnetic field magnitude is reproduced well too, but the jumps in the density, the solar wind speed and the temperature are somewhat overestimated. Note that due to the data gap in the ACE data, the value of the average temperature ($2.5 \times 10^4$ K) in the upstream region has been taken from the plot provided by the IP shocks data base of the University of Helsinki,<fnref n="1" /> based on the WIND data. A similar value is estimated by the CfA interplanetary shock list.<fnref n="2" />

Like for the previous event, we fitted the observed intensity-time profiles provided by the ACE/EPAM instrument and the SEPEM reference data (Jiggens et al. 2012). The proton enhancements observed in this event by GOES/HEPAD in the high-energy channels ranging from 330 to 700 MeV are weak. From Fig. 9.9, it can be seen that already at 166.3 MeV the background-subtracted enhancement is lower than the background level itself. For this reason, instead of fitting the observed

<footnote n="1">¹http://ipshocks.fi/database.</footnote>

<footnote n="2">²https://www.cfa.harvard.edu/shocks.</footnote>
<figure><img src="image_80.png" /><figcaption>**Fig. 9.8** Same as in Fig. 9.2, but for the 2012 January 23 event</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>**Fig. 9.9** Left panel: Examples of the observed proton intensity-time profiles (in red) with superposed fits/synthetic profiles (light blue lines) for the 2012 January 23 SEP event. The short green horizontal lines indicate the pre-event background intensities that have been subtracted from the measured intensities. The meaning of the other lines is the same as in Fig. 9.3. Right panel: Temporal evolution of the proton injection rate G at high energies. Black lines are polynomial fits to the simulated G functions</figcaption></figure>
intensities, we computed synthetic intensity-time profiles for the GOES/HEPAD energy channels by extrapolating the particle injection rate $G(t, E)$ from the highest energies in the SEPEM reference dataset. Such a profile is shown in Fig. 9.9 at 605 MeV. Figure 9.9, right panel shows the proton injection rate resulting from the SaP simulations of this SEP event.

### 9.3.2.2 Simulation of Proton Acceleration at the Shock

To conduct CSA simulations of proton acceleration at the shock in the same fashion as it was done for the 2012 May 17 event, we used white-light (WL) images of the corona from SOHO and STEREO A and B and the Potential Field Source Surface (PFSS) modelling of the magnetic field. The WL images were fitted using a spherical representation for the CME. This allowed us to obtain the magnetic field, the shock speed and the shock-normal angle along different field lines (Fig. 9.10). As concerns the plasma density, we used the following representation:

$$n(r) = n_2 \left(\frac{r_{\oplus}}{r}\right)^2 + n_8 \left(\frac{R_{\odot}}{r}\right)^8, \qquad (9.17)$$

with $n_2 = 10 \text{ cm}^{-3}$ and $n_8 = 8 \times 10^8 \text{ cm}^{-3}$. Since the plasma density was not constrained using the real observations, we performed a CSA simulation only for one magnetic field line. The simulation showed a typical buildup of the proton energy with time and formation of a power-law energy spectrum with a roll-over, similar to the spectra presented in Fig. 9.6. The maximum proton energy at the shock attained the value ~1 GeV at $t = 1130$ s after the start of the simulation, but then was decreasing and reached ~500 MeV at $t = 2000$ s.

<figure><img src="image_82.png" /><figcaption>Fig. 9.10 The shock-normal cosine $\mu_s$, the shock speed along the field line $V_s$ and the magnetic field $B$ with the fits superimposed for a selected magnetic field line</figcaption></figure>
<figure><img src="image_83.png" /><figcaption>**Fig. 9.11** Time-integrated energy spectra of protons precipitated at the Sun, resulting from DSP simulations using the results of SaP (**left panel**) and CSA simulations (**right panel**) for the 2012 January 23 event. The spectra are shown by *blue and green lines* with the corresponding DSP model parameters indicated. The integration time is 6 h. Also shown are the time-integrated spectrum of protons injected by the shock to the downstream (**black dashed line**) and the proton spectra obtained from the *Fermi/LAT observation* (**red line**)</figcaption></figure>

### 9.3.2.3 Modelling of the Proton Transport Back to the Sun

Similar to the DSP modelling for the 2012 May 17 event, we assumed a radial magnetic field in the shock's downstream and $A_s(z) = (R_\odot + z)^2\Omega_0$ was taken for the particle source at the shock. The model parameters were taken to be $z_{s0} = 0.6 R_\odot (3.1 R_\odot)$ in the case of CSA(SaP) input, $V_s = 1450 \text{ km s}^{-1}$, $r_c = 3.7$, $u_{sw} = 414 \text{ km s}^{-1}$, $\Omega_0 = 1 \text{ sr}$, and $a_{\text{fill}} = 0.1$.

Figure 9.11 shows the simulated spectra of protons absorbed at the Sun together with the shock-injected spectrum and the proton spectrum derived from the *Fermi/LAT* observation of >100 MeV γ-photons. Similar to the simulation for the 2012 May 17 event, in the CSA+DSP simulation, there is a substantial excess of shock-injected protons as compared to the number of interacted protons derived from the observation. In contrast, there is a lack of shock-injected protons against the observed number in the SaP+DSP simulation.

## 9.4 Discussion and Conclusions

We have modelled particle acceleration at coronal shocks driven by CMEs and proton transport from the shock to both the Sun and the far upstream region (towards the 1-AU observer). The purpose of our study is to find out, whether shock-accelerated protons streaming back from the shock could be responsible for the long-duration γ-ray events observed by *Fermi/LAT*. We simulated the shock propagation from the Sun to 1 AU using a two-dimensional MHD model. We also employed the empirical models of the CME-driven pressure front propagation,
which allowed us to assess the early evolution of the shock in a system that does not possess the symmetries assumed in the MHD model.

Our results show that the efficiency of particle acceleration crucially depends on the modelled properties of the shock in the corona. Conditions on different field lines vary very much and while the shock on some field lines is able to produce a relativistic particle event, it fails to do so on others. The most important factor governing the acceleration efficiency in our study was the Alfvénic Mach number of the shock: the higher the Alfvénic Mach number, the more likely the shock to accelerate protons to relativistic energies (Afanasiev et al. 2017).

In our simulations, we focused on two events that differ in one important aspect: one of them (2012 May 17) is a GLE and the other one is not. Both of the events are associated with long-duration γ-ray events, which might seem contradictory, but actually is not. In the light of our simulations, there are two possible explanations for this. Firstly, as the particle acceleration efficiency at the shock varies a lot from one flux tube to another, the 1-AU proton event does not necessarily correspond to the best acceleration conditions on the shock surface. Thus, we may well observe a long-duration (pion-decay) γ-ray event due to shock-accelerated protons without a clear increase observed at 1 AU at energies required for pion production at the Sun. The Earth-based observer sees but a small fraction of the complete picture. The other, more subtle explanation deals with the strength of the turbulence in the foreshock region of the coronal shock. The CSA simulations show that the foreshock is extremely turbulent near the Sun, traps a large fraction of particles (almost all) in its vicinity and allows only a minor fraction to escape. While the SaP model contains the possibility to use enhanced foreshock turbulence, it is tuned to reproduce the observations when the shock is detected in situ. Therefore, the source function deduced from the 1-AU observation represents more the fraction of particles that can escape upstream than the fraction that can be transported downstream from the shock. The large discrepancy between the CSA and SaP modelled spectra of precipitated particles is, thus, partly explained by this effect, as well. Furthermore, these explanations shed light on the tendency to get lower numbers of precipitated >300 MeV protons in the SaP+DSP simulations (especially for the 2012 January 23 event), as compared to the observations.

The transport model we employ for the downstream region has several important simplifications in it. Firstly, it employs a shock completely opaque to protons and, thus, allows downstream-advected particles to reside in the region between the shock and the Sun for as long as they get precipitated. The only loss process we employ is adiabatic cooling of the distribution, when the region between the shock and the Sun expands. On the other hand, we do not include any downstream re-acceleration processes, which could also be important and would act in the opposite direction, helping particles to overcome adiabatic energy losses. Such processes include downstream stochastic acceleration (see, e.g., Afanasiev et al. 2014) and compressive acceleration in the highly-compressed regions close to the CME core (Kozarev et al. 2013). Therefore, we do not regard our model to be overly optimistic about the prospect of letting shock-accelerated protons precipitate over large time scales.
Regarding the CSA model, one important reservation has to be made: the model makes a simplification to the quasi-linear resonance condition between particles and scattering Alfvén waves, $k_{\text{res}} = 2\pi f_{\text{res}}/V = B/(R\mu)$, neglecting the dependence of resonant wave number $k_{\text{res}}$ on particle pitch-angle cosine $\mu$ ($R$ is the particle rigidity). This simplification, while without proper physical justification, allows one to build the code using the assumption of isotropic scattering, which speeds up the running times easily by a factor of ten over times obtained when treating the resonance condition in full. We have, however, evaluated the effect of this simplification using a local model (Afanasiev et al. 2015), and shown that in parallel coronal shocks the difference between the two models yields about a factor of 2 in the roll-over rigidity obtained from the model. On the other hand, the spectrum in CSA also cuts off much more rapidly than in the model employing the complete resonance condition, so we do not regard this to be a very serious problem in the performed modelling. A more complete global model to be developed in future, however, should take the full resonance condition into account, also since it affects the foreshock spatial structure as well (Afanasiev et al. 2015). We will undertake the development of such a model in future projects.

Another transport process missing from the CSA model is diffusion perpendicular to the mean magnetic field. This process can be implemented in a Monte Carlo simulation, but its inclusion will require to incorporate at least one more spatial dimension in the model. Therefore, also the requirement for particle statistics in CSA will be tremendously increased to avoid statistical noise in the result, as the number of spatial cells in the model will have to be increased by a factor 30–300, depending on the coarseness of the grid in the perpendicular direction. This is still beyond the reach of the present computers with reasonable running times of the code. Fortunately, as perpendicular diffusion cannot occur due to slab-mode waves, the enhanced Alfvénic turbulence in the upstream region is not strengthening the perpendicular diffusion of the particles. However, one would expect the downstream plasma to have much more isotropic turbulence which, then, would lead to the migration of particles from one flux tube to the other while they are on the downstream side of the shock. For an opaque shock, like we have assumed, this is not affecting the acceleration of ions at the shock too much since their fate (scattered back to the shock or transported to the far downstream region with no return) would be decided (almost) instantly, giving the particles very little time to diffuse perpendicular to the mean field. Thus, the first step to take the perpendicular transport into account would have to be performed on the downstream side of the shock.

We note that the total number of >500 MeV protons as simulated by the CSA code is several orders of magnitude larger than the observational value in both of the simulated events. At first, this might seem to be problematic. However, in addition to the caveats about the resonance condition and the downstream transport modelling discussed above, we point out that the CSA model is set up quite favourably for efficient particle acceleration: we use a seed particle population with a relatively hard suprathelial distribution ($\kappa = 2$) in the model, which guarantees an efficient injection to the acceleration process at all obliquity angles of the shock. Using a
thermal population, only, would limit the injection efficiency of the shock especially at obliquity angles greater than ~20° quite substantially: according to Battarbee et al. (2013), the injection efficiency in an oblique coronal shock would go down by an order of magnitude if we would use a steeper seed population with $κ = 15$. Given all the possible ways to limit the number of precipitating protons in our model, we would regard the result of getting more than enough high-energy protons precipitating at the Sun to be a supporting rather than a countering result for the shock acceleration scenario.

The final shortcoming in the shock models we have employed is their assumption of the open topology of the upstream magnetic fields. Especially if the γ-ray event occurs during a period of closely spaced CMEs, the second one may drive a shock through a set of large closed loop-like or flux-rope structures, which would alter both the shock acceleration properties and, more importantly, the ability of the particles to escape upstream from the shock. In fact, developing codes capable of modelling particle acceleration and transport in more complicated upstream fields than the radial/Archimedean-spiral fields could be listed as one of the most urgent things needing improvement on the way towards physical space-weather modelling capabilities.

One of the most difficult things to estimate is the size of the source of near-relativistic protons in the event. While the 3-D modelling of the shock front can be performed in a relatively accurate and detailed manner, high-resolution density and magnetic-field structure of the corona are crucial for the correct determination of the shock properties and, thus, the total number of interacting protons in the event. Therefore, even a fully global 3-D model of coronal shock evolution and particle acceleration might not capture every detail affecting the total number of relativistic protons in the CME system. In this work, we resorted to estimating numbers based on the filling factor of field lines being capable of facilitating proton acceleration to relativistic energies based on an evaluation of shock properties on a large set of field lines. We believe that such statistical method to estimate the total number of interacting protons is the most efficient way to address the problem.

In conclusion, while a number of simplifications have been introduced in the modelling performed, we have still demonstrated that coronal shock acceleration and subsequent diffusive downstream particle transport is a viable option to explain pion-decay γ-ray events from the Sun observed by Fermi/LAT. More elaborated simulation models are needed (especially for the particle transport back to the Sun) but our results serve as a motivation by indicating that the end result of this vast modelling effort can be positive.

## References

* Ackermann, M., Ajello, M., Albert, A., et al.: Astrophys. J. 787, 15 (2014). doi:10.1088/0004-637X/787/1/15

* Afanasiev, A., Vainio, R., Kocharov, L.: Astrophys. J. 790, 36 (2014). doi:10.1088/0004-637X/790/1/36
Afanasiev, A., Battarbee, M., Vainio, R.: Astron. Astrophys. **584**, A81 (2015). doi:10.1051/0004-6361/201526750

Afanasiev, A., Vainio, R., Rouillard, A., et al.: Modelling of proton acceleration in application to a ground level enhancement. Astron. Astrophys. (2017, Submitted)

Aran, A., Lario, D., Sanahuja, B., et al.: Astron. Astrophys. **469**, 1123 (2007). doi:10.1051/0004-6361/20077233

Atwood, W.B., Abdo, A.A., Ackermann, M., et al.: Astrophys. J. **697**, 1071 (2009). doi:10.1088/0004-637X/697/2/1071

Battarbee, M.: Acceleration of solar energetic particles in coronal shocks through self-generated turbulence. Ph.D. thesis, University of Turku, Finland (2013). http://urn.fi/URN:ISBN:978-951-29-5574-9

Battarbee, M., Vainio, R., Laitinen, T., Hietala, H.: Astron. Astrophys. **558**, A110 (2013). doi:10.1051/0004-6361/201321348

Jiggens, P.T.A., Gabriel, S.B., Heynderickx, D., et al.: IEEE Trans. Nucl. Sci. **59**, 1066 (2012). doi:110.1109/TNS.2012.2198242

Kocharov, L.G., Torsti, J., Vainio, R., Kovaltsov, G.A.: Sol. Phys. **165**, 205 (1996). doi:10.1007/BF00149100

Kozarev, K.A., Evans, R.M., Schwadron, N.A., et al.: Astrophys. J. **778**, 43 (2013). doi:10.1088/0004-637X/778/1/43

Pomoell, J., Aran, A., Jacobs, C., et al.: J. Space Weather Space Clim. **5**(27), A12 (2015). doi:10.1051/swsc/2015015

Rouillard, A.P., Plotnikov, I., Pinto, R.F., et al.: Astrophys. J. **833**, 45 (2016). doi:10.3847/1538-4357/833/1/45

Vainio, R.: Astron. Astrophys. **406**, 735 (2003). doi:10.1051/0004-6361:20030822

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4.0/, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

<figure><img src="image_84.png" /></figure>
**Chapter 10**

# Inversion Methodology of Ground Level Enhancements

B. Heber, N. Agueda, R. Bütikofer, D. Galsdorf, K. Herbst, P. Kühl, J. Labrenz, and R. Vainio

**Abstract** While it is believed that the acceleration of Solar Energetic Particles (SEPs) is powered by the release of magnetic energy at the Sun, the nature, and location of the acceleration are uncertain, i.e. the origin of the highest energy particles is heavily debated. Information about the highest energy SEPs relies on observations by ground-based Neutron Monitors (NMs). SEPs with energies above 500 MeV entering the Earth's atmosphere will lead to an increase of the intensities recorded by NMs on the ground, also known as Ground Level Event or Ground Level Enhancement (GLE). A Fokker-Planck equation well describes the interplanetary transport of near relativistic electrons and protons. An NM is an integral counter defined by its yield function. From the observations of the NM network, the additional solar cosmic ray characteristics (intensity, spectrum, and anisotropy) in the energy range $\geq$500 MeV can be assessed. If the interplanetary magnetic field outside the Earth magnetosphere is known (see Sect. 10.3.2) a computation chain can be set up in order to calculate the count rate increase of an NM for a delta injection at the Sun along the magnetic field line that connects the Sun with the Earth (Sect. 10.3.3). By this computations, we define a set of

B. Heber (✉) • D. Galsdorf • K. Herbst • P. Kühl • J. Labrenz

Christian-Albrechts-Universität zu Kiel, Kiel, Germany

e-mail: heber@physik.uni-kiel.de; galsdorf@physik.uni-kiel.de; herbst@physik.uni-kiel.de;

kuehl@physik.uni-kiel.de; labrenz@physik.uni-kiel.de

N. Agueda

University of Barcelona, Barcelona, Spain

e-mail: agueda@fqa.ub.edu

R. Bütikofer

University of Bern, Physikalisches Institut, Sidlerstrasse 5, CH-3012 Bern, Switzerland

High Altitude Research Stations Jungfraujoch and Gornergrat, Sidlerstrasse 5, CH-3012 Bern, Switzerland

e-mail: rolf.buetikofer@space.unibe.ch

R. Vainio

Department of Physics and Astronomy, University of Turku, Turku, Finland

e-mail: rami.vainio@utu.fi
Green's functions that can be fitted to an observed GLE to determine the injection time profile. If the latter is compared to remote sensing measurements like radio observations, conclusions of the most probable acceleration process can be drawn.

## 10.1 Introduction

The Earth is constantly bombarded by high-energy particles, also known as cosmic rays. Those who are not deflected by the geomagnetic field, as discussed in Chap. 5, enter the atmosphere and undergo interactions with atoms and molecules as well as with the nuclei of the atmosphere. Low-energy cosmic rays (<500 MeV) are absorbed in the atmosphere, i.e. no secondary particles can be observed at sea level from primary particles as described in Chap. 6. Most of the primary particles are of galactic origin and are known as Galactic Cosmic Rays (GCR). In solar eruptive events, such as solar flares and Coronal Mass Ejections (CMEs), protons and other ions can be accelerated to high energies (>30 MeV). The acceleration mechanisms are thought to be related to magnetic reconnection in solar flares (Aschwanden 2012) and the shock waves generated by CME (Cliver et al. 2004).

According to the review by Reames (1999) impulsive events are small, electron rich, tend to be enriched in heavy ions and are $ ^{3}\text{He}-\text{rich} ($$ ^{3}\text{He}/^{4}\text{He} \approx 1 $). They are associated with small magnetic loops in the lower corona with heights smaller than $10^4$ km and have ionization states typical of solar flare temperatures, i.e., $10^7$ K (Mason 2007). The emission is generated at low altitudes, and fast drift (type III) radio emission, which reflects the electron escape into the interplanetary medium (Klein and Trottet 2001), is observed as well.

Gradual SEP events, in contrast, are believed to originate high in the corona from a CME-driven shock whose extent is consistent with the observations of such SEP over relatively wide longitudinal ranges (Rouillard et al. 2011; Dresing et al. 2014). The association of CMEs with type II emission (see Chaps. 1 and 2 for details) also confirms the conclusion that gradual SEP events are accelerated by large CME-driven shocks (Reames 1999). These events produce some of the highest intensity events observed at Earth at energies up to several GeV (Mewaldt et al. 2012; Kühl et al. 2016).

While this two class picture appears to describe the origin of <50 MeV ions with some success, the origin of SEPs that cause GLEs is less studied due to the lack of detailed observations (Moraal and McCracken 2012). It was found that the GLE spectra tend to be slightly harder than non-GLE spectra and that they are consistent with double power laws (Mewaldt et al. 2012; Kühl et al. 2016). Also, the authors found that the composition of GLEs tends to have higher Fe/O ratios, enrichments in $ ^{3}\text{He} $ (Wiedenbeck 2013) and highly-ionized charge-states of heavy elements such as Fe. This lead to the conclusion that GLE ions may be accelerated by CME-driven shocks, with quasi-perpendicular shock geometry and the presence of suprathermal ions from previous flares playing a key role (Tylka et al. 2005). Moraal and Caballero-Lopez (2014) found different scenarios for different GLEs
having a prompt component from the impulsive phase and a gradual one indicating
shock acceleration. To gain further insight the interplanetary and magnetospheric
transport of high energy charged particles needs modeling (Bieber et al. 2004).
While a Fokker-Planck equation well describes the interplanetary transport of near
relativistic electrons (Dröge 2003; Dröge et al. 2010) the processes that need to be
included at high energies have not been fully explored (Dalla et al. 2013, 2015).

Gradual particle events tend to be larger, typically associated with large 10⁵ km X-ray emitting structures, last much longer than impulsive events, are proton-rich without significant enrichment in ³He/⁴He, electron poor, and have elemental abundances and charge states representative of solar coronal or solar wind material and temperatures, i.e., 10⁶ K (Klecker 2013).

However, the relative roles of both components and how we can discriminate
them remains a key problem in solar and solar-terrestrial physics, especially
regarding the diverse interest in GLEs. There is a practical interest in GLEs owing
to their significance for space weather. Solar cosmic rays can damage spacecraft
electronic components and are a significant radiation hazard to astronauts. To
quantify these risks, the full particle distribution in energy and pitch angle as a
function of time needs to be determined from the NM observations. A method to
derive the “physical quantities” is based on forward-modeling of SEP transport
(Sects. 10.3.1 and 10.3.2) in interplanetary space and the Earth’s magnetosphere (see
also Chaps. 4 and 5) by utilizing a power law spectrum in rigidity at the injection
point of the Sun (see Chap. 3) and the response/yield function of the NM (see
Sect. 10.3.2 and Chap. 6). The forward modeling is utilized in Sect. 10.4 to derive an
inversion methodology that is applied to observations in Sect. 10.5. To validate our
model chain, results are compared to spacecraft measurements that are described in
Sect. 10.2.

## 10.2 Space and Ground Based Measurements of GLEs

SEP events, where protons are accelerated to energies above 500 MeV, occur a few
times per solar cycle. Protons with such energies penetrate the Earth's atmosphere
and produce secondary particle showers which can increase the intensities recorded
by NMs on the ground. Such intense SEP events are also known as Ground
Level Events or GLEs. Initially designed by Simpson (1948), NMs are used for
precise monitoring of spectral and directional variations in the cosmic-ray flux. The
detection of a GLE event by an NM on average occurs a few times per solar cycle.
The first GLE was observed on February 28, 1942 and the first GLE observed by
NMs was the one on February 23, 1956 (see gle.oulu.fi).

Since 1942 until the end of 2015 a total of about 70 GLEs have been observed, i.e. ~one GLE per year. Each GLE has its typical characteristics (amplitude, spectrum, duration, spatial distribution of flux, etc.). During a GLE the measurements of the ground-based NM network show an increase in the count rate within typically a few minutes and decreasing intensities to background levels within hours. In some
cases, GLEs show a double-peaked time structure, with an initial fast rise and an anisotropic particle population—called “prompt component”—followed by a more gradual and less anisotropic “delayed component” (McCracken et al. 2012).

SEP events causing GLEs are restricted to events that accelerate ions to energies above ~500 MeV/nucleon. Experiments on board spacecraft located close to Earth such as the Geostationary Operational Environmental Satellites (GOES), the Solar and Heliospheric Observatory (SOHO), the Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics (PAMELA) and the Alpha Magnetic Spectrometer 02 (AMS-02) aboard the International Space Station (ISS) can observe protons below and above this threshold in differential energy channels. These data serve to fill an observational gap between the few hundred MeV and less, typically monitored by other satellite instruments on board e.g. the Advanced Composition Explorer (ACE) (Gold et al. 1998) and Wind (Lin et al. 1995). This latter spacecraft allows us to determine the chemical and isotope abundance (Mewaldt et al. 2007; Nitta et al. 2015), the charge states (Mason et al. 1995; Kartavykh et al. 2007) and to study the lepton component i.e. electrons (Dröge 2000; Dresing et al. 2012, 2014; Agueda et al. 2014) during SEP events. As discussed in Chap. 4 it is important to know not only the intensity-time profiles but also the evolution of the PAD. The pitch angle is defined as the angle between the magnetic field and the velocity vector $\beta$. A particle instrument with a limited opening angle and a pointing direction $\mathbf{n}$ in interplanetary space is sensitive to a pitch angle range around $\cos(\vartheta) = -\frac{B\cdot n}{B}$ with $\mathbf{B}$ the interplanetary magnetic field vector. Spacecraft like ACE, Wind, and the International Monitoring Platform (IMP) 8 outside the Earth’s magnetosphere provide us with measurements of the interplanetary magnetic field vector $\mathbf{B}$ near the first Lagrangian point L1. Contrary to IMP 8 neither GOES nor SOHO measures $\mathbf{B}$. The OMNI dataset<fnref n="1" /> maintained by the US National Space Science Data Center provides us with the magnetic field, and plasma data sets from ACE, Wind and the IMP-8 shifted to the Earth’s bow shock nose. For an observer in the magnetosphere, (PAMELA and AMS-02) the viewing direction $\mathbf{n}$ is the asymptotic direction that has to be calculated as described in Chap. 5.

### 10.2.1 $dE/dx-dE/dx$-Method

This method is utilized by EPHIN aboard SOHO and is based on the energy loss in two detectors. Figure 10.1 right illustrates the measurement principle of the EPHIN instrument (Müller-Mellin et al. 1995). Plotting the energy loss in two

<footnote n="1">¹http://omniweb.gsfc.nasa.gov/owmin.html.</footnote>
<figure><img src="image_85.png" /><figcaption>**Fig. 10.1** The left and right panels display a sketch of the Electron Proton Helium INstrument (EPHIN) aboard SOHO (Müller-Mellin et al. 1995) and the $dE/dx-dE/dx$ measurement principle applied to it, respectively</figcaption></figure>

adjacent detectors against each other, the mean energy losses of H and He follow the characteristic tracks. This is used to identify the particle species in certain areas of the two-dimensional pulse height plane. In these regions, the energy loss in both detectors is attributed to the incoming energy of the particle. Both Goddard Medium Energy (GME) experiment and the Medium Energy Detector (MED) (Meyer and Evenson 1978) utilize this technique.

EPHIN is a multi-element array of solid state detectors with anticoincidence to measure energy spectra of electrons in the range 250 keV to >8.7 MeV, and of hydrogen and helium isotopes in the range 4 MeV/n to >53 MeV/n. The instrument is sketched in Fig. 10.1 left and consists of a stack of silicon solid state detectors (A-F) surrounded by an anticoincidence (scintillator, G). The method to derive energy spectra for penetrating particles with energies above 50 MeV/nucleon is described in detail by Kühl et al. (2015). Since relativistic protons and electrons have the same energy loss $dE/dx$ in matter, electrons with energies above 10 MeV are leading to too high fluxes at energies above ≈ 700 MeV for protons when utilizing the $dE/dx - dE/dx$-method.

The second instrument utilizing this method is the GME aboard IMP-8 that was launched by NASA in 1973 into a 35 $R_E$ geocentric orbit with a 12 days period. The spacecraft was in the solar wind for 7–8 days of every 12-day orbit, where it measured the magnetic fields, plasma, and energetic charged particles (e.g., cosmic rays). The spacecraft spin axis was normal to the ecliptic plane, and the spin rate was 23 rpm. PAD information was obtained in eight angular sectors.

In contrast to EPHIN, the MED design consisted of three pulse-height analyzed CsI detectors with thicknesses of 1 mm, 2 cm and 1 mm, respectively, and a cylindrical plastic scintillator anticoincidence shield. Penetrating particles have an energy above 80 MeV and a geometry factor of nearly 5.0 cm$^2$sr (for details see http://spdf.gsfc.nasa.gov/imp8_GME/GME_instrument.html).
### 10.2.2 $dE/dx$ - C

This method is used for penetrating particles. This method is based on charged particles completely penetrating the Solid State Detector (SSD) and a Cerenkov detector (C), which is placed underneath. If the particles penetrate C faster than with the speed of light in this medium, i.e., with a speed of $\beta > \frac{1}{n}$ where n is the refractive index of the material, they produce a measurable light flash (Cerenkov radiation). Plotting the Cerenkov detector signal against the energy-loss by ionization, $\Delta E$ in the SSD results in characteristic curves for protons and helium, clearly separated with their slopes depending on particle speed. Thus, the method allows an identification of the penetrating particles and gives their energy above a threshold speed. This method is utilized by the HEPAD assembly (Rinehart 1978) aboard a series of GOES satellites in geosynchronous orbit maintained by National Oceanic and Atmospheric Administration (NOAA), which provides differential fluxes in three channels between 350 MeV and 700 MeV, and integral flux above 700 MeV (Sauer 1993). HEPAD consists of two 500 µm thick silicon detectors and a silica radiator. The two surface barrier silicon detectors have an area of 3 cm<fnref n="2" /> and define an effective acceptance aperture of ~24° half-angle.

Charged particles, which completely penetrate the two semiconductor detectors, can be studied with the help of the Cerenkov-detector (McDonald 1956; Linsley 1955), which is placed underneath the semiconductor detectors. The use of a different detector is necessary since, with energy-loss-measurements alone, the energy determination is only possible in a narrow energy range above the maximum energy for stopping particles. The official archive for GOES Energetic Particle data, including proton flux data, can be found at the National Geophysical Data Center.<fnref n="2" />

### 10.2.3 Magnet Spectrometer

In addition we utilize data from the Russian-Italian PAMELA mission (Picozza et al. 2007) and the AMS-02 aboard the ISS (Aguilar et al. 2013; Kounine 2012). PAMELA as well as AMS-02 are magnetic spectrometers in Low Earth Orbit, providing extremely high-quality observations of electrons and ions, which we have used for validation and cross-calibration purposes. Because the orbits of PAMELA and AMS are within the Earth's magnetic field, the two spacecraft do not have a 100% duty cycle for observing low-energy cosmic rays (Adriani et al. 2011). Therefore our approach is to use their observations to validate measurements e.g. from EPHIN.

<footnote n="2">²http://www.ngdc.noaa.gov/stp/satellite/goes/.</footnote>
## 10.3 Forward Modeling from the Sun to the Observer at Ground

The project HESPERIA gathered experts in different fields to take into account interconnections between the solar, heliospheric and NM communities and to advance our knowledge of GLEs further. The group developed a model chain allowing to infer the solar release time profile of relativistic SEPs and their interplanetary transport parameters directly from NM observations. The process chain starts with the propagation of SEPs from the Sun to the Earth. Utilizing the local interplanetary field in Geocentric Solar Ecliptic (GSE) coordinates, the PAD outside the magnetosphere can be converted to an angular distribution in GSE coordinates. Computing the energetic particle transport in the geomagnetic field and taking into account the NM yield function the count rate variation of all NMs that have measured the event can be predicted. Although the physics of the underlying processes is discussed in detail in Chap. 4 (interplanetary transport), Chap. 5 (transport through the magnetosphere) and Chap. 6 (ground-based measurements by NMs) we will show in a compact way how the different models need to be interlinked.

### 10.3.1 Interplanetary Particle Transport: From the Sun to the Magnetosphere

Here we give a summary of Chap. 4 in which we model the particle transport from the Sun ($r_i = 2 r_{Sun}$) in an unperturbed solar wind with constant velocity $v$. The Interplanetary Magnetic Field (IMF) can be described as a smooth average field, represented by an Archimedean spiral, with superposed magnetic fluctuations. The quantitative treatment of the evolution of the particle's phase space density, $f(z, \mu, t)$, can be described by the focused transport equation (Roelof 1969):

$$ \frac{\partial f}{\partial t} + \beta c \mu \frac{\partial f}{\partial z} + \frac{1 - \mu^2}{2L} \beta c \frac{\partial f}{\partial \mu} - \frac{\partial}{\partial \mu} \left( D_{\mu\mu} \frac{\partial f}{\partial \mu} \right) = q(z, \mu, t) \quad (10.1) $$

Here $z$ is the distance along the magnetic field line that depends on the solar wind velocity $v$, $\mu = \cos \alpha$, is the cosine of the particle pitch angle, $\alpha$, and $t$ is the time. The systematic force is characterized by the focusing length, $L(z) = -B(z)/(\partial B/\partial z)$, in the diverging magnetic field $B$, while the stochastic forces are described by the pitch-angle diffusion coefficient $D_{\mu\mu}(\mu)$. As discussed in detail in Chap. 4 the pitch-angle diffusion coefficient has the same form as in Agueda et al. (2008).

Another approximation was introduced by Hasselmann and Wibberenz (1970). If we take the particle radial mean free path, $\lambda_r$, to be spatially constant, then the mean free path parallel to the IMF line is given by $\lambda_{||} = \lambda_r \sec^2 \psi$, where $\psi$ is the angle
between the field line and the radial direction. As shown in Chap. 4 it is sufficient to compute a database of “Green’s-functions” for particles moving with the speed of light ($\beta = 1$). The database consists of 30 different values of the radial mean free path, from $\lambda_r = 0.1$ to $\lambda_r = 2.0$ AU and for solar wind velocities ranging from 300 to 700 km/s. The intensity spectra at the source are given by the solar spectrum that is typically a power-law with $N(R) \propto R^{-\gamma}$, where $\gamma$ is the spectral index of the solar source. Results of this modeling have been discussed in Chap. 4 and (Agueda 2008). There the results could be directly compared to spacecraft measurements. However, observations in the magnetosphere and ground-based measurements need to take into account the geomagnetic filter as described in Chap. 5 and the shielding by the atmosphere and the specific response of the NM (Chap. 6).

> If we assume for each energy interval a number of delta injections in time at the source region the phase space density $f(z, \mu, t)$ can be calculated for any time $t$ at 1 AU corresponding to a distance $z(v)$ depending on the solar wind velocity $v$ for every $\mu$ if the mean free path $\lambda_r$, the solar wind speed $v$ and the direction of the interplanetary magnetic field are known. Acceleration theories (see Chap. 3), as well as measurements, suggests that the source spectrum is given by a power law described by the spectral index $\gamma$. Thus three free parameters describe the temporal evolution of a SEP event caused by a single $\delta$-injection at the Sun.

### 10.3.2 From the Interplanetary Particle Distribution to Neutron Monitor Measurements - Magneto- and Atmospheric Transport of Charged Energetic Particles

The transport of cosmic ray particles in the geomagnetic field and in the Earth’s atmosphere are described in detail in the Chaps. 5 and 6. Here we refer to the objects that are relevant to the investigations of NM data during solar cosmic ray events. The asymptotic viewing directions for each NM station are often computed only for primary cosmic ray particles that penetrate into the Earth’s atmosphere from a vertical direction. The contribution of obliquely incident particles is often neglected in GLE analysis as their contribution to the counting rate is assumed to be small because of the soft spectrum of the solar cosmic rays.

The minimum rigidity that a charged energetic particle must have to reach a location within the magnetosphere and from a given direction of incidence is expressed by the geomagnetic cutoff rigidity $R_C$. This cutoff rigidity varies from a minimum at the magnetic poles ($R_C \approx 0$ GV) to a value of about 15 GV in equatorial regions. The asymptotic direction for a NM, for a given direction of incidence into the Earth’s atmosphere above the location of the NM and a selected
particle rigidity, is defined as the direction of motion of the primary cosmic ray particle before penetrating the Earth's magnetosphere. In contrast to previous work by e.g. Bieber et al. (2004), neither the shape of the PAD nor the incoming direction is a free parameter in the High Energy Solar Particle Events forecasting and Analysis (HESPERIA) approach. The PAD of SEPs outside the Earth's magnetosphere is calculated as described in Sect. 10.3.1. To assign for each rigidity the pitch angle range of a NM we need to know the Interplanetary Magnetic Field (IMF) direction close to the Earth magnetopause and the asymptotic direction in GSE coordinates. In HESPERIA, we make use of the 5-min averaged OMNI data to prepare a set of magnetic field directions outside the magnetosphere, to be used as the directions of the symmetry axis of the directional relativistic proton distribution. The field to which the distribution tends to become gyro-tropic is not the momentarily measured field at the particle position but rather the field it averages over. Although the timescales of the gyro motion are relatively small, the spatial ones are not. A first order estimation, $r_L/(u_{SW}t) \sim 1$, leads to $t = 1670$ s for the averaging time of the fluctuating field for a 1-GV particle in a 5 nT field and a solar wind speed of 400 km/s. However, it should be noted that particles would need a rigidity-dependent averaging time of the field direction. As an example, Fig. 10.2 displays on the left the rigidity-dependent vertical asymptotic directions for a stationary observer in Kiel from 1:30 to 3:00 UT on May 17, 2012. The colored triangles show the direction of the interplanetary field for each period. From these two directions, the pitch angle for each rigidity and time can be calculated. The right panel of Fig. 10.2 displays the corresponding results.

<figure><img src="image_86.png" /><figcaption>**Fig. 10.2 Left:** Calculations of the asymptotic direction at different rigidities during the onset of the May 17, 2012 SEP event as function of time for an observer at Kiel. The triangles indicate the direction of the interplanetary magnetic field. The right panel shows the pitch angle coverage for the same period. For details see text</figcaption></figure>
If we assume that only protons are accelerated in a SEP event to sufficient rigidities $R$ that NMs can measure, then the count rate $N$ caused by solar cosmic ray protons at the NM station $x$ and at time $t$ can be expressed by:

$$N_x(t) = \int_{R_c}^{R_u} I(R, \alpha, t) \cdot S_x(R, h) \cdot dR, \qquad (10.2)$$

where $R_c$ and $R_u$ are the cutoff and upper ($R_u = 20$ GV) rigidities, $\alpha$ the pitch angle of the incoming proton, $S_x(R, h)$ the yield function of the NM station $x$ at the atmospheric depth, $h$ (for details see Chap. 6). Here Eq. (10.2) is approximated by the sum of discrete values covering the full rigidity and pitch angle range. The upper value of $R_u = 20$ GV was chosen because there is no observational evidence for SEP protons with rigidities above 20 GV.

Given the IMF vector **B** in GSE coordinates and the asymptotic directions **n** in the same coordinate system, we compute the pitch angle, $\alpha$ coverage, for each NM. Together with the yield function described in Chap. 6 the count rate increase for each NM is computed.

### 10.3.3 Combined Greens-Function

As part of HESPERIA, the two models described in the previous two sections were combined. The first step is to compute the additional count rate, $N$, caused by SEPs for a given NM station with cutoff rigidity $R_c$ at the time $t_1$, based on the computed rigidity and cosine of the pitch angle ($\mu$) dependent intensity $I_1(\mu, t_1, r)$ of relativistic protons at time $t_1$. This spectrum $I_1$ results from an energy spectrum $I_0 \propto R^{-\gamma}$ injected as a $\delta$-function at a time $t_0$ before $t_1$. The two parameters determining the injection function are the power-law index $\gamma$ and the number of particles injected at a rigidity $R = 1$ GV. The particle transport in the IMF is described by Eq. (10.1) using $\lambda_0$ as the only free parameter. The path length of the particles is determined by the solar wind speed measured at 1 AU (see Sect. 10.3.1). The IMF outside the Earth's magnetosphere together with the asymptotic direction calculated for $t_1$ gives the pitch angle coverage for the NM. Within HESPERIA the Green's function $G_\beta(t, z, \mu, v)$ are derived from the ones $G_c(t, z, \mu, c)$ that are available from SEPServer (http://server.sepserver.eu) for relativistic protons ($\beta = 1$). For near relativistic particles the energy loss in the IMF can be neglected and therefore the Green's function for a lower velocity $\beta$ is given by:

$$G_{\beta}\left(\frac{1}{\beta}t, z, \mu, \beta\right) = \beta G_{c}(t, z, \mu, c). \qquad (10.3)$$
We consider near relativistic protons with rigidities from 0.5 to 20 GV (kinetic energies from 0.12 to 19 GeV) and following a power law in rigidity $I(R) = I_0 \cdot R^{-\gamma}$ with the spectral index $\gamma$.

The proton mean-free path is taken to increase with rigidity following the standard model. The resulting pitch angle dependent energy spectra at 1 AU are transported through the magnetosphere and atmosphere as described above.

## 10.4 Inversion Methodology

The previous two sections are describing the forward modeling of energetic particles from the Sun through the inner helio-, magneto- and atmosphere. Up to now, the inversion has only been attempted in two separate steps:

* Inversion of NM data: The standard analysis of GLE events is based on solving an inverse problem where data by the worldwide network of NMs is used to determine the spectral and angular characteristics of SEPs near Earth but outside the magnetosphere causing the GLE.

* Inversion of spacecraft data: The injection time profile of SEPs at the Sun and the characteristics (mean free path) of their transport in the interplanetary medium is inverted from 1 AU spacecraft measurements.

In what follows we review first the two inversion approaches used so far and then describe the HESPERIA approach.

### 10.4.1 Inversion of Spacecraft Data to the Sun

Numerical simulations of the propagation of SEPs along the IMF are a useful tool to understand SEP events and their sources. We currently have a good theoretical understanding of the transport processes that affect SEPs in the interplanetary medium (see Chap. 4). In Sect. 10.3.1 we discussed a model that simulates the processes undergone by SEPs during their propagation from their source to the observer with the critical parameters summarized in the gray box on page 186.

The approach introduced by Agueda et al. (2008) is to utilize the computed response of the “system” to an impulsive (delta) injection at the Sun, i.e. the Green’s function of particle transport. A convolution of some delta injections allows us to compute different pitch angle dependent proton intensity time profiles that are used as input to the second step in the chain (see Sect. 10.3.2).

Given a system impulse response, $g(t)$, and the input injection profile $q(t)$, the output, $I(t)$, is the convolution of $g(t)$ and $q(t)$:

$$I(t) = \int_{0}^{\infty} g(t')q(t-t')dt' \qquad (10.4)$$
<figure><img src="image_87.png" /><figcaption>**Fig. 10.3** Convolution of different δ-injections close to the Sun on April 15, 2001 (GLE60) from 13:45 h to 14:15 h (left panel) with a set of δ injection γ = 5.5 at the Sun, and a mean free path λ = 0.1 AU in interplanetary space leads to intensity increases at the Apatity, Athens, McMurdo and South Pole NM (right panel)</figcaption></figure>

i.e., *I(t)* is the sum of responses resulting from a series of impulses at the Sun, weighted and shifted in time according to *q(t)*. For simplicity, here we assumed that the response is only time-dependent, but the same holds if one includes the energetic and directional dependence in the Green's function.

Particle intensities measured in the heliosphere as a function of time, energy and direction, are obtained as a temporal convolution of the source function (particle injection profile) and the Green's function of particle transport at the spacecraft location. The origin of SEP events can be unfolded by solving the inverse problem (deconvolving the in-situ measurements). The measurements are used to infer the values of the model parameters. It is a deductive approach, and it has the advantage that a systematic exploration of the parameters' space is possible. A simulated GLE time profile with a proton injection from 13:45 h to 14:15 h at the Sun as displayed in the left panel of Fig 10.3 leads to a proton event observable close to Earth by intensity increases at Apatity, Athens, McMurdo and the South Pole displayed in the right panel. A convolution of a δ-injection with γ = 5.5, and a mean free path λ₀ = 0.1 AU has been used.

Let's consider an arbitrary function $q(t)$—to be determined—that represents the injection function of SEPs close to the Sun. The modeled directional intensities, $M_j^k$, resulting from a series of impulse solar injections can be written as

$$M_j^k(t; \lambda) = \int_{T_1}^{T_2} dt' g_j^k(t, t'; \lambda) q(t'), \quad (10.5)$$

where $g_j^k(t, t'; \lambda)$ represents the impulse response in a given direction $j$ and energy interval $k$, at a given time $t$, when particle injection took place at time $t'$ and assuming an interplanetary mean free path $\lambda$. The duration of the injection function, $t' \in [T_1, T_2]$, is determined by the SEP event time interval selected for fitting, $t \in [t_1, t_2]$, that is, $T_1 = t_1 - \Delta t$ and $T_2 = t_2 - \Delta t$, where $\Delta t$ is the transit time of the first arriving particles at the spacecraft location, assuming a given value of the scattering
mean free path. The number of time points in the event time interval selected for fitting is equal to $n_t = (t_2 - t_1)/\delta t + 2$, where $\delta t$ is the time resolution of the data.

Taking discrete values of time, we have

$$M_j^k(t_h; \lambda) = \sum_{l=1}^{n_t} g_j^k(t_h, t'_l; \lambda) q(t'_l) \quad (10.6)$$

where $j = 1, 2, \dots, n_s$ are numbers representing the directions (sectors or bins) observed by the telescope, $k = 1, 2, \dots, n_c$ numbers for the energy channels and $h = 1, 2, \dots, n_t$ numbers of the time intervals.

Equation (10.6) can be written as

$$M_j^k(t_h; \lambda) \equiv M_i(\lambda) = \sum_{l=1}^{n_t} g_{il}(\lambda) q_l = (\mathbf{g} \cdot \mathbf{q})_i \quad (10.7)$$

where $i = k + (j - 1) \cdot n_t = 1, 2, \dots, n_T$ numbers the total number of observational points and $n_T = n_t n_s$ gives the total number of observational points in all sectors; $\mathbf{g}$ is an $n_T \times n_t$ matrix with $(\mathbf{g})_{il} = g_{il}(\lambda)$.

In matrix form,

$$\begin{pmatrix} M_1 \\ M_2 \\ \vdots \\ M_{n_T} \end{pmatrix} = \begin{pmatrix} g_{11} & g_{12} & \cdots & g_{1n_t} \\ g_{21} & g_{22} & \cdots & g_{2n_t} \\ \vdots & \vdots & \ddots & \vdots \\ g_{n_T 1} & g_{n_T 2} & \cdots & g_{n_T n_t} \end{pmatrix} \cdot \begin{pmatrix} q_1 \\ q_2 \\ \vdots \\ q_{n_t} \end{pmatrix} \quad (10.8)$$

The goal is to compare the modeled intensities with the observations. Let $J_i$ be the observations (background subtracted). We want to derive the $n_t$-vector $\mathbf{q}$ that minimizes the length of the vector $\mathbf{J} - \mathbf{M}$, that means minimizing the value of

$$\|\mathbf{J} - \mathbf{M}\| = \|\mathbf{J} - \mathbf{g} \cdot \mathbf{q}\|, \quad (10.9)$$

subject to the constraint that $q_l \ge 0 \quad \forall l = 1, 2, \dots, n_t$. Thus, the best-fit $\mathbf{q} = (q_1, q_2, \dots, q_{n_t})$ corresponds to a combination of delta-function injection amplitudes. To obtain the best-fit values, we use the NNLS method developed by Lawson and Hanson (1974), which always converges to a solution.

Note, that if each energy channel is fitted separately, the total number of observational points considered in the fit is $n_t n_s$ and the total number of independent fitting parameters (injection amplitudes) is $n_t$. Since $n_t n_s \gg n_t$, the number of degrees of freedom is clearly much larger than the number of model parameters and the inversion problem is well constrained. If one instead neglects the directional information in the data and uses the modeled omnidirectional intensities to fit the omnidirectional event, then the number of observational points and the number of independent fitting parameters is equal, and the problem is not well constrained (multiple injections and transport scenarios can provide an explanation for the data).
The goodness of fit describes how well the model predictions fit a set of observations. One way to evaluate the goodness of the fit, in case the measurement errors are known, is to construct a weighted sum of squared residuals (for details see Agueda 2008). The $\chi^2$ estimator does not work very well for SEP events because during impulsive events the maximum intensities can be several orders of magnitude higher than the intensities observed during the decay phase, thus emphasizing the peak period. Therefore a better goodness-of-fit estimator is provided by the sum of the squared logarithmic differences between the observational and the modeled data. This estimator gives an equal weight of all relative residuals instead of just emphasizing the goodness of fit at the time of maximum. By evaluating the goodness of the fit under different interplanetary transport conditions (different values of $\lambda$), one can objectively discern the “best fit” scenario ($\lambda$-value and associated injection profile) by minimizing the values of the goodness-of-fit estimator.

### 10.4.2 Inversion of NM Data to the Border of the Earth's Magnetosphere

In general, the standard analysis of GLE events is based on an inverse problem, where data by the worldwide network of NMs is used to determine the spectral and angular characteristics of SEPs near Earth causing GLEs for selected times (Shea and Smart 1982; Mishev et al. 2014). Analysis of the characteristics of the primary solar particles causing GLEs from ground-based data records is a serious challenge (Bütikofer and Flückiger 2013). Data from stations with different cutoff rigidities (geomagnetic latitudes) provide information necessary to determine the spectral characteristics. Responses of stations over a wide range of geographical locations are required to determine the axis of symmetry. Data are fitted to directional distributions that are rotationally symmetric about one direction in space. This, in principle, is close but not exactly the direction of the instantaneous magnetic field measured close to the Earth but outside the Earth's magnetosphere (Bieber et al. 2013). Therefore, axis-symmetry is assumed, but the direction of the axis of symmetry is optimized to fit the data.

The PADs of relativistic solar protons in space are assumed to follow a given functionality. A variety of functions has been used in the literature. These include a linear form, an exponential plus a constant, a parabola, two Gaussians, and two exponentials plus constant. The latter three are expected to provide better fits to bidirectional fluxes, if present (Bieber et al. 2013; Mishev et al. 2014).

Similarly, the spectra of relativistic solar protons during a GLE are assumed to follow a power law in rigidity, or energy with extensions that describe the softening of the spectrum at higher energies by multiplying a power law with an exponential cutoff.
The parameters of the rigidity spectrum, the PAD and the direction of symmetry are determined by minimizing the sum of squared differences between the modeled relative change in the NM count rate at the station x

($\frac{\Delta N_x}{N_x}$) mod. and

the corresponding observed relative count rate change

($\frac{\Delta N_x}{N_x}$) obs.

$$F = \sum_{x=1}^{m} \left[ \left( \frac{\Delta N_x}{N_x} \right)_{\text{mod.}} - \left( \frac{\Delta N_x}{N_x} \right)_{\text{obs.}} \right]^2 \quad (10.10)$$

for the set of selected NM data. The Levenberg-Marquardt algorithm (LMA) (Marquardt 1963) provides a numerical solution to the problem of minimizing a nonlinear function over the space of parameters of the function. The goodness of the fit is commonly expressed by a weighted sum of squared residuals or by computing a correlation coefficient, $\rho$ between measurements and the model.

### 10.4.3 The HESPERIA Approach

For the first time, models of the transport of SEPs in the interplanetary medium, the Earth's magnetosphere and atmosphere, and the response of NMs are linked to each other. The first goal is to compute the expected additional count rate, $N$, caused by SEP for a given NM station, based on the intensities of the primary cosmic rays near Earth, $I(R, \alpha, t)$, as function of rigidity, $R$, pitch angle, $\alpha$, and time, $t$. As detailed above the interplanetary transport is described by Green's functions representing characteristics of the SEP interplanetary transport conditions. To ascribe the magnetospheric transport, the magnetic field direction outside the magnetosphere is computed from interplanetary measurements, and the asymptotic directions are calculated utilizing the PLANETOCOSMICS code in GSE coordinates. The latter can be found at the HESPERIA webpage during each GLE in the past and for the NMs of the worldwide network.

Since the parameters are the same than the ones used by Agueda et al. (2008) the differences lie in the indirect measurement of the pitch angle dependent count rate and the integration over wide energy ranges: For all NM stations, the counting rate increases can then be computed for a series of delta injections from the Sun and for selected times. Note that the hardness of the source spectra described by $\gamma$ is assumed to be constant in the HESPERIA approach for all $\delta$-injections. The amplitudes of the source components, for a given scenario, are inferred by fitting the NM observations with the modeled NM counting rate increases. The NNLS algorithm described in Sect. 10.4.1 is used to determine the best set of parameters. Regularized inversion approaches will be explored for refinements. The goodness of the fit will be evaluated by computing a weighted sum of the squared residuals.
The result of the inversion problem is a detailed time profile of the injection process at the Sun. The shape of this profile is presumably determined by details of the acceleration process and possible transport processes in the corona (see Chap. 3 for details).

In contrast to a classical approach with a total number of at least 10 fit parameters to derive the injection function at the Sun (see Sect. 10.4.2) the HESPERIA approach relies on several well-documented assumptions (see Chaps. 3, 4, and 5) reducing the number of free parameters. Assuming for each rigidity a number of delta injection at the source region the phase space density $f(z, \mu, t)$ close to Earth can be calculated for any given time $t$ corresponding to the distance $z(v)$ that depends on the measured solar wind velocity $v$. The mean free path $\lambda_0$, the number of particles, the injection time profile and its rigidity spectrum described by a power law with index $\gamma$ are utilized to derive $f(z, \mu, t)$. The cosmic ray particle trajectories through the Earth's magnetosphere are computed for selected rigidities $R_c, R_c + 0.1$ GV ... 20 GV as a function of GSE coordinates of the NM station, for a selected time. Applying the Tsyganenko 1989 (Tsy89) model and the yield function for each NM we can ascribe the SEP time profiles for the NM network by a minimum of free parameters.

## 10.5 Results and Validation

Since the HESPERIA approach depends on the knowledge of the IMF outside the Earth's magnetosphere that is provided by OMNI since November 4, 1973, our investigation is restricted to GLEs with numbers larger than number 26 occurring on April 30, 1976. To validate the results, pitch angle dependent intensity time profiles of protons with energies above 500 MeV are needed. Here we utilize SOHO/EPHIN and Wind 3DP measurements. The latter are needed in order to estimate the pitch angle coverage of SOHO/EPHIN.

Since SOHO was launched in December 1995, we restrict our validation to GLEs with numbers $\ge 55$ (GLE55 occurred on November 6, 1997). Kühl et al. (2015, 2016) showed that EPHIN is capable of measuring the proton spectra in the range from 100 MeV to above 1 GeV. Since no magnetic field measurements are available on SOHO, we compare the intensity-time-profiles of near relativistic electrons measured by EPHIN and Wind 3DP with each other. In Fig. 10.4 different colored curves show the pitch angle dependent time-profiles for Wind 3DP with energies between 230 and 646 keV. Three thick black lines show the EPHIN measurements with energies between 0.25 and 0.7 MeV multiplied by different factors (0.05, 0.08, and 0.1 as solid, dashed, and dashed-dotted lines) to take into account cross
<figure><img src="image_88.png" /><figcaption>**Fig. 10.4** *From top to bottom*: Hourly averaged intensity time profiles of 234–646 keV electrons measured by Wind 3DP. The different colors give the pitch angle coverage for the WIND observations shown in the *bottom panels*. Also the *solid, dashed and dashed-dotted lines* represent the 0.25–0.7 MeV EPHIN electron measurements multiplied by a factor 0.05, 0.08, and 0.1 respectively</figcaption></figure>

calibration. Only the dashed-dotted line agrees with all Wind sectors when the flux becomes isotropic. Thus we assume that the curve that agrees best with the dashed-dotted line ascribes the viewing direction of EPHIN. In our example, the purple curve describing the measurements at a pitch angle of 120° fulfills this criteria best.

Figure 10.5 is an extension of Fig. 10.3 that only showed predicted intensity time profiles close to Earth and by different NMs. It compares the results computed with the GLE inversion software with measurements of other facilities as radio telescopes, NMs and particle detectors in space. The parameters used in the prediction were derived by the inversion of the NM data for GLE 60. The middle panel of Fig. 10.5 compares the simulated and measured intensity time profiles of the Apatity, Athens, McMurdo and South Pole NMs with each other showing a good agreement between measurements and the model. The injection profile close to the Sun is shown in the left panel together with the micro wave profile measured by the Radio Solar Telescope Network (RSTN) showing a good agreement between the particle component injected into interplanetary space and towards the lower corona. Note, that the injection profile is in very good agreement with the one derived by Bieber et al. (2004). For details on micro waves and their importance, the reader is referred to Chap. 2.

The middle bottom panel displays the energy spectrum between 100 MeV and 1 GeV predicted by the model (red curve) and measured by EPHIN (black curve). Note, that the prediction was scaled down by a factor of ~20. From the figure, it is evident that the model predicts the same spectral index when taking into account the contribution of electrons at 900 and 700 MeV, respectively (for details see Kühl et al. 2017).
<figure><img src="image_89.png" /><figcaption>**Fig. 10.5** Predicted δ-injections (upper left panel), NM increases (upper right panel) and differential intensity spectrum (lower panel) for GLE60 (see also Fig. 10.3) in comparison to actual measurements (for details see text)</figcaption></figure>

**10.6 Concluding Remarks**

A new approach has been presented here that allows computing the injection function of SEP close to the Sun based on the data of the worldwide NM network during a GLE. This injection function is described by a power law in rigidity $R$ with two parameters that are the intensity at $R_0 = 1$ GV and the spectral index $\gamma$. For the interplanetary transport, we utilize a 1-dimensional model with the mean free path $\lambda_r$ as free parameter. The solar wind speed taken from interplanetary measurements determines the length of the magnetic field line connecting the Earth back to the Sun. The IMF direction is taken from the OMNI data set using appropriate accumulation periods. The transport through the Earth's magnetosphere is computed with the PLANETOCOSMICS software using the Tsy89 model for the outer Earth's magnetic field. The NM yield function describes the passage of the cosmic ray particles through the atmosphere and the detection of the secondary nucleon component by the NMs. Using this model chain GLE 60 could be reproduced very well. A discrepancy between the prediction and the measurements in interplanetary space exists that needs further investigations. Especially the whole number of observed GLEs should be analyzed and validated by using GOES HEPAD data that
go back to the 1970s. Further improvements like utilizing more sophisticated models of the magnetosphere and the interplanetary transport, as well as different NM yield functions, should be taken into account in order to improve the model presented here.

## References

Adriani, O., Barbarino, G.C., Bazilevskaya, G.A., Bellotti, R., Boezio, M., Bogomolov, E.A., Bonechi, L., Bongi, M., Bonvicini, V., Borisov, S., et al.: Astrophys. J. 742, 102 (2011). doi:10.1088/0004-637X/742/2/102

Agueda, N.: Near-relativistic electron events. Monte Carlo simulations of solar injection and interplanetary transport. Ph.D. thesis, Dep. Astronomia i Meteorologia University of Barcelona Martí i Franquès 1 08028 Barcelona (2008)

Agueda, N., Vainio, R., Lario, D., Sanahuja, B.: Astrophys. J. 675, 1601 (2008). doi:10.1086/527527

Aguilar, M., Alberti, G., Alpat, B., Alvino, A., Ambrosi, G., Andeen, K., Anderhub, H., Arruda, L., Azzarello, P., Bachlechner, A., et al.: Phys. Rev. Lett. 110(14), 141102 (2013). doi:10.1103/PhysRevLett.110.141102

Agueda, N., Klein, K.L., Vilmer, N., Rodríguez-Gasén, R., Malandraki, O.E., Papaioannou, A., Subirà, M., Sanahuja, B., Valtonen, E., Dröge, W., et al.: Astron. Astrophys. 570, A5 (2014). doi:10.1051/0004-6361/201423549

Aschwanden, M.J.: Space Sci. Rev. 171, 3 (2012). doi:10.1007/s11214-011-9865-x

Bieber, J.W., Evenson, P., Dröge, W., Pyle, R., Ruffolo, D., Rujiwarodom, M., Tooprakai, P., Khumlumlert, T.: Astrophys. J. Lett. 601, L103 (2004). doi:10.1086/381801

Bieber, J.W., Clem, J., Evenson, P., Pyle, R., Sáiz, A., Ruffolo, D.: Astrophys. J. 771, 92 (2013). doi:10.1088/0004-637X/771/2/92

Bütikofer, R., Flückiger, E.O.: J. Phys. Conf. Ser. 409(1), 012166 (2013). doi:10.1088/1742-6596/409/1/012166

Cliver, E.W., Kahler, S.W., Reames, D.V.: Astrophys. J. 605, 902 (2004). doi:10.1086/382651

Dalla, S., Marsh, M.S., Kelly, J., Laitinen, T.: J. Geophys. Res. (Space Phys.) 118, 5979 (2013). doi:10.1002/jgra.50589

Dalla, S., Marsh, M.S., Laitinen, T.: Astrophys. J. 808, 62 (2015). doi:10.1088/0004-637X/808/1/62

Dresing, R., Gómez-Herrero, N., Klassen, A., Heber, B., Kartavykh, Y., Dröge, W.: Sol. Phys. 281, 281 (2012). doi:10.1007/s11207-012-0049-y

Dresing, N., Gómez-Herrero, R., Heber, B., Klassen, A., Malandraki, O., Dröge, W., Kartavykh, Y.: Astron. Astrophys. 567, A27 (2014). doi:10.1051/0004-6361/201423789

Dröge, W.: Astrophys. J. 537, 1073 (2000). doi:10.1086/309080

Dröge, W.: Astrophys. J. 589, 1027 (2003). doi:10.1086/374812

Dröge, W., Kartavykh, Y.Y., Klecker, B., Kovaltsov, G.A.: Astrophys. J. 709, 912 (2010). doi:10.1088/0004-637X/709/2/912

Gold, R.E., Krimigis, S.M., Hawkins, S.E.I., Haggerty, D.K., Lohr, D.A., Fiore, E., Armstrong, T.P., Holland, G., Lanzerotti, L.J.: Space Sci. Rev. 86(1), 541 (1998)

Hasselmann, K., Wibberenz, G.: Astrophys. J. 162, 1049 (1970)

Kartavykh, Y.Y., Dröge, W., Klecker, B., Mason, G.M., Möbius, E., Popecki, M., Krucker, S.: Astrophys. J. 671, 947 (2007). doi:10.1086/522687
Klecker, B.: J. Phys. Conf. Ser. **409**(1), 012015 (2013). doi:10.1088/1742-6596/409/1/012015

Klein, K.L., Trottet, G.: Space Sci. Rev. **95**, 215 (2001)

Kounine, A.: Int. J. Modern Phys. E **21**(08), 1230005 (2012). doi:10.1142/S0218301312300056.
http://www.worldscientific.com/doi/abs/10.1142/S0218301312300056

Kühl, P., Banjac, S., Dresing, N., Goméz-Herrero, R., Heber, B., Klassen, A., Terasa, C.: Astron. Astrophys. **576**, A120 (2015). doi:10.1051/0004-6361/201424874

Kühl, P., Dresing, N., Heber, B., Klassen, A.: Solar energetic particle events with protons above 500 MeV between 1995 and 2015 measured with SOHO/EPHIN. Sol. Phys. **292**(10), 10 (2017). arXiv:1611.03289. doi:10.1007/s11207-016-1033-8. http://adsabs.harvard.edu/abs/2017SoPh..292...10K

Kühl, P., Dresing, N., Heber, B., Klassen, A.: Sol. Phys. **292**, 10 (2017). doi:10.1007/s11207-016-1033-8

Lawson, C.L., Hanson, R.J.: Solving Least Squares Problems. SIAM, Philadelphia (1974)

Lin, R.P., Anderson, K.A., Ashford, S., Carlson, C., Curtis, D., Ergun, R., Larson, D., McFadden, J., McCarthy, M., Parks, G.K., et al.: Space Sci. Rev. **71**, 125 (1995). doi:10.1007/BF00751328

Linsley, J.: Phys. Rev. **97**, 1292 (1955). doi:10.1103/PhysRev.97.1292

Marquardt, D.W.: J. Soc. Ind. Appl. Math. **11**(2), 431 (1963). doi:10.1137/0111030

Mason, G.M.: Space Sci. Rev. **130**, 231 (2007). doi:10.1007/s11214-007-9156-8

Mason, G.M., Mazur, J.E., Looper, M.D., Mewaldt, R.A.: Astrophys. J. **452**, 901 (1995). doi:10.1086/176358

McCracken, K.G., Moraal, H., Shea, M.A.: Astrophys. J. **761**, 101 (2012). doi:10.1088/0004-637X/761/2/101

McDonald, F.B.: Phys. Rev. **104**, 1723 (1956). doi:10.1103/PhysRev.104.1723

Mewaldt, R.A., Cohen, C.M.S., Mason, G.M., Cummings, A.C., Desai, M.I., Leske, R.A., Raines, J., Stone, E.C., Wiedenbeck, M.E., von Rosenvinge, T.T., et al.: Space Sci. Rev. **130**, 207 (2007). doi:10.1007/s11214-007-9187-1

Mewaldt, R.A., Looper, M.D., Cohen, C.M.S., Haggerty, D.K., Labrador, A.W., Leske, R.A., Mason, G.M., Mazur, J.E., von Rosenvinge, T.T.: Space Sci. Rev. **171**, 97 (2012). doi:10.1007/s11214-012-9884-2

Meyer, P., Evenson, P.: IEEE Trans. Geosci. Electron. **16**, 180 (1978)

Mishev, A.L., Kocharov, L.G., Usoskin, I.G.: J. Geophys. Res. **119**, 670 (2014). doi:10.1002/2013JA019253

Moraal, H., Caballero-Lopez, R.A.: Astrophys. J. **790**, 154 (2014). doi:10.1088/0004-637X/790/2/154

Moraal, H., McCracken, K.G.: Space Sci. Rev. **171**, 85 (2012). doi:10.1007/s11214-011-9742-7

Müller-Mellin, R., Kunow, H., Fleißner, V., Pehlke, E., Rode, E., Röschmann, N., Scharnberg, C., Sierks, H., Rusznyak, P., McKenna-Lawlor, S., et al.: Sol. Phys. **162**(1), 483 (1995)

Nitta, N.V., Mason, G.M., Wang, L., Cohen, C.M.S., Wiedenbeck, M.E.: Astrophys. J. **806**, 235 (2015). doi:10.1088/0004-637X/806/2/235

Picozza, P., Galper, A.M., Castellini, G., Adriani, O., Altamura, F., Ambriola, M., Barbarino, G.C., Basili, A., Bazilevskaja, G.A., Bencardino, R., et al.: arXiv.org (4), 296 (2006)

Reames, D.V.: Space Sci. Rev. **90**, 413 (1999). doi:10.1023/A:1005105831781

Rinehart, M.C.: Nucl. Inst. Methods **154**, 303 (1978). doi:10.1016/0029-554X(78)90414-7

Roelof, E.C.: Propagation of solar cosmic rays in the interplanetary magnetic field. In: Ögelman, H., Wayland, J.R. (Eds.) Lectures in High-Energy Astrophysics, p. 111 (1969). http://adsabs.harvard.edu/abs/1969lhea.conf..111R. Provided by the SAO/NASA Astrophysics Data System

Rouillard, A.P., Odstrecil, D., Sheeley, N.R., Tylka, A., Vourlidas, A., Mason, G., Wu, C.C., Savani, N.P., Wood, B.E., Ng, C.K., et al.: Astrophys. J. **735**, 7 (2011). doi:10.1088/0004-637X/735/1/7

Sauer, H.H.: Int. Cosmic Ray Conf. **3**, 250 (1993)

Shea, M.A., Smart, D.F.: Space Sci. Rev. **32**, 251 (1982). doi:10.1007/BF00225188
Simpson, J.A.: Phys. Rev. **73**, 1389 (1948). doi:10.1103/PhysRev.73.1389

Tylka, A.J., Cohen, C.M.S., Dietrich, W.F., Lee, M.A., Maclennan, C.G., Mewaldt, R.A., Ng, C.K., Reames, D.V.: Astrophys. J. **625**, 474 (2005). doi:10.1086/429384

Wiedenbeck, M.E., Mason, G.M., Cohen, C.M.S., Nitta, N.V., Gómez-Herrero, R., Haggerty, D.K.: Astrophys. J. **762**, 54 (2013). doi:10.1088/0004-637X/762/1/54

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

<figure><img src="image_90.png" /></figure>
# Index

## A

Abundance ratios, 8, 46, 149, 150, 151, 152

Asymptotic directions, 79, 80, 91, 187, 188, 193

Atmospheric effects, 100, 105–106

## B

Bremsstrahlung, 27, 28, 31, 32, 34, 36, 40, 115, 135, 136, 137, 139, 140

## C

Charged particles, 17, 18, 36, 37, 38, 40, 41, 50, 51, 63, 67–69, 79, 80, 81, 89, 181, 183, 184

Charged particle transport, 63–76

Chromosphere, 29, 32, 35, 135, 145, 157, 158

CME-driven shock, 3, 5, 14, 16, 46, 144, 148, 152, 157, 158, 180

Collapsing magnetic trap, 57–58

Combined Greens-function, 188–189

Compressional acceleration, 46, 57–58

Convolution, 189, 190

Coronal mass ejection, 1, 27, 37–42, 45, 63, 96, 117, 180

Coronal shock acceleration, 9, 159–161, 176

Cosmic ray particle transport, 79–93

Cutoff rigidities, 80, 91–93, 192

Cyclotron resonance, 51, 52

Cyclotron waves, 52, 59

## D

δ-injections, 186, 189, 190, 193, 196

Differential intensity spectrum, 196

Diffusive shock acceleration (DSA), 15, 41, 53–58

Diffusive transport, 71–72, 161

Discontinuity, 91, 92, 159

Downstream propagation (DSP) model, 159, 161–164, 169, 170, 173

Dynamo, 79

## E

Earth's atmosphere, 95, 97–100, 103, 105, 106, 181, 186

Earth's Magnetic field, 79–83, 86–88, 90, 92, 97, 100, 107, 184, 196

Earth's magnetosphere, 19, 66, 79–93, 114, 182, 187, 188, 192, 193, 194, 196

Electric field acceleration, 47–49

Electodynamic coupling, 29

Electromagnetic fields, 45, 47, 48, 69, 80

Electromagnetic waves, 33, 52, 69–70

Electron beams, 33–34, 36

Electron intensities, 120, 121, 122

Energetic particles, 1–23, 41, 63–64, 69, 98, 114, 134, 159, 185–189

Energy spectra, 4, 11, 12, 13, 19, 22, 118–120, 149, 150, 167, 168, 169, 173, 183, 189

Eruptive activity, 27–41
F

First adiabatic invariant, 67–69

Focused transport, 70–71, 185

Forecasting tools, 23, 113–130

Forward modeling, 181, 185–189

G

Gamma-ray, 28, 31, 34, 38, 133–151, 157–173

Gamma(γ)-ray emission, 27, 34, 134–138, 140–145, 149, 169

Geomagnetic field, 80, 82, 84, 85, 87, 88, 89, 92, 95, 96, 104, 180, 186

Geomagnetosphere, 79, 80, 81, 82, 83, 87, 88, 89, 91, 92, 93, 96

Gradual SEP events, 2, 3, 5–11, 46, 60, 64, 180

Ground-based measurements, 95–109

Ground level enhancement, 11–12, 97, 114, 124, 158, 179–196

Gyrosynchrotron, 31–33, 39, 115

Gyrosynchrotron radiation, 31–33, 115

H

Hard X-ray, 28, 30, 31, 32, 35, 36, 38, 115, 135, 137, 138, 139, 140, 141, 142, 143, 148, 151, 152

Hazard assessment, 20–22

Heliosphere, 6, 9, 12–16, 38, 40, 64, 71, 190

Historical microwave data, 115–118

Human health, 19–20

Human health effects, 19–20

I

Impulsive flare, 3, 29, 40, 134, 137, 138, 139, 140

Impulsive SEP events, 2, 3, 151

Intensity time-profiles, 13, 171, 194

Interplanetary magnetic field, 3, 14, 64–66, 69, 87, 89, 166, 182, 185, 186, 187

Interplanetary medium, 63–76, 180, 189, 193

Interplanetary particle distribution, 186–188

Interplanetary transport, 72, 77, 181, 185, 192, 193

Inversion methodology, 179–196

Inversion model, 72, 75, 76

Ionosphere, 19, 83, 86, 87, 114

L

Langmuir waves, 33, 34

Large gradual SEP events, 5–11, 60

Long duration γ-ray events, 136, 141–143, 145, 158, 173, 174

Lorentz force, 36, 38, 47, 49, 67, 68, 80–81, 97

M

Magnetic flux rope, 38, 39

Magnetic polarities, 39

Magnetic structuring, 28

Magnetohydrodynamic, 2, 41, 56, 159

Magnetospheric electric currents, 86–87

Mean free path, 57, 71, 74, 99, 102, 162, 185, 186, 189, 190, 191, 194, 196

Microwave flux density, 114, 116, 117, 119, 120, 129

Mitigation procedures, 20, 22–23

Mitigation strategies, 17, 22

Multi-spacecraft observations, 12–14

N

Neutron monitor data, 105, 108

Neutron monitor detector, 100–106

Neutron monitors, 2, 22, 91, 95–109

Neutron monitor station, 96, 97, 98, 100, 104, 106, 107, 108, 109

Neutron monitor yield function, 103–104

Nuclear gamma-ray emission, 27

Numerical techniques, 72

P

Parallel diffusion coefficient, 9

Particle acceleration, 2, 10, 14–15, 27–37, 40–42, 45–60, 115, 129, 152, 169, 173, 174, 175, 176

Particle cascade, 99–100

Particle effects, 18

Particle transport equations, 69–70

Peak proton intensity, 6, 11

Perpendicular diffusion coefficient, 9

Photon energies, 31, 134–136, 141

Photosphere, 28, 29, 36, 38, 39, 41, 145

Pion-decay, 34, 134, 135, 136, 137, 139, 140, 142, 151, 152, 158, 174, 176

Plasma emission, 31, 33–34, 39

Propagation, 13, 19, 54, 63, 64, 68, 74, 79, 89–91, 114, 159–165, 170, 173, 185, 189

Proton acceleration, 158, 167–168, 172–173, 176

Proton differential intensity, 73

Proton fluxes, 116, 117, 146

Proton injection rate, 165, 166, 167, 171, 172
Proton intensity, 5, 6, 11, 12, 116, 117, 120, 171, 189

Proton spectra, 4, 119, 129, 139, 164, 169, 173, 194

Proton transport, 159, 161, 169–170, 173

R

Real-time predictions, 23, 130

Real-time tools, 113

Relativistic gamma-ray events, 133–152

Remote sensing, 16, 22, 180

Reservoir region, 12, 13

Resonant wave acceleration, 50–53

S

SEP effects, 17–19

Shock-accelerated particles, 9, 56, 71, 75, 141

Shock acceleration, 7, 9, 11, 15, 41, 46, 47, 53–57, 152, 158–161, 176, 181

Shock drift acceleration, 40, 53

Shock modelling, 168

Shock passage, 12, 13, 73, 74, 165

Shock surfing acceleration, 53

Shock waves, 2, 7, 27, 28, 36, 37, 40–42, 54, 70, 141, 144–145, 158, 180

Simulation model, 159, 161, 176

Snowstorm effect, 18

Soft X-ray, 7, 29, 30, 31, 114, 129, 140, 141, 142, 143, 146, 147, 149, 151, 152

Solar activity, 27, 42, 66, 87, 100, 106, 124

Solar atmosphere, 28, 31, 32, 34, 40, 115, 134, 137, 141, 144, 145, 151, 152, 158

Solar corona, 16, 28, 33, 36, 37, 39, 42, 45, 53, 143, 181

Solar cosmic ray, 96, 97, 105, 106, 181, 186

Solar cycle 23, 11, 12

Solar energetic particle events, 27–42, 63–64, 72–76, 100, 114, 134, 137, 145–151

Solar energetic particles, 1–23, 45, 124, 134

Solar flare, 1, 2, 3, 7, 16, 22, 29–36, 40, 53, 63, 135, 137, 138, 158, 180

Source plasma temperature, 12

Space weather, 1–23, 96, 118, 130, 176, 181

Spectrometer, 96, 98, 102, 107–108, 134, 182, 184

Stochastic acceleration, 46, 58–60, 158, 174

Suprathermal seed populations, 8

Sustained emission events, 137–141

T

Tsyganenko model, 88, 89

V

Validation, 184, 194–196