Theory of Electron Injection at Oblique Shock of Finite Thickness

A natural frame of reference to analyze a magnetized oblique shock wave is the so-called normal incident frame (NIF), in which the shock is at rest and the upstream flow is parallel to the shock normal. We define the shock speed as the upstream flow speed U _s in the NIF. The magnetic obliquity θ_Bn is defined as the angle between the upstream magnetic field and the shock normal direction. In the following, we will always work in HTF, in which the plasma flow and the magnetic field are parallel to each other. Since the motional electric field vanishes in this special frame, it is one of the most convenient frames to analyze the transport of charged particles around a collisionless oblique shock. Such a frame may be found by frame transformation in the direction transverse to the shock with a speed of ${U}_{s}\tan {\theta }_{{Bn}}$. The transformation speed must obviously be less than the speed of light ${U}_{s}\tan {\theta }_{{Bn}}\lt c$ (where c is the speed of light). Such a shock is called subluminal. On the other hand, if ${U}_{s}\tan {\theta }_{{Bn}}\geqslant c$, we cannot find HTF, and the shock is called superluminal. In this paper, we will consider only subluminal shocks because most of the nonrelativistic oblique shocks of practical interest are subluminal. For instance, θ_Bn < 89° is sufficient for a young SNR shock of a shock speed of U _s = 3000 km s⁻¹ to be subluminal. Appendix A shows that the transport equation conventionally written in an arbitrary frame can be cast into a simpler form in HTF. As in the case of single-particle dynamics, the explicit use of HTF makes it easy to understand the physics involved.

In the following, we will denote the flow velocity in HTF by V . It is, by definition, always parallel to the local magnetic field, i.e., V = V b , where b is the magnetic field unit vector. The flow velocity in NIF is denoted by U .

We investigate the electron injection at a plane oblique shock of finite thickness of order L _s by considering the steady-state solution (i.e., ∂/∂t = 0) of the diffusion–convection equation

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial f}{\partial t}+V\cos \theta \,\displaystyle \frac{\partial f}{\partial x}+\displaystyle \frac{1}{3}\left(\displaystyle \frac{\partial \mathrm{ln}B}{\partial x}-\displaystyle \frac{\partial \mathrm{ln}V}{\partial x}\right)V\cos \theta \displaystyle \frac{\partial f}{\partial \mathrm{ln}p}\\ \quad =\,\displaystyle \frac{\partial }{\partial x}\left(\kappa {\cos }^{2}\theta \,\displaystyle \frac{\partial f}{\partial x}\right),\end{array}\end{eqnarray} \tag{ 1 }$$

where f = f(x, p) is the isotropic part of the distribution function as a function of position x and plasma-frame momentum p. The meanings of other quantities are as follows: B = B(x) is the magnetic field strength, θ = θ(x) is the angle between the magnetic field and x-axis, V = V(x) is the field-aligned flow speed, and κ = κ(x, p) is the diffusion coefficient along the local magnetic field. Note that Equation (1) is equivalent to Equation (A8) except that we have dropped the subscript for f₀, as we will consider only the isotropic part throughout this paper (see Appendix A for details). We assume that the shock normal is parallel to the x-axis and that the flow velocity and the magnetic field strength change smoothly from the upstream values V₁, B₁, θ₁( =θ_Bn ) at x → − ∞ to the downstream values V₂, B₂, θ₂ at x → +∞. Note that the upstream flow speed in HTF is given by ${V}_{1}={U}_{s}/\cos {\theta }_{{Bn}}$.

We find it convenient to use the following form to estimate the spectral index:

$$\begin{eqnarray}&&{\int }_{{X}_{1}}^{{X}_{2}}\displaystyle \frac{\partial {V}_{x}}{\partial x}\left[f+\displaystyle \frac{1}{3}\displaystyle \frac{\partial f}{\partial \mathrm{ln}p}\right]{dx}={\left[{V}_{x}f-{\kappa }_{{xx}}\displaystyle \frac{\partial f}{\partial x}\right]}_{{X}_{1}}^{{X}_{2}},\end{eqnarray} \tag{ 2 }$$

which is obtained by integrating Equation (1) with respect to x on an arbitrary interval [X₁, X₂]. Note that we have introduced the notations ${V}_{x}=V\cos \theta$ and ${\kappa }_{{xx}}=\kappa {\cos }^{2}\theta$.

It is easy to find from Equation (1) that the characteristic scale length is given by ${l}_{\mathrm{diff}}={\kappa }_{{xx}}/{V}_{x}=\kappa {\cos }^{2}\theta /{U}_{{\rm{s}}}$, which we call the diffusion length. More specifically, the spatial profile is given by $f(x)\propto \exp \left(x/{l}_{\mathrm{diff}}\right)$ for a system that is approximately homogeneous. We thus find that the particle acceleration at a shock of finite thickness can be classified by the ratio between the diffusion length and the shock thickness l_diff/L _s .

We first look at the conventional DSA as an obvious example. If the diffusion length for a particular particle energy is much larger than the thickness of the shock l_diff/L _s ≫ 1, the shock as seen by the particle is essentially a discontinuity in the plasma flow:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {V}_{x}}{\partial x}\approx \left({V}_{2}\cos {\theta }_{2}-{V}_{1}\cos {\theta }_{1}\right)\,\delta (x).\end{eqnarray} \tag{ 3 }$$

We adopt the boundary condition such that f = 0 at x = X₁ and ∂f/∂x = 0 at x = X₂. Carrying out the integral in Equation (2) in the limit of infinite integration interval X₁ → −∞ and X₂ → +∞ , we have

$$\begin{eqnarray}&&\begin{array}{l}\left({V}_{1}\cos {\theta }_{1}-{V}_{2}\cos {\theta }_{2}\right)\left[{f}_{2}+\displaystyle \frac{1}{3}\displaystyle \frac{\partial {f}_{2}}{\partial \mathrm{ln}p}\right]\\ \quad =\,{V}_{2}\cos {\theta }_{2}{f}_{2},\end{array}\end{eqnarray} \tag{ 4 }$$

where f₂ = f(x = X₂, p). It becomes immediately clear that the spectrum becomes a power law f₂(p) ∝ p^−q , with its index given by

$$\begin{eqnarray}&&q=\displaystyle \frac{3{V}_{1}\cos {\theta }_{1}}{{V}_{2}\cos {\theta }_{2}-{V}_{1}\cos {\theta }_{1}}=\displaystyle \frac{3r}{r-1},\end{eqnarray} \tag{ 5 }$$

where $r={V}_{1}\cos {\theta }_{1}/{V}_{2}\cos {\theta }_{2}={V}_{x,1}/{V}_{x,2}$ is the compression ratio. It is well known that the spectral index of particles accelerated at an oblique shock is determined solely by the compression ratio and independent of the magnetic obliquity θ_Bn even though the energy gain clearly includes both the first-order Fermi acceleration and SDA effects (e.g., Drury 1983).

Let us now consider the general case in which the effect of a finite shock thickness needs to be taken into account. As in the case of the discontinuous shock solution, we may anticipate that the particle density increases within the shock from the far upstream value and then becomes constant in the downstream. We thus use the same boundary condition as DSA. By assuming the solution of the form f(x, p) = g(x)p^−q , we estimate the spectral index by

$$\begin{eqnarray}&&q=3\left[1-{\left({\int }_{{X}_{1}}^{{X}_{2}}\frac{g(x)}{{g}_{2}}\frac{\partial }{\partial x}\left(\frac{{V}_{x}}{{V}_{x,2}}\right){dx}\right)}^{-1}\right],\end{eqnarray} \tag{ 6 }$$

where g₂ = g(x = X₂). Note that the integration interval [X₁, X₂] is arbitrary as long as it contains the entire particle acceleration region. It is easy to see that the canonical DSA spectral index may be recovered for a discontinuous flow profile. For more general cases, the evaluation of the integral requires a specific form of g(x) that will only be known by solving Equation (1). Nevertheless, this formula is useful for comparisons with in situ observations and kinetic simulations because in both cases the profiles of g(x) and V_x (x) are already known.

It is instructive to proceed further to obtain an analytic estimate for the spectral index without solving Equation (1). By looking at fully numerical solutions (see Section 3), we have found that the exponential profile $g(x)\propto \exp (x/{l}_{\mathrm{diff}})$ is approximately correct even within the shock transition layer ∂V_x /∂x ≠ 0. Since the condition ${l}_{\mathrm{diff}}/{L}_{s}=\kappa \cos {\theta }_{{Bn}}^{2}/{U}_{s}\lesssim 1$ favors quasi-perpendicular shocks, we here only consider particle acceleration at a finite-thickness quasi-perpendicular shock, which we call SSDA. In this case, we may ignore the particle energy changes due to the flow compression (V₁ ≈ V₂) and the kink in the magnetic field line (θ₁ ≈ θ₂) at quasi-perpendicular shocks. We then have

$$\begin{eqnarray}&&\displaystyle \frac{\partial {V}_{x}}{\partial x}\approx -{V}_{1}\cos {\theta }_{1}\displaystyle \frac{\partial \mathrm{ln}B}{\partial x}.\end{eqnarray} \tag{ 7 }$$

If we further adopt an average magnetic field gradient $\partial \mathrm{ln}B/\partial x\approx \left\langle \partial \mathrm{ln}B/\partial x\right\rangle \approx \mathrm{const}.$ to evaluate the integral, we finally obtain the following estimate for the spectral index:

$$\begin{eqnarray}&&q\approx 3\left[1+{\left({l}_{\mathrm{diff}}\left\langle \displaystyle \frac{\partial \mathrm{ln}B}{\partial x}\right\rangle \right)}^{-1}\right].\end{eqnarray} \tag{ 8 }$$

It is determined by the ratio between the diffusion length and the magnetic field transition scale length ${\left({\rm{\partial }}{\rm{l}}{\rm{n}}B/{\rm{\partial }}x\right)}^{-1}\,\sim {L}_{s}$. Since, in general, the diffusion coefficient is energy dependent, the spectrum of the accelerated particles will not be described by a single power law. A power-law spectrum will be obtained only when the diffusion length l_diff (and thus the diffusion coefficient κ) is independent of energy.

It is easy to see that the spectral index for SSDA is steeper than the DSA prediction. Since, by definition, the diffusion length should be smaller than the shock thickness l_diff ≲ L _s , it should have a lower limit q ≳ 6. Even the hardest spectrum predicted for SSDA is softer than q = 4 predicted by DSA. The spectrum becomes even softer as the diffusion length decreases.

The steep spectrum for particles with small diffusion lengths is a natural consequence of efficient particle confinement. In other words, the convective escape time of particles τ_esc ≈ r(l_diff/U _s ) from the acceleration region becomes shorter as the confinement efficiency increases, while the acceleration timescale τ_acc ≈ 3L _s /U _s remains unchanged. This is in clear contrast to DSA, in which τ_acc ≈ 3r/(r − 1)(l_diff/U _s ) has the same dependence on l_diff. This is indeed the reason why DSA gives the nearly universal power law that is independent of the diffusion coefficient: q = τ_acc/τ_esc + 3 = 3r/(r − 1) (see, e.g., Drury 1991).

