Acceleration of solar wind particles by traveling interplanetary shocks

TL;DR

This study models how interplanetary shocks accelerate solar wind particles, revealing that shock speed, magnetic turbulence, and prior solar wind heating significantly influence the energy and intensity of accelerated particles.

Contribution

It introduces a simulation framework combining $ppa$-distribution functions and shock parameters to analyze solar wind particle acceleration, highlighting the roles of shock speed and turbulence.

Findings

Fast shocks produce higher energy particles.

Prior solar wind heating enhances particle acceleration.

Maximum energies limited by shock transit times.

Abstract

The acceleration of thermal solar wind protons at spherical interplanetary shocks driven by coronal mass ejections is investigated. The solar wind velocity distribution is represented using $\kappa$-functions, which are transformed in response to simulated shock transitions in the fixed-frame flow speed, plasma number density, and temperature. These heated solar wind distributions are specified as source spectra at the shock from which particles with sufficient energy can be injected into the diffusive shock acceleration process. It is shown that for shock-accelerated spectra to display the classically expected power-law indices associated with the compression ratio, diffusion length scales must exceed the width of the compression region. The maximum attainable energies of shock-accelerated spectra are found to be limited by the transit times of interplanetary shocks, while spectra may be accelerated to higher energies in the presence of higher levels of magnetic turbulence or at faster-moving shocks. Indeed, simulations suggest fast-moving shocks are more likely to produce very high-energy particles, while strong shocks, associated with harder shock-accelerated spectra, are linked to higher intensities of energetic particles. The prior heating of the solar wind distribution is found to complement shock acceleration in reproducing the intensities of typical energetic storm particle events, especially where injection energies are high. Moreover, simulations of $\sim$0.2 to 1 MeV proton intensities are presented that naturally reproduce the observed flat energy spectra prior to shock passages. Energetic particles accelerated from the solar wind, aided by its prior heating, are shown to contribute substantially to intensities during energetic storm particle events.

Tables

Table 1. Table 1: Reference configuration of the parameters needed to define the fixed-frame flow speed V s ​ w ′ superscript subscript 𝑉 𝑠 𝑤 ′ V_{sw}^{\prime} , number density n 𝑛 n , temperature T 𝑇 T and κ 𝜅 \kappa -function as discussed in Section 2.1 .

$V_{sw}^{\prime }$:	$V_1^{\prime }$	600	km s-1
	$V_{sh}$	2400	km s-1
	$s$	3.0
	$L$	0.005	AU
$n$:	$n_0$	6.0	cm-3
$T$:	$T_e$	$5\times {10}^5$	K
	${\gamma }_c$	3.5/3
$f_{\kappa }$:	$\kappa$	2.25

Table 2. Table 2: Parameters used to fit Eq. 25 to the corresponding model solutions shown in Figure 4 . The λ 0 subscript 𝜆 0 \lambda_{0} -values are constants used to scale the MFPs in Eq. 12 . E tr subscript 𝐸 tr E_{\text{tr}} and E cut subscript 𝐸 cut E_{\text{cut}} are energies at which spectra transition to power-law indices associated with the full shock compression and roll over into exponential decreases, respectively. γ a subscript 𝛾 𝑎 \gamma_{a} is the spectral index for E < E tr 𝐸 subscript 𝐸 tr E<E_{\text{tr}} . Note that E tr subscript 𝐸 tr E_{\text{tr}} -values are omitted for the two largest λ 0 subscript 𝜆 0 \lambda_{0} -values, because γ a = γ b = − 1.25 subscript 𝛾 𝑎 subscript 𝛾 𝑏 1.25 \gamma_{a}=\gamma_{b}=-1.25 , which reduces Eq. 25 to the form of Eq. 23 .

${\lambda }_0$ (AU)	${\gamma }_a$	$E_{\text{tr}}$ (MeV)	$E_{\text{cut}}$ (MeV)
0.005	$-2.5$	1.5	30.0
0.01	$-2.4$	0.5	11.0
0.02	$-2.0$	0.2	5.0
0.035	$-1.25$	0.08	2.5
0.06	$-1.25$	-	1.5
0.1	$-1.25$	-	0.9

Table 3. Table 3: Configurations used throughout Section 5 , in which 1: the shock speed V s ​ h subscript 𝑉 𝑠 ℎ V_{sh} is varied, 2: the compression ratio s 𝑠 s is varied, and 3: V s ​ h subscript 𝑉 𝑠 ℎ V_{sh} and s 𝑠 s are varied together such that V 2 ′ / V 1 ′ superscript subscript 𝑉 2 ′ superscript subscript 𝑉 1 ′ V_{2}^{\prime}/V_{1}^{\prime} remains constant. Here, V 1 ′ = superscript subscript 𝑉 1 ′ absent V_{1}^{\prime}= 600 km s -1 and V 2 ′ = ( V s ​ h ​ ( s − 1 ) + V 1 ′ ) / s superscript subscript 𝑉 2 ′ subscript 𝑉 𝑠 ℎ 𝑠 1 superscript subscript 𝑉 1 ′ 𝑠 V_{2}^{\prime}=(V_{sh}(s-1)+V_{1}^{\prime})/s are the fixed-frame flow speeds, respectively up- and downstream of the shock. V 2 = ( V s ​ h − V 1 ′ ) / s subscript 𝑉 2 subscript 𝑉 𝑠 ℎ superscript subscript 𝑉 1 ′ 𝑠 V_{2}=(V_{sh}-V_{1}^{\prime})/s is the downstream flow speed in the shock frame. The last line contains the reference configuration.

	$V_{sh}$	$s$	$V_2$	$V_2^{\prime }$	$V_2^{\prime }/V_1^{\prime }$
	(km s-1)		(km s-1)	(km s-1)
1:	3000	3	800	2000	3.67
	2800	3	733	2067	3.44
	2600	3	667	1933	3.22
	2280	3	560	1720	2.87
	2200	3	533	1667	2.78
2:	2400	2	900	1500	2.5
	2400	2.2	818	1581	2.64
	2400	2.5	720	1680	2.8
	2400	3.5	514	1886	3.14
	2400	4	450	1950	3.25
3:	3000	2	1200	1800	3
	2800	2.2	1000	1800	3
	2600	2.5	800	1800	3
	2280	3.5	480	1800	3
	2200	4	400	1800	3
ref:	2400	3	600	1800	3.0

Table 4. Table 4: Parameters used to fit Eq. 23 to the solutions of Figure 11 and corresponding to the configurations of shock speed V s ​ h subscript 𝑉 𝑠 ℎ V_{sh} and compression ratio s 𝑠 s introduced in Table 3 . Here, γ ​ ( s ) 𝛾 𝑠 \gamma(s) is the spectral index associated with s 𝑠 s as given by Eq. 24 and E cut subscript 𝐸 cut E_{\text{cut}} is the energy at which spectra roll over into exponential decreases. The last line contains the reference configuration.

	$V_{sh}$ (km s-1)	$s$	$\gamma \left(s\right)$	$E_{\text{cut}}$ (MeV)
1:	3000	3	-1.25	2.0
	2800	3	-1.25	1.9
	2600	3	-1.25	1.8
	2280	3	-1.25	1.4
	2200	3	-1.25	1.35
2:	2400	2	-2.0	1.4
	2400	2.2	-1.75	1.5
	2400	2.5	-1.5	1.5
	2400	3.5	-1.1	1.7
	2400	4	-1.0	1.7
3:	3000	2	-2.0	1.65
	2800	2.2	-1.75	1.6
	2600	2.5	-1.5	1.6
	2280	3.5	-1.1	1.65
	2200	4	-1.0	1.55
ref:	2400	3	-1.25	1.5

Table 5. Table 5: Parameters used to fit Eq. 23 to the solutions of Figure 16 for different injection energies E inj subscript 𝐸 inj E_{\text{inj}} . All quantities are as described for Table 2 .

$E_{\text{inj}}$ (keV)	${\gamma }_a$	$E_{\text{tr}}$ (keV)	$E_{\text{cut}}$ (MeV)
10	$-2.0$	30	0.7
30	$-1.8$	30	1.0
60	$-1.25$	-	1.5

Equations

$$f_{\kappa }=A_{\kappa }{\left(1+\frac{v^2}{\kappa v_{\kappa }^2}\right)}^{-\kappa -1}\text{,}$$

$$v_{\kappa }^2=\left(2\kappa -3\right)\frac{k_{b}T}{\kappa m_p}\text{,}$$

$$A_{\kappa }=\frac{n}{{\left(\pi \kappa v_{\kappa }^2\right)}^{3/2}}\frac{\Gamma \left(\kappa +1\right)}{\Gamma \left(\kappa -\frac{1}{2}\right)}$$

$$V_{sw}^{\prime }$$

$$-\frac{1}{2s}\left(\left(s-1\right)\left(V_{sh}-V_1^{\prime }\right)\tanh \left(\frac{r-r_{sh}}{L}\right)\right)\text{,}$$

$$r_c={\tanh }^{-1}\left(\frac{2s{V}_c-V_1^{\prime }(s+1)-V_{sh}(s-1)}{(1-s)(V_{sh}-V_1^{\prime })}\right)L\text{,}$$

$$n=n_0\frac{V_{sh}-V_1^{\prime }}{V_{sh}-V_{sw}^{\prime }}{\left(\frac{r_0}{r}\right)}^2\text{,}$$

$$d{\mathcal{E}}_{sw}=\frac{1}{2}{m}_p\left(n_0{\left(V_1^{\prime }-V_{sh}\right)}^2-n{\left(V_{sw}^{\prime }-V_{sh}\right)}^2\right)\text{,}$$

$$d{T}_{sw}=\frac{2}{{\mathcal{N}}_{df}}\frac{d{\mathcal{E}}_{sw}}{n{k}_b}\text{,}$$

$$T_0=T_e{\left(r/r_0\right)}^{2\left(1-{\gamma }_c\right)}$$

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

$$+\frac{1}{3r^2}\frac{d\left(r^2{V}_{sw}^{\prime }\right)}{dr}E\xi (E)\frac{\partial f}{\partial E}+Q\text{,}$$

$${\kappa }_{rr}=\frac{v}{3}{\lambda }_{rr}\text{,}$$

$${\lambda }_{rr}={\lambda }_0{\left(\frac{R}{R_0}\right)}^{1/3}$$

