Preacceleration in the Electron Foreshock. I. Electron Acoustic Waves

Astrophysical shocks occur ubiquitously in our universe, with detections of nonthermal X-ray emission from sources such as supernova remnants implying the presence of accelerated electrons. Their primary acceleration mechanism is widely accepted to be diffusive shock acceleration (DSA) (Bell 1978), yet alone, it does not provide a complete picture. This is because of the injection problem: For DSA to work, the Larmor radii of the particles it accelerates must be comparable to the finite shock width, a criterion not met by thermal electrons. Therefore, they require some preacceleration, the exact physical nature of which is crucial to fully comprehend many astrophysical phenomena.

The study of electron-scale phenomena, which may be responsible for preacceleration, requires a method capable of resolving these small kinetic scales. Particle-in-cell (PIC) simulations fulfill this criterion and have previously proved invaluable in helping to establish the contribution to electron preacceleration in perpendicular shocks from mechanisms such as shock-surfing acceleration, magnetic reconnection, and stochastic Fermi acceleration (Bohdan et al. 2017, 2019a, 2019b, 2020a, 2020b).

The obliquity angle, θ_Bn, defined as the angle subtended between the shock normal and upstream magnetic field, can heavily influence the behavior of the shock, with DSA the dominant particle energization process for quasi-parallel shocks and shock-drift acceleration (SDA) increasingly important at larger subluminal magnetic field inclination angles (Ellison et al. 1995; Sironi & Spitkovsky 2009). In comparison to perpendicular shocks (where θ_Bn = 90°), which typically have a well-defined shock transition with a narrow foot region, oblique shocks permit the escape (or reflection) of particles along the upstream magnetic field lines and back into the shock upstream, which creates extended regions known as the electron and ion foreshocks, depending on the particle species. (e.g., Burgess 1995; Fitzenreiter 1995; Treumann 2009).

In previous studies of oblique shocks, the mechanism of electron reflection to the upstream region has been shown to be mirror reflection (Honda & Honda 2005). More energetic electrons are favorably reflected back upstream, with these electrons usually pre-energized by shock-surfing acceleration (Amano & Hoshino 2007). Notably, 1D PIC simulations have demonstrated that such reflected electrons are capable of generating electrostatic waves in the shock foot, where they may undergo SDA before being injected downstream (Xu et al. 2020; Kumar & Reville 2021).

As any electrons upstream of an oblique shock must first encounter the foreshock and any turbulence present there before interacting with the shock, understanding the properties of this region and how they influence the upstream plasma are of utmost importance when addressing electron preacceleration and the shock injection problem.

In this paper, we investigate the impact of changing the orientation of the upstream magnetic field on the electron foreshock. The structure of the manuscript is as follows: First, we outline the initial mechanism responsible for reflecting electrons. Then, we show that these reflected electrons generate electron acoustic waves and examine their properties. Finally, we look at what happens to the energy associated with these waves, and how they affect electrons in the upstream plasma.

The work presented here makes use of code adapted from TRISTAN (Buneman 1993), which allows us to study fully relativistic shocks while tracking 2-spatial and all 3-velocity components (2D3V) of the individual particles in our plasma (Niemiec et al. 2008). Particle positions are updated using the Vay solver (Vay 2008). When applied to perpendicular shocks, the data from these simulations have proven to accurately reproduce the fundamental physical processes present in larger and computationally more expensive full 3D simulations (e.g., Bohdan et al. 2017; Matsumoto et al. 2017) and measured in situ (Bohdan et al. 2020a, 2021). This version of the code also permits the tracing of individual particles, where we can track the components of their velocities and local fields, allowing for a detailed analysis of the underlying physical processes.

Previous works investigating perpendicular shocks successfully employed a flow–flow method in which two counterstreaming plasma slabs collide and generate two shocks (Wieland et al. 2016; Bohdan et al. 2017, 2019a, 2019b, 2020a, 2020b, 2021). Relative to perpendicular shocks, an investigation of an oblique shock requires a larger simulation box, and therefore a greater computational expense, to ensure reflected electrons are captured, so in this work, we instead implement a reflecting wall setup (injection method). In this configuration, a single shock is produced by reflecting an incoming plasma beam from a conducting wall located at x = 0 in our simulation box, where the velocity transformation v_x → −v_x is performed (v_y and v_z remain unaltered) (Quest 1985; Burgess et al. 1989). To prevent an electromagnetic transient from being induced by the large values of ∇ × E and the corresponding ∂ B /∂t, we taper the fields near the reflection boundary following Wieland et al. (2016), where a drift velocity is also assigned to ions in these regions to compensate for the now nonzero ∇ × B . After reflection, this drift velocity is removed from the ions because the drifts at a real contact discontinuity compensate each other. By adding further domains, fresh plasma is regularly injected into the upstream from the right boundary of the simulation box with the same initial parameters as outlined above. One can visualize the injection of new particles into a new domain as dividing it into small subregions with one of each particle species injected at a random, but identical, position within each subregion. The sizes of these subregions are selected such that the desired number of particles per cell, n_ppc, is returned.

Excepting the upstream magnetic field configuration, all of the simulations presented here use the same initial parameters, which are expressed in Heaviside–Lorentz units. First, we initialize upstream plasma with a bulk velocity of ${{\boldsymbol{v}}}_{{ups}}=-0.2\,c\,\hat{{\boldsymbol{x}}}$ in the simulation plane, which corresponds to a shock velocity of v_sh ≃ 0.264c in the upstream rest frame (or ${v}_{\mathrm{sh}}^{\prime} \simeq 0.0673\,c$ in the simulation frame). This upstream plasma is initialized with a uniform large-scale magnetic field with a field orientation described by the obliquity angle, θ_Bn, and angle relative to the simulation plane, ϕ, such that ${{\boldsymbol{B}}}_{{\bf{0}}}={B}_{0}(\cos {\theta }_{\mathrm{Bn}},\sin {\theta }_{\mathrm{Bn}}\cos \phi ,\sin {\theta }_{\mathrm{Bn}}\sin \phi )$. This gives rise to a motional electric field, which is computed as E = − v _ups × B ₀ . Such a magnetic field configuration results in particles that effectively have three degrees of freedom, thus the adiabatic index is γ = 5/3. The shock compression ratio depends on this quantity, and for all simulations here has the value r ≈ 4.