It is important to point out that SSDA is a fast process instead. Since the particle acceleration via SSDA operates locally within the shock transition layer, there is no need for particles bouncing back and forth across the shock, thereby saving time. Note that high-frequency waves, which play the role of scattering of low-energy electrons, are most intense in the close vicinity of the shock (e.g., Hull et al. 2012; Wilson et al. 2014; Oka et al. 2017; Amano et al. 2020). Therefore, it does not require high-frequency wave activity far from the shock, where the waves will be subject to strong collisionless damping. Considering that the diffusion length may, in general, be considered as a monotonically increasing function of energy, one may expect that SSDA and DSA tend to be dominant mechanisms for particle acceleration at low energy and high energy, respectively.

We should mention that the spectral index given by Equation (8) is similar to, but different from, the estimate given by Katou & Amano (2019) (in which the spectral index is defined for N(p) = 4π p² f(p) ∝ p^−q+2 for nonrelativistic particles). More specifically, l_diff in Equation (8) is replaced by the typical shock thickness L _s in Katou & Amano (2019). This difference arises because they estimated the typical escape time by τ_esc ≈ L _s /U _s , i.e., the convection timescale over the shock thickness L _s , rather than the diffusion length l_diff. The erroneous expression should be considered appropriate only when the diffusion length is energy independent and comparable to the shock thickness. Interestingly enough, in situ measurements at Earth's bow shock suggested that both of these conditions were reasonably satisfied in the energy range where the spectrum is described by a power law (Amano et al. 2020), which is fully consistent with the theory.

As we have seen in previous subsections, the acceleration of particles for a given diffusion coefficient may be understood in terms of DSA or SSDA. So far, we have developed the general discussion without specifying the shock thickness, particle species, and particle energy. From now on, we consider the specific application to the electron injection problem.

It is well known that the transition layer thickness of supercritical collisionless shocks is of the order of the ion gyroradius U _s /Ω_ci, where Ω_ci is the ion cyclotron frequency in the upstream (e.g., Livesey et al. 1984). Note that the thickness always refers to the length scale of the shock transition layer including the foot, ramp, and overshoot in this paper. We thus define the thickness by L _s = η U _s /Ω_ci, with η being a numerical factor of order unity. The ratio between the diffusion length and the shock thickness is thus given by

$$\begin{eqnarray}&&\frac{{l}_{\mathrm{diff}}}{{L}_{s}}=\frac{1}{6\eta }{\left(\frac{{D}_{\mu \mu }}{{{\rm{\Omega }}}_{\mathrm{ci}}}\right)}^{-1}{\left(\frac{v}{{U}_{s}/\cos {\theta }_{{Bn}}}\right)}^{2}.\end{eqnarray} \tag{ 9 }$$

To evaluate the quantity, we may consider that the minimum velocity of the particles being accelerated corresponds to that expected by adiabatic SDA $v\gtrsim {U}_{s}/\cos {\theta }_{{Bn}}$. Since the scattering rate will be given at most by the Bohm limit, ions with D_{μ μ} /Ω_ci ≲ 1 should not be able to satisfy the condition unless the shock is excessively thicker than normal (η ≫ 1). On the other hand, the scattering of electrons can be much faster D_{μ μ} /Ω_ci ≲ m_i /m_e , where m_i /m_e is the ion-to-electron mass ratio. Therefore, low-energy electrons may be subject to particle acceleration via SSDA within the shock transition layer.

Consider a low-energy electron that is being accelerated by SSDA. Because of the upper limit on the scattering rate, there should be the maximum energy ${E}_{\max ,\ \mathrm{SSDA}}$ that can be achievable through SSDA. This may be determined by the condition l_diff/L _s ≃ 1, yielding

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{\max ,\ \mathrm{SSDA}}}{{E}_{0}}\simeq 6\eta \left(\displaystyle \frac{{D}_{\mu \mu }}{{{\rm{\Omega }}}_{\mathrm{ci}}}\right),\end{eqnarray} \tag{ 10 }$$

where ${E}_{0}={m}_{e}{({U}_{s}/\cos {\theta }_{{Bn}})}^{2}/2$. Note that this expression is appropriate only for nonrelativistic electrons. It is clear that electrons beyond this energy cannot be trapped inside the shock transition layer but may potentially be accelerated even further via DSA in a much larger spatial extent. In other words, one may expect that a steep spectrum will develop at low energy $E\lesssim {E}_{\max ,\ \mathrm{SSDA}}$, which may then be smoothly connected to a harder spectrum at high energy $E\gtrsim {E}_{\max ,\ \mathrm{SSDA}}$, as predicted by the standard DSA theory. This is the injection scheme that we consider in this study.

If we consider the particle acceleration problem in the infinitely large system, this injection must always happen as long as the scattering rate is finite. The idealized situation is certainly not satisfied in reality. One has to estimate the scattering efficiency to understand whether or not the proposed scheme provides the solution to the injection problem in realistic systems of finite size or finite age.

Let us now consider particle acceleration in a system of finite spatial extent. In DSA theory, the typical scale size of the particle acceleration region is given by the diffusion length l_diff, which is, in general, energy dependent. If the diffusion length is comparable to or larger than the system size, no particle acceleration should be expected. Alternatively, one may think of a time-dependent particle acceleration problem. We see from ${\tau }_{\mathrm{acc}}\propto {l}_{\mathrm{diff}}\propto {({D}_{\mu \mu }/{{\rm{\Omega }}}_{\mathrm{ci}})}^{-1}$ that the acceleration time becomes longer if the scattering is weaker. Again, particle acceleration will not take place for those particles having a typical acceleration time comparable to or longer than the age of the system. In both cases, the energy dependence of the scattering rate D_{μ μ} determines whether or not DSA is capable of accelerating the particles of a given energy. Therefore, the injection from SSDA to DSA requires that the particles escaped out from SSDA (i.e., those with energies $E\sim {E}_{\max ,\ \mathrm{SSDA}}$) experience strong enough scattering outside the shock transition layer.

If one considers the standard DSA theory in which the particle scattering outside the shock transition layer is provided by cyclotron resonance with low-frequency magnetostatic turbulence with a power-law spectrum, the scattering rate D_{μ μ} will be larger at higher energies. In other words, the resonance condition kr _g ∼ 1 (where r _g denotes the particle gyroradius) indicates that the resonant wavelength becomes longer with increasing particle energy. Therefore, higher-energy particles will interact with more intense fluctuations at longer wavelengths. Since thermal and slightly suprathermal electrons have gyroradii much smaller than the typical dissipation length scale of the turbulence, there will be essentially no power available for their scattering. This is the severest obstacle prohibiting the injection of low-energy electrons into DSA. It is clear that the injection threshold energy depends crucially on the assumption of the fluctuation power spectrum both upstream and downstream of the shock.

First, we consider preexisting turbulence in the unshocked medium or that is generated by the shock-accelerated ions. In the former case, a power-law wave spectrum will develop in between the energy injection scale at a long wavelength and the dissipation scale. One may naively estimate that the dissipation scale is determined by the cyclotron damping of thermal ions, which becomes strong at the ion inertial length kc/ω_pi ∼ 1 in a low-beta plasma or the thermal ion gyroradius ${kc}/{\omega }_{\mathrm{pi}}\sim {\beta }_{i}^{-1/2}$ in a high-beta plasma, respectively. Note that ω_pi and β _i denote the ion plasma frequency and plasma beta, respectively. In either case, it is reasonable to assume that there will be finite wave power for scattering at a wavelength longer than the ion scale. We may apply the same argument for ion-generated turbulence. The primary scale at which specularly reflected ions ($v\sim {U}_{s}/\cos {\theta }_{{Bn}}$) generate large-amplitude fluctuations may be estimated as ${kc}/{\omega }_{\mathrm{pi}}\sim \cos \theta /{M}_{{\rm{A}}}$. While this appears at a very long wavelength, the development of turbulence will eventually generate fluctuations down to the ion scale. Although we do not exclude the possibility that the spectrum develops into even smaller scales, the integrated power at the electron scale will drop off by orders of magnitude. It is thus reasonable to assume that the turbulence spectrum will have a break at the critical wavenumber ${k}_{* }c/{\omega }_{\mathrm{pi}}\sim \min (1,{\beta }_{i}^{-1/2})$. Electrons should have sufficiently high energies so that they interact with fluctuations with wavelengths longer than the dissipation scale.

Another possibility is the self-generation of turbulence by the shock-accelerated electrons themselves. The energetic electrons that diffuse ahead of the shock may self-generate left-hand-polarized Alfvén waves (i.e., Alfvén ion cyclotron waves) propagating parallel to the ambient magnetic field via the cyclotron resonance. Note that this is exactly the same process that has been considered for the generation of right-hand-polarized waves by shock-accelerated ions. The self-generation by the energetic electron streaming will be possible only when the instability overcomes the ion cyclotron damping that becomes strong beyond the critical wavenumber k_* c/ω_pi. We then arrive at the same conclusion that the electrons should have sufficiently high energies for injection such that they can resonantly interact with long-wavelength fluctuations below the critical wavenumber ${k}_{* }c/{\omega }_{\mathrm{pi}}\sim \min (1,{\beta }_{i}^{-1/2})$.

We should nevertheless also mention that there are possibilities to self-generate (right-hand-polarized) higher-frequency whistler waves (either parallel or oblique propagation) by lower-energy electrons (Levinson 1992, 1996; Amano & Hoshino 2010). These models consider that DSA may be activated for low-energy electrons if the self-generation of high-frequency waves is realized. However, the growth rates expected for these instabilities can be very high, and the wave growth and resulting scattering will happen within the shock transition layer. Therefore, we think that these waves will play a role of scattering for SSDA rather than DSA operating outside the shock transition layer. We will discuss this issue in more detail later in Section 4.

In summary, we conclude that the threshold energy required for the electron injection is determined by the resonance condition with the wave at the critical wavenumber k_* c/ω_pi. Since the corresponding wave frequency is sufficiently low, we can use the magnetostatic approximation ω = kv − Ω_ce/γ ≈ 0 (where Ω_ce is the electron cyclotron frequency). We then have the condition

$$\begin{eqnarray}&&\displaystyle \frac{\gamma v}{c}\gtrsim {\left(\displaystyle \frac{{k}_{* }c}{{\omega }_{\mathrm{pi}}}\right)}^{-1}\left(\displaystyle \frac{{m}_{i}}{{m}_{e}}\right)\left(\displaystyle \frac{{V}_{{\rm{A}}}}{c}\right)\end{eqnarray} \tag{ 11 }$$

as the threshold velocity for electron injection, where V_A is the Alfvén speed. We see that the threshold energy is determined only by the Alfvén speed for a given critical wavenumber k_*. The corresponding threshold energy as a function of the Alfvén speed is shown in Figure 1 for three different parameters of k_* c/ω_pi = 1, 0.1, and 0.01. For typical Alfvén speeds of V_A ∼ 100 km s⁻¹ in the interstellar medium and the solar wind at 1 au, we may estimate the injection threshold energy as 0.1–1 MeV by assuming that the resonant scattering starts to occur below a reasonably small wavenumber k_* c/ω_pi ∼ 0.1–1.