$$-\frac{\partial f}{\partial t}={\mu }_r\frac{\partial f}{\partial r}+{\mu }_E\frac{\partial f}{\partial E}+\frac{1}{2}{\sigma }_r^2\frac{{\partial }^{2}f}{\partial r^2}\text{,}$$

$$dr$$

$$dE$$

$$dr$$

$$dE$$

$$r_{sh}=r_{\text{obs}}-V_{sh}t$$

$$\Delta t=\text{min}\left\{\frac{0.1L}{\text{max}\left(\left|\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{\kappa }_{rr}\right)\right|,V_{sw}^{\prime }\right)},\frac{0.1L^2}{{\kappa }_{rr}}\right\}\text{.}$$

$${\alpha }_i(t-\Delta t)={\alpha }_i(t)+Q\Delta t\text{.}$$

$$f_{\text{obs}}\left(r_{\text{obs}},E_{\text{obs}}\right)=\frac{1}{N_p}\sum _{i=1}^{N_p}{\alpha }_i$$

$$Q=\frac{f_{\kappa }}{m_p^3}\left|\frac{d{V}_{sw}^{\prime }}{dr}\right|H\left(\frac{E-E_{\text{inj}}}{E_{\text{inj}}}\right)\text{.}$$

$$j_{\text{DSA}}=j_0{\left(\frac{E}{E_0}\right)}^{{\gamma }_s}{\mathrm{e}}^{-{\left(E/E_{\text{cut}}\right)}^2}\text{,}$$

$${\gamma }_s=\frac{s+2}{2-2s}$$

$$j_{\text{DSA}}=j_0{\left(\frac{E}{E_0}\right)}^{{\gamma }_a}{\left(\frac{E^{\xi }+E_{\text{tr}}^{\xi }}{E_0^{\xi }+E_{\text{tr}}^{\xi }}\right)}^{\left({\gamma }_b-{\gamma }_a\right)/\xi }{\mathrm{e}}^{-{\left(E/E_{\text{cut}}\right)}^2}\text{.}$$

$${\left(\frac{{\kappa }_{rr}}{V_{sw}}\right)}_{\left(r=r_{sh}\right)}$$

$$={\kappa }_{rr}\left(\frac{1}{V_1}+\frac{1}{V_2}\right)\text{,}$$

$${\tau }_a=\frac{3}{V_1-V_2}\underset{p_{\text{inj}}}{\overset{p}{\int }}{\kappa }_{rr}\left(\frac{1}{V_1}+\frac{1}{V_2}\right)\frac{d{p}^{\prime }}{p^{\prime }}\text{,}$$

$$I_{E_l}=\underset{E_l}{\overset{\infty }{\int }}j(E)dE\text{.}$$

$$d\mathcal{N}=f\left(\vec{x},\vec{p},t\right)d^{3}x{d}^{3}p\text{.}$$

$$n=\int f\left(\vec{x},\vec{p},t\right)d^{3}p\text{.}$$

Full text

acceleration of particles — solar wind — Sun: heliosphere — shocks

\AuthorCallLimit

=3

1 Introduction

Energetic particle enhancements observed at Earth are often associated with the passage of interplanetary (IP) shocks, driven by e.g. coronal mass ejections (CMEs), during energetic storm particle (ESP) events (Bryant et al., 1962; Lario & Decker, 2002; Ho et al., 2009; Huttunen-Heikinmaa & Valtonen, 2009; Mäkelä et al., 2011). These enhancements are thought to occur as a result of particle acceleration at the shock, for which diffusive shock acceleration (DSA, Axford et al., 1977; Krymskii, 1977; Bell, 1978; Blandford & Ostriker, 1978) continues to be regarded a viable mechanism, often coupled with self-generated turbulence in the form of magnetohydromagnetic waves (Lee, 1983). See also the relevant discussions of Desai & Giacalone (2016) on ESPs and DSA. More recently, attention has been drawn to the notable flattening of proton energy spectra in observations directly preceding the passage of IP shocks at Earth (Lario et al., 2018). This is attributed to follow as a result of the propagation of accelerated particles from the shock to the observer, but is not fully explained. Moreover, while CME shock properties such as speed and compression ratio have been linked to the efficiency of particle acceleration (Lario et al., 2005b; Giacalone, 2012), shock obliquity (Ellison et al., 1995) and the nature of the seed particles (Desai et al., 2006) have also been identified as important factors.

It has been suggested that IP shocks can accelerate particles from thermal energies (e.g. Giacalone, 2005). Indeed, for quasi-parallel shocks, where particles can repeatedly cross the shock front along magnetic field lines, injection energies are assumed to be small, and Maxwellian-like distributions are proposed to provide adequate seed particles for DSA (Giacalone et al., 1992; Neergaard Parker & Zank, 2012). By contrast, injection energies are generally assumed to be high for perpendicular shocks, because the perpendicular diffusion coefficient is small for particles interacting resonantly with micro-scale turbulence. However, the injection energy is considerably lowered when particles experience perpendicular diffusion along intermediate-scale meandering magnetic field lines. The injection speed at perpendicular shocks can become comparable to those at parallel shocks, and DSA thereby more effective, if the amplitude of field line meandering is large enough (Giacalone, 2005). On the other hand, the shock front itself may be rippled: while an idealized shock propagating through an approximately radial magnetic field may vary from being quasi-perpendicular at the CME flanks to quasi-parallel near the nose, realistically, local geometries may resemble any obliquity (see Klein & Dalla, 2017).

Aside from conducive shock geometry, injection into the DSA process is facilitated by the formation of non-thermal seed particles (Neergard Parker et al., 2014; Zank, 2017), which has also recently been investigated using kinetic hybrid simulations (Caprioli et al., 2015; Sundberg et al., 2016). Suprathermal tails are often observed in solar wind (SW) velocity distributions (Collier et al., 1996; Maksimovic et al., 1997; Chotoo et al., 2000; Qureshi et al., 2003) and are considered conducive features for the injection of particles into DSA (Desai et al., 2006; Kang et al., 2014). To parameterize these distributions, Vasyliunas (1968) introduced the Kappa($\kappa$)- distribution function, which characterizes both the Maxwellian core and the suprathermal tail at higher energies. Refer to Pierrard & Lazar (2010) and Livadiotis & McComas (2013) for complete reviews on the theory and applications of $\kappa$-functions. These functions characterise SW distributions using only three parameters, namely, the $\kappa$-parameter, which relates to the spectral index of the high-energy tail, and the equivalent temperature and number density of the plasma or particle species considered (Formisano et al., 1973; Chateau & Meyer-Vernet, 1991). Of course, quantities such as the plasma density and temperature are observed to change during the passage of an IP shock (e.g. Lario & Decker, 2002), which also transforms the $\kappa$-function and thereby the properties of the potential DSA seed population it represents.

In this study, two broad scientific questions are addressed, namely, how the prior heating of the SW distribution affects acceleration at IP shocks, and what role these shock-accelerated particles play in producing the spectral features observed during ESP events. The acceleration of SW particles is modelled by solving a set of stochastic differential equations (SDEs; introduced in Section 3), equivalent to the Parker (1965) transport equation, in spherical symmetry. The applications are limited to halo CME shocks expanding radially at constant speeds. In this spherically symmetrical scenario, the shock normal is radially aligned so that the shock obliquity can be approximated using the Parker (1958) spiral angle.

The manner in which the shock passage affects the energy distribution of SW particles is investigated both before and after injection into the DSA process: It is firstly considered how the initial SW distribution transforms in response to changes in the plasma properties (Section 2), and later how the resulting distribution affects the subsequent DSA process as a seed population (Section 6). The classical spectral characteristics of DSA for travelling shocks, and the comparative efficiency of DSA at fast and strong shocks, are revisited in Sections 4 and 5, respectively. In addition to reproducing these more typical DSA features, the model is also applied to investigate the spectral flattening reported by Lario et al. (2018) ahead of IP shock passages (Section 6.2), as well as to identify the original spatial distributions and energies of seed particles (Section 6.3) that would eventually contribute to intensities during ESP events.

2 Modelling CME shock-induced changes in the solar wind

To describe both the thermal and suprathermal velocity distributions observed in the SW, the $\kappa$-function is implemented in terms of particle speed $v$ as

[TABLE]

where $v_{\kappa}^{2}$ is the generalised thermal speed given by

[TABLE]

with $k_{B}$ the Boltzmann constant and $m_{p}$ the proton mass, and where

[TABLE]

is the normalization constant obtained when setting $\int f_{\kappa}d^{3}v=n$, with $n$ taken as the SW number density and where $\Gamma$ is the Gamma function. This normalization is discussed in some detail in Appendix A.3.

The $\kappa$-parameter is related to the power-law index of the suprathermal tail. Note that Eq. 1 is defined for $\kappa>3/2$, and that if $\kappa\rightarrow\infty$, it reduces to a Maxwellian. For the correct interpretation of Eq. 2 and the associated temperature $T$, it is instructive to consider the discussion by Hellberg et al. (2009): $v_{\kappa}$ is introduced by Vasyliunas (1968) as the most probable particle speed, associated with the non-relativistic kinetic energy of $E_{\kappa}=m_{p}v_{\kappa}^{2}/2$. Evaluating the second moment of $f_{\kappa}$ yields a total energy of $NE_{N}=\frac{3}{2}NE_{\kappa}\kappa/(\kappa-3/2)$, where $E_{N}$ and $N$ are the mean energy per particle and the total number of particles, respectively. Eq. 2 then follows upon the introduction of the plasma temperature $T$ (originally by Formisano et al., 1973) through the invocation of the equipartition theorem, $E_{N}=\frac{3}{2}k_{B}T$, for a monatomic gas. Although this temperature definition and the foregoing assumption of the equipartition of energy are not strictly valid for non-Maxwellian distributions, their use in this manner has become standard practice and is generally considered appropriate (see Hellberg et al., 2009, and the references therein).

2.1 Changes in plasma properties across the shock