We assume that upstream ions and electrons are initially in thermal equilibrium with each other, such that ${k}_{B}{T}_{i}={k}_{B}{T}_{e}={m}_{i}{v}_{i,\mathrm{th}}^{2}/2={m}_{e}{v}_{e,{th}}^{2}/2$, where k_B is the Boltzmann constant, T denotes temperature, and v_th is the most probable particle speed in the upstream reference frame. The subscripts i and e represent ions and electrons, respectively, with their masses characterized by m_i /m_e = 64. This corresponds to numerical values of v_th,e = 0.0713c and v_th,i = 0.00089c and allows us to determine the sound speed as ${c}_{{\rm{s}}}={(\gamma {k}_{B}{T}_{{\rm{i}}}/{m}_{{\rm{i}}})}^{1/2}=0.00814c$. These values give the ratio of electron plasma to magnetic pressures (plasma beta) as β = 1.

With n_ppc = n_i = n_e = 20 particles per cell as the number densities of ions and electrons, respectively, the Alfvén velocity is defined as ${v}_{{\rm{A}}}={B}_{0}c/\sqrt{{n}_{\mathrm{ppc}}({m}_{e}+{m}_{i})}$ = 0.0125c, where μ₀ is the permeability of free space. This sets the Alfvénic Mach number as M_A = v_sh/v_A ≈ 21, which we define by convention in the upstream reference frame.

Our simulations all share a ratio of electron plasma frequency, ${\omega }_{\mathrm{pe}}=\sqrt{{q}_{e}^{2}{n}_{e}/({\epsilon }_{0}{m}_{e})}$, to electron cyclotron frequency, Ω_e = q_e B₀/m_e , of ω_pe /Ω_e = 9.9, where ₀ is the vacuum permittivity. The simulation runtime is characterized by the ion cyclotron time, ${{\rm{\Omega }}}_{\mathrm{ci}}^{-1}={m}_{i}/({q}_{i}{B}_{0})$, where each simulation is run for a total of ${T}_{\mathrm{tot}}=7.8\,{{\rm{\Omega }}}_{\mathrm{ci}}^{-1}$ to enable a comparison with the results of Matsumoto et al. (2017). To adequately resolve all of the important timescales in our simulation, we use a simulation time step of $\delta t=1/40\,{\omega }_{\mathrm{pe}}^{-1}$.

In these simulations, we study the rates and implications of the reflection of incident electrons at the shock. Accordingly, the transverse size of our simulation box must be large enough to encapsulate electron-scale phenomena. From our setup, the gyroradius of thermal electrons measures r_g,e = v_e,th/Ω_ce =14Δ, where Δ represents one cell in our simulation grid while the electron skin length measures λ_e = c/ω_pe = 20Δ. For ions, these quantities are both higher by a factor of $\sqrt{{m}_{i}/{m}_{e}}$. The transverse size of our simulation box measures 768Δ (about 55r_ge or 4.8λ_i ). This width is therefore large enough to probe electron physics, yet small enough to save computational expense and allow resources to be allocated to the necessary extension of our simulation box with respect to time so that it contains the reflected electrons.

In this paper, we present a number of simulations of varying obliquity angles, θ_Bn, between the shock normal and upstream magnetic field. In some cases, we also varied the field inclination to the simulation plane, ϕ. For clarity, these angles are defined in Figure 1, with a full list of angles used provided in Table 1. We fix all of the other parameters outlined above in order to allow a physical comparison of the electron foreshock region in each case. The inclusion of θ_Bn = 743 allows us to directly compare our results with full 3D simulations of Matsumoto et al. (2017), compared to which our parameters are identical. To investigate the effect of the in-plane angle ϕ and to allow a more rigorous comparison with the 3D case, we run further simulations with ϕ = 90° for our two largest values of θ_Bn. For these parameters, the angle at which the shock becomes superluminal is given by ${\theta }_{\mathrm{crit}}={\tan }^{-1}(c/{v}_{\mathrm{sh}})\sim 75.2^\circ$, thus our range of angles allows a substantial investigation into how the properties of the electron foreshock depend on the initial upstream magnetic field configuration.

Zoom In Zoom Out Reset image size

Figure 2 shows the shock structure for the in-plane simulation with θ_Bn = 45° at the final time step of the simulation, where $t\,=\,7.8\,{{\rm{\Omega }}}_{\mathrm{ci}}^{-1}$. Here the shock is fully formed with an overshoot at x ∼ 35 λ_i , and panel (a) shows a map of the electron density relative to the upstream density. In the downstream region toward the left of the figure, we measure a compression ratio close to the analytically expected value of ρ/ρ₀ ≈ 3.97, with the same true of the magnetic field (panel (e)). From this panel, we also see a smooth transition of the magnetic field profile from downstream to upstream, where it remains at around the initial upstream value where ${(B/{B}_{0})}^{2}=1$.

Zoom In Zoom Out Reset image size

In contrast to the magnetic field, panel (f) shows that the electric field steadily decreases with increasing distance from the shock. Panels (b) and (c) show that waves are present in the upstream region where ${(E/{E}_{0})}^{2}\gt 1$. These waves have no magnetic component and are purely electrostatic and in this in-plane simulation lie in the xy plane by definition. Visually, we note that a higher density of large-amplitude waves occupies the regions with greater ${(E/{E}_{0})}^{2}$ (panel (f)) and that the physical length scales on which they occur are of order λ_se and comparable to electron scales, immediately implying that they may be related to electrons. Panel (d) shows the phase-space momentum distribution for all electrons in the selected region of the simulation. There is a clear excess of electrons with positive Γβ_x with a number density that decreases farther from the shock. The regions in which these electrons are present correspond to the region occupied by electrostatic waves, the nature of which we discuss in the following section. Electrons reflected at the shock are able to carry energy back upstream, and understanding what happens to this energy is the object of this paper.