Zoom In Zoom Out Reset image size

It is worth noting that the often-quoted condition that the particle gyroradius is larger than the shock thickness r _g ≳ L _s is not necessarily relevant for determining the injection threshold. To understand this, let us introduce the mean free path normal to the shock surface ${\lambda }_{\mathrm{mfp}}=v\cos \theta /2{D}_{\mu \mu }={r}_{g}{({D}_{\mu \mu }/{{\rm{\Omega }}}_{\mathrm{ce}})}^{-1}\cos \theta /2$. For low-energy electrons to which SSDA applies, we have ${\lambda }_{\mathrm{mfp}}/{l}_{\mathrm{diff}}=3{U}_{s}/v\cos \theta \ll 1$ because $v\gg {U}_{s}/\cos \theta$ for the accelerated particles at a nonrelativistic shock. As the electron energy increases, the dominant acceleration mechanism will change from SSDA to DSA, which nevertheless happens without violating the condition λ_mfp/L _s ≪ l_diff/L _s ∼ 1. In other words, the particle acceleration for small gyroradius electrons can continue under the diffusion approximation. As the energy increases even further, λ_mfp/L _s ≳ 1 is reached, for which the electron interaction with the shock is nearly scatter-free. In this case, substantial anisotropy may develop near the shock front because large pitch-angle particles coming from the upstream will be reflected back by the shock via adiabatic SDA, while those coming from the downstream will always be transmitted to the upstream. The strong anisotropy will nevertheless be relaxed a few mean free paths away from the shock, and the overall particle acceleration may be described again by the diffusion approximation. It is well known that the adiabatic reflection effect is integrated into DSA and will not impose any restrictions on the particle energy. Therefore, if the transition between SSDA and DSA happens under the strong scattering regime such that L _s ≫ λ_mfp ≳ r_g , the diffusion approximation is always applicable. As is clear from the above discussion, the validity of the diffusion approximation will be violated if the HTF shock speed becomes relativistic, ${U}_{s}/\cos \theta \sim c$, which thus requires a separate study. Related discussion on particle acceleration at an oblique shock can be found in, e.g., Drury (1983).

Let us consider the condition that an electron of a given momentum p can be accelerated by SSDA. We rewrite the maximum energy estimate (10) as the condition required for the Alfvén Mach number in HTF ${M}_{{\rm{A}}}^{\mathrm{HTF}}={M}_{{\rm{A}}}/\cos {\theta }_{{Bn}}$ in the following form:

$$\begin{eqnarray}&&{M}_{{\rm{A}}}^{\mathrm{HTF}}\gtrsim {\left(6\eta \right)}^{-1/2}{\left(\frac{{D}_{\mu \mu }(p)}{{{\rm{\Omega }}}_{\mathrm{ce}}}\right)}^{-1/2}\left(\frac{p}{{p}_{{\rm{A}},e}}\right),\end{eqnarray} \tag{ 12 }$$

where ${p}_{{\rm{A}},e}={m}_{e}{V}_{{\rm{A}},e}={({m}_{i}{m}_{e})}^{1/2}{V}_{{\rm{A}}}$ is the electron momentum associated with the electron Alfvén speed. Note that we have expressed the momentum dependence of the scattering rate D_{μ μ} (p) explicitly. Particle acceleration by SSDA continues as long as this condition is satisfied. Combining this with the injection threshold given by Equation (11), we estimate the condition for the electron injection:

$$\begin{eqnarray}&&{M}_{{\rm{A}}}^{\mathrm{HTF}}\gtrsim {\left(6\eta \right)}^{-1/2}{\left(\displaystyle \frac{{D}_{\mu \mu }(p)}{{{\rm{\Omega }}}_{\mathrm{ce}}}\right)}^{-1/2}{\left(\displaystyle \frac{{m}_{i}}{{m}_{e}}\right)}^{1/2}{\left(\displaystyle \frac{{k}_{* }c}{{\omega }_{\mathrm{pi}}}\right)}^{-1},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\approx 55{\left(\displaystyle \frac{\eta }{1}\right)}^{-1/2}{\left(\displaystyle \frac{{D}_{\mu \mu }(p)/{{\rm{\Omega }}}_{\mathrm{ce}}}{0.1}\right)}^{-1/2}{\left(\displaystyle \frac{{k}_{* }c/{\omega }_{\mathrm{pi}}}{1}\right)}^{-1},\end{eqnarray} \tag{ 14 }$$

where the scattering rate D_{μ μ} (p) must be evaluated within the shock transition layer and for the threshold momentum corresponding to Equation (11).

The above conditions for the Mach number ${M}_{{\rm{A}}}^{\mathrm{HTF}}$ suggest that there will be three different regimes of electron acceleration for a given D_{μ μ} (p), if it is independent of ${M}_{{\rm{A}}}^{\mathrm{HTF}}$. First, if ${M}_{{\rm{A}}}^{\mathrm{HTF}}$ is too small and is not able to satisfy Equation (12) for slightly above the thermal momentum, electron acceleration via SSDA will not take place. Second, at higher ${M}_{{\rm{A}}}^{\mathrm{HTF}}$ shocks that satisfy Equation (12) but not Equation (13), SSDA will be able to accelerate electrons within the shock transition layer. However, the highest-energy electrons escaping out from the system will not experience further energization. Finally, if Equation (13) is satisfied, those escaping electrons have sufficiently high energies and will be subject to further particle acceleration via the standard DSA.

In reality, the scattering rate will be dependent on the shock parameters. Indeed, higher ${M}_{{\rm{A}}}^{\mathrm{HTF}}$ shocks may have more efficient scattering rates because of self-generated whistler-mode waves (Levinson 1992, 1996; Amano & Hoshino 2010). Therefore, the problem will be nonlinear and requires in-depth analyses. Nevertheless, it is true that quasi-parallel shocks need to have very high Mach numbers for the electron injection, M_A ≳ 17, even with the most optimistic Bohm scattering rate D_{μ μ} /Ω_ce ∼ 1. Nearly perpendicular shocks will be more favorable sites for the electron injection, as even moderately strong shocks with M_A ∼ 5 will be sufficient for injection if θ_Bn ∼ 85° and D_{μ μ} /Ω_ce ∼ 0.1. It is nevertheless important to keep in mind that the required Mach number is dependent on the dissipation scale k_* c/ω_pi surrounding the shock. The electron injection is regulated by the rate of scattering both inside and outside the shock transition layer.

In this section, we will discuss fully numerical solutions to the diffusion–convection Equation (1) in the steady state to check the validity of the discussion developed so far. Note that we do not intend to compare the results with specific observations or kinetic simulations. Therefore, we here adopt a simple analytic functional form of the shock profile. We assume that the magnetic field always lies in the x-z plane (i.e., the coplanarity plane) and use the profile for B_z (x) given by

$$\begin{eqnarray}&&{B}_{z}(x)=\displaystyle \frac{{B}_{z,2}-{B}_{z,1}}{2}\tanh \left(\displaystyle \frac{x}{{L}_{s}}\right)+\displaystyle \frac{{B}_{z,2}+{B}_{z,1}}{2}.\end{eqnarray} \tag{ 15 }$$

The shock compression ratio is given by r = B_z,2/B_z,1. Since a stationary MHD shock satisfies ${V}_{x}{B}_{z}=\mathrm{const}.$ across the shock, V_x (x) is readily obtained as

$$\begin{eqnarray}&&{V}_{x}(x)={V}_{x,1}\displaystyle \frac{{B}_{z,1}}{{B}_{z}(x)}.\end{eqnarray} \tag{ 16 }$$

The divergence of the plasma flow is thus given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {V}_{x}}{\partial x}=-{V}_{x}(x)\displaystyle \frac{\partial \mathrm{ln}{B}_{z}}{\partial x}.\end{eqnarray} \tag{ 17 }$$

Note that the conditions ${B}_{x}(x)=B\cos \theta ={B}_{1}\cos {\theta }_{{Bn}}=\mathrm{const}.$ and ${V}_{z}(x)=V\sin \theta ={V}_{1}\sin {\theta }_{{Bn}}=\mathrm{const}.$ always guarantee V × B = 0.

To represent a finite system size, we introduce the so-called free escaping boundary in the upstream of the shock. In other words, the upstream boundary condition is taken such that f(x) = 0 at x = −L_FEB, indicating that particle acceleration will not take place for those particles with l_diff/L_FEB ≳ 1. For the downstream boundary condition at x = +∞ , we use ∂f/∂x = 0 throughout in this study.

We wish to obtain the solution in the momentum range p₁ ≤ p ≤ p₂. Since ∇ · V ≤ 0 is satisfied throughout the region of interest, there will be no energy loss of the particles. This indicates that the boundary condition at the highest momentum p = p₂ is not needed. On the other hand, we have to specify the boundary condition at the lowest momentum p = p₁, for which we take the fixed boundary

$$\begin{eqnarray}&&f(x,{p}_{1})={f}_{\mathrm{inj}}\exp \left(-\displaystyle \frac{{x}^{2}}{2{\sigma }^{2}}\right),\end{eqnarray} \tag{ 18 }$$

where f_inj is constant. This indicates that there is an injection into the system that is spatially localized at around x = 0. Although we always use σ = L _s in this study, this will not affect the solution as long as we consider the momentum range sufficiently far away from p₁. Note that, since the problem is linear, f_inj can be chosen arbitrarily, which merely determines the normalization of the numerical solution.

The problem definition will be completed once a model for κ = κ(x, p) is given. Since the shock transition layer is the region where the main shock dissipation is taking place, the fluctuation level may be substantially different from the regions outside. Therefore, we will consider the solutions in these regions separately, which are then connected with each other through appropriate boundary conditions. More specifically, we divide the entire domain into the following three regions: (1) upstream − L_FEB ≤ x ≤ X₁, (2) shock transition layer X₁ ≤ x ≤ X₂, (3) downstream X₂ ≤ x. By taking X₁ and X₂ a few times L _s away from x = 0, we may use the analytic solution for a homogeneous medium in both the upstream and downstream. Taking into account the boundary conditions at x = −L_FEB and x = +∞ , we have