An objective of this study is to consider how the $\kappa$-function changes during the passage of an IP shock, and to implement this transformed distribution as a seed-particle spectrum for DSA. To this end, while the $\kappa$-index is assumed constant across the shock, it is necessary to model the change in $n$ and $T$, upon which $f_{\kappa}$ depends (see also Livadiotis, 2015). For consistency, transitions across the shock are modelled to correspond to that of the flow speed as would be observed by e.g. a spacecraft in Earth’s orbit. The SW is consequently modelled to transition across the shock between up- and downstream flow speeds in the spacecraft (or fixed) frame, that is, $V_{1}^{\prime}$ and $V_{2}^{\prime}=\left(V_{sh}\left(s-1\right)+V_{1}^{\prime}\right)/s$, respectively, according to

[TABLE]

where $s$ is the shock compression ratio (or ratio of up- and downstream flow velocities in the shock frame), $V_{sh}$ is the shock speed, $r-r_{sh}$ is the radial position relative to the position of the shock $r_{sh}$, and $L$ is a characteristic length used to specify the broadness of this transition. Figure 1 illustrates this transition for the reference parameters listed in Table 1. Should $V_{sh}=0$, Eq. 2.1 reduces to an expression (used by e.g. le Roux et al., 1996) for a stationary shock.

Note furthermore that $L$ is not a direct measure of the shock width. Its extent can however be approximated from Eq. 2.1 as $\triangle x_{sh}\approx 2r_{c}$, with

[TABLE]

where a fraction $\mathcal{C}\in[0,0.5]$ is chosen such that $V_{c}=V_{1}^{\prime}+\mathcal{C}(V_{2}^{\prime}-V_{1}^{\prime})/2$ is the flow speed at some distance $r_{c}$ from the shock. Typically, to approximate the shock width, $\mathcal{C}=0.01$ so that $V_{c}\approx V_{1}^{\prime}$. It is illustrated in Figure 1 how the lengths $\triangle x_{sh}$, $L$, and the actual shock width compare.

With the flow speed transition given, the change in number density $n$ across the shock follows from the continuity equation. Of course, the total factor by which the number density jumps across the shock must be equal to the compression ratio. In terms of $V_{sw}^{\prime}$ it follows, in spherical coordinates, that

[TABLE]

where $n_{0}$ is the number density at some reference position such as the Earth ($r_{0}=$ 1 AU), and where the flow speeds are transformed back to the shock frame, where their ratio is equal to $s$. This transition is also shown in Figure 1.

To estimate the jump in temperature across the shock, instead of implementing the relevant hydromagnetic Rankine-Hugoniot jump condition, a simpler approach capturing the essential physics is used. The magnetic energy is expected to be much smaller than the thermal energy in the flow-dominated region considered in this study. It is therefore assumed that the change in kinetic energy density $d\mathcal{E}_{sw}$ of the bulk solar wind flow, as a function of position through the shock, is converted to thermal energy. That is,

[TABLE]

which is related through the equipartition theorem to the change in temperature $dT_{sw}$ by

[TABLE]

where $\mathcal{N}_{df}=2/(\gamma_{c}-1)$ is the number of degrees of freedom and $\gamma_{c}$ is the ratio of specific heats. $dT_{sw}$ is subsequently added to the unshocked temperature value, that is, $T=T_{0}+dT$, where

[TABLE]

and $T_{e}$ is the plasma temperature at the Earth. For a monatomic gas for which $\gamma_{c}=5/3$ it follows that $T_{0}\propto r^{-4/3}$. It is assumed that the temperature gain during the shock passage is not significant enough to disassociate molecules, and hence $\gamma_{c}$ remains constant as for a calorically ideal gas. Furthermore, given the multi-species composition of the SW plasma, $\gamma_{c}$ may deviate from the monatomic value of $5/3$. It is noted that while the hydromagnetic treatment of the shock jump conditions may break down in the transition region itself, it approximates the up- and downstream quantities adequately. Since DSA requires particles to interact with the shock on length scales exceeding its width (e.g. Jones & Ellison, 1991; Krülls & Achterberg, 1994), this description is sufficient for the current study. See also Section 4.1.

The shock transitions of $V_{sw}^{\prime}$, $n$ and $T$ are displayed in Figure 2 as viewed by an observer at the position of the Earth. The parameters used to define these quantities are listed in Table 1 and are mostly informed by observations of CMEs during the 2003 Halloween epoch (Skoug et al., 2004; Richardson et al., 2005; Gopalswamy et al., 2005; Lario et al., 2005a; Wu et al., 2005). In particular, SW flow speeds in excess of 1800 km s*-1* (Skoug et al., 2004) and shock speeds of $\sim$2400 km s*-1* (Gopalswamy et al., 2005; Wu et al., 2005) are noted. Since density measurements during some larger events were uncertain (e.g. Skoug et al., 2004; Wu et al., 2005), the compression ratio for such events can be calculated using the ratios of shock-frame flow velocities. For example, $s=(V_{sh}-V_{1}^{\prime})/(V_{sh}-V_{2}^{\prime})=3$ for $V_{sh}=$ 2400 km s*-1*, $V_{1}^{\prime}=$ 600 km s*-1* and $V_{2}^{\prime}=$ 1800 km s*-1*.

2.2 Evolution of the $\kappa$-function during the shock passage

Aside from $\kappa$, which is kept fixed, $f_{\kappa}$ is shown above to depend on parameters $n$ and $T$, both of which change during the shock passage. Additionally, to illustrate how the SW distribution shifts when the flow speed is shocked to higher values, it is necessary to express the particle speed in $f_{\kappa}$ relative to flow speed in the spacecraft frame, that is, $v:\rightarrow v-V_{sw}^{\prime}$ (see also Leubner, 2004; Kong et al., 2017). The SW distribution is therefore expected to be transformed by the shock in at least three ways, depending on how $V_{sw}^{\prime}$, $n$ and $T$ change. To illustrate this transformation, the SW distribution, as represented by $f_{\kappa}$ and viewed by an observer at Earth, is shown at different stages during the shock passage on the right-hand side of Figure 2 for parameters as specified in Table 1. Bear in mind these distributions have not yet been injected into the DSA process, and that these transformations occur solely because the shock had heated the SW plasma. Also note that for consistency with the units in which observations are typically presented, $f_{\kappa}$ is converted throughout this study into units of differential intensity (as discussed in Appendix A.2) and is subsequently denoted as $j_{\kappa}$.

With the shock still 0.02 AU away from the Earth, Figure 2 shows that the SW distribution at Earth is as of yet unchanged, since none of the plasma quantities have been shocked at this point. When the shock passes Earth, that is, when $r_{sh}=r_{0}$, the temperature is shown to have increased substantially, accompanied by more modest increases in the flow speed and density. The most obvious changes the SW distribution incurred at this point is that it notably broadened as a result of the temperature increase, and that the thermal peak shifted to higher energies as a result of the increased flow speed. Note that because $f_{\kappa}$ is normalised to the number density, the conservation of particles demands that the peak intensity of the distribution decreases when it broadens. Considering, finally, how the distribution changes when the shock has moved 0.02 AU beyond the Earth, it shifts further, to higher energies, as the flow speed attains its full downstream value, while increasing overall intensities due to the increase in number density. At this point the temperature had already plateaued and hence the distribution shows no appreciable broadening from when the shock passed Earth’s position.

These SW distributions are specified as source spectra for DSA in Section 3.2.

3 Modelling energetic particle acceleration at CME shocks

The events of interest in this study are large ESP events with small associated anisotropies. These events are typically associated with fast-moving shocks driven by halo CMEs (Mäkelä et al., 2011), that is, those propagating radially outward and approximately centred on the solar disk. Therefore, to describe the transport and DSA of energetic particles associated with such events, it sufficient to solve the Parker (1965) transport equation (TPE) for a single spatial dimension. See also the motivation offered by Giacalone (2015) in this regard. Hence, in radial coordinates, the TPE is written as

[TABLE]

with $\xi(E)=(E+2E_{p})/(E+E_{p})$, where $E_{p}$ is the proton rest-mass energy. The TPE contains the relevant particle transport processes such as SW convection, spatial diffusion, and energy changes due to transport in regions of compressing or expanding SW flows. The last-mentioned process implicitly simulates DSA for regions such as shocks with negative SW velocity divergences. $Q$ represents a particle source function. Note that the TPE is solved for the pitch-angle averaged, omni-directional distribution function $f=f_{0}\left(r,p,t\right)$, which is only a function of position, scalar momentum and time; see Appendix A. As before, when presented, the distribution function is converted to units of differential intensity, that is, $j=p^{2}f$.

The diffusion coefficient considered is the effective radial diffusion coefficient $\kappa_{rr}$, which relates to the equivalent mean free path (MFP) $\lambda_{rr}$ according to

[TABLE]

where $v$ is the particle speed and where $\lambda_{rr}$ is given by

[TABLE]

where $R=\sqrt{E^{2}+2EE_{p}}$ is the particle rigidity and where $\lambda_{0}$ is a reference MFP defined at $R_{0}\equiv$ 1 GeV. Given that the parallel diffusion coefficient is much larger than the perpendicular coefficient, and magnetic field lines are largely radial for the fast flow speeds considered, diffusion along field lines is assumed to dominate. The rigidity dependence implemented here is hence chosen to emulate that predicted for parallel diffusion by quasi-linear theory (Jokipii, 1966) for a Kolmogorov-distributed turbulence power spectrum. Kallenrode et al. (1992) reports no pronounced variation in MFPs with radial distance between the Sun and 1 AU for diffusive events. As such, a radial dependence is omitted in Eq. 12. Furthermore, radial MFPs of 0.02 to 0.15 AU are reported for 0.3 to 0.8 MeV electrons (Kallenrode et al., 1992), and the ratio of MFPs for 18 MeV protons to 1 MeV electrons as $1.6\pm 0.9$ (Kallenrode, 1993). Taking the lower limits of the aforementioned quantities into account, the corresponding lower limit of $\lambda_{0}$ (defined at 1 GeV) for protons is estimated from Eq. 12 as $\sim$0.05 AU. $\lambda_{0}$ is varied in Section 4 and chosen as 0.06 AU elsewhere, using the above lower limit and the results of the aforementioned section as guidelines.

3.1 The numerical model: SDEs

The TPE specified in multiple computational dimensions typically requires numerical methods to solve. However, instead of solving Eq. 3 using finite-difference methods (e.g. Giacalone, 2015), the transport and DSA of energetic particles are simulated here by solving an equivalent set of SDEs (see also Krülls & Achterberg, 1994; Marcowith & Kirk, 1999; Zhang, 2000). The aspects of the SDE approach that are important for this study are detailed below. Refer to Strauss & Effenberger (2017) for a comprehensive review.