Our oblique-shock simulations reveal the presence of electrostatic waves that are spatially coincident with the reflected electrons. To determine their nature, we note that following Bohdan et al. (2022), their wavelength and growth rate indicate that they are excited by the electron acoustic instability. In contrast to the Buneman instability (Buneman 1958), which is driven directly at the shock by the ion/electron interaction, the electron acoustic instability results from the velocity difference between the cold incoming and the hot reflected electrons along the magnetic field direction (Gary 1987). Simulations with periodic boundary conditions presented in Bohdan et al. (2022) confirm that reflected ions play no role in generating these electrostatic waves, so we can immediately rule out electron cyclotron drift instability (EDCI), which is excited by reflected ions interacting with incoming electrons (e.g., Muschietti & Lembège 2006; Yang et al. 2020), as the instability that drives them. Similarly, we can rule out synchrotron maser instability (Hoshino & Arons 1991) as we do not see the characteristic electron-ring structure in their phase-space distributions (e.g., Plotnikov & Sironi 2019; Babul & Sironi 2020).

As they are driven by reflected electrons, these electrostatic waves are predominantly propagating along the background magnetic field. Because the mechanism driving this electron/electron instability does not involve the more inertial ions, it grows much faster and saturates already far ahead of the shock front. We refer to these waves as electron acoustic waves (EAWs) and discuss their properties, and those of the reflected electrons that drive them, before outlining the electron reflection mechanism in the following sections.

4.1. Electron Spectra

To illustrate the properties of the reflected electrons, we calculate spectra and reflection rates over a region of the simulation box. To define this region, we first estimate the location of the shock foot, x_foot, given as the location where the mean ion density averaged across the transverse direction of the simulation box first reaches ρ/ρ₀ = 1.05, as measured from upstream to downstream. By defining x_foot based on ion density, we ensure that it remains independent from the electron-induced electrostatic waves that we analyze in this work. It is well established that electrons can undergo shock-surfing acceleration along the leading edge of the shock (Bohdan et al. 2019a), after which they can either enter the downstream or be reflected back upstream (Amano & Hoshino 2007). We wish to omit these electrons from our spectra and reflection rate calculations, hence they are calculated in the simulation frame for a region where the location in the simulation box, x, satisfies x_foot(t) + 2λ_si < x ≤ x_foot(t) + 5λ_si, where λ_si is the ion skin depth, which is resolved by 160 cells. This range is farther from the shock than shock-surfing electrons and so ensures they are not included in our calculations. Similarly, if this region is too large, we risk including regions of the simulation box that contain no reflected electrons because they have not yet had time to travel this far. This region is shown for Ω_ci t = 7.8 for θ_Bn = 45° in panels (e) and (f) of Figure 2.

4.1.1. In-plane Simulations

Electron spectra for in-plane simulations after Ω_ci t = 7.8 are shown in Figure 3. Two features immediately stand out. First, we see that for lower values of θ_Bn, more electrons are reflected. The second feature is that the peak value of Γ − 1 for these reflected electrons becomes higher for larger θ_Bn. The physical interpretation of these features is interconnected. Previous results have verified that upstream electrons arriving at the shock may only be successfully reflected should their velocity, v_e , satisfy ${v}_{e}\cos {\theta }_{\mathrm{Bn}}\geqslant {v}_{\mathrm{sh}}$, where all quantities are in the shock frame (Amano & Hoshino 2007; Matsumoto et al. 2017). As reflected electrons travel on average in the direction of the upstream magnetic field lines, this result is somewhat intuitive as they require a significant velocity component in the direction parallel to the shock propagation; otherwise, they cannot outrun it.

Zoom In Zoom Out Reset image size

It follows that electrons require larger velocities for larger values of θ_Bn. Furthermore, as less energetic electrons are relatively more abundant due to the lower amount of work done on them to achieve such energies, more of them fulfill the reflection criterion when θ_Bn is smaller, which in turn leads to a higher reflection rate. Similarly, at larger θ_Bn the number of electrons eligible for reflection is reduced, hence we measure a smaller reflection rate overall. Figure 3 shows that the peak value of Γ − 1 typically increases with θ_Bn, which further supports this explanation.

One important quantity we assess here is the electron reflection rate, which we define as R = N_e,ref/(N_e,ref + N_e,ups), where N_e,ref is the number of reflected electrons and N_e,ups is the number of upstream electrons in the region over which spectra are calculated. This was obtained over the regions for which our spectra were calculated and separated the low-energy and high-energy components at each time step at the value of Γ − 1 for which d(N_e (Γ − 1))/d(Γ − 1) = 0. Electrons with Γ − 1 below and above this value were labeled as being upstream or reflected, respectively. The reflection rate as calculated by this method is shown in Figure 4. As is consistent with earlier arguments, we see that simulations with lower values of θ_Bn have larger electron reflection rates, with these reaching a high of ∼5% for θ_Bn = 30° and falling to only ∼0.05% for θ_Bn = 63°. We do not observe electron reflection for θ_Bn = 743. We further note that the onset of electron reflection typically begins at later times for larger values of θ_Bn. This is consistent with them requiring a larger time to be preaccelerated so that they have a large-enough velocity to outrun the shock (Amano & Hoshino 2007).

Zoom In Zoom Out Reset image size

4.1.2. Out-of-plane Simulations