$$\begin{eqnarray}f(x,p)=\left\{\begin{array}{ll}{C}^{{u}}(p)\left[\exp \left(\frac{{U}_{s}x}{{\kappa }_{{xx}}^{u}}\right)-\exp \left(-\frac{{U}_{s}{L}_{\mathrm{FEB}}}{{\kappa }_{{xx}}^{u}}\right)\,\right] & -{L}_{\mathrm{FEB}}\leqslant x\leqslant {X}_{1}\\ {C}^{d}(p) & {X}_{2}\leqslant x\end{array}\right.\end{eqnarray} \tag{ 19 }$$

for the solutions in the upstream and downstream regions, respectively. Note that ${\kappa }_{{xx}}^{u}={\kappa }^{u}(p){\cos }^{2}{\theta }_{1}$ represents the xx component of the momentum-dependent diffusion coefficient in the upstream. Because of the constraint imposed by the boundary condition at x = +∞ , the problem is independent of the downstream diffusion coefficient. We will denote the diffusion coefficient inside the shock transition layer by κ without the superscript u, which may or may not be different from κ ^u . The momentum-dependent unknown constants C ^u (p), C ^d (p) must be determined by connecting these solutions to the one inside the shock transition layer.

In general, it is difficult to obtain a simple analytic solution in the region X₁ ≤ x ≤ X₂. Therefore, we expand the solution using orthogonal basis functions and use the pseudospectral method. We demand the continuity of the density f and the flux κ_xx ∂f/∂x as the boundary condition at x = X₁ and x = X₂ for connecting the internal and external solutions. Note that the first-order derivative ∂f/∂x may be discontinuous if the diffusion coefficient changes discontinuously across the boundary. Ultimately, the problem is reduced to solving a large system of linear equations. More detailed descriptions of numerical techniques may be found in Appendix B.

The following parameters will always be fixed for the numerical solutions discussed below: r = 4, B₁ = 1, U _s = 1, L _s = 1, X₂ = −X₁ = 4L _s . The orders of expansion for the orthogonal basis in the x- and y-directions will be denoted by N_x and N_y , respectively.

We first consider the simplest case in which the diffusion coefficient is constant κ = κ ^u = κ₀ in both space and momentum. Parameters used for the results discussed in this subsection are as follows: N_x = 32, N_y = 64, p₂/p₁ = 10², L_FEB/L _s = 10⁵. Note that the results do not depend on the free escaping boundary, as it is taken sufficiently far away from the shock.

Figure 2 shows the downstream momentum spectra for θ_Bn = 60° obtained with several different values of l_diff/L _s . Note that we always show the spectrum multiplied by p⁴ so that the standard DSA limit will be represented as a flat spectrum. We confirm that, as soon as the momentum becomes a few times larger than the injection momentum, a power-law spectrum develops. The spectrum is much steeper for l_diff/L _s ≲ 1 and approaches gradually to q = 4 with increasing l_diff/L _s , which is consistent with the theoretical prediction. We note that a similar conclusion was also obtained with an analytic approach by Drury et al. (1982).

Zoom In Zoom Out Reset image size

The dependence on l_diff/L _s is more clearly seen in Figure 3. The spectral indices evaluated directly from the numerical solutions by calculating numerical derivatives $\partial \mathrm{ln}f/\partial \mathrm{ln}p$ at p/p₁ = 5 are shown with the filled circles. We have also calculated the spectral indices by integrating Equation (6) using the numerical solutions, which are shown with the solid lines. It is seen that they both agree very well, indicating that Equation (6) is indeed useful to calculate a theoretical spectrum from known profiles of f(x) and V_x (x).

Zoom In Zoom Out Reset image size

As is clearly seen, the spectral index approaches the DSA prediction at l_diff/L _s ≫ 1 and steepens substantially for l_diff/L _s ≲ 1. While this trend itself is always true, the actual slope is dependent on θ_Bn as well. This is because the gradual kink in the magnetic field line changes the diffusive particle flux. Although the particle diffusion always occurs along the local field line, the decrease in θ from the downstream to the upstream side increases the particle flux toward the negative x (i.e., upstream) direction when the magnetic curvature is nonnegligible. This makes the effective diffusion length longer at small θ_Bn . Since $\theta \approx {\theta }_{{Bn}}\approx \mathrm{const}$ at quasi-perpendicular shocks, this effect is negligible and steeper spectra will result for larger θ_Bn .

We now consider the diffusion coefficient given by $\kappa (p)={\kappa }^{u}(p)={\kappa }_{1}\left(p/{p}_{1}\right)$, where κ₁ is constant. It is often used as a phenomenological diffusion model for CR transport, as it implies that ${\lambda }_{\mathrm{mfp}}/{r}_{g}=\mathrm{const}.$ for relativistic particles. We here employ the model simply to show qualitative characteristics expected for a diffusion coefficient monotonically increasing with momentum. We use N_x = 32, N_y = 128, p₂/p₁ = 10⁴, θ_Bn = 80°, and κ₁ = 10. This choice makes the diffusion length for the lowest momentum p = p₁ smaller than the shock thickness ${l}_{\mathrm{diff}}/{L}_{s}={\kappa }_{1}{\cos }^{2}{\theta }_{{Bn}}/({U}_{s}{L}_{s})\approx 0.30$. The diffusion length increases with momentum and eventually becomes much larger than the thickness l_diff/L _s ≫ 1. Therefore, we expect that SSDA should be dominant for low-energy particles and the acceleration regime gradually shifts to DSA at higher energies.

The top panel of Figure 4 shows the downstream spectra obtained with L_FEB/L _s = 10¹, 10², 10³, and 10⁴. The middle and bottom panels respectively show l_diff/L _s as a function of momentum and the local spectral slope for L_FEB/L _s = 10⁴ evaluated by direct numerical derivative. We find that the spectrum at low energy p/p₁ ≲ 10 is very steep. It is clear that the spectrum is not a pure power law and has a concave shape. This is qualitatively consistent with the increasing trend of l_diff/L _s as a function of momentum. We see from the middle panel that l_diff/L _s ≳ 10 is reached for p/p₁ ≳ 30, where the slope becomes close to the canonical DSA prediction as seen in the bottom panel. The high-energy part of the spectrum will be cut off if the particles diffusing ahead of the shock reach the free escaping boundary l_diff/L _s ∼ 1. Comparison between the top and middle panels confirms that the idea is consistent with the numerical solutions.

Zoom In Zoom Out Reset image size

We think that the spectrum shown in Figure 4 qualitatively represents the entire history of particle acceleration: injection by SSDA into DSA at low energy, further acceleration by DSA, and escape of the highest-energy particles from the system, resulting in the high-energy cutoff. With a more complicated but "realistic" momentum- and space-dependent diffusion model, we investigate the possibility of injection more quantitatively in the next subsection.

For the numerical solutions discussed so far, we have not specified either a physical scale length for the shock thickness or a timescale of particle pitch-angle scattering. For a more quantitative discussion of electron injection, we here introduce specific physical scales based on available knowledge of the collisionless shock dynamics.

It is physically more transparent to define the scattering rate normalized to the gyrofrequency Ω_ce/γ (including the relativistic effect) and convert it to the diffusion coefficient. We thus introduce the following functional form for the pitch-angle scattering rate:

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{\mu \mu }(p)}{{{\rm{\Omega }}}_{\mathrm{ce}}/\gamma }={\tilde{D}}_{\mu \mu ,0}{\left[1+{\left(\displaystyle \frac{p}{{p}_{* }}\right)}^{-2\nu }\right]}^{-1/2},\end{eqnarray} \tag{ 20 }$$

where ${\tilde{D}}_{\mu \mu ,0}$ is a dimensionless constant. For nonrelativistic energies and with ν > 0, it smoothly connects D_{μ μ} (p) ∝ p^ν for p ≪ p_* and ${D}_{\mu \mu }(p)\approx \mathrm{const}.$ for p ≫ p_*. For relativistic energies, it approaches D_{μ μ} (p) ∝ p⁻¹, which gives a constant scattering rate when normalized to the relativistic gyrofrequency Ω_ce/γ. It is reasonable to consider the scattering rate monotonically increasing with momentum because higher-energy particles can resonantly interact with lower-frequency and larger-amplitude waves. On the other hand, it is natural to assume that the normalized scattering rate saturates to a constant at high energies because it is likely to be limited by the Bohm scattering rate D_{μ μ} ∼ Ω_ce/γ. Equation (20) qualitatively represents such a momentum dependence.

Since it is not an easy task to obtain an accurate and physically grounded analytic model for the scattering rate within the shock transition layer, we assume a phenomenological model based on recent in situ observations (Amano et al. 2020). The particle intensity profiles measured at Earth's bow shock indicate that D_{μ μ} ∝ p² at energies below a few keV, while it saturates to a constant at higher energies. This motivates us to choose ν = 2, ${\tilde{D}}_{\mu \mu ,0}=3\times {10}^{-2}$, and p_*/(m_e c) = 1 × 10⁻¹ (which corresponds to ∼2.5 keV) as fiducial parameters. It is interesting to note that D_{μ μ} ∝ p² at lower energies indicates that the diffusion length is energy independent. Therefore, we expect that the spectrum in this energy range is described by a single power law, which is consistent with the measured particle spectrum. The saturation of D_{μ μ} at higher energies indicates that the diffusion length monotonically increases with energy, leading to the escape of higher-energy particles from the shock transition layer.

The particles with diffusion lengths longer than the shock thickness primarily interact with turbulence in the upstream and downstream regions. As we have discussed in Section 2.6, the fluctuation power available for particle scattering will be very small at short wavelength k ≳ k_*. Therefore, we assume that the scattering rate is finite only for those particles with momenta larger than the threshold $p\gtrsim {p}_{* }^{u}$ but is essentially negligible otherwise. Note that we denote quantities defined in the upstream by the superscript u. For simplicity, we use the same functional form given by Equation (20) with ν ^u = 6, ${\tilde{D}}_{\mu \mu }^{u}=3\times {10}^{-3}$, and ${p}_{* }^{u}/{({m}_{e}c)=({k}_{* }c/{\omega }_{\mathrm{pi}})}^{-1}({m}_{i}/{m}_{e})$, where the last relationship between ${p}_{* }^{u}$ and k_* comes from Equation (11).

In the following, unless otherwise noted, we use the following parameters representative of a moderately high Mach number and quasi-perpendicular Earth's bow shock: M_A = 10, θ_Bn = 85°, ${p}_{1}={m}_{e}{U}_{{\rm{s}}}/\cos {\theta }_{{Bn}}$, p₂/p₁ = 10³, V_A/c =2 × 10⁻⁴, η = 0.5. Note that we determine the lower bound of momentum by assuming that particle acceleration starts from the shock speed in HTF $v\sim {U}_{s}/\cos {\theta }_{{Bn}}$. Figure 5 shows the scattering rates as functions of energy for the shock transition layer (blue) and the upstream (red) as defined above. The threshold scattering rate defined by l_diff/L _s = 1 is also shown with the black dashed line. Note that we have assumed k_* c/ω_pi = 1. In this case, since particles with energies E/m_e c² ≲ 5 × 10⁻² (i.e., ≲25 keV) have the diffusion length smaller than the shock thickness, their major acceleration mechanism will be SSDA. In particular, a pure power-law spectrum is expected for E/m_e c² ≲ 5 × 10⁻³ (i.e., ≲5 keV) for which D_{μ μ} ∝ p² (${l}_{\mathrm{diff}}/{L}_{s}\approx \mathrm{const}.$) As the maximum energy expected for SSDA ${E}_{\max ,\ \mathrm{SSDA}}/{m}_{e}{c}^{2}\sim 5\times {10}^{-2}$ is comparable to the threshold energy associated with ${p}_{* }^{u}$ (which corresponds to E/m_e c² ∼ 10⁻¹), higher-energy particles escaping out from the shock transition layer will be injected into the particle acceleration cycle through DSA.

Zoom In Zoom Out Reset image size