Eq. 3 is conveniently written in the form of the time-backward Kolmogorov equation

[TABLE]

from which the SDEs can be cast into the form of the Itô equation (see e.g. Zhang, 1999)

[TABLE]

where $dW=\sqrt{dt}\ \Lambda(t)$ represents the Wiener process and $\Lambda(t)$ is a simulated Gaussian-distributed pseudo-random number. Note that the forward time $t^{\prime}$ as used in Eq. 3 is related to the backward time $t$ according to $t=t_{T}-t^{\prime}$, where $t_{T}$ is the total simulation time. It thus follows that $\partial/\partial t^{\prime}=-\partial/\partial t$ in Eq. 13. Substituting the coefficients corresponding to $\mu_{r}$, $\mu_{E}$, and $\sigma_{r}$ from Eq. 3 into Eqs. 14 and 15 yields

[TABLE]

which are equivalent to the TPE in Eq. 3.

These SDEs are integrated in a time-backwards fashion: Starting from an observational point $(r_{\text{obs}},E_{\text{obs}})$ at which the value of the distribution function is sought, Eqs. 16 and 17 are solved iteratively using the Euler-Maruyama numerical scheme (Maruyama, 1955) for a finite time step $\Delta t$. The coordinates $r$ and $E$ are updated upon each iteration until $t^{\prime}=0$, or equivalently, until $t=t_{T}$. When simulating ESP events, $t_{T}$ is chosen as the transit time of the shock travelling (in a time-backwards fashion) from $r_{\text{obs}}$ to the inner modulation boundary near the Sun ($r_{\text{min}}=5R_{\odot}$), that is, $t_{T}=(r_{\text{obs}}-r_{\text{min}})/V_{sh}$. The shock position $r_{sh}$ is therefore updated in step with $r$ and $E$, according to

[TABLE]

where $0<t\leq t_{T}$. A constant shock speed is assumed, since the mean acceleration of CME-driven shocks associated with ESP events is reportedly close to zero (Mäkelä et al., 2011). Note that the time is not incremented by the same amount during each iteration. It is instead specified to vary, depending on the dominant transport process at the current position, to maintain a fixed step length, that is,

[TABLE]

This has at least two advantages: Firstly, limiting the step length to some fraction of the length scale $L$, which is associated with the shock, ensures that the shock structure is properly resolved. Secondly, scaling $\Delta t$ in this manner saves computation time, since instead of applying a very small fixed time step, $\Delta t$ will only be small when it is required. As a result of this variable time step, computation times are typically longer for larger values of $\kappa_{rr}$, such as at higher energies, and for narrower shocks (smaller $L$-values).

The above time-integration is repeated $N_{p}$ (typically, $\sim 10^{6}$) times for each observational point $(r_{\text{obs}},E_{\text{obs}})$, thereby tracing out $N_{p}$ trajectories in $r$ and $E$ of phase-space density elements, conventionally referred to as pseudo particles. The source function $Q$ in Eq. 3 is handled in the SDE approach as a correction term (Strauss & Effenberger, 2017). Expressing it as a rate of contribution to the distribution function allows the contribution per pseudo particle (or its amplitude) to be calculated iteratively along the integration trajectory. That is, $Q=d\alpha/dt\implies\Delta\alpha=Q\Delta t$, from which it follows that

[TABLE]

The source contribution is then tallied for $N_{p}$ pseudo particles such that

[TABLE]

gives the value of the distribution function at the observational point at time $t^{\prime}=t_{T}\ (t=0)$. Section 3.2 discusses how the source function itself is specified. Furthermore, physical particle distributions can also be obtained by tallying the amplitude-weighted flux contributions of pseudo particles within appropriately sectored spatial and energy intervals (or bins) at any time during the simulation, divided by the number of pseudo particles counted within each sector. This method is used in Section 6.3 to trace likely acceleration sites and seed-particle energies in a time-backwards fashion.

Finally, a reflective inner boundary is implemented such that if $r<r_{\text{min}}\implies r=2r_{\text{min}}-r$, while at the outer boundary, if $r>r_{\text{max}}=1.4$ AU, the time-integration routine is interrupted and that pseudo-particle’s contribution is discarded, defining an absorbing boundary condition.

3.2 Modelling the shock source function

As per Eq. 20, when a pseudo-particle is traced back to a point where it interacts with the shock, its amplitude is attributed a value that depends on the product of the source function $Q$ and the time $\Delta t$ it spends interacting with the shock. Indeed, considering Eq. 3, $Q$ should have units corresponding to that of the distribution function per unit time. Otherwise, the specification of $Q$ is largely arbitrary: e.g. Giacalone (2015) employs delta functions, others (e.g. le Roux et al., 1996) specify distribution functions to represent seed populations, while Malkov & Völk (1995) also consider shock-heated seed populations. Drawing on aspects of the aforementioned examples, the following source function is proposed:

[TABLE]

Here, $f_{k}$ is the $\kappa$-function representing the SW distribution. As discussed in Appendix A.3, the $m_{p}^{3}$ factor is included for dimensional consistency between $f_{\kappa}$ and the distribution function in the TPE. Since $f_{\kappa}$ is already normalised to the number density, scaling the SDE solutions is not necessary. The dimensionless Heaviside function $H\left((E-E_{\text{inj}})/E_{\text{inj}}\right)$ ensures that only the contributions of particles with energies larger than the injection energy $E_{\text{inj}}$ are included. Note that $E_{\text{inj}}=$ 60 keV unless stated otherwise; see also Section 6.1. The absolute value of the SW velocity gradient $\left|dV_{sw}^{\prime}/dr\right|$ is included to provide a spatial region with a width and position that is self-consistently associated with the shock and that has a reasonable probability of being sampled by pseudo particles. Figure 3 shows how $\left|dV_{sw}^{\prime}/dr\right|$ is representative of the width of the shock. While it is sufficient for a physical particle to pass between the upstream and downstream media in order to gain energy, using this time-backward SDE method, a finite-width region has to be sampled by pseudo particles to register flux emanating from the shock.

3.3 Advantages and limitations of the SDE approach

The SDE approach has been successfully applied in many instances to simulate space particle transport (Zhang, 1999; Pei et al., 2010; Strauss et al., 2011, 2013; Moloto et al., 2019) and DSA in particular (Krülls & Achterberg, 1994; Marcowith & Kirk, 1999; Zhang, 2000; Zuo et al., 2011; Hu et al., 2017). The time-backward approach is favoured here due to its efficiency, since every simulated pseudo-particle contributes to intensities at the desired observational point. Additionally, in a similar fashion to how the boundary interactions of pseudo particles are used to trace the most probable points of entry into the heliosphere for cosmic rays (e.g. Strauss et al., 2011), this backward tracing of phase-space trajectories is similarly utilised in this study to map probable acceleration sites and seed-particle energies; see Section 6.3. Since these pseudo-particle trajectories are solved entirely independent of each other, this approach is also conducive to the utilisation of parallel computing platforms.

However, this mutual independence of the SDE solutions also limits applications to the test-particle case. Non-linear effects of shock acceleration such as self-generated turbulence (Lee, 1983; le Roux & Arthur, 2017) or particle mediation of shock structures (Mostafavi et al., 2017) can therefore not be considered. Moreover, any process that requires the calculation of particle intensity gradients becomes computationally expensive (e.g. Moloto et al., 2019). Nevertheless, the SDE approach remains suitable to investigate the intricacies of classical DSA, without the limitations imposed by a numerical grid, such as instabilities involving the large gradients that are typically encountered at shocks. It is shown in Appendix B that the SDEs produce appreciably similar results to a finite-difference approach in reproducing ESP events for the same parameter set. Furthermore, as opposed to finite-differences schemes, the SDE model introduced in this study can seamlessly be expanded to a larger number of computational dimensions (e.g. Pei et al., 2010) as required to study particle acceleration at more complicated shock geometries.

4 Spectral features and diffusion dependence of shock-accelerated particles.

The most distinctive characteristics of shock-accelerated particles are observed in their energy spectra, which are considered here at the time the shock passes an observer, e.g. a spacecraft, near the Earth. Implementing the model configuration discussed in Sections 2 and 3, and scaling the diffusion coefficient by varying the value of $\lambda_{0}$ in Eq. 12, the dependence of shock-accelerated spectral features on this transport process is investigated. The resultant spectra are shown in Figure 4.

Qualitatively, the typical DSA-associated features are evident: A power-law distribution at lower energies transitioning to an exponential-like decrease at higher energies (Ellison & Ramaty, 1985). This distribution can be described using a simple function in the form of

[TABLE]

where $j_{0}$ is a differential intensity defined at some reference energy $E_{0}$, $E_{\text{cut}}$ is the cut-off (or roll-over) energy above which the distribution begins to decrease exponentially, and where

[TABLE]

is the spectral index (for $E\ll E_{p}$) associated with the compression ratio $s$ of the shock (see also Ellison et al., 1990). For $s=3$, as specified for these solutions, it is therefore expected that $\gamma_{s}=-1.25$. However, not all of the spectra in Figure 4 appear to follow a $E^{-1.25}$ power law. This is accentuated in Figure 5, which shows spectral indices for solutions with smaller $\lambda_{0}$-values actually varying with energy even before the exponential decreases ensue. To describe this behaviour, it is useful to generalise Eq. 23 as follows:

[TABLE]

This describes a function transitioning between power-law indices $\gamma_{a}$ and $\gamma_{b}$ about $E_{\text{tr}}$, with $\xi$ specifying the smoothness of this transition, and which rolls over into an exponential decrease above $E_{\text{cut}}$. Note that if $\gamma_{a}=\gamma_{b}$, this expression simplifies to the form of Eq. 23. Setting $\gamma_{b}=\gamma_{s}$, where $\gamma_{s}=-1.25$ for $s=3$, and choosing $\xi=0.8$, the function given in Eq. 25 is fitted to the solutions in Figure 4 with parameters as presented in Table 2.

Both $\gamma_{a}$ and $E_{\text{cut}}$ respond to varying the value of $\lambda_{0}$, affecting both the hardness of the spectra and the energies up to which they are accelerated. The behaviour of these two sets of features and the underlying physics are discussed separately in the subsections below.

