Solar energetic particle propagation in wave turbulence and the possibility of wave generation

TL;DR

This paper models the propagation of solar energetic electrons in turbulent magnetic fields, incorporating dynamical effects and wave growth, to better understand particle transport and turbulence interactions near the Sun.

Contribution

It introduces a physics-first simulation approach with dynamical perpendicular diffusion and wave growth effects, advancing the understanding of particle-turbulence interactions.

Findings

Reasonable qualitative agreement with observations

Wave growth significantly alters turbulence levels near the Sun

Inclusion of dynamical effects improves propagation modeling

Abstract

A complete theory for the complex interaction between solar energetic particles and the turbulent interplanetary magnetic field remains elusive. In this work we aim to contribute towards such a theory by modelling the propagation of solar energetic particle electrons in plasma wave turbulence. We specify a background turbulence spectrum, as constrained through observations, calculate the transport coefficients from first principles, and simulate the propagation of these electrons in the inner heliosphere. We have also, for the first time, included dynamical effects into the perpendicular diffusion coefficient. We show that such a "physics-first" approach can lead to reasonable results, when compared qualitatively to observations. In addition, we include the effect of wave growth/damping due to streaming electrons and show that these particles can significantly alter the turbulence levels close to the Sun for the largest events.

Equations

$$\vec{B}=\overrightarrow{B_0}+\vec{b}(x,y,z).$$

$$\vec{b}(x,y,z)={\vec{b}}_{\mathrm{slab}}(z)+{\vec{b}}_{2D}(x,y)$$

$$\delta B^2=\delta B_{\mathrm{slab}}^2+\delta B_{2D}^2.$$

$$b_{\mathrm{slab},x}(k_{||},t)=b_{\mathrm{slab},x}(k_{||})e^{-i\omega t+{\gamma }^{\prime }t}.$$

$$\delta B_{\mathrm{slab}}^2(t)$$

$${\gamma }^{\prime }={\tau }^{-1}+\gamma .$$

$${\tau }^{-1}=\alpha \left|\omega \right|=\alpha \left|V_A{k}_{||}\right|,$$

$$\frac{1}{g^{\mathrm{slab}}(k_{||},t)}\frac{\partial g^{\mathrm{slab}}(k_{||},t)}{\partial t}\approx 2\gamma (k_{||},t).$$

$$\delta B_{\mathrm{slab}}^2=2\pi {\int }_0^{\infty }\sum _{j=\pm 1}\sum _{n=\pm 1}{g}_{n,j}^{\mathrm{slab}}(k_{||})d{k}_{||}.$$

$${\sigma }_c\left(k_{||}\right)=\frac{{\sum }_{n=\pm 1}{g}_{n,j=+1}^{\mathrm{slab}}-{\sum }_{n=\pm 1}{g}_{n,j=-1}^{\mathrm{slab}}}{{\sum }_{n=\pm 1}{g}_{n,j=+1}^{\mathrm{slab}}+{\sum }_{n=\pm 1}{g}_{n,j=-1}^{\mathrm{slab}}},$$

$${\sigma }_m\left(k_{||}\right)=\frac{{\sum }_{j=\pm 1}{g}_{n=+1,j}^{\mathrm{slab}}-{\sum }_{j=\pm 1}{g}_{n=-1,j}^{\mathrm{slab}}}{{\sum }_{j=\pm 1}{g}_{n=+1,j}^{\mathrm{slab}}+{\sum }_{j=\pm 1}{g}_{n=-1,j}^{\mathrm{slab}}}.$$

$$b_{2D,x}(k_{||},t)=b_{2D,x}(k_{||})e^{-i(V_A{k}_{\perp })t},$$

$$\frac{1}{g^{2D}(k_{\perp },t)}\frac{\partial g^{2D}(k_{\perp },t)}{\partial t}\approx 0,$$

$$\delta B_{2D}^2$$

$$k_d^{\mathrm{LH},\mathrm{RH}}\approx \frac{|\Omega {|}_{p,e}}{V_A+3v_{\mathrm{th}}^{p,e}},$$

$${\omega }^{\mathrm{LH},\mathrm{RH}}\in |\Omega {|}_{p,e}\pm 3k_{||}{v}_{\mathrm{th}}^{p,e}$$

$$T_{p,e}=T_{p,e}^0{\left(\frac{r}{r_0}\right)}^{-{\chi }_{p,e}}$$

$${\sigma }_c=\exp \left[-\frac{r}{r_0}\right],$$

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

$$L^{-1}=\nabla \cdot \hat{b},$$

$$A(r,\varphi ,t)=3\frac{{\int }_{-1}^{+1}\mu fd\mu }{{\int }_{-1}^{+1}fd\mu },$$

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