Figure 6 confirms the electron injection scenario, in which numerical solutions for k_* c/ω_pi = 0.1, 0.32, and 1.0 are shown with L_FEB/L _s = 10⁵. The spectra in the top panel clearly show that the low-energy steeper power-law spectrum E/m_e c² ≲10⁻² is smoothly connected to the flat DSA spectrum at high energies for k_* c/ω_pi = 1.0. With decreasing k_* c/ω_pi, however, the injection becomes less efficient. The middle panel shows l_diff/L _s for the upstream and the shock transition layer as functions of energy. We see that l_diff/L_FEB ≪ 1 for k_* c/ω_pi = 1.0 at the energy where l_diff/L _s ∼ 1. Therefore, DSA is applicable at this energy. On the other hand, since l_diff/L_FEB ∼ 1 for k_* c/ω_pi = 0.32 at the same energy, a substantial fraction of particles will be lost from the system. This explains the partial cutoff before the spectrum becomes harder. With k_* c/ω_pi = 0.1, the cutoff is almost complete as expected from l_diff/L_FEB ≫ 1. In the bottom panel, local spectral slopes calculated by numerical derivatives are shown as functions of energy. In all three cases, the low-energy spectral indices are given by q ∼ 6. At sufficiently high energies, on the other hand, the spectral indices (for k_* c/ω_pi = 0.32, 1.0) are consistent with q = 4 as predicted by DSA. Note that the highest-energy cutoff is due to the free escaping boundary.

Zoom In Zoom Out Reset image size

We have seen that the upstream turbulence property as characterized by k_* is important to control the electron injection. Similarly, the capability of electron injection is affected strongly by M_A and θ_Bn as predicted by Equation (13). To see this, the numerical solutions obtained with θ_Bn = 75°, 80°, and 85° and a fixed value of L_FEB/L _s = 10³ are compared in Figure 7 in a similar format. Note that we have employed the scaling ${p}_{* }\propto {m}_{e}{U}_{{\rm{s}}}/\cos {\theta }_{{Bn}}$, which is the same as p₁ and p₂, assuming that the velocity scales with the shock speed in HTF ${U}_{s}/\cos {\theta }_{{Bn}}$. On the other hand, ${p}_{* }^{u}$ does not depend on the shock parameters, and a fixed value of k_* c/ω_pi = 1.0 is used for all three cases. The low-energy parts of spectra in three cases are very similar but shifted in energy because of the energy scaling $\propto 1/{\cos }^{2}{\theta }_{{Bn}}$. With increasing shock obliquity, the maximum energy of SSDA increases and eventually reaches the injection threshold. This results in the flat spectrum above E/m_e c² ∼ 3 × 10⁻² seen only for θ_Bn = 85°. As in the case of Figure 6, the characteristics of the spectra can be reasonably understood by looking at the energy dependence of l_diff/L _s shown in the middle panel. More specifically, l_diff/L_FEB ≪ 1 is satisfied only for θ_Bn = 85° at the theoretical maximum energy for SSDA determined by l_diff/L _s ∼ 1.

Zoom In Zoom Out Reset image size

The result shown here is consistent with the discussion presented in Section 2.7, in particular, Equation (13). We note that the model for the scattering rate adopted here has a relatively smooth transition below and above ${p}_{* }^{u}$. Also, the choice of ν ^u = 6 makes the condition for electron injection more or less sensitive to L_FEB. The sensitivity on L_FEB will be eliminated if one chooses a model that has a more abrupt transition of energy dependences across ${p}_{* }^{u}$. In any case, the efficiency of particle acceleration within the shock transition layer alone is not able to determine the injection condition. Depending on the efficiency and the energy dependence of the scattering rate outside the shock transition layer, the injection may or may not be achieved in a system of a given scale size L_FEB.

Although we have concentrated only on the downstream spectrum so far, it is also important to look at the spatial dependence of the spectrum for comparison with observations and kinetic simulations. Figure 8 shows the profiles of magnetic field B_z (x), as well as the distribution functions f(x) at specific energies for the solution obtained with θ_Bn = 85°. Note that f(x) is normalized by the downstream value for each energy. We see that the profiles at low energies are nearly identical to each other, which is a natural consequence of ${l}_{\mathrm{diff}}/{L}_{s}\approx \mathrm{const}$. The characteristic length scale becomes larger with increasing energy. The profile at the highest energy is nearly flat, meaning that l_diff/L _s ≫ 1. This is consistent with the fact that the particle acceleration at this energy is essentially the same as DSA at a discontinuous shock. Recall that the solution within the shock transition layer shown in Figure 8 is connected to the external (analytic) solutions with the boundary conditions. Therefore, the abrupt drop down in density at x/L _s = −4 seen at low energies is not a boundary effect but due to the physical escape beyond the free escaping boundary.

Zoom In Zoom Out Reset image size

Figure 9 shows the energy spectra and spectral slopes measured at different spatial positions x/L _s = −2, 0, and +2. At low energy E/m_e c² ≲ 3 × 10⁻³, the spectral slopes are nearly independent of x, while the density becomes substantially larger as it approaches the downstream. This is again due to the constant l_diff/L _s . The immediate upstream spectrum beyond this energy is harder than the downstream because of the monotonically increasing diffusion length as a function of energy. The spatial dependence of the spectrum provides yet another clue to test the theory against observations and kinetic simulations.

Zoom In Zoom Out Reset image size

In the proposed scenario of electron injection, the pitch-angle scattering over a wide range of energy plays a crucial role. The "realistic" model that we discussed in Section 3.4 adopts a phenomenological model based on in situ measurements of a single bow shock crossing event. The universality of the result is thus not guaranteed at present. In particular, since we know that the collisionless shock dynamics is largely controlled by parameters such as M_A and θ_Bn , it is quite likely that the wave intensity may also vary substantially depending on these parameters. We should consider the generation of plasma waves over a broad range of frequencies corresponding to the broad range of particle energy. Based on in situ measurements, the waves of particular interest may be subdivided into three characteristic frequencies: (1) high-frequency quasi-parallel whistlers 0.1 ≲ ω/Ω_ce ≲ 0.5 (e.g., Zhang 1999; Hull et al. 2012; Oka et al. 2017); (2) low-frequency oblique whistlers ω ∼ ω_LH, where ω_LH is the lower-hybrid frequency (Hull et al. 2012; Oka et al. 2019; Hull et al. 2020); (3) Alfvén ion cyclotron (AIC) waves or the rippling mode ω ≲ Ω_ci (e.g., Winske & Quest 1988; Johlander et al. 2016). More details on these wave characteristics will be discussed in Appendix C.

Among these, the high-frequency whistler waves are the key ingredient for triggering electron acceleration in the first place. In particular, because the wave intensity is much smaller at higher frequency in general, the condition l_diff/L _s ≲ 1 is more difficult to satisfy at low energy than high energy. Interestingly, the spatial profiles of nonthermal electrons indicate that ${l}_{\mathrm{diff}}/{L}_{{\rm{s}}}\approx \mathrm{const}.$ is satisfied for the energy range where a power-law spectrum is observed. This is quite consistent with the theoretical prediction. It might thus be possible to speculate that the generation of both high- and low-frequency whistler waves and the resulting electron confinement occur in a self-sustaining manner. In other words, if the weakly scattered accelerated electrons tend to escape upstream, the electron streaming becomes more intense and may destabilize waves that scatter themselves to reduce the diffusion length to be comparable to the shock thickness. Unless there exists such a self-sustaining mechanism, it is hard to believe that ${l}_{\mathrm{diff}}/{L}_{{\rm{s}}}\approx \mathrm{const}.$ is satisfied merely by chance.

The self-generation of waves by the back-streaming low-energy electrons that can scatter themselves has been discussed theoretically. Levinson (1992, 1996) proposed that upstream-directed oblique and low-frequency whistler-mode waves may be generated by low-energy electrons that obey the diffusion–convection equation. The instability condition for an oblique shock may be written as ${M}_{{\rm{A}}}^{\mathrm{HTF}}\gtrsim {\left({m}_{i}/{m}_{e}\right)}^{1/2}{\beta }_{e}^{-1/2}$, where β _e is the electron plasma beta. Note that the β _e dependence appears if the thermal electrons are responsible for the wave excitation. For the self-generation by electrons of momentum p, one may replace ${\beta }_{e}^{-1/2}$ by p_A,e /p. This wave generation process is very similar to the low-frequency whistler excitation mechanism driven by an electron heat flux (e.g., Krafft & Volokitin 2010; Roberg-Clark et al. 2018; Verscharen et al. 2019). The oblique whistler waves may be destabilized by fast electrons with a large p/p_A,e resonantly interacting with the wave via the anomalous cyclotron resonance. The consequence of the wave generation is to suppress the upstream-directed electron streaming (i.e., heat flux) through the wave emission in the same direction, which is exactly what is needed for SSDA. The excitation of oblique waves by the shock-reflected electrons (sometimes referred to as electron firehose instability) and the electron scattering associated with the waves observed in PIC simulations may be related to the self-generation mechanism discussed here (e.g., Matsukiyo et al. 2011; Guo et al. 2014a, 2014b).

On the other hand, Amano & Hoshino (2010) considered another wave excitation mechanism, in which quasi-parallel high-frequency whistler-mode waves are excited. The instability is driven by an upstream-directed beam of electrons with a loss cone that naturally arises if the beam is produced by magnetic mirror reflection at the shock (or by classical SDA). The electron cyclotron resonance around the loss cone destabilizes the waves propagating downstream. This mechanism may also provide means to induce scattering of upstream-directed electrons. The condition for the instability is roughly given by ${M}_{{\rm{A}}}^{\mathrm{HTF}}\gtrsim {({m}_{i}/{m}_{e})}^{1/2}{\beta }_{e}^{1/2}$ (Amano & Hoshino 2010). The β _e dependence here appears to overcome the cyclotron damping by thermal electrons. Noting that we have ignored numerical factors of order unity in the threshold Mach numbers, we think both of them give very similar conditions at β _e ∼ 1. Furthermore, the threshold Mach numbers are consistent with statistical analysis of Earth's bow shock observations (where β _e ∼ 1 on average), in which harder electron spectra were preferentially found when ${M}_{{\rm{A}}}^{\mathrm{HTF}}\gtrsim {({m}_{i}/{m}_{e})}^{1/2}$ (Oka et al. 2006). It is thus possible to conjecture that, once the Mach number goes beyond the threshold, quasi-parallel high-frequency whistlers are generated first to trigger SSDA at the lowest energy, and moderately accelerated electrons then excite oblique low-frequency whistlers to realize the continuous particle acceleration over a wide range of energy. We should mention that the self-generation mechanisms originally considered the wave generation in the upstream region to trigger DSA. In reality, however, the wave growth time may be sufficiently short to amplify the wave within the shock transition layer substantially. We thus expect that they should play a crucial role in the acceleration of low-energy electrons through SSDA.