4.1 Fractional compression sampling: On shock widths and diffusion length scales

Classically, the spectral index associated with a DSA-produced spectrum is a function of the shock compression ratio alone. The behaviour observed in Figure 4, where shock-accelerated spectra become softer, displaying smaller spectral indices, for smaller diffusion coefficients, is therefore not theoretically expected. Note, however, the solutions for $\lambda_{0}\geq 0.035$ AU do follow the theoretically predicted power law of $E^{-1.25}$ for $s=3$. Figure 5 shows the spectral indices of the solutions presented in Figure 4. Here, spectral indices are also shown to be equal to $-1.25$ for large $\lambda_{0}$-values at low energies. Indices become progressively smaller for smaller $\lambda_{0}$-values, but increase toward higher energies. Consider that particles with $\lambda_{0}\lesssim 0.035$ AU may be sampling only a fraction of the total compression of the shock. The spectra harden toward higher energies, because the MFPs themselves increase with energy and progressively greater fractions of the total compression are sampled.

This effect is illustrated in Figure 6. Recall from Eq. 11 that $\kappa_{rr}\propto\lambda_{rr}$, which is scaled using $\lambda_{0}$. Since $\kappa_{rr}$ does not have dimensions of length, when comparing to other length scales it is useful to define the diffusion length scale $\kappa_{rr}/V_{sw}$, which is often expressed as the sum of the down- and upstream values at the shock (e.g. Steenberg & Moraal, 1999). That is,

[TABLE]

where $(\kappa_{rr})_{1}=(\kappa_{rr})_{2}=\kappa_{rr}$, since $\kappa_{rr}$ is assumed not to change across the shock, with subscripts 1 and 2 denoting up- and downstream values, respectively. Here, $V_{1}=V_{sh}-V_{1}^{\prime}$ and $V_{2}=(V_{sh}-V_{1}^{\prime})/s$ are the up- and downstream flow speeds in the shock frame. This quantity is plotted as a function of energy in Figure 6, along with the shock width, calculated using Eq. 5 for $L=0.005$ AU. Note that for $\lambda_{0}=0.06$ and $0.1$ AU, for which spectra are aligned with $E^{-1.25}$, the diffusion length scales are greater than the shock width for all energies in the considered domain. Those particles therefore sample the full compression ratio, and hence their energy distributions display the power law associated with $s=3$. A similar effect is reported for the heliospheric termination shock (e.g. Arthur & le Roux, 2013). For $\lambda_{0}<0.06$ AU, particles sample fractional compression ratios up to the energies where their respective diffusion length scales begin to exceed the shock width. Indeed, the energies at which $\kappa_{rr}/V_{sw}$ and the shock width intersect provide good estimates for $E_{\text{tr}}$ in Table 2.

Since the prediction of DSA given in Eq. 24 is only observed where $\kappa_{rr}/V_{sw}$ is larger than the width of the shock, these results are in agreement with the length scale hierarchy of $\kappa_{rr}/V_{sw}\gg\lambda_{rr}\gg\Delta x_{sh}$ required for classical DSA to be valid (Blandford & Ostriker, 1978; Jones & Ellison, 1991). Given the $L$-dependence of the variable time step (Eq. 19), to limit computation times, the shock width in this study is chosen much broader ($L=0.005$ AU) than typical IP shock widths ($L\sim 10^{-6}$ AU; e.g. Sapunova et al., 2017). Nevertheless, the results of the model should remain valid as long as the aforementioned length scale hierarchy is observed.

4.2 Finite-time acceleration and the termination of shock-accelerated spectra

The second set of features considered entails the highest energies attained by shock-accelerated spectra before intensities begin to decrease exponentially. Figure 4 illustrates that spectra for smaller $\lambda_{0}$-values are accelerated to higher energies before terminating. These energies, represented by $E_{\text{cut}}$, are listed in Table 2. Since they notably respond to varying $\lambda_{0}$, it seems reasonable to expect that this spectral transition also occurs due to $\kappa_{rr}/V_{sw}$ attaining some characteristic length. Shock-accelerated spectra have previously been reported to roll over due to diffusion length scales becoming comparable to system sizes or some related length of the shock geometry (Ellison & Ramaty, 1985; Steenberg & Moraal, 1999). Indeed, Figure 6 shows the diffusion length scales corresponding to each of the cut-off energies of the solutions in Figure 4 are reasonably similar, distributed between 0.1 and 0.4 AU. However, these do not resemble any obvious length scale in either the physical system or within the numerical model. Simulations where the position of modulation boundaries were varied yielded negligible effects on energy spectra, suggesting that system size is not the limiting factor in this instance.

Another instructive quantity that is related to the diffusion length scale is the acceleration time (Drury, 1983; Ellison et al., 1990), which is the time required to accelerate particles from a particular energy (or equivalent momentum) to another at a planar shock. This is expressed in terms of the diffusion length scale as

[TABLE]

where $V_{1}$ and $V_{2}$ are as defined for Eq. 4.1 and $p_{\text{inj}}$ is the momentum equivalent to $E_{\text{inj}}$. Note that momentum is generally related to kinetic energy according to $p=(1/c)\sqrt{E^{2}+2EE_{p}}$ (for protons), with $E_{p}$ the proton rest-mass energy. Eq. 27 is used in Figure 7 to approximate the time required to accelerate spectra up to different energies for each of the length scale profiles shown in Figure 6. Also shown in Figure 7 as symbols are the times required to accelerate spectra up to the observed cut-off energies listed in Table 2 for each respective $\kappa_{rr}/V_{sw}$-profile. Similar to the diffusion length scales corresponding to the cut-off energies, the associated acceleration times also have similar values. These $\tau_{a}$-values also compare well to the total duration of the simulation, which is equal to the time taken by the shock to travel from near to the Sun to Earth at 2400 km s*-1*, that is, $\sim$17.3 hours. This suggests that for each of the spectra presented in Figure 4, $E_{\text{cut}}$ is the highest energy that could be attained in the available time (see also Channok et al., 2005). At $E>E_{\text{cut}}$, the rate of escape of particles from the shock begins to exceed the acceleration rate, which leads to the observed exponential intensity decreases. In this context, smaller $\lambda_{0}$-values, and consequently smaller diffusion coefficients, serve to better confine particles near the shock. This, in turn, reduces the time required to accelerate particles up to a particular energy, or stated differently, allows particles to be accelerated to higher energies within the available time frame.

4.3 Spectral features and acceleration efficiency

The preceding discussions reveal that the characteristics of shock-accelerated spectra are sensitive to the value of diffusion length scales in two opposing ways: Diffusion length scales should exceed shock widths for spectra to display DSA-predicted power-law indices, but accelerated spectra terminate at lower energies for larger diffusion length scales. The top panel of Figure 8 illustrates this dichotomy. The energies above which $\kappa_{rr}/V_{sw}$ exceeds the shock width, that is, $E_{\text{tr}}$, become smaller with increasing $\lambda_{0}$-values, implying greater overall spectral hardening, while similarly decreasing $E_{\text{cut}}$-values imply that spectra terminate at lower energies. Harder spectra yield larger intensities of energetic particles, but so do higher cut-off energies. Yet, with regards to diffusion properties, they are attained in opposite ways. Neither of these therefore necessarily provide meaningful measures of acceleration efficiency on their own.

To evaluate which combination of spectral features yield greater intensities of energetic particles, it is useful to consider particle fluences. These are calculated by integrating differential intensities over energy as follows:

[TABLE]

This integral converges, since $j$ fortuitously always decreases as $E\rightarrow\infty$. Eq. 28 is evaluated for each of the spectra shown in Figure 4 with two different lower limits, namely $E_{l}=E_{\text{inj}}$ and 1 MeV. The resulting energy-integrated intensities are shown in the bottom panel of Figure 8. The total intensity for all particles with $E>E_{\text{inj}}$ peaks at $\lambda_{0}\sim$ 0.035 AU, where the power-law segments of spectra have reached, or nearly reached, their maximum hardness. $E_{\text{cut}}$ becomes more important when considering higher-energy particles. For intensities integrated upward from 1 MeV, fluences peak between $\lambda_{0}=$ 0.02 and 0.035 AU and notably decrease for larger $\lambda_{0}$-values due to spectra rolling over at lower energies.

It is revealed that both spectral hardness and the maximum energies spectra attain contribute appreciably to the number of energetic particles the DSA process is able to produce. Of course, the energies up to which spectra extend become increasingly important for the intensities of higher-energy particles.

5 Strong shocks versus fast shocks: Which is the more efficient particle accelerator?

Fast-moving CME shocks with large compression ratios are reported to be more efficient at accelerating energetic particles (Lario et al., 2005b; Mäkelä et al., 2011; Giacalone, 2012). How conducive each of these shock properties are to producing large numbers of high-energy particles, especially as opposed to each other, raises an interesting subject for investigation. Both the shock speed $V_{sh}$ and the compression ratio $s$ affect how particle distributions evolve, and often in contradicting (or even self-contradicting) ways. Prior to injection, both properties affect the heating of the SW, while the hardness of shock-accelerated spectra is already shown to depend on the compression ratio. To explore and compare the effects of these two parameters, they are varied in three different ways: $V_{sh}$ and $s$ are each varied separately with the other remaining fixed, and they are varied together such that the factor by which the fixed-frame flow speed jumps across the shock remains constant. For ease of reference, these three sets of parameter configurations are shown in Table 3, along with the corresponding downstream flow speeds in both the fixed frame and the shock frame, denoted $V_{2}^{\prime}$ and $V_{2}$, respectively.