Here we describe results from the simulations with a fixed θ_Bn = 63° but varying ϕ, such that ϕ = 0°, 45°, 90°. We focus on these rather than the simulations with θ_Bn = 743 because, for this value of θ_Bn, independent of ϕ, we see no electron reflection. Figure 5 shows the spectra from the simulations with θ_Bn = 63° after Ω_ci t = 7.8 in the simulation frame. We see significantly more reflection for the out-of-plane case where ϕ = 90° than for ϕ = 0° and ϕ = 45°. The corresponding reflection rates were about 0.03% for ϕ = 0° and ϕ = 45°, but near 4% for ϕ = 90°. This immediately suggests an electron reflection mechanism that is more appropriately captured by out-of-plane simulations.

Zoom In Zoom Out Reset image size

4.2. Energy Content of EAWs

In the wider context of astrophysical shocks, it is important to address the questions of how much kinetic energy these electrons carry back upstream, and whether it can drive any turbulence that affects the upstream plasma. The kinetic energy associated with a reflected electron is given by E_K = (Γ − 1)m_e c², and we have previously shown in Figure 3 that for larger θ_Bn, individual reflected electrons carry on average a higher E_K but there are fewer of them in total.

A more informative quantity is the energy density of reflected electrons, ${U}_{\mathrm{ref}}=({m}_{e}/2){{Rn}}_{\mathrm{ppc}}{\overline{{v}_{\mathrm{ref}}}}^{2}$, which is defined using the reflection rate, R, and an average velocity of reflected electrons, $\overline{{v}_{\mathrm{ref}}}$. We normalize this quantity relative to the upstream electron energy density, ${U}_{\mathrm{ups}}=({m}_{e}/2){n}_{\mathrm{ppc}}{v}_{\mathrm{ups}}^{2}$, with the results shown for all simulations in Figure 6, along with the density ratios of the reflected and upstream electron beams.

Zoom In Zoom Out Reset image size

We see immediately that it is the higher reflection rate that dominates the total energy density of the reflected electron beam as opposed to the mean electron kinetic energy, as shocks with lower θ_Bn have larger reflected energy densities, and their electron foreshocks comprise larger numbers of lower energy electrons relative to those at larger θ_Bn. We therefore expect that it is more likely for lower-θ_Bn shocks to have additional features in the electron foreshock region. Figure 6 also shows shocks with θ_Bn = 63° and varying plane angle ϕ (see solid lines). Here, the energy density ratios are comparable for ϕ = 0° and ϕ = 45°, but significantly higher for the out-of-plane case with ϕ = 90°.

We now assess what happens to the energy content carried upstream by reflected electrons and focus on any turbulence it can drive in this region. This is especially important, as upstream electrons must interact with the electron foreshock before encountering the shock, so any mechanisms that can either lead to preacceleration or potentially compromise their transition to the shock downstream are necessary to understand the global picture of shock physics.

Our in-plane simulations reveal the presence of EAWs far upstream of the shock, over the same spatial extent as where reflected electrons have propagated. Though these waves are also present in the out-of-plane shocks, they are less well captured as a result of the reflected electrons streaming on average along the upstream magnetic field lines, and hence wavevectors parallel to the streaming are not allowed (Bohdan et al. 2022). For this reason, in this subsection, our analysis focuses on our in-plane results.

To determine the energy content of the EAWs and demonstrate that they are powered by reflected electrons, we perform a two-dimensional discrete Fourier transform (DFT) across the region of our simulation box where they are present. Following Vafin et al. (2018), the dimensionless 2D Fourier amplitudes are calculated as

$$\begin{eqnarray}&&\begin{array}{l}\widetilde{E}({k}_{m,n})=\displaystyle \frac{e}{{m}_{e}c{\omega }_{p}}\displaystyle \frac{1}{{N}_{x}{N}_{y}}\cdot \\ \quad \times \displaystyle \sum _{i=0}^{{N}_{x}-1}\displaystyle \sum _{j=0}^{{N}_{y}-1}{{\boldsymbol{e}}}_{{\boldsymbol{k}}}\cdot {\boldsymbol{E}}({x}_{i},{y}_{j})\exp \left(2\pi i\left[\displaystyle \frac{{mi}}{{N}_{x}}+\displaystyle \frac{{nj}}{{N}_{y}}\right]\right),\end{array}\end{eqnarray} \tag{ 1 }$$

where m and n denote the respective indices of the x and y wavevector components, and e _k is in the direction of the unit wavevector. Using Parseval's theorem, we can calculate the energy density associated with the electron acoustic waves by summing up $| \widetilde{E}({k}_{m,n}){| }^{2}$ and dividing by the volume.

The top panel of Figure 7 shows the DFT of the x and y components of the electric field in an upstream region containing EAWs for θ_Bn = 30°. The wavelength of the fastest-growing mode here is consistent with the linear theory of electron acoustic instability. A Gaussian momentum distribution of reflected electrons excites electron acoustic waves with the growth rate of Bohdan et al. (2022):

$$\begin{eqnarray}&&\displaystyle \frac{\gamma ({k}_{\parallel })}{{\omega }_{\mathrm{pe}}}=-\sqrt{\displaystyle \frac{\pi }{2}}\displaystyle \frac{{n}_{\mathrm{ref}}}{{n}_{\mathrm{ups}}}\displaystyle \frac{{m}_{e}^{2}{\omega }_{\mathrm{pe}}^{2}}{{p}_{\mathrm{th}}^{2}}\ \displaystyle \frac{g({k}_{\parallel })}{{k}_{\parallel }^{2}}\exp \left(-{\displaystyle \frac{g({k}_{\parallel })}{2}}^{2}\right).\end{eqnarray} \tag{ 2 }$$

Here n_ref and n_ups denote the densities of reflected and incoming upstream electrons, respectively, and

$$\begin{eqnarray*}&&g(k)=({p}_{\mathrm{res}}(k)-{p}_{0})/{p}_{\mathrm{th}}\end{eqnarray*}$$

depends on the resonant momentum p_res(k), which follows from the resonance condition