There are certainly possibilities other than the self-generation mechanism discussed above. The high-frequency whistlers may also be generated from adiabatically heated electrons with T_⊥/T_∥ > 1 that will emit waves in both parallel and antiparallel directions. However, in situ observations indicate that quasi-parallel whistler waves sometimes exhibit unidirectional propagation. This may be due to the quasi-static cross-shock potential that accelerates the bulk of electrons toward the downstream along the local magnetic field line, which makes the electron distribution asymmetric. This will result in the emission of whistlers in the upstream direction (Tokar et al. 1984), which, therefore, cannot scatter the upstream-directed electrons. The low-frequency wave generation by modified two-stream instability (MTSI) favors ${M}_{{\rm{A}}}^{\mathrm{HTF}}\lesssim {({m}_{i}/{m}_{e})}^{1/2}$ (Matsukiyo & Scholer 2003, 2006), which is the opposite of the self-generation model. Although we almost always find intense wave activities around ω ∼ ω_LH in observations, the wave excitation mechanism might be different depending on the shock parameters. Another candidate is lower-hybrid drift instability (LHDI), which may be almost always unstable. Nevertheless, since LHDI is driven by the intense cross-field current at the shock ramp, it may not be a favorable candidate for the low-frequency waves often found in the foot region where the magnetic field gradient is relatively small.

Although our discussion in this paper focuses almost exclusively on the acceleration of suprathermal electrons, it is closely related to the heating of thermal electrons. It is known that the electron heating at supercritical collisionless shocks is often less efficient compared to ion heating. The reason for this has been considered that electrons behave adiabatically in the shock structure determined by the ion dynamics (Goodrich & Scudder 1984). In this scenario, the major contribution to electron heating is provided by the cross-shock electrostatic potential and associated microscopic plasma instabilities (Thomsen et al. 1987b; Scudder et al. 1986; Schwartz et al. 1988). More specifically, the field-aligned component of the electric field in HTF accelerates the bulk of electrons toward the downstream. The accelerated electrons may drive electrostatic instabilities, which results in the formation of the so-called flat-top velocity distribution in the downstream (Feldman et al. 1982, 1983). The distribution is nearly flat up to the energy corresponding to the shock potential in HTF, which is connected to a power-law spectrum at high energy. Indeed, such a flat distribution, also known as a self-similar distribution, is known to be generated as a result of ion-acoustic turbulence heating (Dum 1978a, 1978b). We think that the SSDA theory is responsible for the power-law part of the spectrum beyond the potential energy, which is typically of the order of ∼ 100 eV at Earth's bow shock. This indicates that the electron heating provides the seed population for SSDA. Modeling the electron heating is thus crucial for the overall normalization of the electron spectrum.

Strictly speaking, our model under the diffusion approximation is not able to deal with the low-energy thermal population with the plasma-frame velocity below the HTF shock speed. Such particles can diffuse in the plasma rest frame. However, as the flow velocity is faster than the individual particle velocities, they cannot propagate back upstream against the fast plasma flow. The diffusion approximation is clearly inadequate for them. Instead, we should adopt the pitch-angle diffusion Equation (A6) with the added effect of the parallel electric field associated with the cross-shock potential and turbulent components. This allows us to model the entire electron spectrum with a given source spectrum provided as the boundary condition in the upstream. Although the computational cost for solving the transport equation increases substantially, we do not need to assume any artificial injection into the system in this approach. This is certainly an important subject for future investigation.

It is important to note that recent high-resolution in situ spacecraft measurements suggested that the heating cannot be understood by the simple classical picture. In other words, reflected-ion-driven or current-driven waves can contribute significantly to the overall heating (e.g., Wilson et al. 2010, 2014; Goodrich et al. 2018). In fact, our injection model implies that even low-energy electrons should somehow behave nonadiabatically. Understanding the roles of large-amplitude waves often observed with in situ spacecraft is important not only for the heating itself but also for the nonthermal particle production efficiency. Indeed, the electron heating associated with the cross-shock potential alone may not always be sufficient. For the classical SDA to operate in initiating SSDA (and potentially subsequent DSA), the particle must have a sufficiently large pitch angle defined in HTF. The condition becomes progressively more difficult to satisfy as ${M}_{{\rm{A}}}^{\mathrm{HTF}}$ increases. Since the electron thermal velocity in the solar wind is much larger than the solar wind speed itself, it is not a problem at planetary bow shocks. On the other hand, additional electron heating processes may be needed to trigger SDA at very high Mach number young SNR shocks. It is known that the electrostatic Buneman instability driven by the reflected ion beam may be destabilized in this regime. It was shown that strong energization of electrons via the shock-surfing acceleration occurring when they first enter into the shock transition layer may be able to trigger the subsequent SDA (Amano & Hoshino 2007; Matsumoto et al. 2017). This suggests that the SSDA may or may not be initiated depending on the heating efficiency through microscopic plasma instabilities.

In constructing the simplified theoretical model, we have explicitly made use of HTF, assuming that the shock is stationary. However, it is well known that the collisionless shock exhibits intrinsic nonstationarity. The cyclic self-reformation (e.g., Leroy et al. 1982; Lembege & Savoini 1992), shock-surface rippling (e.g., Winske & Quest 1988), and any other types of nonstationarity violate the existence of HTF. It is also noted that a fast-mode MHD shock naturally produces a magnetic field component out of the coplanarity plane within its internal structure (Thomsen et al. 1987a). In such a case, it is impossible to find a frame in which the electric field vanishes throughout the structure even if it is stationary. Considering the highly dynamic and complicated nature of collisionless shocks, one may naturally question the robustness of the electron injection mechanism that assumes the stationarity and the existence of HTF.

We emphasize that the existence of HTF is not an essential ingredient in our model. Our motivation to adopt the specific frame of reference is to simplify the discussion as much as possible, but the application of the particle acceleration mechanism itself should not be limited to a stationary shock. It is known that the particle energization through the classical SDA may be equally understood as the magnetic gradient drift in the direction (antiparallel direction for electrons) of the motional electric field in NIF (Krauss-Varban & Wu 1989). The use of HTF substantially simplifies the analysis but is not essential. Since the energy gain mechanism of SSDA is the same as SDA, the existence of HTF is not a strict requirement. It does, however, require that the shock is at least locally subluminal. This is because we have assumed that the diffusive particle transport is limited only along the magnetic field line.

It is perhaps important to note that the accelerated electrons will be confined within the shock transition layer for a substantial period of time. It is known that, for the classical SDA, the typical timescale for the electron interaction with the shock structure is comparable to the ion gyroperiod (Krauss-Varban et al. 1989). In the SSDA mechanism, more efficient confinement via diffusion within the shock transition layer indicates that the interaction time will be much longer. Since the timescale of nonstationarity is typically of the order of the ion gyroperiod, energetic electrons will fully experience the nonstationary shock. It is possible to consider that a nonvanishing electric field associated with the nonstationarity contributes an additional energy gain, but at the same time the particle confinement efficiency might somewhat deteriorate. While the effect of nonstationarity has not yet been analyzed in detail, there have already been reports on the SSDA signatures with fully kinetic simulations where the shocks are clearly not stationary (Matsumoto et al. 2017; Kobzar et al. 2021; see Section 4.4 for details). This suggests that SSDA can survive even in the nonstationary shock structure.

We conclude that there is always the particle energization mechanism via the magnetic gradient drift, as the magnetic compression is, by definition, an indispensable property of a fast-mode shock. Nonstationary behaviors of the shock may quantitatively affect the particle acceleration efficiency, but not qualitatively. The applicability of the model is rather limited by the magnetic field orientation, which is, however, practically not an issue in nonrelativistic shocks.

Fully kinetic PIC simulation has been a powerful tool to investigate the dynamics of collisionless shocks. Recent simulations of nonrelativistic oblique shocks have already found signatures of SSDA in the accelerated particle trajectories (Matsumoto et al. 2017; Kobzar et al. 2021; Ha et al. 2021). More specifically, efficient pitch-angle scattering in momentum space during the energization of particles is qualitatively consistent with the theoretical prediction. On the other hand, the dependence of the spectral index on the shock parameters is not understood. In addition, the spatial dependence of the spectrum has not been carefully analyzed. The theoretical model presented in this paper provides the basis for more quantitative analyses of kinetic simulation results.

Nevertheless, we have to point out the artifacts of numerical simulations. First of all, PIC simulations typically adopt the shock speed one or two orders of magnitude higher than reality. Therefore, even for a modest shock obliquity, the shock may become very close to superluminal. Use of a relativistic HTF shock speed ${U}_{s}/\cos {\theta }_{{Bn}}\sim c$ violates the assumption of weak anisotropy made in the present theory. One has to check whether the pitch-angle anisotropy for a particular energy range of interest is not substantial before applying the theory. For the nearly isotropic energy range, Equation (8) may again be used for calculating the theoretical spectral index.

For more general cases, the effect of pitch-angle anisotropy must be taken into account explicitly by solving Equation (A6). For instance, we will be able to model energetic particles escaping upstream in a nearly scatter-free manner, which are often seen in simulations of oblique shocks. If the upstream scattering becomes stronger, the escaping particles will be scattered back into the shock to gain further energy in a stochastic manner. Note that such a signature has already been identified in published PIC simulation results and referred to as "multiple SDA" (e.g., Guo et al. 2014a, 2014b; Ha et al. 2021). In this process, electrons accelerated first at the shock by SDA return back to the shock to gain further energy by another SDA. The particle energy gain thus proceeds without going into the downstream region. As we have seen in this paper, if the pitch-angle anisotropy is sufficiently weak, this process has already been integrated into DSA at an oblique shock. Even though the accelerated particles do not travel deep into the downstream region, the multiple SDA is the signature that is consistent with the particle acceleration by DSA.

Although we have considered only steady-state solutions in this study, time-dependent solutions will also be useful. This allows us to compare the theoretical spectrum with the simulated energy spectrum that may not have reached a steady state. Theoretically, DSA at an oblique shock will accelerate particles with both SDA and the first-order Fermi process. Typical trajectory analysis, however, tends to pick up for those particles accelerated preferentially by multiple SDA because the first-order Fermi, which requires particles to diffuse into the downstream, is a relatively slow process at a quasi-perpendicular shock. We think that a time-dependent theoretical model will provide a useful guideline for the analysis of particle trajectories obtained by kinetic simulations.

Concerning plasma waves and instabilities in and around collisionless shocks, kinetic simulations have been used extensively over the decades. High-frequency quasi-parallel whistlers, low-frequency oblique whistlers, and AIC waves have all been identified by PIC simulations in the past (e.g., Matsukiyo & Matsumoto 2015; Tran & Sironi 2020; Kobzar et al. 2021). Since they have substantially different frequencies, scale separation by using a large ion-to-electron mass ratio m_i /m_e is important. It should be noted that a PIC simulation of a quasi-perpendicular shock with the realistic mass ratio by Matsukiyo & Matsumoto (2015) demonstrated the appearance of all three kinds of waves at the same time. It is interesting to investigate under what conditions those waves of different scales are excited and their roles on particle heating and acceleration.

The use of a large mass ratio is also desirable from the viewpoint of separation of energy scales. Assuming that SSDA starts from the momentum typically gained by the classical SDA ${p}_{1}\sim {m}_{e}{U}_{s}/\cos {\theta }_{{Bn}}$, we estimate the ratio between the injection threshold p_* and the initial momentum p₁ as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{p}_{* }}{{p}_{1}}\sim {\left(\displaystyle \frac{{k}_{* }c}{{\omega }_{\mathrm{pi}}}\right)}^{-1}\left(\displaystyle \frac{{m}_{i}}{{m}_{e}}\right){\left({M}_{{\rm{A}}}^{\mathrm{HTF}}\right)}^{-1}.\end{eqnarray} \tag{ 21 }$$