Varying $V_{sh}$ and $s$ affects the shock transitions of the SW flow speed, number density, and temperature, and thereby also influences how the energy distribution of SW particles changes during the passage of the shock. Each of the aforementioned parameters’ shock transitions are shown in Figure 9 for the configurations listed in Table 3. The fixed-frame flow speed, that is, the flow speed as viewed by an observer stationed at Earth, increases with both $V_{sh}$ and $s$ when each is varied separately. The results, as shown in the top row of frames in Figure 9, yield larger jumps in flow speed should either $V_{sh}$ or $s$ be increased. As intended, when $V_{sh}$ and $s$ are varied together, the downstream flow speed remains constant in the fixed frame. Hence the factor by which the flow speed increases is the same for these instances. The corresponding shock transitions of the number density and temperature in Figure 9 follow from Eqs. 6 and 8. Figure 10 shows how the SW energy distributions change in response to these shock transitions. It appears the SW energy distribution is heated by comparable amounts for varying shock speeds and strengths. However, $j_{\kappa}$ appears marginally more sensitive to changes in temperature, which in turn is affected most appreciably by changes in $V_{sh}$. Therefore, it can be argued that $V_{sh}$ contributes more towards heating the SW distribution, as opposed to $s$, considering the large shifts in energy and broadening of $j_{\kappa}$ that follows from varying it.

The heated SW distributions discussed above essentially serve as input spectra, from which particles are injected into the DSA process at the shock for $E>E_{\text{inj}}$. The shock-accelerated spectra obtained by solving the SDEs of Section 3.1 for each of the configurations in Table 3 are presented in Figure 11. The analyses of these spectra are analogous to those of spectra presented in Figure 4, while the discussions that follow draw on concepts introduced in Section 4. Accordingly, functions are fitted to the SDE solutions and the parameters tabulated for each configuration in Table 4. The value of $\lambda_{0}$ is chosen as 0.06 AU throughout this section. With reference to Figure 6, this implies the diffusion length scale exceeds the shock width for all considered energies and that the full compression ratio is sampled. It is therefore sufficient to fit the single power-law function of Eq. 23, with spectral indices as predicted by DSA.

The shock-accelerated spectra presented in Figure 11 share the same general features as those encountered before: a power-law segment extending from near the injection energy up to where it rolls over into an exponential decrease. The following discussion considers how these features change as a result of varying the shock speed and compression ratio. Firstly, varying $V_{sh}$, while keeping the compression ratio fixed at $s=3$, produces spectra as shown in the left column of Figure 11. They are power-law distributed with a spectral index of $-1.25$, as illustrated in the accompanying frame below. This is the index expected from DSA for $s=3$ (see Eq. 24). The sizeable differences in the intensities of these three spectra are due in large part to the heated SW spectra, which themselves are notably affected by varying $V_{sh}$. When varying $s$ instead and keeping the shock speed fixed at $V_{sh}=$ 2400 km s*-1* as shown in the right-most column of Figure 11, intensities differ mostly due to changes in the power-law indices of shock-accelerated spectra. These indices, shown in the accompanying frame below as $-2$, $-1.25$ and $-1$, correspond to those expected from Eq. 24 for $s=2$, 3, and 4, respectively.

The middle column of Figure 11 illustrates the case where both $V_{sh}$ and $s$ are varied so that $V_{2}^{\prime}/V_{1}^{\prime}=3$. It is worth drawing attention to the fact that the spectral indices these spectra display are not the index associated with $V_{2}^{\prime}/V_{1}^{\prime}=3$, even though this would be the jump factor in the flow speed observed by spacecraft during the passage of the shock. The power-law indices of all shock-accelerated spectra are consistently those associated with $s$. The observed factor by which flow speeds change across a shock should not be confused with the actual compression ratio, which is the factor by which either the density or the flow speeds in the shock frame change across the shock. This is important when calculating the power-law index of the DSA-predicted spectrum for comparison against an observed energy spectrum.

The spectra shown in the left and right-hand sides of Figure 11 extend to higher energies for both faster and stronger shocks, where $V_{sh}$ and $s$ are each varied separately with the other held constant. Spectra extend up to similar energies where $V_{sh}$ and $s$ are varied together. This is also illustrated in Figure 12, where the cut-off energies of the aforementioned spectra are plotted against both $V_{sh}$ and $s$, showing almost no reponse to changing shock conditions. Similar to Section 4.2, the acceleration times corresponding to each of the spectra in Figure 11 are calculated and compared to the simulation times, as shown in Figure 13. As before, the maximum attainable energies for shock-accelerated spectra are determined predominantly by the shock transit time.

5.1 Acceleration efficiency: Strong versus fast shocks

The preceding sections analyse the collated effects of varying $V_{sh}$ and $s$ during the complete acceleration process, including the heating of the SW energy distribution and the subsequent shock acceleration. It can be argued that the shock speed has the largest effect, since varying $V_{sh}$ yields pronounced changes in both the heated SW distributions and in reducing acceleration times, thereby extending spectra to higher energies. On the other hand, the most notable effect of varying $s$ is the changes in the power-law indices displayed by shock-accelerated spectra, which also affects energetic particle intensities considerably.

Implementing the same technique used in Section 4.3, particle fluences are calculated using Eq. 28. The resulting energy-integrated intensities are shown in the top frames of Figure 14 as functions of both $V_{sh}$ and $s$. To constrain the intensities resulting solely as a result of DSA and not including those of the input spectra, fluences are also calculated for the heated SW distributions ($I_{\text{kappa}}$) and subtracted from those calculated for the shock-accelerated spectra ($I_{\text{total}}$). These two additional sets of energy-integrated intensities, namely $I_{\text{kappa}}$ and $I_{\text{DSA}}=I_{\text{total}}-I_{\text{kappa}}$, are also shown in Figure 14. Note that $I_{\text{kappa}}$ is roughly an order of magnitude smaller than $I_{\text{total}}$ on average, with comparably weak dependences on $V_{sh}$ and $s$. $I_{\text{DSA}}$ therefore closely resembles $I_{\text{total}}$. Considering particle fluences as functions of $V_{sh}$ and $s$, where each is varied alone, it is apparent that both faster and stronger shocks yield greater numbers of energetic particles, whether considered for $E>E_{inj}$ or for $E>$ 1 MeV. Since compression ratios for shocks in the SW do not typically exceed 4, and there is no hard limit on the speed shocks can attain, faster shocks arguably have greater potential as particle accelerators.

It is insightful to consider fluences for the configurations where $V_{sh}$ and $s$ are varied together, since this is essentially comparing the acceleration efficiencies of strong slower-moving shocks against weaker fast-moving shocks. For these cases, fluences remain fairly constant for all combinations of $V_{sh}$ and $s$ for $E>E_{\text{inj}}$. However, if only the number of higher-energy particles is considered, that is, for $E>$ 1 MeV, the results are more interesting: As a function of compression ratio, the fluences increase as shocks become stronger, despite simultaneously becoming slower. As a function of the shock speed, these fluences become smaller for faster shocks for which the accompanying compression ratios are smaller. This implies that strong slower-moving shocks are able to produce a greater number of energetic particles than weak fast-moving shocks, at least within the parameter ranges considered and for shocks with $V_{2}^{\prime}/V_{1}^{\prime}=3$.

6 From the solar wind to energetic storm particles

In this section, the development of particles injected from the suprathermal SW into energetic particles associated with ESP events is considered in greater detail. The underlying question of how SW particles can be accelerated to the much higher energies at which ESPs are observed is visually represented in Figure 15. The intensities observed during the Halloween ESP event of 2003 October 29 are shown as an example of typical energy spectra observed during such events. Note this particular spectrum is reproduced well by the power law associated with $s=4$ (e.g. Giacalone, 2015).

The acceleration process is initiated with the heating of the SW energy distribution, where the $\kappa$-function describing it is allowed to change in response to changes in flow speed, density and temperature across the shock as shown in Section 2. The importance of this initial heating cannot be understated: consider, firstly, a spectrum shock-accelerated from the original, undisturbed SW distribution. It would likely not extend up to the energies of observed ESPs because of the limited time available for acceleration. Secondly, assuming the injection threshold for DSA is much larger than thermal energies, the larger intensities associated with the heated SW could allow shock-accelerated intensities to reproduce those of typical ESPs. Indeed, the bottom panel of Figure 15 illustrates that the heated SW distribution increases intensities of potential DSA seed particles by up to a factor of $10^{4}$, while increasing suprathermal tail intensities a hundredfold at least. Interestingly, observations suggest that peak fluxes during ESP events display some dependence on CME sheath temperature (Dayeh et al., 2018), while the broadening of SW distributions (which influence eventual intensities of shock-accelerated particles) is predominantly dictated by temperature changes across the shock.

6.1 Injection from the heated SW distribution

In the simplest terms, the injection threshold for DSA can be considered the energy at which particles should at least propagate upwind in order to make repeated shock crossings. This injection energy can be inferred by matching theoretically expected shock-accelerated spectra with observations (e.g. Neergaard Parker & Zank, 2012). It can also be estimated using the argument that the particle anisotropy must be small in order for DSA to be valid (Giacalone & Jokipii, 1999; Zank et al., 2006). For the typical values implemented in this study ($V_{1}^{\prime}=$ 600 km s*-1*, $s=3$), injection energies are in the order of a few keV for parallel shocks, for which injection from an unheated Maxwellian seed particle distribution would indeed be possible (Neergaard Parker & Zank, 2012). However, for quasi-perpendicular shocks, and in the absence of magnetic field-line wandering, it is estimated that $E_{\text{inj}}\gtrsim$ 0.1 MeV (see also Li, 2017). The prior heating of the SW distribution becomes especially useful in this situation. See also the simplified method used by Hu et al. (2017) to estimate injection energies.

The halo-CME shocks considered in this study are mostly quasi-parallel if a Parker (1958) spiral is assumed for the magnetic field. Here, $E_{\text{inj}}$ is varied between 10 and 60 keV, spanning an energy domain extending from the thermal peak to the suprathermal tail of the heated SW distribution, and bounded by the aforementioned estimates for parallel and quasi-perpendicular shocks. The effects of these varied injection thresholds on shock-accelerated spectra are shown in Figure 16. As before, Eq. 25 is fitted for each case and the parameters listed in Table 5. The reference configuration as specified in Table 3 is used, with $\lambda_{0}=$ 0.06 AU, $V_{sh}=$ 2400 km s*-1*, and $s=3$. The standard features are visible: the spectra are distributed according to $E^{-1.25}$ as expected for $s=3$ and rolls over exponentially at higher energies. Spectra appear to terminate at similar energies for the injection thresholds considered, as this is presumably governed by the acceleration time.