$$\begin{eqnarray}&&\displaystyle \frac{{p}_{\mathrm{res}}}{{\left[{m}_{e}^{2}{c}^{2}+{p}_{\mathrm{res}}^{2}\right]}^{1/2}}=\displaystyle \frac{{\omega }_{\mathrm{pe}}}{{kc}},\end{eqnarray} \tag{ 3 }$$

where the reflected electron beam has a momentum distribution centered around p₀ with a standard deviation p_th.

Zoom In Zoom Out Reset image size

We analytically determine the expected value of k_∥ by measuring the required quantities from our simulation data. Using the reflected electrons that occupy the region upstream of the shock shown in the central two panels of Figure 7, we obtain values of p₀ = 0.92 m_e c, p_th = 0.309 m_e c, and n_ref = 0.01 n_ups. For these parameters, Equation (2) yields a peak growth rate of 0.024 ω_pe at k_∥ λ_se = 1.76, large enough for significant wave growth within small parts of the simulation time, $T\approx 5000\,{\omega }_{\mathrm{pe}}^{-1}$.

The top panel of Figure 7 shows the Fourier transform of the electric field in this same region, allowing for a direct comparison of k_∥ obtained from our simulations to be made to the analytically expected value. The dashed silver line indicates the location of k_∥, which is in the direction of the reflected electron beam along the direction of the upstream magnetic field. From this, we calculate the observed scale of the electrostatic waves as ${k}_{\parallel }={k}_{x}\cos {\theta }_{\mathrm{Bn}}+{k}_{y}\sin {\theta }_{\mathrm{Bn}}\sim 1.79$. The analytical value assumes the reflected electron beam follows a Gaussian distribution. As previously discussed, electrons with higher energies are more likely to be reflected back upstream, hence the true distribution is not normally distributed in momentum space. Additionally, Figure 4 shows that we observe temporal fluctuations in the electron reflection rate, which from Equation (2) we expect to manifest proportionally in the growth rate of the EAWs, meaning this quantity will not necessarily be constant across a given upstream region of our simulation box. Despite these discrepancies between theoretical calculations and simulation results we nevertheless achieve a very good match in k_∥, strongly supporting the notion that these waves are indeed generated by the electron acoustic instability.

Furthermore, we note that the most energetic reflected electrons are found at the leading edge of the foreshock. This time-of-flight effect arises naturally as these more energetic and therefore faster electrons are more capable of outrunning the shock. That they have larger streaming speeds means that the associated electron acoustic waves have larger wavelengths farther out from the shock, while the lower number density of these electrons reduces the total power contained in these waves.

The total energy contained in the EAWs across the upstream region of our simulations is computed by using Parseval's theorem applied to a window that we slide upstream in our simulation domain. This window is square in shape with a size of 40λ_se × 40λ_se, which ensures the electron acoustic waves are appropriately sampled. As we are not interested in the downstream region here, the first window begins at the same point upstream as was used to define our spectra calculation. Likewise, we utilize the same methods to measure spectra and quantify the kinetic energies of the reflected and upstream components from them.

Figure 8 plots the normalized energy density contained in these electrostatic waves as calculated from summing up the Fourier power in each 40λ_se × 40λ_se window as a function of the kinetic energy of the contained reflected electrons, again normalized to that upstream. It is clear that when an extended foreshock is present, a strong correlation exists between the two, further indicating that these waves are indeed driven by reflected electrons. Here, the floor at around U_ES/U₀ ∼ 7 × 10⁻⁴ is likely a result of electrostatic noise, making a threshold reflection rate for the generation of EAWs difficult to ascertain.

Zoom In Zoom Out Reset image size

4.3. Electron Reflection Mechanisms

As EAWs are excited by a velocity difference between cold incoming and hot reflected electrons, it is important to know what causes the reflection of the hot component. It has long been established that ions that are reflected by the electrostatic potential at the leading edge of the shock can generate electrostatic waves upon their interaction with incoming electrons (Hoshino & Shimada 2002). Note that electrons are oppositely charged to ions and so cannot be reflected by the shock potential, and instead would be accelerated downstream. The requirement to excite these so-called Buneman waves is that the relative velocity difference between the two-particle species, Δv, must exceed the thermal speed of the upstream electrons (Buneman 1958). They have a wavelength of (e.g., Amano & Hoshino 2009b) ${\lambda }_{B}=(2\pi /c){\rm{\Delta }}v{\lambda }_{\mathrm{se}}\approx (2\pi /c)\tfrac{3}{2}{v}_{\mathrm{sh}}{\lambda }_{\mathrm{se}}$, where λ_se is the electron skin length.

Figure 9 shows the early-time reflection of a typical electron within our simulation. Its trajectory is shown by the black line, with the color maps illustrating the strength of the x (leftmost) and y (center) components of the electric field. Here, the electron, which is initially traveling with a velocity x component directed from right to left (− x direction) interacts with this Buneman wave and turns back. We note also from the trajectory that this electron has a relatively large y momentum, as without this it would simply move with the bulk flow from right to left. This detail is important, as a comparatively large transverse momentum allows the electron to relatively freely enter the potential well associated with the electrostatic wave. Figure 10 confirms that the more energetic electrons injected into the upstream plasma are more favorably reflected in this manner. Once trapped, the electron can undergo shock-surfing acceleration where it is accelerated by the motional electric field or be accelerated by the electrostatic field itself (see appendix of Bohdan et al. 2019a). To escape the potential well, the trapping force (qE_BI, where E_BI is the electric field associated with a Buneman wave) must be exceeded by the Lorentz force (q_e v_e B₀/c) (Hoshino & Shimada 2002; Matsumoto et al. 2012).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Importantly, the involvement of Buneman waves in the initial reflection process is supported by the fact our simulations show a larger electron reflection rate when the plane angle ϕ = 90°. For an out-of-plane magnetic field configuration, the energization of electrons by Buneman waves relative to the in-plane case is more efficient (Bohdan et al. 2017, 2019a), thus the probability that they gain enough energy to outrun the shock and escape upstream is increased. This is because the gyromotions of the ions driving the Buneman waves are appropriately captured in the simulation plane for ϕ = 90°, but this is not the case for other angles. Accordingly, electron shock-surfing on Buneman waves becomes more efficient, though this alone cannot explain why they are reflected.