For instance, at Earth's bow shock with M_A = 8 and θ_Bn = 85°, we have p_*/p₁ ≈ 20. This indicates that the electron energy must be increased by a factor of ≈ 400 for efficient injection. On the other hand, the required energy gain becomes much smaller if a small mass ratio is used instead. For instance, with the same M_A and θ_Bn , a reduced mass ratio of m_i /m_e = 100 often adopted in PIC simulations will make it possible for electrons to interact with ion-scale fluctuations such as the rippling mode even by a single SDA. Therefore, we understand that a small mass ratio artificially makes the electron acceleration more efficient. For a realistic mass ratio, the condition p_*/p₁ ∼ 1 requires ${M}_{{\rm{A}}}^{\mathrm{HTF}}\gtrsim {10}^{3}$ or M_A ≳ 10² if $\cos {\theta }_{{Bn}}\sim 0.1$. We have recently shown that the dynamics of such high Mach number shocks will likely be dominated by the Weibel instability rather than the rippling (Nishigai & Amano 2021). In this case, whistler waves may no longer be relevant for particle acceleration because electrons can directly interact with the ion-scale turbulence generated by the Weibel instability. The electrons scattered by the large-amplitude Weibel-generated turbulence may be accelerated to ultrarelativistic energies within the shock transition layer (Matsumoto et al. 2017). This suggests that the SSDA applies to a wide range of shock parameters.

While the basis of the theory has already been confirmed by using recent in situ observations at Earth's bow shock (Amano et al. 2020), the findings in this paper will allow an even more detailed comparison with observations. In particular, the consistency between the scattering rate as a function of energy and the detailed features in the energy spectrum needs to be investigated using the observed profiles of magnetic field and flow velocity. Equation (6) will be particularly useful for comparison between the theoretical and observed spectral slopes for each energy channel of the particle instruments without solving the diffusion–convection equation.

One thing to keep in mind when analyzing the bow shock data is the finite curvature of the shock. Since the accelerated electrons may travel long distances along the magnetic field, the curvature of the shock has a substantial effect on the particle acceleration, especially at high energy. Krauss-Varban & Burgess (1991) estimated the typical energy beyond which the curvature effect becomes nonnegligible for the classical SDA by comparing a ballistic travel distance and the planer scale length. On the other hand, since we consider electrons suffering from pitch-angle scattering, we instead estimate the energy by comparing the diffusion length parallel to the field line ${l}_{\parallel }={l}_{\mathrm{diff}}/\cos {\theta }_{{Bn}}$ and the scale size. We assume that the variation of the shock obliquity δ as seen by electrons during the particle acceleration is negligible when ${l}_{\parallel }\lesssim {R}_{c}\sin \delta \approx {R}_{c}\delta$, where R _c is the curvature radius. This leads to the following condition:

$$\begin{eqnarray}\begin{array}{rcl}E & \lesssim & 31\,\mathrm{keV}\left(\frac{{D}_{\mu \mu }/{{\rm{\Omega }}}_{\mathrm{ce}}}{0.1}\right)\left(\frac{\delta }{\pi /180\,\mathrm{rad}}\right)\left(\frac{{U}_{s}}{400\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)\\ & & \times \left(\frac{B}{10\,\mathrm{nT}}\right)\left(\frac{{R}_{c}}{20{R}_{{\rm{E}}}}\right){\left(\frac{\cos {\theta }_{{Bn}}}{\cos 85^\circ }\right)}^{-1},\end{array}\end{eqnarray} \tag{ 22 }$$

where R_E is the Earth radius. This indicates that, for a typical curvature radius at the bow shock R _c ∼ 20R_E, the scale size of the shock will have a nonnegligible effect for electron acceleration beyond a few tens of keV. Note that the diffusive confinement of the particles makes the upper limit much larger than the estimate by Krauss-Varban & Burgess (1991). The finite curvature effect may roughly correspond to the presence of the free escaping boundary at a distance of ${L}_{\mathrm{FEB}}\,\sim {R}_{c}\delta \cos {\theta }_{{Bn}}$ away from the shock. Therefore, under the normal circumstance, Earth's bow shock will not be the place where the electron injection into DSA occurs efficiently. In any case, the high-energy part of the spectrum must be analyzed carefully by taking into account the finite curvature effect.

Investigating the spectrum at high energy around and beyond the transition between SSDA and DSA may be possible at shocks of larger scale sizes such as interplanetary shocks and the bow shocks at Saturn and Jupiter. While interplanetary shocks are favorable in scale size, they are typically much weaker in terms of Mach number. On the other hand, the Saturnian bow shock appears to be a very high Mach number and has a relatively large scale. In one of the shock crossing events observed by the Cassini spacecraft, Masters et al. (2013, 2017) found that the highest energy range (≳100 keV) of the spectrum was harder than the lower energy. We may be able to interpret the hardening of the energy spectrum with the present theory. Although such events are very rare, careful reanalysis might be able to check the consistency with the theory.

Detailed analysis of plasma waves is also desirable. The most important parameter is obviously the intensity of the waves. The condition l_diff/L _s ≲ 1 required for efficient particle acceleration by SSDA indicates that the wave intensity must be larger than a frequency-dependent threshold (Katou & Amano 2019; Amano et al. 2020). The self-generation hypothesis mentioned above predicts that the absolute wave intensity will increase as a function of ${M}_{{\rm{A}}}^{\mathrm{HTF}}$. In particular, an abrupt transition across ${M}_{{\rm{A}}}^{\mathrm{HTF}}\sim {({m}_{i}/{m}_{e})}^{1/2}$ might be realized. On the other hand, even if the wave intensity is kept more or less constant by some other wave excitation mechanisms, the relative wave intensity normalized to the theoretical threshold will increase as a function of ${M}_{{\rm{A}}}^{\mathrm{HTF}}$ because the threshold is proportional to ${({M}_{{\rm{A}}}^{\mathrm{HTF}})}^{-2}$ if quasi-linear theory applies. Both of the above scenarios will explain the observed correlation between the electron spectral index and ${M}_{{\rm{A}}}^{\mathrm{HTF}}$ found by Oka et al. (2006). Nevertheless, differentiating the two will provide crucial information for modeling the electron injection.

In addition to the intensity, the wave propagation property is also important. Since quasi-parallel whistlers can scatter low-energy electrons that are moving opposite to the wave propagation direction, their propagation directions must be bidirectional for efficient electron confinement. We have indeed confirmed the bidirectional characteristics of quasi-parallel whistlers in some of the bow shock crossings, but not always. In contrast, the oblique low-frequency whistlers may scatter electrons moving in both directions because of the anomalous cyclotron resonance in addition to the normal cyclotron resonance. If the waves are driven by the upstream-directed electrons, they should always propagate upstream. Since MTSI may generate both upstream- and downstream-propagating waves, one may distinguish the wave excitation mechanisms using the wave characteristics. While one such attempt has been reported recently (Hull et al. 2020), the dependence on the shock parameters has not been known yet. It is interesting to investigate the possibility that the wave property changes systematically depending on the shock parameters as predicted by the self-generation hypothesis.

Laboratory experiment using high-power lasers has now become an important tool to study collisionless shocks physics (e.g., Huntington et al. 2015; Schaeffer et al. 2017, 2019). Recent development in experimental technologies enables one to perform experiments with various pre-shock plasma magnetizations. In particular, Fiuza et al. (2020) reported nonthermal electron acceleration up to ∼500 keV in Weibel-mediated collisionless shocks reproduced by laser experiments conducted at the National Ignition Facility. A supporting PIC simulation representing the experiment demonstrated a power-law tail in the electron energy spectrum. While they showed that the particle energy gain is consistent with the first-order Fermi mechanism, the spectral index ${dN}/{dE}\propto {E}^{-3}$ or f(p) ∝ p⁻⁷ is clearly steeper than the standard DSA prediction but is more consistent with SSDA. The accelerated particle trajectories indicate that they are confined within the shock transition layer during the energy gain, which implies that l_diff/L _s ≲ 1 is satisfied. Note that it may not be appropriate to refer to the particle acceleration mechanism as SSDA for this application because the shock seems to be essentially unmagnetized, in which case SDA arising from the ambient magnetic field compression is absent. Nevertheless, magnetic field compression is not a necessary ingredient for the present particle acceleration mechanism. As long as the scattering is strong enough to confine the particles within the acceleration region, the particle acceleration by the first-order Fermi mechanism in a smooth shock profile may produce a power-law-like spectrum that is steeper than DSA.

Theoretically, we expect that more strongly magnetized conditions may result in more efficient electron acceleration because of the additional energy gain by SDA. On the other hand, if the magnetization is too strong, the Weibel instability will be suppressed, and the shock will behave essentially like the magnetized Earth's bow shock (Nishigai & Amano 2021). It would be interesting to investigate the transition from the classical magnetized shock to a more violent Weibel-dominated shock and how this affects the overall electron acceleration efficiency.

The injection model, preferably calibrated with in situ observations and kinetic simulations, may apply to predict the amount of CR electrons in astrophysical nonthermal emission sources. However, the biggest obstacle for the astrophysical application is the resolution of observations. The ideal assumption that the upstream medium is homogeneous will not be satisfied in practice. It is rather more realistic to assume that observations reflect a superposition of different upstream conditions. In other words, the efficiency of injection will be determined by averaging over the variation of local shock parameters introduced by upstream turbulence interacting with the shock. We note that the particle acceleration mechanism discussed in this paper is dependent only on the local parameters as long as the shock is locally planar. Therefore, if the statistics of the local shock parameters are given, we will be able to predict the absolute flux of CR electrons as a function of ${M}_{{\rm{A}}}^{\mathrm{HTF}}$ using the present injection model.

Observed radio and X-ray synchrotron spectra may be compared directly with the theoretical prediction. While synchrotron photons are emitted from highly relativistic electrons (GeV for radio and TeV for X-ray), we can consistently construct a theoretical spectrum from subrelativistic to relativistic energies to calculate a theoretical synchrotron spectrum. Another potential probe is nonthermal bremsstrahlung emission in the hard X-ray band from nonrelativistic suprathermal electrons. Since the hard X-rays are directly emitted typically from the electrons in the energy range where SSDA plays the dominant role, such observations will provide an important clue for testing the electron injection model. For instance, recent NuSTAR observation of W49B indicates a power-law-like photon spectrum that was interpreted as nonthermal bremsstrahlung emission (Tanaka et al. 2018). Future X-ray observations of SNRs may provide more detailed information on the electron spectrum in the nonrelativistic energy range, which will allow us to perform a direct comparison with the theoretical model.

Let us start with the relativistic Vlasov equation describing the evolution of phase-space density f( x , p , t):

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+\displaystyle \frac{{\boldsymbol{p}}}{\gamma m}\cdot {\rm{\nabla }}f+q\left({\boldsymbol{E}}+\displaystyle \frac{{\boldsymbol{p}}}{\gamma {mc}}\times {\boldsymbol{B}}\right)\cdot \displaystyle \frac{\partial f}{\partial {\boldsymbol{p}}}=0,\end{eqnarray} \tag{ A1 }$$

for a particle species with charge q and rest mass m. The momentum is defined by p = γ m v , with the particle three-velocity v and Lorentz factor $\gamma =1/\sqrt{1-{v}^{2}/{c}^{2}}$. The pitch-angle diffusion equation, also known as the focused transport equation, may be obtained by taking the gyrophase average of the Vlasov equation and later introducing a phenomenological pitch-angle diffusion term on the right-hand side, assuming that the gyroradius is small compared to the scale size of the inohomogeneity of the electromagnetic field (e.g., Skilling 1975; Isenberg 1997). We will follow the same procedure in HTF. In other words, we assume that the plasma flow V is everywhere parallel to the local magnetic field: V = V b , where b = B /B denotes the magnetic field unit vector. In this case, the convection electric field E = − V × B /c is always zero. Ignoring also the field-aligned electric field E · b = 0, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+\displaystyle \frac{p\mu }{\gamma m}{\boldsymbol{b}}\cdot {\rm{\nabla }}f+\displaystyle \frac{1-{\mu }^{2}}{2}\displaystyle \frac{p}{\gamma m}{\rm{\nabla }}\cdot {\boldsymbol{b}}\displaystyle \frac{\partial f}{\partial \mu }=0,\end{eqnarray} \tag{ A2 }$$

where p = ∣ p ∣ and μ denotes pitch-angle cosine. The absence of a term proportional to ∂f/∂p clearly indicates that the particle energy is a conserved quantity. Rewriting ${\rm{\nabla }}\cdot {\boldsymbol{b}}\,=-({\boldsymbol{b}}\cdot {\rm{\nabla }})\mathrm{ln}B$, we see that the transport described by Equation (A2) is nothing more than the adiabatic magnetic mirroring in an inhomogeneous magnetic field.

In Equation (A2), both the configuration and velocity space coordinates are represented in HTF. However, it is customary and more convenient to use the momentum coordinates defined in the local plasma rest frame, in which momentum space diffusion via wave−particle interaction can be written in the simplest form. By using the Lorentz transform from HTF (p, μ) to the local plasma rest-frame coordinates $(p^{\prime} ,\mu ^{\prime} )$,

$$\begin{eqnarray}&&\gamma c={\rm{\Gamma }}\left(\gamma ^{\prime} c+p^{\prime} \mu ^{\prime} V\right),\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}&&p\mu ={\rm{\Gamma }}\left(p^{\prime} \mu ^{\prime} +\gamma ^{\prime} V\right),\end{eqnarray} \tag{ A4 }$$