It can be seen from Figure 16 that injecting particles from the heated SW distribution yields high intensities in the energy domain typically associated with ESPs. In the case of $E_{\text{inj}}=$ 60 keV, the shock-accelerated intensities are quite similar to those observed for the ESP event included in Figure 16 as an example. Should the injection speed be defined as the minimum possible speed a particle needs to stay ahead of the shock, whilst following a Parker field line with a 45*∘* angle ahead of the shock nose, the injection speed can be written from focused transport theory as $v_{\text{inj}}=\lvert V_{1}^{\prime}-V_{sh}\rvert/\cos\left(45^{\circ}\right)$. For a fast halo-CME shock with $V_{sh}=2400$ km s*-1* and a more typical upstream flow speed of 400 km s*-1*, the injection speed at 1 AU in the upstream flow frame is $v_{\text{inj}}=$ 2828 km s*-1*, which corresponds to an injection energy of $E_{\text{inj}}\sim 50$ keV. This is reasonably similar to the 60 keV injection energy required to fit the presented observations using the heated $\kappa$-distribution as a source.

Note the presented simulations are not intended to reproduce observations exactly, nor is it posited that the spectrum observed during the presented ESP event is the unambiguous result of particles accelerated from the SW. The aforementioned results do however demonstrate that the prior heating of the SW plasma and energy distribution complements the shock acceleration process well and might even be necessary when considering shocks with large injection energies.

6.2 The evolution of shock-accelerated distributions during the shock passage

Particles of sufficient energy can be injected into the DSA process at any time during a CME shock’s passage between the Sun and the Earth. However, ESP events entail the local enhancement of particle intensities as viewed by spacecraft at Earth. The largest enhancement is naturally expected when the shock reaches Earth. Solving the time-backward SDEs for the reference configuration and different starting positions of the shock, the energy spectra and intensity profiles of shock-accelerated particles are simulated for the approach and aftermath of the shock’s passage at Earth. These are shown in Figures 17 and 18, respectively.

The time-evolution of energy spectra during the approach of the shock is considered first. From Figure 17, the most recognizable aspect is of course the power-law form of the spectrum at the time of the shock’s arrival at Earth, that is, where $r_{sh}=$ 1 AU, followed by a cut-off at higher energies. However, this spectrum appears notably different when the shock is just 0.01 AU away, which for $V_{sh}=2400$ km s*-1* is a few minutes before its arrival at Earth. The differences, which become more stark for larger distances between the shock and the Earth, include lower overall intensities and a downturn in the spectrum at low energies. The progressively lower intensities follow merely as a result of the source of energetic particles (i.e. the shock) being further away. The low-energy downturns are more severe when the shock is far from the Earth, but manifest as flattened spectra for smaller distances from the Earth. This likely follows because shock-accelerated particles are adiabatically cooled in the expanding SW while they propagate ahead of the shock toward Earth. These steep spectra with positive power-law indices at lower energies are known spectral characteristics of adiabatic energy losses (e.g. Moraal & Potgieter, 1982; Strauss et al., 2011).

A particularly interesting feature visible in Figure 17 is the flattening of spectra during the approach of the shock. Similar flattening has recently been reported in observations of proton spectra between 50 keV and 1 MeV prior to the passage of an IP shock at Earth (Lario et al., 2018). It is conceivable that this is the result of the competing effects of shock acceleration and adiabatic cooling of particles. Whereas shock acceleration tends to distribute particles according to $E^{\gamma_{s}}$, where $\gamma_{s}=-1.25$ for $s=3$ according to Eq. 24, cooling tends to force protons into a characteristic $E^{+1}$ spectrum at Earth. When the shock is nearer to Earth, the shock-accelerated component dominates, whereas particles transported from the shock while it is still further upwind have been cooled to a greater extent. Where these effects balance, the spectrum flattens. As a result of progressively stronger cooling, Figure 17 shows the flattened segments narrow and move to higher energies for larger distances between the shock and the Earth. Figure 18 shows intensity profiles similar to how the observations of Lario et al. (2018) are presented. The flat energy spectra can be discerned from the coinciding intensity profiles of 0.25 to 1 MeV particles. This energy range can be broadened in the simulations by specifying a lower injection energy or assuming a stronger shock, e.g. $s=4$, for which the shock-accelerated spectrum will be less steep. These flat segments are visible for shock positions up to at least 0.15 AU away from the Earth, that is, for nearly three hours prior to its arrival. Note the shock speed of 2400 km s*-1* implemented in these simulations is very fast. If a more typical shock speed of e.g. 1000 km s*-1* (e.g. Mäkelä et al., 2011) is assumed, this effect would be visible for a longer time prior to the shock’s arrival at Earth.

The intensity profiles shown in Figure 18 largely resemble that of typical ESP events: a large onset of particle intensities over a relatively short time is visible before the arrival of the shock, followed by a more gradual decline of intensities after it has passed. Note that the peaks of the simulated profiles, corresponding to the arrival of the IP shock at Earth, are not as sharp as observed intensities tend to show. These peaks are often associated with large anisotropies (e.g. Lario et al., 2005a), which can be more appropriately modelled using focused transport models (Zuo et al., 2011; le Roux & Webb, 2012). The TPE in this study is solved for an omni-directional distribution function (as discussed in Appendix A) and describes only the near-isotropic particle component.

6.3 Seed-particle energies and initial positions

It only remains to be investigated where and from which energies SW particles contributing to any given observational point are accelerated. To do this, the binning technique discussed in Section 3.1 is implemented. Firstly, pseudo particles are traced time-backwards from the observational point $(r_{\text{obs}},E_{\text{obs}})=(\text{1 AU, 1 MeV})$, whereafter the average particle amplitudes are calculated for each $(r,E)$-bin at $t^{\prime}=0$. As before, the simulation is run for the reference parameters listed in Table 3. The result is shown as a colour-scaled plot in Figure 19.

This plot essentially maps the relative contributions of particles from different initial positions and energies to the intensities of 1 MeV particles at Earth’s position at the time of the shock’s passage. As such, they reveal a great deal of insight into the probable original energies and locations of seed particles. Figure 19 shows the undisturbed distribution (that is, before the departure of the shock from near to the Sun) of SW particles that would eventually contribute to 1 MeV intensities at Earth. It is illustrated that the largest contributions are from particles initially located upwind from Earth and accelerated from energies ranging from tens to a few hundred keV. The lower limit of this range corresponds to the 60 keV injection energy implemented in this study. The contributions of particles either cooled to 1 MeV from higher energies or convected to Earth without incurring energy changes are minuscule by comparison. During the shock transit, the distribution of contributing particles will move closer to the observational point, until, at the time of the shock passage at Earth, only that point will be populated on the plot.

The much larger relative contribution of particles accelerated to $E_{\text{obs}}$ from lower energies illustrates that DSA at the travelling IP shock is the chief contributor to intensities during the simulated ESP event.

7 Discussion and conclusions

In this study, the acceleration of SW particles at halo CME-driven IP shocks is investigated. The thermal and suprathermal components of the SW velocity (or equivalent energy) distributions are collectively described using $\kappa$-functions. Furthermore, these SW distributions are transformed in response to simulated shock transitions in plasma properties upon which they depend, such as the flow speed, number density, and temperature. These transformed (or heated) distributions are consequently specified as source spectra, from which particles with sufficient energy can be injected into the DSA process at the shock. To model this acceleration process, SDEs equivalent to the Parker TPE are solved in a time-backwards fashion. Using the combined approach of the pre-injection heating of the SW distribution and DSA, simulations reveal a number of noteworthy results with regards to the particle acceleration at IP shocks and particularly of SW particles.

The simulations are shown to produce the classical spectral features of DSA. However, in cases where diffusion length scales are small relative to the shock width, shock-accelerated spectra do not display the spectral indices associated with the full compression ratio. Such situations may arise in the case of high levels of magnetic turbulence, which can decrease diffusion length scales, or in the case of broader compressions in the SW (e.g. Giacalone et al., 2002). At any rate, the softer power-law distributions that are attained as a result notably reduces the overall intensities of shock-accelerated particles. Furthermore, it is shown that the highest attainable energies of shock-accelerated spectra are limited by the transit time of the shock. Reduced diffusion length scales may serve to better confine particles near the shock, thereby allowing spectra to be accelerated to higher energies within the available time. Intensities of higher-energy particles ($E>$ 1 MeV) are shown to be particularly sensitive to this time limit.

The dependence of the acceleration process on shock properties such as its speed and compression ratio is also investigated. Fast shocks contribute appreciably to the heating of the SW distribution. While fast-moving shocks have shorter associated transit times, allowing less time for particle acceleration, larger shock speeds have a net positive effect on particle acceleration. Through their contribution to larger flow speeds they reduce diffusion length scales, improving particle confinement at the shock. They also reduce the time required for acceleration: since larger shock speeds imply larger differences of flow speeds across the shock, a particle scattered across it experiences a larger mean energy gain per crossing; the scattering centres can also be thought of converging on the shock at a higher rate for fast shocks.

The compression ratio, as measure of shock strength, also affects the magnitude of flow speed transitions across the shock, providing the actual factor by which both the shock-frame flow speed and number density change across the shock. In particular, larger compression ratios increase intensities of the SW distribution across the shock, since the $\kappa$-function is normalised to the number density. While it has a more modest effect on acceleration times than the shock speed, its effect on overall particle intensities is significant, since the spectral indices of the shock-accelerated spectra depend directly on the compression ratio. When strong and fast shocks are compared as particle accelerators, it is found that strong slower-moving shocks produce larger numbers of energetic particles than weak fast-moving shocks. The compression ratio, being directly associated with the steepness of accelerated energy distributions, is therefore identified to be the greater limiting factor between these two shock properties.

With regards to simulating ESP events, the prior heating of the SW distribution during the shock passage is found to complement the DSA process well: it provides greater intensities of potential seed particles for DSA, especially where large injection energies are considered, and allows shock-accelerated spectra to achieve large enough intensities at sufficiently high energies to reproduce typical ESP events. This result is consistent with observations reporting larger peak particle fluxes during ESP events for warmer CME sheath temperatures (Dayeh et al., 2018). Furthermore, simulations of shock-accelerated intensities at Earth reveal significant flattening of spectra forming ahead of the shock during its approach. This is found to result from shock-accelerated particles experiencing adiabatic cooling in the expanding SW while propagating toward Earth ahead of the shock. This provides a potential explanation for similar features recently reported in observations (Lario et al., 2018).