To address this, we consider the direction in which work is done on individual electrons. This can be separated into components either in the perpendicular or parallel direction to the local magnetic field. One can visualize the former as bringing about a change in the gyroradius, r_ge = p_⊥/(q_e B) while solely increasing the parallel velocity leads to a decrease in pitch angle, $\alpha =\arctan \left({v}_{\perp }/{v}_{\parallel }\right)$.

Perpendicular work naturally arises from acceleration due to the motional upstream electric field, defined as E = − v _ups × B ₀ , and thus is oriented perpendicularly to the initial upstream magnetic field in the electron foreshock that we are interested in here. It essentially constitutes gyration in the rest frame of the magnetic field. One such mechanism that increases perpendicular momentum is shock-surfing acceleration (e.g., Shapiro & Üçer 2003; Amano & Hoshino 2009a), which relies on the formation of Buneman waves described earlier. Upstream electrons can then become trapped in the potential wells of these waves and are accelerated perpendicularly by the motional electric field (e.g., Shimada & Hoshino 2000; Hoshino & Shimada 2002).

Conversely, for electrons to be accelerated parallel to the local magnetic field, the plasma must be somehow disturbed, otherwise E and B remain perpendicular to each other. Supposing, hypothetically, an electron was instantaneously reflected back upstream with the transformation v_x → −v_x and an oblique upstream magnetic field orientation described by $B=({B}_{0}\cos {\theta }_{\mathrm{Bn}},{B}_{0}\sin {\theta }_{\mathrm{Bn}})$. Here, we can see intuitively that the sign of p_∥ reverses, while we expect ∣p_⊥∣ to remain constant. An example of this is magnetic mirroring (e.g., Ball & Melrose 2001), in which a particle encountering a shock is either reflected back upstream or transmitted downstream, depending on the pitch angle. A signature property of this process is the conservation of the first adiabatic invariant, μ, defined by

$$\begin{eqnarray}&&\mu =\displaystyle \frac{{m}_{e}{v}_{\perp }^{2}}{2B}=\displaystyle \frac{{m}_{e}{v}^{2}{\sin }^{2}\alpha }{2B}\end{eqnarray} \tag{ 4 }$$

for pitch angle α. For μ to remain constant, the kinetic energy of the particle must scale with the magnetic field strength. The condition for reflection can be derived by Lorentz transforming into the de Hoffmann–Teller frame in which the fluid velocities are parallel to the magnetic field lines on both sides of the shock (de Hoffmann & Teller 1950). As there are no electric fields (e.g., Toptygin 1980) and noting that v² is the same before and after reflection, we find ${\sin }^{2}{\alpha }_{u}/{B}_{u}={\sin }^{2}{\alpha }_{d}/{B}_{d}$ for a particle crossing from upstream (subscript u) to downstream (subscript d). Following Ball & Melrose (2001), we see that if a particle had ${\sin }^{2}{\alpha }_{u}\gt {B}_{u}/{B}_{d}$, it would have an undefined pitch angle downstream. Instead, it is reflected, with the reflection condition in this frame being α_u > α_c for ${\sin }^{2}{\alpha }_{c}={B}_{u}/{B}_{\max }$, with the subscript c here denoting the critical value and B_max the peak of the magnetic field along the particle trajectory.

Figure 11 shows some of the properties of one of the earliest reflected electrons in the simulation with θ_Bn = 63°. The green line in the lower panel shows the trajectory of the electron, with its location along the x-axis of the simulation box plotted against time. At around Ω_ci t ∼ 1.75, the electron encounters Buneman waves in the shock foot and begins to shock-surf, which is indicated both by the sudden increase of the first adiabatic invariant and the large initial increase in ∣p_⊥∣ relative to ∣p_∥∣.

Zoom In Zoom Out Reset image size

At around Ω_ci t ∼ 2.5, the electron is reflected from the shock with the first adiabatic moment remaining roughly constant, indicating magnetic mirroring. This two-step reflection process is that described by Amano & Hoshino (2007), who showed that electrons are initially energized in the perpendicular direction by shock-surfing acceleration in the shock foot and are then further accelerated and reflected by magnetic mirroring (or SDA). The initial perpendicular acceleration from SSA helps to maintain a large pitch angle, which results in a higher probability of the electron fulfilling the requirements for reflection instead of transmission through the shock.

Once enough electrons have been reflected in this manner, their velocity difference relative to the upstream electrons can excite the aforementioned electron acoustic instability and generate EAWs in the electron foreshock.

Having discussed the general properties of the electron foreshock in the previous section, we discuss what happens upstream following the production of EAWs. To achieve this, we take advantage of the particle-tracing feature in our code, periodically recording the locations, velocity, and local-field components for each particle. In this section, we once again focus on our in-plane simulations for the reason that, by definition, the electron acoustic waves generated by the reflected electrons reside in the simulation plane, and so are well captured when ϕ = 0.

For each simulation, we trace a sample of 10,000 electrons. To ensure that we can trace the entire evolutionary history of these particles beginning with their injection into the simulation upstream until after their reflection, we select particles that already reside in the foreshock (so have already been reflected) at Ω_ci t = 5. Our code can then be rerun at earlier time steps and record data for these particles from their injection into the upstream up until the end of the simulation. We select foreshock electrons by imposing a number of criteria. First, we select particles more than 5λ_i further upstream of the shock-foot location, x_foot. This prevents electrons that are currently undergoing shock-surfing acceleration, but may end up downstream, from being selected. Second, we ensure that the magnitude of the x component of their velocity in the upstream rest frame exceeds the shock velocity, v_e,x > u_sh. This guarantees they can, unless they are decelerated, outrun the shock and that at some point they must have been reflected. Finally, we also impose the condition that their velocity y-component must be positive, v_y > 0, as reflected electrons should be tied to the upstream magnetic field, which in these simulations always lies in the first quadrant of the x–y plane. The individual and unique IDs of all electrons fulfilling these criteria are output, and a random selection is made to ensure different behaviors are captured and limit any selection bias.