where ${\rm{\Gamma }}=1/\sqrt{1-{V}^{2}/{c}^{2}}$, Equation (A2) may be rewritten as follows:

$$\begin{eqnarray}&&\begin{array}{l}\left(1+\displaystyle \frac{v^{\prime} \mu ^{\prime} V}{{c}^{2}}\right)\displaystyle \frac{\partial f}{\partial t}+\left(V+v^{\prime} \mu ^{\prime} \right)\displaystyle \frac{\partial f}{\partial s}\\ \quad -\,\,\displaystyle \frac{1-\mu {{\prime} }^{2}}{2}\left[(v^{\prime} +V\mu ^{\prime} )\displaystyle \frac{\partial \mathrm{ln}B}{\partial s}\right.\\ \quad \left.+\,2{{\rm{\Gamma }}}^{2}\left(\mu ^{\prime} +\displaystyle \frac{V}{v^{\prime} }\right)\displaystyle \frac{\partial V}{\partial s}\right]\displaystyle \frac{\partial f}{\partial \mu ^{\prime} }\\ \quad +\,\left[\displaystyle \frac{1-\mu {{\prime} }^{2}}{2}V\displaystyle \frac{\partial \mathrm{ln}B}{\partial s}-{{\rm{\Gamma }}}^{2}\right.\\ \quad \times \,\left.\left(\mu ^{\prime} +\displaystyle \frac{V}{v^{\prime} }\right)\mu ^{\prime} \displaystyle \frac{\partial V}{\partial s}\right]\displaystyle \frac{\partial f}{\partial \mathrm{ln}p^{\prime} }\\ \quad =\,\,{\rm{\Gamma }}\displaystyle \frac{\partial }{\partial \mu ^{\prime} }\left[(1-\mu {{\prime} }^{2}){D}_{\mu \mu }\displaystyle \frac{\partial }{\partial \mu ^{\prime} }f\right].\end{array}\end{eqnarray} \tag{ A5 }$$

Note that ∂/∂s = b · ∇ represents the spatial derivative along the local magnetic field line. On the right-hand side, we have introduced a pitch-angle scattering of particles. The scattering conserves energy in the plasma rest frame, and the coefficient D_{μ μ} is defined as the rate of scattering measured in that frame. The Lorentz factor Γ on the right-hand side indicates the relativistic time dilation effect seen in the laboratory frame.

Henceforth, we will always use the momentum coordinates defined in the local plasma rest frame and drop primes unless otherwise noted. We define x as the coordinate parallel to the normal vector of a plane shock surface, with the positive direction pointing toward the downstream. The angle between the magnetic field and the x-axis is denoted by θ = θ(x) (thus θ(x) = θ_Bn at x → −∞ ). We may now rewrite Equation (A5) in the following form:

$$\begin{eqnarray}&&\begin{array}{l}\left(1+\displaystyle \frac{v\mu V}{{c}^{2}}\right)\displaystyle \frac{\partial f}{\partial t}+\left(V+v\mu \right)\cos \theta \displaystyle \frac{\partial f}{\partial x}\\ \quad -\,\displaystyle \frac{1-{\mu }^{2}}{2}\left[(v+V\mu )\displaystyle \frac{\partial \mathrm{ln}B}{\partial x}\right.\\ \quad \left.+\,2{{\rm{\Gamma }}}^{2}\left(\mu +\displaystyle \frac{V}{v}\right)\displaystyle \frac{\partial V}{\partial x}\right]\cos \theta \displaystyle \frac{\partial f}{\partial \mu }\\ \quad +\,\left[\displaystyle \frac{1-{\mu }^{2}}{2}V\displaystyle \frac{\partial \mathrm{ln}B}{\partial x}-{{\rm{\Gamma }}}^{2}\right.\\ \quad \left.\times \,\left(\mu +\displaystyle \frac{V}{v}\right)\mu \displaystyle \frac{\partial V}{\partial x}\right]\cos \theta \displaystyle \frac{\partial f}{\partial \mathrm{ln}p}\\ \quad =\,{\rm{\Gamma }}\displaystyle \frac{\partial }{\partial \mu }\left[(1-{\mu }^{2}){D}_{\mu \mu }\displaystyle \frac{\partial }{\partial \mu }f\right],\end{array}\end{eqnarray} \tag{ A6 }$$

where we have used $\partial /\partial s=\cos \theta \,\partial /\partial x$. This equation may be used to investigate the particle acceleration at plane oblique shocks for arbitrary shock speeds and particle energies, as long as the shock is subluminal. We note that the transport equation used by Katou & Amano (2019) may be obtained if both the particle velocity and the shock speed are nonrelativistic (v ≪ c and V ≪ c).

If the scattering is so strong that the pitch-angle distribution is always kept nearly isotropic, one can further reduce the transport equation to a simpler form. The weak-anisotropy assumption requires that individual particle velocities are much larger than the flow speed V²/v² ≪ 1, because otherwise substantial anisotropy will develop in the upstream region. Therefore, a necessary (but not sufficient) condition for the weak-anisotropy approximation is that the shock speed in HTF is nonrelativistic ${U}_{s}/\cos {\theta }_{{Bn}}\ll c$. While the anisotropy of electrons at suprathermal energies ( ≳200 eV) was shown to be weak by in situ observations at Earth's bow shock (Amano et al. 2020), its validation, in general, requires the knowledge of D_{μ μ} even for nonrelativistic shock speeds. We here simply adopt this approximation, but the application should be limited to beyond an unspecified threshold energy where the anisotropy becomes weak. The general case with finite pitch-angle anisotropy will be left for future investigation.

Under the assumption of weak anisotropy, we may expand the pitch-angle distribution in terms of the Legendre polynomials:

$$\begin{eqnarray}&&f(x,p,\mu )=\sum _{n=0}^{N}{f}_{n}(x,p){P}_{n}(\mu ),\end{eqnarray} \tag{ A7 }$$

where P_n (μ) is the nth-order Legendre polynomial with argument μ. We truncate the expansion to the first order, assuming ∣f₀∣ ≫ ∣f₁∣. This leads to the following diffusion–convection equation for the isotropic part of the distribution function f₀(x, p) (see, e.g., Zank 2014),

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {f}_{0}}{\partial t}+V\cos \theta \,\displaystyle \frac{\partial {f}_{0}}{\partial x}+\displaystyle \frac{1}{3}\left(\displaystyle \frac{\partial \mathrm{ln}B}{\partial x}-\displaystyle \frac{\partial \mathrm{ln}V}{\partial x}\right)V\\ \quad \times \cos \theta \displaystyle \frac{\partial {f}_{0}}{\partial \mathrm{ln}p}=\displaystyle \frac{\partial }{\partial x}\left(\kappa {\cos }^{2}\theta \,\displaystyle \frac{\partial {f}_{0}}{\partial x}\right).\end{array}\end{eqnarray} \tag{ A8 }$$

The diffusion coefficient along the ambient magnetic field is given by κ = v²/6D_{μ μ} , where we have assumed that D_{μ μ} is independent of μ. It is instructive to note

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial x}\left(V\cos \theta \right)=-V\cos \theta \left(\displaystyle \frac{\partial \mathrm{ln}B}{\partial x}-\displaystyle \frac{\partial \mathrm{ln}V}{\partial x}\right),\end{eqnarray} \tag{ A9 }$$

which indicates

$$\begin{eqnarray}&&\displaystyle \frac{1}{p}\displaystyle \frac{{dp}}{{dt}}=-\displaystyle \frac{1}{3}{\rm{\nabla }}\cdot {\boldsymbol{V}}.\end{eqnarray} \tag{ A10 }$$

In other words, the particle momentum gain is proportional to the divergence of the plasma flow. We thus find that Equation (A8) is equivalent to the standard diffusion–convection equation written in HTF.

Katou & Amano (2019) showed that the particle energy gain at a nearly perpendicular shock with $\cos {\theta }_{{Bn}}\ll 1$ is dominated by the gradient-B term. This may easily be understood by recognizing that the compression of the flow is only in the direction normal to the shock surface, which is only a small fraction of the HTF shock speed. Therefore, we understand that $| \partial \mathrm{ln}B/\partial x| \gg | \partial \mathrm{ln}V/\partial x|$ is satisfied at a quasi-perpendicular shock. Physically, the energy gain represented by $\partial \mathrm{ln}B/\partial x$ is due to SDA, or the gradient-B drift in the direction of the convection electric field. On the other hand, $\partial \mathrm{ln}V/\partial x$ represents the first-order Fermi mechanism by converging scattering centers.