Finally, taking advantage of the time-backward tracing of phase-space density elements, with their flux contributions weighted according to particle amplitude, it is revealed that most SW particles contributing to intensities during a simulated ESP event are transported to Earth from upwind and are accelerated from energies ranging from tens to a few hundred keV. Due to the outsized fraction of particles accelerated to 1 MeV from lower energies, it can be concluded that DSA is indeed a chief contributor to intensities at Earth during ESP events.

It is ultimately concluded that with the combination of prior heating and shock acceleration, energetic particles can be accelerated from the suprathermal SW at IP shocks, and that these particles may make appreciable contributions to intensities during observed ESP events. Also, the reproduction of the observed flat energy spectra ahead of shock passages (Lario et al., 2018), resulting from the competing processes of DSA and adiabatic cooling, illustrates the advantage of studying particle acceleration in association with more general transport processes. It is furthermore demonstrated that the SDE approach, while shown in Appendix B to yield comparable results to finite-difference methods, can be used to reveal unique physical insights with regards to particle transport and acceleration at IP shocks. Along with its computational advantages, the SDE model developed in this study can be expanded to include more spatial dimensions, providing a means to explore more complicated shock geometries. Further research using this SDE approach may also include more refined characterisations of transport coefficients, the solution of focused transport equations to study large-anisotropy events, and the reproduction of various shock-related particle events, such as the acceleration of particles associated with co-rotating interactions regions.

PLP and RDS acknowledges the financial support of the South African National Research Foundation (NRF). Opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the NRF.

Appendix A Counting particles

A.1 The distribution function and number density

Consider a number of particles $d\mathcal{N}$ in a small volume $d^{3}x$ at a position $\vec{x}=(x,y,z)$ and with momentum in $d^{3}p$ around $\vec{p}=(p_{x},p_{y},p_{z})$ as illustrated in Figure 20 (left). The distribution function $f\left(\vec{x},\vec{p},t\right)$ is introduced such that

[TABLE]

This can be related to the ordinary number density $n$ by integrating over all momentum space,

[TABLE]

The above volume integral can be re-written by considering a momentum coordinate system as shown in Figure 20 (right). Assuming the gyro-centre at the origin and the magnetic field line about which the particle gyrates along $p_{z}$, that is $\vec{B}=B\hat{z}$, the azimuthal angle $\phi$ represents the particle’s gyrophase and the polar angle $\theta$ its pitch angle.

This allows the momentum-space volume element to be expressed in spherical coordinates as $d^{3}p=p^{2}\ d\Omega\ dp$, where $d\Omega=\sin\theta\ d\theta\ d\phi$ is a solid-angle element about $\vec{p}$ and $p=|\ \vec{p}\ |$ is the particle’s scalar momentum. Noting that $\vec{p}$ can hence be expressed in terms of the coordinates $(p,\ \phi,\ \theta)$, Eq. A2 can be rewritten as

[TABLE]

where $\int_{\Omega}$ represents the integral over all solid angles. However, since particle detectors typically cannot resolve statistically significant numbers of particles incident from only one particular direction, the average number of particles per unit solid angle in momentum space is considered. Accordingly, the omni-directional distribution function (that is, averaged over all directions in $\vec{p}$-space for a fixed $p$) is defined as

[TABLE]

since $\int_{\Omega}d\Omega=\int\limits_{0}^{2\pi}d\phi\int\limits_{0}^{\pi}\sin\theta d\theta=4\pi$. Note that $f_{0}\left(\vec{x},p,t\right)$ is the distribution function solved for in the Parker (1965) transport equation. Note furthermore from Eq. A4 that $f\left(\vec{x},\vec{p},t\right)$ can be taken out of the integral over solid angle if it is assumed to be gyro- and isotropic, that is, independent of both $\phi$ and $\theta$. In this case it follows that $f\left(\vec{x},\vec{p},t\right)=f_{0}\left(\vec{x},p,t\right)$. From Eqs. A.1 and A4 the number density can be expressed as

[TABLE]

where $U_{p}\left(\vec{x},p,t\right)=4\pi p^{2}\ f_{0}\left(\vec{x},p,t\right)$ is known as the differential number density. Here, $U_{p}$ represents the number density between two spherical shells in $\vec{p}$-space with radii of $p$ and $p+dp$, respectively. While spacecraft cannot resolve volume densities, $U_{p}$ does prove useful to relate quantities they measure to those introduced above.

A.2 Toward the differential intensity

Spacecraft essentially measure the flux of particles observed for a particular viewing direction within a particular momentum (or energy) band. This is known as the differential intensity $j_{p}$, and is the number of particles detected per momentum interval per unit time per unit area of observed space per unit solid angle. Dimensionally, this can be attained by the product of the differential number density and the speed at which particles move toward the detector through a surface perpendicular to their motion, averaged over all solid angles. That is,

[TABLE]

which, in turn, relates the differential intensity to the distribution function. It is furthermore useful to express the differential intensity in terms of kinetic energy to compare with spacecraft observations. The conversion of $j_{p}$, expressed per momentum interval $dp$, to $j_{E}$, expressed per interval of kinetic energy $dE$, is carried out by noting that $j_{p}\ dp=j_{E}\ dE$ due to the conservation of particles. It hence follows that

[TABLE]

Note that since $p^{2}=\left(E_{T}^{2}-E_{0}^{2}\right)/c^{2}$, where $E_{T}=E+E_{0}$ with $E_{T}$, $E$ and $E_{0}$ denoting the total, kinetic and rest-mass energies, respectively, it follows that

[TABLE]

where $\gamma$ is the Lorentz factor and $m_{0}$ is the rest mass. Then, using Eq. A8, $j_{E}$ can be related to the distribution function as follows

[TABLE]

A.3 On kappa velocity distribution functions

The normalization constant $A_{\kappa}$ of the standard $\kappa$-function introduced in Eq. 1 is obtained by setting its integral over all phase space equal to the number density of the solar wind. Assuming an omni-directional (or three-dimensional) $\kappa$-function, averaged as per Eq. A4, the expression given in Eq. A5 can be used to calculate the normalization constant. Note, however, because $f_{\kappa}$ is often expressed as a function of velocity instead of momentum, the analogous expression in velocity coordinates, namely

[TABLE]

is used to calculate the number density. Bear in mind that for use in the transport equation, e.g. when specifying $f_{s}$ in the source function, the conversion $f_{s}=f_{\kappa}/m_{p}^{3}$ applies, because $\int 4\pi p^{2}\ f_{s}\ dp=\int 4\pi(m_{p}v)^{2}\ f_{s}\ d(m_{p}v)\equiv\int 4\pi v^{2}\ f_{\kappa}\ dv$. Here, only the proton rest mass $m_{p}$ is used, since this study is concerned with solar wind protons with kinetic energies typically much smaller than their rest-mass energy. A non-relativistic description is therefore applicable. Note furthermore that Eq. A11 is also known as the zeroth-order velocity moment of $f_{\kappa}$, and that it can be cast into a more general form for the $n^{\text{th}}$-order moment, namely

[TABLE]

which, upon substituting Eq. 1, becomes

[TABLE]

from which the normalization constant $A_{\kappa}$ can be calculated. The integral in Eq. A13 is evaluated by applying the transformation $v^{2}/(\kappa v_{\kappa}^{2}):\rightarrow\eta$, which allows the expression to be rewritten as

[TABLE]

Note that the integral is now in the form of the $\beta$-function, given by

[TABLE]

Recognizing that $x=(n+3)/2$, $y=\kappa-(n+1)/2$, and $x+y=\kappa+1$, the integral in Eq. A14 can be solved by invoking the identity given in Eq. A15, so that

[TABLE]

from which the normalization constant can be determined as

[TABLE]

by setting $\langle v^{n=0}\rangle=n_{sw}$, yielding the same expression given in Eq. 3. The standard kappa function is hence given by

[TABLE]

Appendix B Equivalence of SDE and finite-difference numerical methods: Application to an event during the 2003 Halloween epoch

In a similar study, where DSA is also assumed as an acceleration mechanism in simulating ESP events, Giacalone (2015) utilizes a finite-difference numerical scheme to solve the Parker (1965) TPE in a single spatial dimension for a strong fast-moving shock ($V_{sh}=$ 1900 km s*-1*, $s=4$). In that study, the energetic particle event, or Halloween event, of 2003 October 29 is considered as an application for which observed energy spectra and temporal profiles are reproduced. To demonstrate their equivalence as numerical methods, the SDE approach of this study is implemented to reproduce particle intensities for the same event using a similar parameter configuration to that of Giacalone (2015). The source function in that study is bound to a numerical grid, which does not exist in the SDE approach. For the purposes of this application, the source function is instead specified at the shock as a very soft power law in momentum ($p^{-7}$), where simulated intensities are retroactively normalised to observed intensities. Furthermore, the diffusion coefficient is implemented as specified in that study, that is, of order $10^{19}$ cm2 s*-1* near the Sun and with a momentum dependence of $\sim p^{1.5}$ (refer to Giacalone, 2015, for further details).

The resultant energy spectra and time profiles of energetic particle intensities are shown in Figure 21. Overall, the finite-difference and SDE solutions are similar. Recall that the SDEs must be solved for a large number of pseudo particles for each energy, position or point in time at which intensities are sought. The resolution of SDE solutions therefore merely depends on the number of times the simulation is repeated and how the observational points are distributed, as opposed to the fixed grid resolution of finite-difference schemes. For example, in Figure 21, a greater number of more tightly spaced observational points are chosen in energy than in time, as shown in the left and right-hand side panels, respectively. Note that the SDE solutions tend to undershoot the finite-difference solutions where steeper gradients exist. This follows because there is no passage of information between pseudo particles and their associated observational points as would exist between adjacent points in finite-difference schemes. The SDEs typically yield less gradual intensity gradients as a result. This is useful to bear in mind when fitting functions to shock-accelerated spectra simulated using these two different approaches, since the high-energy cut-off is likely to be more abrupt for the SDEs than for finite-difference solutions.

Note, as an aside, that simulated intensities of either method reproduce observations better near the shock than away from it. This follows since the only source of particles in the simulations is situated on the shock itself, while the observed intensities include contributions from other sources. In any event, the comparison of the simulations generated using the two different numerical methods demonstrate the equivalence of their results.