We illustrate our main results in Figures 12 and 13, which summarize the results for the simulations with θ_Bn = 55° and θ_Bn = 45°, respectively. In both figures, the top panel shows the density profile in the simulation box averaged over the transverse direction as a color map that is plotted as a function of time and the distance to the shock, x_sh, which is calculated assuming it moves with its analytical velocity value of ${v}_{\mathrm{sh}}^{\prime} \approx 0.067\,c$ in the simulation frame. To show the effect of the electrostatic waves present in x–y (simulation) plane of the electron foreshock, the second panel shows ${\mathrm{log}}_{10}({E}_{x}^{2}+{E}_{y}^{2})$. In both figures we see a clear shock transition, which is sharper in the density profile for θ_Bn = 45° relative to θ_Bn = 55°, and likewise stronger in ${\mathrm{log}}_{10}({E}_{x}^{2}+{E}_{y}^{2})$, as a consequence of the higher reflection rate and larger total energy content carried by the reflected electrons (see Figure 6).

Zoom In Zoom Out Reset image size

These upper panels have the trajectories of selected reflected electrons superimposed. Together, they encapsulate the range of physical behaviors that electrons experience when traveling through the electron foreshock and encountering the EAWs. The lower panels illustrate the momentum history of the selected electrons, separated into the components parallel, p_∥, and perpendicular, p_⊥, to the upstream magnetic field. The color and line style of each electron trajectory in the upper panel matches the lower panel corresponding to it, with each electron assigned a label to the right of each lower panel. It is immediately apparent that a wide range of possible trajectories, with accompanying changes in particle momentum, is possible. The physical origin of each case is discussed in turn.

We see from Figures 12 and 13 that some of the selected electrons exhibit sharp increases in p_∥ while p_⊥ remains approximately constant, which is indicative of them being reflected by magnetic mirroring. This effect is most profound for electrons 1₅₅ and 4₅₅, between Ω_ci t = 4−5, and is also apparent for electron 2₄₅ at Ω_ci t ≈ 1.5. While the momentum history of these particles is comparable around the time of mirroring, they initially also undergo some shock-surfing, experiencing the reflection mechanism described in Amano & Hoshino (2007) (see Section 4.3). As already discussed, we find that this two-step acceleration plays a role at early times and helps to reflect the initial foreshock electrons.

Zoom In Zoom Out Reset image size

Figures 12 and 13 reveal further changes in the particle momenta after reflection. This effect is particularly noticeable in the parallel direction; as such, the work required cannot be done by the motional electric field and instead comes from the EAWs in the electron foreshock oriented, on average, along the direction of the upstream magnetic field (because reflected electrons travel, on average, in this direction). We note that the values of p_∥ and p_⊥ typically do not remain constant after the initial reflection, ultimately leading a range of possible outcomes within the Ω_ci t = 7.8 time span of our simulations. We attribute this behavior to electrons becoming trapped by an electrostatic wave and being convected either back toward or farther away from the shock. For electrons 1₅₅, 2₅₅, and 2₄₅ these associated momentum changes appear to be relatively inconsequential, with the parallel particle momenta jumping around a mean value with respect to time. For electron 2₅₅, a reduction in p_∥ can be seen at late times in the top two panels of Figure 12, with the gradient of the blue line decreasing, allowing the relative distance between this electron and the shock to become smaller. Electron 3₅₅ initially undergoes steady parallel acceleration before shock-surfing and ultimately obtains enough energy to be reflected. During this stage, the parallel component of its velocity matches the phase velocity of the EAWs (Bohdan et al. 2022), supporting the notion that it is Landau trapped. Eventually, while this electron is spatially separated from the shock at Ω_ci t ≈ 3−6, it is advected back toward the shock, losing p_∥ at Ω_ci t ≈ 5. The relative percentage momentum change is significantly smaller for faster moving electrons such as electron 1₅₅, which can be explained by limited Landau damping as the velocity is not in resonance with the phase velocity of the EAWs.

Multiple reflection events are not the only possible outcomes for foreshock-scattered electrons. Electron 4₅₅ ends up in the downstream region some four gyrotimes later than the initial reflection at Ω_ci t ≈ 4. Here, the electron re-interacts with the shock with significantly more energy than it was injected into the upstream with, which can result in further and more efficient acceleration via stochastic Fermi acceleration (Bohdan et al. 2019b) or stochastic SDA (Matsumoto et al. 2017; Katou & Amano 2019). This is important, as the extra energy they gain from the earlier acceleration processes could allow them to more easily pass into the downstream than compared to if they had thermal energies.

Another interesting feature of Figure 13 is the differing trajectories of electrons 1₄₅ and 3₄₅. The location history of these electrons shows that they do not get within 50λ_i of the shock discontinuity yet are still turned back upstream. We note that trajectories like this did also occur for θ_Bn = 55°, yet, due to the less energetic foreshock, were somewhat rarer, occurring for only 70 of our 10,000 traced electrons in comparison to ∼15% for θ_Bn = 45°. For magnetic mirroring to be the mechanism responsible for reflection, we need a converging magnetic field, which is not present in the electrostatic foreshock.

We further demonstrate the ability of the EAWs to turn electrons back upstream in Figure 14, which shows the time evolution of p_∥ between Ω_ci t = 3 and Ω_ci t = 4 in panel (a). This interval was selected as it is during this period that the electron is turned away from the shock, despite always residing in the upstream. The quantities here are measured in the upstream rest frame. The accompanying panels display the corresponding time evolution of (b) E_∥, (c) d p_∥/d t, (d) B, and (e) μ. Panel (f) illustrates a strong correlation between dp_∥/dt and E_∥, suggesting that this field component, which is associated with the EAWs, is responsible for the change in p_∥. This is further reinforced by the fact that the local magnetic field experienced by the electron is consistent with the upstream value and shows no correlation with d p_∥/d t, and that the first adiabatic invariant is not constant and in fact steadily increases during the reflection, ruling out magnetic mirroring.