$$f(r$$

$$D_{\mu \mu }=\frac{\pi v^2(1-{\mu }^2)}{B_0^2{r}_L^2}{\int }_0^{\infty }d{k}_{||}\sum _{n=\pm 1}\sum _{j=\pm 1}{\mathcal{R}}_{n,j}{g}_{n,j}^{\mathrm{slab}},$$

$${\mathcal{R}}_{n,j}=\frac{{\tau }^{-1}}{{\left({\tau }^{-1}\right)}^2+{\left(v\mu j{k}_{||}-\omega +n\Omega \right)}^2},$$

$${\lambda }_{||}=\frac{3v}{8}{\int }_{-1}^{+1}\frac{{\left(1-{\mu }^2\right)}^2}{D_{\mu \mu }}d\mu .$$

$${\mathcal{R}}_{n,j}^{\mathrm{ND}}\approx \pi \delta \left(v\mu j{k}_{||}-\omega +n\Omega \right),$$

$$D_{\mu \mu }^{\mathrm{ND}}=\frac{{\pi }^{2}v(1-{\mu }^2)}{B_0^2{r}_L^2|\mu |}\sum _{n=\pm 1}\sum _{j=\pm 1}{g}_{slab}^{n,j}(k_{||}^{\mathrm{res}}),$$

$$k_{||}^{\mathrm{res}}=\frac{n\Omega }{V_A-v\mu j}.$$

$$D_{\perp }^{\mathrm{dyn}}\approx ⟨l_{\perp }⟩\left(v|\mu |+V_A\right){\frac{\left(\delta B_{2D}^2\right)}{B_0}}^{1/2},$$

Full text

cosmic rays — diffusion — Sun: heliosphere, particle emission — turbulence

\AuthorCallLimit

=3 \fullcollaborationNameThe Friends of AASTeX Collaboration

1 Introduction

During transient solar phenomena, such as solar flares and/or coronal mass ejections, high energy solar energetic particles (SEPs) are accelerated (e.g. Reames, 2013, 2015), where after these particles propagate in the turbulent heliospheric magnetic field to reach Earth and the inner heliosphere. A complete description of the interaction (i.e. scattering) between particles and these turbulent fields remains elusive. One the one hand, more observations (of both SEP particles, but also of the turbulent magnetic fields) near the Sun are needed, while, on the other hand, more detailed modelling of SEP propagation is needed on a more fundamental level. We advocate a so-called “physics-first” approach in modelling SEPs: specify the background turbulence quantities as informed by observations and/or simulations and calculate the relevant transport quantities (i.e. diffusion coefficients) by using a realistic turbulence spectrum as input. Finally, a comparison between modelled and SEP intensities can be performed. To some extent this is a continuation of our effort described in Strauss et al. (2017). This approach is supplemented by empirical modelling where the transport coefficients are treated as free parameters and fitted to particle observations (e.g. Dröge et al., 2016).

In this study we perform further refinements of our model and again apply it to the transport of impulsively accelerated, near-relativistic ($\sim 100$ keV) electrons. We therefore focus on so-called electron rich events, where the SEP proton component can, for the most part, be neglected. We calculate the transport coefficients for dynamical turbulence and also include previous neglected effects such as those that arise from a non-vanishing magnetic and cross helicity. We also attempt to quantify the background turbulence in more detail, specifically the wave-number at which the dissipation of turbulence starts being significant. In addition to the background (or rather, omni-present) solar wind turbulence, we also investigate whether SEP electrons can modify the background turbulence via the streaming instability: can a SEP event amplify/damp solar wind turbulence? Although this has been modelled and observed before for proton events (e.g. Ng et al., 1999; Desai et al., 2012), we are unaware of any such simulations for electrons. Moreover, such a study is additionally motivated by the recent findings of Agueda & Lario (2016) that the derived mean-free-path for some events may depend on the intensity of the SEP events; an indication that more particles can lead to enhanced levels of scattering and hence a decrease in the mean-free-path.

2 Solar wind fluctuations

Magnetic turbulence is usually described by Reynolds decomposing the magnetic field, $\vec{B}$, into a locally uniform background field, $\vec{B_{0}}$, and a random turbulent component, $\vec{b}(x,y,z)$, such that

[TABLE]

Averaging over long time periods, yields $\langle\vec{B}\rangle=\vec{B_{0}}$, while $\delta B^{2}=\langle\vec{b}^{2}\rangle$ is the variance of the fluctuating component. Furthermore, we assume the turbulence to be transversal, i.e. $\vec{B_{0}}\cdot\vec{b}=0$ (this is equivalent to assuming a strong guide field, i.e. if $\vec{B_{0}}=B_{0}\mathbf{z}$, then the $\mathbf{z}$ component of $\vec{b}$ is negligible, $\vec{b}\cdot\mathbf{z}\ll B_{0}$) and that the energy contained in the turbulence is much less than that of the background field, $\delta B^{2}\ll B_{0}^{2}$.

2.1 Nature of the assumed turbulence

Following Shalchi (2009), we assume the fluctuating field consists of a slab (with fluctuations directed along the mean field, $k_{||}$) and 2D (with fluctuations directed perpendicular to the mean field, $k_{\perp}$) component (see also Matthaeus et al., 1995), so that the fluctuations can be described as

[TABLE]

For the total variance of the fluctuations, it follows that

[TABLE]

We work in terms of the so-called wave approach, where the slab fluctuations are considered to consist of a spectrum of (magneto-hydrodynamic) MHD waves with a wave frequency of $\omega$ and a growth/damping rate of $\gamma^{\prime}$ (so-called plasma wave turbulence; Schlickeiser, 2002). In this model, we can simply specify any component of the turbulence, in wave-number space, as

[TABLE]

The strength of the background slab component is then calculated as

[TABLE]

where $g(k_{||})$ is the so-called turbulence power spectrum and $\langle\cdot\rangle$ indicates an appropriate averaging procedure.

For this work, we consider only non-dispersive, parallel propagating and circularly polarized Alfvén waves, where $\omega=V_{A}k_{||}$, and $V_{A}$ is the Alfvén speed. This assumption is valid for near-relativistic electrons ($\sim 100$ keV) under consideration here: in Fig. 1, $k_{\mathrm{res}}\sim r_{L}^{-1}$ ($r_{L}$ being the maximal Larmor radius) indicates the maximum approximate wave-number at which these electrons will resonate, indicating that this is indeed in the Alfvén branch of each node.

It takes some care to interpret the growth/damping rate present in Eq. 4. We assume this rate to consist of two contributions; the first is the growth/damping rate due to the streaming and non-thermal SEPs, given by $\gamma$, where, of course, $\gamma=0$ when there are no streaming particles. The dynamical character of the “background turbulence” (that is, the quasi-stationary turbulence that is present in the undisturbed solar wind before the SEP event occurs), however, contributes to fluctuations changing at a dynamical timescale of $\tau$. Hence, we use

[TABLE]

Following Bieber et al. (1994), we choose the dynamical timescale of the turbulence to be

[TABLE]

where $\alpha=1$ is chosen throughout for maximal dynamical effects. The SEPs are therefore introduced into a fully evolved turbulent medium in an equilibrium state. This leads to time-dependent changes in the level of the slab component to be solely due to wave growth/damping by non-thermal SEPs, i.e.

[TABLE]

Assuming that the gyrotropic slab wave field includes forward ($j=+1$) and backward ($j=-1$) propagating waves with either a right-handed ($n=-1$) or left-handed polarization ($n=+1$), Eq. 2.1 reduces to

[TABLE]

The ratios of the different wave modes determine the cross helicity (ratio of forward to backward propagating waves),

[TABLE]

and the magnetic helicity (ratio of left hand to right hand polarized waves)

[TABLE]

For the 2D component we again include dynamical effect by assuming these are caused by perpendicular propagating Alfvén waves, so that we may write

[TABLE]

where, similar to the slab component, the 2D fluctuations are in a dynamic, but quasi-steady state,

[TABLE]

with the total strength thereof calculated as

[TABLE]

2.2 Onset and importance of the dissipation range

In the inner heliosphere, low energy cosmic ray electrons resonate with the slab turbulence spectra near the so-called dissipation range (Engelbrecht & Burger, 2015). This is where the cascading waves start being heavily damped by e.g. a cyclotron resonance with the thermal plasma particles (e.g. Woodham et al., 2018). The assumed form of this part of the spectrum is therefore a crucial modelling component to correctly simulate the propagation of these electrons.

Our assumed slab turbulence spectrum is shown in Fig. 3. This form is exactly the same as used by Strauss et al. (2017), but with a modified value of $k_{d}$, the wavenumber corresponding to the onset of the turbulence dissipation range, based on the model of Engelbrecht & Strauss (2018), given by

[TABLE]

where we assume that left handed waves are damped by gyroresonant interactions with thermal protons and right handed waves interact with thermal electrons. The regions where this interaction is believed to be very effective, are contained within

[TABLE]

for LH and RH waves separately (see also the simulations by Schreiner et al., 2017) and indicated on Fig. 1 as the grey bands. The wavenumber where the grey bands intersect the non-dispersive wavemodes, are believed to be where wave damping becomes significant, and is given by Eq. 15.

To use the expressions above, we also need an approximation for the electron and proton temperatures in the solar wind. These are assumed polytropic,

[TABLE]

with $r_{0}=1$ AU and $T^{0}_{p,e}$ a normalization value at Earth. Following Cranmer et al. (2009), we use $T^{0}_{p}=30^{5}$ K and $T^{0}_{e}=10^{5}$ K. The proton temperature is assumed to decrease adiabatically, $\chi_{p}=4/3$, while the electrons decrease is assumed to be more isothermal, $\chi_{e}=1/3$ (e.g. Sittler & Scudder, 1980; Phillips et al., 1995). The resulting temperature profiles are shown in the left panel of Fig. 2.

Using these assumed temperature profiles, the resulting dissipation range onset are shown in the right panel of Fig. 2 as a function of radial distance. Note that the dissipation onset that corresponds to LH polarized waves (and damped by protons) is always much smaller than the resulting quantity for RH polarized waves. This difference results in different forms of the slab turbulence spectrum; see Fig. 3 where the LH-polarized inertial range (coloured blue in the figures) is much more narrow than the RH-polarized inertial range (the combined blue and red region). This will result in a non-zero magnetic helicity (shown in the top panel of Fig. 4) with electrons scattered, in this region, predominantly by RH polarized waves.

2.3 Magnetic and cross helicities

The magnetic helicity is calculated from the turbulence spectra discussed in the previous paragraph and shown in the top panel of Fig.4. At high wave-numbers this value becomes non-zero and indicates, for our assumptions, an excess of RH polarized slab waves.

We model the background cross-helicity by a simple exponential function,

[TABLE]

which is independent of wavenumber and with $r_{0}=1$ AU. This assumed form is shown in the bottom panel of Fig. 4 along with the observations of Roberts et al. (1987a, b) (blue dots) and Breech et al. (2005) (green dots). This form is motivated by the stimulations of e.g. Adhikari et al. (2015) that shows how the highly anisotropic slab wave-field at the Alfvén radius (only forward moving Afvén waves can escape this surface and be convected by the solar wind into the inner heliosphere) can become increasingly isotropic near Earth’s radius due to e.g. wave-wave interactions.

Although seemingly insignificant, the effects of a non-zero magnetic and cross-helicity has important implication for particle transport as discussed later on. The upcoming missions, Parker Solar Probe and Solar Orbiter, will be able to measure these quantities in detail.

3 Numerical transport model

The propagation of energetic electrons is described by the so-called focussed transport equation (e.g. Skilling, 1971), given by

[TABLE]

and solved by means of the numerical approach outlined by Strauss & Fichtner (2015) to yield the gyro-tropic particle distribution function $f$. In Eq. 19, $\hat{b}$ is a unit vector pointing along the mean heliospheric magnetic field, $v$ is the particle speed, $\mu$ is the cosine of the pitch-angle, $\mathbf{D}^{(x)}_{\perp}$ contains the perpendicular diffusion coefficients and is specified in spherical coordinates (radial distance, $r$, and azimuthal angle, $\phi$), $D_{\mu\mu}$ is the pitch-angle diffusion coefficient and the focusing length is calculated as

[TABLE]

for a Parker (1958) heliospheric magnetic field (HMF), normalized to 5 nT at Earth. The first order anisotropy is calculated as

[TABLE]

and the omni-directional distribution function as

[TABLE]

From $F$, the isotropic differential intensity is calculated as $j=p^{2}F$, with $p$ being particle momentum.

As an inner boundary condition, the following mono-energetic isotropic injection function

[TABLE]

is prescribed at the inner boundary, located at $r_{0}=0.05$ AU. Gaussian injection in $\phi$ is assumed with $\phi_{0}=\pi/2$ and $\sigma$ determining the broadness thereof. The model is solved for $E_{0}=100$ keV. We use a value of $\sigma=20^{\circ}$ for the simulations presented here. It is believed that magnetic structures near the Sun (e.g. Klein et al., 2008) can significantly broaden the compact source regions from a solar flare below the Alfvénic radius. We can account for this broadening in our model by specifying $\sigma$, although the simulations of Strauss et al. (2017) have shown that the modelled intensities are relatively insensitive to the choice of this parameter as long as significant levels of perpendicular diffusion is present. A Reid-Axford (Reid, 1964) temporal injection profile is specified with $\tau_{a}=1/10$ hr, $\tau_{e}=1$ hr, and $C$ is a parameter that is adjusted to give the correct intensity at Earth. Other important quantities assumed in the model are: a solar wind number density of 5 particles.cm*-3* at Earth decreasing as $r^{-2}$ and a constant solar wind speed of $V_{sw}=400$ km.s*-1*.

The pitch-angle and perpendicular diffusion coefficients, needed as input to Eq. 19, are discussed and calculated in the next section.

4 Transport coefficients

In this work we try to to use, for as far as possible, a “physics-first” approach to simulation SEP transport. That is, we attempt to specify, as self-consistently as possible, the background turbulence spectra, and calculate the SEP transport coefficients from these assumed forms.

4.1 Pitch-angle diffusion

The plasma wave pitch-angle diffusion coefficient is given by e.g. Schlickeiser (2002) as

[TABLE]

where we have assumed that $V_{A}\ll v$. The resonance functions are given by

[TABLE]

with $\gamma\ll\tau^{-1}$; the broadening of the resonance function is therefore assumed to be entirely due to the dynamic character of the background fluctuations. This is the form of $D_{\mu\mu}$ that will be used for the rest of the study.

We can use this expression to see when electrons will resonate with what type of slab wave (i.e. when is $\mathcal{R}_{n,j}>0$). This is shown in Fig. 5, where the shaded regions indicate possible resonance with either forward (red regions) or backward (blue regions) propagating waves with either polarization state. The vertical dashed lines show the classical resonance gap which is present when either only LH-polarized waves are assumed, or if the waves are assumed static. Note that when a beam of electrons is present (i.e. $\mu>0)$, these particles will strongly resonate with forward propagating LH waves and backward propagating RH polarized waves.

The resulting $D_{\mu\mu}$ is shown in the left panel of Fig. 6 at three different radial positions. Note that this coefficient is not symmetric about $\mu\sim 0$, but enhanced for backwards moving particles. This is a direct result of the non vanishing magnetic and cross-helicities assumed in the model; there is an excess of RH-polarized, forward propagating waves in the model, as therefore backwards propagating particles experience more scattering. Refer again to Fig. 5. Near Earth, $D_{\mu\mu}$ becomes more symmetrical as the cross-helicity decreases to zero.

The parallel mean free path (MFP, $\lambda_{||}$) is calculated from $D_{\mu\mu}$ following the usual definition of Hasselmann & Wibberenz (1968),

[TABLE]

The resulting $\lambda_{||}$ is shown in Fig. 7 as a function of radial distance. Note the strong radial dependence of $\lambda_{||}$ and the large values thereof close to the Sun.

If we, however, consider the case of negligible damping (ND) or dynamical effects, the resonance function reduces to

[TABLE]

and the pitch-angle diffusion coefficient becomes

[TABLE]

where we have introduced the resonant wavenumber as

[TABLE]

This approximate form of the pitch-angle diffusion coefficient, $D^{\mathrm{ND}}_{\mu\mu}$, will be used later on to derive the wave damping/growth rate.

4.2 Perpendicular diffusion

No consistent theory exists to describe perpendicular diffusion, on the pitch-angle level, in the presence of dynamical turbulence. However, we may include the dynamic contribution of propagating Alfvén waves (or rather, particles drifting in the random electric fields induced by the propagating Alfvén waves) by using the approximate result of Strauss et al. (2016)

[TABLE]

which is valid for the limit of $r_{L}\ll\langle l_{\perp}\rangle$, where $\langle l_{\perp}\rangle$ is, approximately, the perpendicular correlation scale. In the magnetostatic limit of $V_{A}=0$, Eq. 30 should, however, converge to the field-line random walk (FLRW) model of Jokipii (1966), where

[TABLE]

with $a^{2}\in\left[0,1\right]$ a parameter determining the probability of particles being stuck to wandering field-lines, and the field-line diffusion coefficient $\kappa_{\mathrm{FL}}$ given by e.g. Qin & Shalchi (2014) as

[TABLE]

Making the rather ad hoc identification,

[TABLE]

which is remarkably similar to the form used by Matthaeus et al. (2007), but given, by these authors, in terms of the ultrascale, we use the following estimate for the dynamical perpendicular diffusion coefficient,

[TABLE]

The value of $a$ is changed in later sections to illustrate the role of perpendicular diffusion in wave amplification, while $\alpha$ was introduced for consistency with Eq. 7.

The isotropic perpendicular mean free path, is calculated as

[TABLE]

Both $D_{\perp}$ and $\lambda_{\perp}$ are shown in Figs. 6 and 7. Note that the inclusion of dynamical effects, although small, leaves to a non-zero value of $D_{\perp}$ near $\mu\sim 0$. This effect may, however, become increasingly dominant at lower energies where $v\sim V_{A}$.

It should also be noted that if $g^{\mathrm{2D}}(k_{\perp})$ is interpreted as the probability of finding a fluctuation at a given wave-number with a given energy (amplitude), then the characteristic length scale (which can be identified as the correlation length) can be evaluated as a root-mean-squared average scale using the second-order moment of the distribution

[TABLE]

which leads to Eq. 33 after some re-arrangement.

5 Wave growth/damping and its numerical implementation

The basic idea behind wave growth/damping is illustrated in Fig. 8: a distribution of particles with speed $v_{1}$ is isotropic with respect to the bulk flow frame (or, fixed frame) denoted by $v_{||}=0$ in the figure, i.e. their pitch-angles are equally populated so that they form a shell distribution (the grey half-circle on the figure). Forward propagating Alfvén waves (with a speed of $v_{||}=V_{A}$) are now introduced into the system and they interact with the particles through gyro-resonance which leads to pitch-angle scattering. The Alfvén waves will now attempt to isotropize the distribution with respect to the wave frame (the dashed grey half-circle on the figure). Note that the energy of particles, as observed in the wave frame, is conserved during the interaction, $v^{\prime}_{1}=v_{1}$. However, the particle energy, as measured in the flow frame, decreased to $v_{2}<v_{1}$ as the particle moved from $a\rightarrow c$, with the scattering event at $b$ (although pitch-angle scattering in this sense will be a continuous process). Because of the conservation of energy, the energy that the particle lost was gained by the Alfvén waves responsible for that scattering event, and hence we may state that forward propagating particles ($v_{||}>0$) that interact with forward moving Alfvén waves ($v_{||}=v_{A}>0$) will loose energy during gyro-resonant interactions, and will in turn grow and/or amplify forward moving Alfvén waves. A similar sketch and argument can be found in Afanasiev & Vainio (2013).

Of course, the opposite interaction is also possible: backward propagating particles that interact with forward moving Alfvén waves will gain energy during gyro-resonant interactions, and will in turn damp the forward moving Alfvén waves. If the particle distribution is isotropic (equal number of forward and backward propagating particles), the net effect will be that the interactions average out, and no net energy change (for either the particles or the waves) will take place.

A more rigorous treatment of the wave growth/damping problem is given in Appendix A and leads to the result

[TABLE]

in terms of the resonant wavenumber ($k_{||}^{\mathrm{res}}$) and pitch angle ($\mu_{\mathrm{res}}$), and the energy contained in the background magnetic field, $\mathcal{B}$. By looking at this result, it is clear that for a isotropic distribution, $\partial f/\partial\mu\approx 0$, we should have $\gamma\approx 0$ and no net wave amplification will occur. Moreover, we should note that the growth rate is generally expected to be low, due to the presence of the factor $V_{A}/v\ll 1$. But, as we will show, wave growth may still reach appreciable levels if enough streaming particles are present, and if the amount of scattering is sufficient (this is quantified by the plasma quantities present in Eq. 37), as the case $D_{\mu\mu}\approx 0$ will lead to $\gamma\approx 0$ (the particles cannot interact with the waves; see e.g. Eq. A6). Lastly, we note that the sign of $\gamma$ is determined by the wave propagation direction, $j$ in accordance with the sketch presented in Fig. 8.

In the numerical model, at each timestep, the numerical derivative $\partial f/\partial\mu$ is calculated, which allows us to calculate the growth/damping rate as given by Eq. 37. The different wave amplitudes are then adjusted accordingly by calculating

[TABLE]

or, approximated in its discretized form as

[TABLE]

where the numerical timestep is $\Delta t=t-t_{0}$. The time-dependent turbulent levels are then used to calculate time-dependent transport coefficients, which are then incorporated in the time-dependent simulations. A self-consistent description is therefore obtained. For the present study, transport effects of the newly generated waves, in both physical and wave-number space, are neglected. This is partly motivated by the fact that the SEP electrons propagate relativistically from the Sun to Earth and that the distribution peaks within minutes, while it is unlikely that any newly generated waves will propagate appreciably far within this timeframe. Although wave attenuation might therefore be important later in an SEP event, and can considerably change the characteristic of the generated turbulence, the effects thereof on SEP propagation can most likely be safely neglected. Wave cascading effects, or more importantly for electrons, wave dissipation effects, were assumed negligible following a similar argument; while it is unclear if these effects will be significant in the assumed timeframes, we also do not currently have a firm grasp of what the wave dissipation rate should be. Although outside the scope of the present study such effect will be included in future work.

What do we expect for a beam (i.e. an excess of forward propagating particles, i.e. $\partial f/\partial\mu>0$) of electrons? We expect to amplify forward moving waves and damp backward moving waves of both polarities. The wavenumber where this occurs (i.e. the resonant wavenumber) however will also depend on the polarity of the waves. Moreover, because forward moving particles resonate more with forward moving LH-polarized waves, we expect these waves to be preferably amplified. Refer again to Fig. 5.

6 Results

6.1 Reference solution

First, we show so-called reference solutions. These are the resulting model intensities without including any wave-growth in the model. It is extremely important to normalize these intensities to the correct levels, as the SEP intensity will directly influence the amount of wave modification. We start by injecting an arbitrary intensity of SEPs into the model, solve the model to obtain the intensity at Earth, and normalize this to observed levels. Note that the normalization is performed only once in order for the results to remain consistent. The temporal profile for this reference solution is shown in the left panel of Fig. 9 for various radial positions; all magnetically connected to the SEP source, i.e. along the same magnetic field line. The maximum intensity at 1 AU is normalized to the maximum value reported by Dresing et al. (2014) for electrons of these energies (indicated on the figures by the red circle). The right panel of the figure shows the corresponding anisotropies. Interesting, the anisotropy near the Sun becomes negative after the initial injection phase. This is due to the large $\lambda_{||}$ value close to the Sun that does not effectively isotropizes the distribution; after the initial injection, SEPs are scattered and isotropized further outwards, where after they propagate back towards the Sun.

The differential intensity of the reference solution shown in Fig. 10, as a contour plot in the ecliptic plane, at $t=2$ hrs after the SEP events was injected into the model. The position of the three “virtual spacecraft”, where the intensities are calculated throughout this work and in Fig. 9, are also shown. Note that the azimuthal extent of this distribution will depend heavily on the effectiveness of perpendicular diffusion. The effect thereof is discussed in a later section.

6.2 Wave growth/damping effects

The left panel of Fig. 11 shows the pre-event energy spectra of the four waves components of the slab turbulence spectrum at a radial position of $r=0.5$ AU. The right panel shows the same spectra after the SEP events passed over the virtual spacecraft. As expected, forward propagating wavemodes were amplified, while backward propagating waves were damped. The forward propagating LH polarized modes is most heavily affected as a beamed electron population will resonate most effectively with these wavemodes. Interestingly, for 100 keV electrons under consideration here, the wave growth and/or damping occurs at the higher end of the inertial range, or possibly near the onset of the dissipation range, of the turbulence spectrum. Newer spacecraft definitely have the possibility to resolve this part of the turbulence spectrum, and results such as shown on the right panel of Fig. 11 should be observable.

To quantity and illustrate where in space most of the wave modification takes place, we calculate the so-called “residual turbulence energy” for each wave as

[TABLE]

where $k^{(1)}_{||}$ and $k^{(2)}_{||}$ is the range of the wavenumber space under consideration (see Fig. 11), $g^{\mathrm{slab}}(k_{||},t,n,j)$ the energy spectra of the different wavemodes and $t_{1}=0$ hrs and $t_{2}=10$ hrs the start and end of the simulation interval. This quantity is shown, for each wavemode separately, as a contour plot in the ecliptic plane of the heliosphere in Fig. 12. To some extent these maps mimic the results of the differential intensity contour (see Fig. 10). This is expected as higher fluxes generally lead to more wave modification. The details are, however, much more complex, with the wave growth dependent on e.g. $D_{\mu\mu}$, the anisotropic part of the SEP distribution, and the resonance functions. As an illustration, note that the most efficient wave growth for LH polarized, forward moving waves is close to the Sun (see panel a), while most effective wave growth for RH-polarized, forward moving waves are closer to Earth orbit.

The pitch-angle diffusion coefficient, as calculated before and after the SEP event, is shown in Fig. 13. The change in $D_{\mu\mu}$ is pitch-angle dependent with the biggest changes occurring at intermediate pitch-angle values, where the level of scattering increased by more than an order of magnitude. This is very large, considering that only the forward propagating waves are amplified. The kink in the calculated $D_{\mu\mu}$ after the event is due to our rather coarse pitch-angle grid used in the model. See also the zig-zags present in the solutions presented in Fig. 11.

To quantify the enhance level of scattering, we calculate $\lambda_{||}$, at different longitudes, after the simulated SEP events. The results are shown in Fig. 14. Note that the different azimuthal positions co-inside with the different virtual spacecraft placed in the computational domain (see again Fig. 10). The effect of wave growth on $\lambda_{||}$ is dramatic; near the Sun $\lambda_{||}$ decreases by at least two orders of magnitude. This decrease, however, depends very strongly on magnetic connectivity to the SEP source with an observer $\sim 45^{\circ}$ away from the source unlikely to see any larger effects. The effects of wave growth also diminishes away from the Sun, with only moderate effects apparent at Earth.

6.3 The effects of perpendicular diffusion

Perpendicular diffusion may also influence the calculated wave growth-rates indirectly by controlling the intensity of SEPs near the Sun: with more effective perpendicular diffusion in the model, it is expected that the SEP particles will spread out more effectively, in terms of longitude, with less particles reaching Earth. However, as we normalize the intensities to Earth values, more particle spreading means that the initial SEP distribution specified near the Sun must be higher. Naively, this argument suggests that more efficient perpendicular diffusion will lead to more efficient wave-growth.

To test this argument, we have performed simulations with various choices of $a\in[0,1]$ leading to different levels of perpendicular diffusion. The resulting peak intensities are normalized to the same value at Earth, and Fig. 15 shows the radial dependence thereof. Note that we calculate the peak intensity along the magnetic field line that connect each point to the SEP source. Referring to Fig. 15, we see that the behaviour is well described by a $\sim 1/r^{2}$ dependence, and that the peak intensity near the source only changes by a factor of $\sim 2$ between the extreme values of $a=0.01$ and $a=1$.

Interestingly, our model results compare well with these of He et al. (2017) who found a $r^{-1.7}$ dependence using a much more involved 3D modelling approach.

7 Discussion

In this paper we have calculated diffusion coefficients for 100 keV electrons in the inner heliosphere based on observed turbulence spectra. Effects such as non-zero magnetic and cross helicities, which are usually neglected, have been included and it was shown shown to have large effects near the Sun. We have also included dynamical effects in the standard FLRW perpendicular diffusion coefficient which may have large effects for low energy particles.

Using these theoretically motivated transport coefficients, we have included, for the first time, wave generation by streaming SEP electrons. Our results show that, for the biggest events, SEP electrons can significantly grow forward propagating waves (especially left-handed polarized waves near the Sun) Our results therefore confirm, to some extent, the observations of Agueda & Lario (2016). These effect diminish towards Earth and are not expected to observable outside of 0.5 AU and when the observer is not magnetically well connected to the flaring region. However, the upcoming Parker Solar Probe and Solar Orbiter spacecraft should be able to observe these amplified levels of turbulence in the coming years.

Our model of wave growth is by no means complete; a full spectrum of particles needs to be injected into interplanetary space, with wave cascading and dissipation processes treated more self consistently. Wave attenuation effects will also have to be included in future. Such a model is however not computationally viable at the moment. We are also not able, at present, to confirm whether a so-called streaming limit (e.g. Reames & Ng, 2010) is present in the electron simulations. However, we expect that the simulated temporal evolution of the SEP electron intensity at e.g. Earth will be greatly affected by wave growth for the largest SEP events given the large changes in $D_{\mu\mu}$. Interestingly, we note that the enhanced turbulence levels (and correspondingly small $\lambda_{||}$) near the Sun can strongly enhance the diffusive shock acceleration rate in this region, possibly boosting the level of energetic electrons reaching Earth when a shock from a coronal mass ejection is present.

We have modelled the peak intensity along the MHF line of optimal magnetic connectivity. Our calculations compare well with previous 3D modelling results and the observations of Lario et al. (2013). Moreover, it was shown that the level of perpendicular diffusion can influence the wave growth rate by changing the SEP intensity near the Sun (as normalized to constant levels at Earth). We again emphasize that perpendicular diffusion is an important transport process, as also illustrated by Dröge et al. (2010), He (2015), Laitinen et al. (2016) and Zhao & Zhang (2018), amongst others, and should be included in SEP transport models.

We are confident that our approach of using theoretically derived coefficients, combined with turbulence quantities, consistent with the limited available observations, can produce realistic SEP intensities in the inner heliosphere. The approach can be refined with more measurement, but is the most effective way to progress towards a predictive model and understanding of the fundamental processes that couple particle scattering to solar wind fluctuations. Moreover, most of these turbulence quantities can be more constrained by the upcoming Parker Solar Probe and Solar Orbiter spacecraft.

This work is based on research supported in part by the National Research Foundation (NRF) of South Africa (grant no. 106049). Opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the NRF. RDS acknowledges partial financial support from the Fulbright Visiting Scholar Program.

Appendix A Calculation of the growth/damping rate

Below we show the detailed calculation of the wave growth/damping rate. The derivation is loosely based on work presented by e.g. Ng & Reames (1994) and Vainio (2003), amongst others.

The change in wave energy during a scattering event can be approximated as $\Delta E_{w}=-jV_{A}p\Delta\mu$, with the total contribution from all particles in the distribution giving the total change in the energy density of the waves

[TABLE]

Assuming pitch-angle scattering to be the only process changing the distribution, we may estimate this change as

[TABLE]

where the second part of the expression follows from the standard definition of the Fokker-Planck equations, resulting in

[TABLE]

and hence that the average change in wave energy density is given by

[TABLE]

Assuming $f$ to be gyro-tropic, we find

[TABLE]

In order to make wave-generation implementable in our numerical model, we have to assume a mono-energetic SEP distribution, i.e. $f(\vec{r},p,\mu,t)=p_{0}f(\vec{r},\mu,t)\delta(p-p_{0})$, so that

[TABLE]

If the amplified waves are Alfvénic, both the energy density of the kinetic and magnetic contributions change (with an equipartition for non-dispersive Alfvén waves), resulting in

[TABLE]

By using Eq. 9, the change of energy density can be evaluated as a change in wave power,

[TABLE]

Combining Eqs. A6 and A8, we get

[TABLE]

and we now need to estimate the pitch-angle diffusion coefficient, $D_{\mu\mu}$. In order to get a tractable expression, we use the negligible damping version, given by Eqs. 27 – 29, to obtain

[TABLE]

As this expression should be true for each independent wavemode, we can use Eq. 8 to find the growth/damping rate for the different wavemodes as

[TABLE]

Manipulating the $\delta$-function gives

[TABLE]

where the resonant pitch-angle is

[TABLE]

and the energy is the background magnetic field is

[TABLE]

Appendix B Plasma dispersion relations

For the sake of completeness, we list the dispersion relations for circularly polarized waves in the cold plasma limit, that is, assuming $T_{\mathrm{e}}=T_{\mathrm{p}}=0$. These relations are given by (see, amongst other, Stix, 1962)

[TABLE]

and

[TABLE]

where $k_{||}$ is the parallel wavenumber, $\omega$ the frequency in the bulk (flow) plasma frame, $c$ the speed of light and $V_{A}$ the Alfvén speed. $\Omega_{p}$ and $|\Omega_{e}|$ are the proton and electron cyclotron frequencies, given as

[TABLE]

where $\Lambda=m_{p}/m_{e}$ is the ratio between proton and electron masses, $e$ the elementary charge and $B$ the mean background magnetic field, while the different plasma frequencies are given as

[TABLE]

with $n_{e}$ the electron number density and $\epsilon_{0}$ the permittivity of free space. Note that we assume quasi-neutrality with $n_{e}\approx n_{p}$. In terms of these frequencies, $V_{A}$ is

[TABLE]

In deriving Eqs. B1 and B2 it was assumed that $\omega\ll\omega_{p}\ll\omega_{e}$. The dispersion relations are shown in Fig. 1 as the solid blue lines.

In the cold plasma limit, the waves are strongly damped near either the proton or electron cyclotron frequency. However, when a warm plasma is considered (such as the solar wind), the thermal motion of the plasma particles, with a speed of

[TABLE]

for either protons or electrons, need to be considered. In Eq. 15 this was done in a rather ad-hoc way by relying on the results of Schreiner et al. (2017). Indeed, this approach was also shown by Engelbrecht & Strauss (2018) to compare well with observed quantities. The difficulty lies is calculating, in a tractable analytical fashion, the damping rate of MHD waves by a warm plasma. Unfortunately, this is only possible under the assumption of a small growth rate, $\gamma_{\omega}\ll\omega$ and leads to

[TABLE]

and needs to be evaluated for each polarization mode separately (for details, see e.g. Chen et al., 2013). The function $D$ is related to the distribution of the particles under consideration, and for our purposes, we additionally assume no temperature anisotropy, so that

[TABLE]

and

[TABLE]

where, for the latter equation, we have furthermore assumed that $v^{p,e}_{\mathrm{th}}\rightarrow 0$. The calculated damping rates, using 1 AU plasma values, are shown in Fig. 1 as the solid red curves.