Zoom In Zoom Out Reset image size

We propose the following overall picture: Initially, electrons undergo shock-surfing acceleration on Buneman waves and are reflected via magnetic mirroring. If enough electrons are reflected with large-enough total energy content, we begin to see electron acoustic waves traveling into the shock foot, with the overall energy density of these waves increasing with the reflection rate, hence is higher closer to the shock, and at later times once the reflection rate begins to saturate. These waves have a k-vector approximately aligned with the upstream magnetic field and can do work parallel to it. On average, electrons encountering these waves will experience a small net force in the direction of the upstream magnetic field, which leads to their momenta being aligned increasingly toward it instead of in the direction of the bulk flow (toward the shock). Analogous to Buneman waves, it is possible for upstream electrons to become trapped in the potential wells of the EAWs if the potential difference is sufficiently large. Acceleration by the motional electric field leads to an increase in p_⊥, as is the case for electron 3₄₅. Furthermore, they can be accelerated by the on-average parallel electric field. Subsequent interactions between these electrons and the EAWs can lead to them being turned farther away from the shock, or convected back toward it, where they may be caught up with the shock and again interact with it.

As with Buneman waves, the EAWs can accelerate upstream electrons and even reflect them back upstream. The reflection of electrons in this manner will again help to drive further upstream turbulence and contribute to further electron reflection and preacceleration. This is indicated in Figure 15, which illustrates that the mean reflection location of an electron with respect to x_foot increases as a function of time. Here, for each reflected electron, the reflection point was defined as the minimum distance from x_foot and only includes the first reflection point for electrons with multiple shock interactions.

Zoom In Zoom Out Reset image size

By using our tracing data in combination with the definition of the shock foot, x_foot, we can estimate how often reflected electrons interact with the shock. We define a single shock interaction when an electron crosses from x > x_foot to x ≤ x_foot. In general, for all simulations, 70%–80% of all reflected electrons have a single interaction, indicating the reflection mechanism is that described in Amano & Hoshino (2007). The percentage of all reflected electrons with two or more interactions was ∼30% for θ_Bn = 63°, decreasing to ∼9% for θ_Bn = 45°. This is because at higher obliquities, electrons require a larger speed to escape the shock and therefore need more preacceleration, such as more cycles of SDA. In comparison, the higher EAW energy density contained in EAWs means that the percentage of reflected electrons turned back by interacting with them was higher at smaller θ_Bn, with ∼15% for θ_Bn = 45° and ∼0.7% for θ_Bn = 55° and negligible for θ_Bn ≥ 55°. For θ_Bn = 30°, the shock takes longer to form (e.g., Gingell et al. 2017) and is in fact dominated by an ion foreshock close to the foot, making a comparison challenging and beyond the scope of this work.

In Bohdan et al. (2022), a similar simulation to those here is presented, with the main differences a wider simulation box, a total runtime of Ω_ci t = 50.4, and an initial upstream magnetic field configuration characterized by θ_Bn = 60° and ϕ = 90°. By integrating the energy density of shock-reflected electrons up to a threshold value where EAWs can be detected, it is estimated that the EAW region will saturate at approximately Ω_ci t ∼ 270. Furthermore, similarly to our θ_Bn = 30° simulation (although with a later onset), substantial ion reflection begins in the Bohdan et al. (2022) simulation at around Ω_ci t ∼ 10, which causes Whistler waves to form immediately ahead of the shock transition. We expect such features to occur in the simulations presented here at later times and aim to investigate preacceleration associated with the Whistlers in future work.

In this work, we have performed PIC simulations that have revealed the presence of electrostatic waves, identified as electron acoustic waves, which are generated when reflected electrons excite the electron acoustic instability. We have investigated the properties of these waves and their influence on the upstream electrons, with the key findings of our work summarized below:

• 1.  
  PIC simulations were completed for the range of angles between the shock normal and upstream magnetic field of ${\theta }_{\mathrm{Bn}}=\left[30,45,55,63,74.3\right]\,\mathrm{degrees}$. This allowed for a detailed study of the rates of incident electrons that are reflected back upstream and an assessment of their effect on the upstream plasma.
• 2.  
  A higher percentage of incoming electrons is reflected back upstream for smaller θ_Bn. About 5% of electrons were reflected for θ_Bn = 30°, but only around 0.05% for θ_Bn = 63°.
• 3.  
  The initial reflection mechanism involves shock-surfing acceleration on Buneman waves followed by magnetic mirroring and is described in Amano & Hoshino (2007). Our results show a higher reflection rate for out-of-plane simulations, verifying the role of Buneman waves in the initial electron reflection.
• 4.  
  Although fewer electrons are reflected at higher θ_Bn, the typical peak energy of these electrons is higher. This is because they require velocities ${v}_{e}\cos {\theta }_{\mathrm{Bn}}\gt {v}_{\mathrm{sh}}$. Despite this, the greater number of electrons reflected in shocks with smaller θ_Bn means overall there is a greater energy density in reflected electrons for these angles.
• 5.  
  Reflected electrons lead to the production of electron acoustic waves in the electron foreshock on spatial scales comparable to Buneman waves.
• 6.  
  These electron acoustic waves are themselves capable of accelerating and reflecting upstream electrons farther away from the shock, preventing them from reaching it and compromising injection for ∼1% of upstream electrons for θ_Bn ≤ 45°.

M.P. acknowledges support by DFG through grant PO 1508/10-1. The numerical simulations were conducted on resources provided by the North-German Supercomputing Alliance (HLRN) under project bbp00033.