Electron Heating in Low-Mach-number Perpendicular Shocks. I. Heating Mechanism

The work done on an isotropic gas with pressure P and volume V is simply ${dW}=-{PdV}$. We shall generalize this expression to the case of an anisotropic gas having pressure perpendicular (parallel) to the magnetic field lines equal to P_⊥ (${P}_{\parallel }$, respectively). Consider a magnetic flux tube with length L, cross-sectional area A, volume V = LA, and field strength B. The magnetic flux through the area A is ${\rm{\Phi }}={BA}$. In response to a compression (or expansion) perpendicular to the magnetic field, the volume will change as

$$\begin{eqnarray}&&{{dV}}_{\perp }=L\,{dA}=L\,d\left(\displaystyle \frac{{\rm{\Phi }}}{B}\right)=-L{\rm{\Phi }}\displaystyle \frac{{dB}}{{B}^{2}}=-{Vd}\mathrm{ln}B,\end{eqnarray} \tag{ 1 }$$

where we have used the fact that, because of flux freezing, Φ is a constant. In contrast, for compression (or expansion) along the field, the volume will change as

$$\begin{eqnarray}&&{{dV}}_{\parallel }=A\,{dL}=A\,d\left(\displaystyle \frac{V}{A}\right)=\displaystyle \frac{{AN}}{{\rm{\Phi }}}\,d\left(\displaystyle \frac{B}{n}\right)=-{Vd}\mathrm{ln}\left(\displaystyle \frac{n}{B}\right),\end{eqnarray} \tag{ 2 }$$

where N is the total number of particles in the volume element, with number density n = N/V. It follows that the work done on an anisotropic gas can be written as

$$\begin{eqnarray}{dW} & = & -{P}_{\perp }{{dV}}_{\perp }-{P}_{\parallel }{{dV}}_{\parallel }\\ & = & {P}_{\perp }{Vd}\mathrm{ln}B+{P}_{\parallel }{Vd}\mathrm{ln}\left(\displaystyle \frac{n}{B}\right).\end{eqnarray} \tag{ 3 }$$

Defining the work done per particle as dw = dW/N, we find that it can be separated into a "perpendicular" component ${{dw}}_{\perp }$ and a "parallel" component ${{dw}}_{\parallel }$ as

$$\begin{eqnarray}&&{dw}=\mathop{\underbrace{{k}_{{\rm{B}}}{T}_{\perp }d\mathrm{ln}B}}\limits_{{{dw}}_{\perp }}+\mathop{\underbrace{{k}_{{\rm{B}}}{T}_{\parallel }d\mathrm{ln}\left(\displaystyle \frac{n}{B}\right)}}\limits_{{{dw}}_{\parallel }}.\end{eqnarray} \tag{ 4 }$$

It follows that ${{dw}}_{\perp }$ will change the internal energy per particle ${u}_{\perp }$ associated with motions perpendicular to the field, while ${{dw}}_{\parallel }$ will affect the energy per particle ${u}_{\parallel }$ associated with parallel motions.

In writing the energy equation for the perpendicular and parallel components, we need to take into account two additional processes: (1) In the presence of pitch-angle scattering, heat can be transferred between the two components (as we show below, this will give rise to an entropy increase). We denote the differential amount of transferred heat as ${{dq}}_{\perp \to \parallel }$, with the convention that ${{dq}}_{\perp \to \parallel }\gt 0$ if heat flows from the perpendicular to the parallel component. (2) Pitch-angle scattering may be caused by self-generated waves (e.g., sourced by the plasma anisotropy), whose energy needs to be provided by the plasma itself. The energy balance relations then read

$$\begin{eqnarray}&&{{du}}_{\perp }={{dw}}_{\perp }-{{dq}}_{\perp \to \parallel }-{{de}}_{w,\perp },\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{du}}_{\parallel }={{dw}}_{\parallel }+{{dq}}_{\perp \to \parallel }-{{de}}_{w,\parallel }\,,\end{eqnarray} \tag{ 6 }$$

where we denote the wave energy per particle coming from the perpendicular (parallel) plasma energy as ${{de}}_{w,\perp }$ (${{de}}_{w,\parallel }$, respectively). By summing the above two equations, we obtain the expected result that the net change of internal energy per particle is equal to the external work minus the energy given to waves:

$$\begin{eqnarray}&&{du}\equiv {{du}}_{\perp }+{{du}}_{\parallel }={dw}-{{de}}_{w\mathrm{tot}},\end{eqnarray} \tag{ 7 }$$

where we denote the total energy per particle transferred to waves as ${{de}}_{w\mathrm{tot}}\equiv {{de}}_{w,\perp }+{{de}}_{w,\parallel }$, including magnetic, electric, and bulk kinetic contributions (in practice, the magnetic term always dominates).

While the total wave energy per particle ${{de}}_{w\mathrm{tot}}$ is easy to extract from our simulations, the two contributions ${{de}}_{w,\perp }$ and ${{de}}_{w,\parallel }$ are hard to separate. We show below that for the entropy calculation it suffices to measure the total energy per particle transferred to waves ${{de}}_{w\mathrm{tot}}$. We also remark that ${{de}}_{w\mathrm{tot}}$ accounts for the differential energy per particle transferred to waves, which might not necessarily equal the differential change in the energy residing in waves, which we shall call ${{de}}_{w}$. More specifically, while for electron-driven waves ${{de}}_{w\mathrm{tot}}={{de}}_{w}$, we will show in Section 6 that proton-generated waves will lose energy by performing work on the electron plasma, so the change in the energy residing in proton waves ${{de}}_{w}$ will be smaller than the differential energy ${{de}}_{w\mathrm{tot}}$ transferred from protons to waves.

For a nonrelativistic bi-Maxwellian plasma with perpendicular temperature ${T}_{\perp }$ and parallel temperature ${T}_{\parallel }$, the entropy per particle (or specific entropy) is

$$\begin{eqnarray}&&s\equiv -\displaystyle \frac{\int {d}^{3}p\,f\mathrm{ln}f}{\int {d}^{3}p\,f}=\mathrm{ln}\left(\displaystyle \frac{{T}_{\perp }{T}_{\parallel }^{1/2}}{n}\right)+C,\end{eqnarray} \tag{ 8 }$$

where $f({\boldsymbol{p}})$ is the phase space distribution and C is a normalization constant. By differentiating,

$$\begin{eqnarray}&&{ds}=\displaystyle \frac{{{dT}}_{\perp }}{{T}_{\perp }}+\displaystyle \frac{1}{2}\displaystyle \frac{{{dT}}_{\parallel }}{{T}_{\parallel }}-\displaystyle \frac{{dn}}{n}.\end{eqnarray} \tag{ 9 }$$

The temperature ${T}_{\perp ,\parallel }$ can be related to the internal energy per particle ${u}_{\perp ,\parallel }$ via the respective adiabatic index ${{\rm{\Gamma }}}_{\perp ,\parallel }$ as

$$\begin{eqnarray}&&{u}_{\perp ,\parallel }=\displaystyle \frac{{k}_{{\rm{B}}}{T}_{\perp ,\parallel }}{{{\rm{\Gamma }}}_{\perp ,\parallel }-1}.\end{eqnarray} \tag{ 10 }$$

For a nonrelativistic gas, ${{\rm{\Gamma }}}_{\perp }=2$ (two degrees of freedom are available in the perpendicular direction), whereas ${{\rm{\Gamma }}}_{\parallel }=3$ (one degree of freedom). The equation above then becomes

$$\begin{eqnarray}&&{ds}=\displaystyle \frac{{{du}}_{\perp }}{{T}_{\perp }}+\displaystyle \frac{{{du}}_{\parallel }}{{T}_{\parallel }}-\displaystyle \frac{{dn}}{n}.\end{eqnarray} \tag{ 11 }$$

Using Equations (5) and (6), we have

$$\begin{eqnarray}&&{ds}={{dq}}_{\perp \to \parallel }\left[\displaystyle \frac{1}{{T}_{\parallel }}-\displaystyle \frac{1}{{T}_{\perp }}\right]-\displaystyle \frac{{{de}}_{w\perp }}{{T}_{\perp }}-\displaystyle \frac{{{de}}_{w\parallel }}{{T}_{\parallel }}\,,\end{eqnarray} \tag{ 12 }$$

which shows that the entropy of the gas can change as a result of heat flowing internally between the parallel and perpendicular components (first term on the right-hand side) or when generating the waves (second term). This can be rewritten in two equivalent forms:

$$\begin{eqnarray}&&{ds}=\left[\displaystyle \frac{1}{2}d\mathrm{ln}\left(\displaystyle \frac{{T}_{\parallel }}{{(n/B)}^{2}}\right)\right]\cdot \left[1-\displaystyle \frac{{T}_{\parallel }}{{T}_{\perp }}\right]-\displaystyle \frac{{{de}}_{w\mathrm{tot}}}{{T}_{\perp }},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{ds}=-\left[d\mathrm{ln}\left(\displaystyle \frac{{T}_{\perp }}{B}\right)\right]\cdot \left[\displaystyle \frac{{T}_{\perp }}{{T}_{\parallel }}-1\right]-\displaystyle \frac{{{de}}_{w\mathrm{tot}}}{{T}_{\parallel }}.\end{eqnarray} \tag{ 14 }$$

As anticipated above, the two separate components ${{de}}_{w,\parallel }$ and ${{de}}_{w,\perp }$ of the wave energy per particle do not explicitly enter the entropy equation.

In Equations (13) and (14), the first term on the right-hand side typically dominates. This clearly demonstrates that two ingredients are required for entropy generation: (1) the presence of a temperature anisotropy and (2) a mechanism to break the adiabatic invariance. Note that the Chew–Goldberger–Low (CGL) double adiabatic theory of Chew et al. (1956) predicts that, for adiabatic perturbations, ${T}_{\perp }\propto B$ and ${T}_{\parallel }\propto {(n/B)}^{2}$, which follow from the conservation of the first and second adiabatic invariants. The form of Equations (13) and (14) is thus easy to understand. In most cases, it is the temperature anisotropy that provides the free energy for generating the waves responsible for breaking the adiabatic invariance.

We conclude with an important remark on the nature of the magnetic field B. This should be meant as a large-scale field, so the particle response to its variation is properly modeled by the CGL approximation. In particular, the field that we have denoted as B must not include the magnetic contribution of the waves that break the particle adiabatic invariance. In practice, B will take into account all of the magnetic contributions at scales much larger than the particle Larmor radius (for the species in question) and at frequencies much lower than the relevant gyration frequency. It follows that proton-generated waves that break the proton adiabatic invariance can still serve as large-scale B fields for the electron energy and entropy equations, as we further discuss in Section 6.

In this subsection, we describe the proton dynamics, with a focus on proton isotropization and thermalization downstream of the shock. Figure 1 shows the profile of various quantities in the shock at time $t=25.6\,{{\rm{\Omega }}}_{{ci}}^{-1\ }$, as a function of the x coordinate relative to the shock location ${x}_{\mathrm{sh}}$, in units of the proton Larmor radius $\,{r}_{\mathrm{Li}}$ defined in Equation (21).

Zoom In Zoom Out Reset image size

Panel (a) shows the y-averaged profile of the proton number density n_i in units of the proton density in the upstream n_i0 (black line). The density compression at the shock reaches ${n}_{i}/{n}_{i0}\sim 3.5$ over a distance of $\sim {r}_{\mathrm{Li}}$, consistent with the expectation that the thickness of a perpendicular shock should be of the order of the proton Larmor radius (Bale et al. 2003; Scholer & Burgess 2006). The density oscillates on a typical length scale of $\sim {r}_{\mathrm{Li}}$ after the overshoot and then relaxes to the Rankine–Hugoniot value of ∼2.8 beyond a distance of $\sim 5\,{r}_{\mathrm{Li}}$ behind the shock.

The density pileup at the shock is related to the electrostatic potential Φ that develops in the shock transition region. This phenomenon has been well studied via hybrid simulations of collisionless shocks (e.g., Leroy et al. 1981, 1982; Leroy 1983). As shown in Figure 1(b), the potential energy $e{\rm{\Phi }}$ reaches $\sim 60 \%$ of the incoming proton energy ${m}_{i}{V}_{1}^{2}/2$. As a result, a significant fraction of the incoming protons are reflected back toward the upstream, leading to a pileup of particles just in front of the shock. The reflected protons can be identified as the ones with positive ${p}_{i,x}$ and ${p}_{i,z}$ ahead of the shock in the phase spaces of Figures 1(c) and (d), respectively. As the reflected protons gyrate in the shock-compressed magnetic field, they gain energy from the upstream motional electric field. Upon their second encounter with the cross-shock potential, the reflected protons now have sufficient energy to penetrate the shock. In the downstream region just behind the shock, the protons keep gyrating in the x–z plane perpendicular to the shock-compressed magnetic field (compare the phase spaces in Figures 1(c) and (d) at $-4\lesssim (x-{x}_{\mathrm{sh}})/{r}_{\mathrm{Li}}\lesssim 0$). The peaks in density seen in Figure 1(a) are then correlated with the locations where the proton gyrophase is such that most protons have small ${p}_{i,x}$ (e.g., at $x-{x}_{\mathrm{sh}}\sim -0.25\,{r}_{\mathrm{Li}}$, $-1.25\,{r}_{\mathrm{Li}}$, and $-2.75\,{r}_{\mathrm{Li}}$). The amplitude of the density oscillations gets smaller as the gyrating reflected protons become more and more phase-mixed with the directly transmitted protons, at $x-{x}_{\mathrm{sh}}\lesssim -5\,{r}_{\mathrm{Li}}$.

Since the postshock protons gyrate in the x–z plane perpendicular to the field, the momentum dispersion along the y direction of the field is expected to be nearly the same on the two sides of the shock (see the ${p}_{i,y}$ phase space in Figure 1(e) near the shock). Farther behind the shock, the dispersion in ${p}_{i,y}$ increases. This can be also quantified with the y-averaged profiles of the proton temperature perpendicular (${T}_{i,\perp }$) and parallel (${T}_{i,\parallel }$) to the background magnetic field, as in Figure 1(f). Here, the jk component of the temperature tensor is defined as ${k}_{{\rm{B}}}{T}_{{jk}}/{m}_{i}{c}^{2}\equiv \langle \gamma ^{\prime} {v}_{j}^{{\prime} }{v}_{k}^{{\prime} }\rangle /{c}^{2}$, where ${v}_{j}^{{\prime} },{v}_{k}^{{\prime} }$ are the particle velocities in the fluid comoving frame and $\gamma ^{\prime}$ is the comoving particle Lorentz factor, and the average is performed over the particle distribution at a given spatial location. As Figure 1(f) shows, the mean proton temperature T_i, defined as⁴

$$\begin{eqnarray}&&{T}_{i}=\displaystyle \frac{2\,{T}_{i,\perp }+{T}_{i,\parallel }}{3},\end{eqnarray} \tag{ 22 }$$

is nearly uniform in the downstream region (green line), but the parallel temperature (orange line)—which is continuous across the shock—increases with distance behind the shock, while the perpendicular temperature (blue line) shoots up at the shock and then experiences a modest decline. This is the same trend shown by the phase spaces in Figures 1(c)–(e).

The decrease in perpendicular temperature and the resulting increase in parallel temperature suggest that protons are being scattered in pitch angle. In fact, in the region $-4\lesssim (x-{x}_{\mathrm{sh}})/{r}_{\mathrm{Li}}\lesssim -1$, where the variation in ${T}_{i,\perp }$ and ${T}_{i,\parallel }$ is most pronounced, strong magnetic waves are observed in Figure 2. Their wavelength is comparable to the proton skin depth, indicating that they are driven by protons (as opposed to electrons). In Figure 2(a), we compare the 1D profiles (averaged over the y direction) of the magnetic fluctuations $\delta {B}_{x}^{2}$, $\delta {B}_{y}^{2}$, and $\delta {B}_{z}^{2}$, normalized to ${B}_{\mathrm{ff}}^{2}$, where ${B}_{\mathrm{ff}}$ is defined as the magnitude of the flux-frozen magnetic field (i.e., ${{\boldsymbol{B}}}_{\mathrm{ff}}\equiv {B}_{0}(n/{n}_{0})\hat{y}$, where n is the y-averaged particle density).⁵ The energy of proton-driven waves peaks at $x\,-{x}_{\mathrm{sh}}\sim -2.5\,{r}_{\mathrm{Li}}$. In Figure 1(a), they are responsible for the excess of magnetic field strength (red curve) above the flux-freezing prediction (which would correspond to the density profile, in black).

The dominant mode at $-4\lesssim (x-{x}_{\mathrm{sh}})/{r}_{\mathrm{Li}}\lesssim -1$ in the x and z directions has a wave vector nearly parallel to the background field (Figures 2(b) and (c)), consistent with the proton cyclotron instability (Kennel 1966; Davidson & Ogden 1975). The waves in $\delta {B}_{y}$ are slightly weaker (compare the green line with the blue and orange curves in Figure 2(a)) and have oblique wave vectors (Figure 2(d)), as expected for the mirror mode (Chandrasekhar et al. 1958; Barnes 1966; Hasegawa 1975; McKean et al. 1993). The presence of mirror modes breaks the flux freezing condition, as shown by the fact that in Figure 1(a) the y-averaged transverse magnetic field profile ${B}_{y}/{B}_{0}$ deviates at $-5\lesssim (x-{x}_{\mathrm{sh}})/{r}_{\mathrm{Li}}\lesssim -2$ from the density profile (in black, which tracks the flux freezing prediction).

Both the proton cyclotron instability and the mirror instability are sourced by proton temperature anisotropy. In fact, since the motion of downstream protons right behind the shock is mostly confined in the x–z plane, a large temperature anisotropy arises, with ${T}_{i,\perp }\gg {T}_{i,\parallel }$ (Figure 1(g)). The anisotropy provides free energy for the growth of proton cyclotron waves and mirror modes, which scatter the protons in pitch angle and reduce their anisotropy back to the upper bound corresponding to marginal stability (Gary et al. 1997; see the red dashed line in Figure 1(g)), which is

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{i,\perp }}{{T}_{i,\parallel }}-1\simeq \displaystyle \frac{1.1}{{\beta }_{i,\parallel }^{0.55}}.\end{eqnarray} \tag{ 23 }$$

Here, ${\beta }_{i,\parallel }=8\pi {n}_{i}{k}_{{\rm{B}}}{T}_{i,\parallel }/{B}^{2}$ is the local value of the proton plasma beta, computed with the parallel proton temperature.

In this subsection, we describe the electron dynamics, with a focus on electron heating in the shock layer and in the downstream region. Due to their opposite charge and much smaller Larmor radius, the dynamics of electrons is drastically different from that of the protons.

Figure 3(a) shows the electron density profile (black line), which strongly resembles that of the protons (black line in Figure 1(a)) and thus ensures approximate charge neutrality. While a small degree of charge separation at the shock is responsible for establishing the electric potential Φ shown in Figure 1(b), the fact that Φ is nearly uniform at $x-{x}_{\mathrm{sh}}\lesssim -5\,{r}_{\mathrm{Li}}$ suggests that charge neutrality is satisfied very well in the far downstream.

Figures 3(b)–(d) shows the electron phase space. Since electrons have opposite charge than protons, they are not reflected back upstream by the cross-shock potential. In fact, unlike for protons, there is no reflected electron population with ${p}_{e,x}\gtrsim 0$ just ahead of the shock (compare Figures 3(b) and 1(c)).

Figure 3(e) shows the temperature profile of electrons, for the perpendicular component ${T}_{e,\perp }$ (blue), the parallel component ${T}_{e,\parallel }$ (orange), and the mean temperature T_e (green), which is defined as

$$\begin{eqnarray}&&{T}_{e}=\displaystyle \frac{2\,{T}_{e,\perp }+{T}_{e,\parallel }}{3}.\end{eqnarray} \tag{ 24 }$$

The profile of perpendicular temperature (blue line) follows closely the density compression (compare with the black line in Figure 3(a)) and starts to rise just ahead of the shock at $x-{x}_{\mathrm{sh}}\sim 0.5\,{r}_{\mathrm{Li}}$. This is consistent with the double adiabatic theory, also known as the CGL theory (Chew et al. 1956), which predicts ${T}_{\perp }\propto B$ (and in flux freezing, $B\propto n$).⁶ The double adiabatic theory applies to electrons, since the density and magnetic field compression occurs on scales much larger than the electron Larmor radius. This is not true for protons, since the shock thickness and the scale length of the downstream oscillations seen in Figure 1(a) are set by the proton Larmor radius $\,{r}_{\mathrm{Li}}$.

The parallel electron temperature (orange line in Figure 3(e)) initially remains unchanged, as the CGL theory predicts ${T}_{\parallel }\propto {(n/B)}^{2}$ and $B\propto n$ as a result of flux freezing (compare the green and black lines in Figure 1(a) in the vicinity of the shock). The increase in perpendicular temperature at the shock, while the parallel temperature stays the same as in the upstream, leads to a strong electron anisotropy, up to ${T}_{e,\perp }/{T}_{e,\parallel }-1\sim 0.6$ (Figure 3(f)). This excites the electron whistler instability, which creates the small-wavelength transverse magnetic waves in $\delta {B}_{x}$ and $\delta {B}_{z}$ seen in the region $x-{x}_{\mathrm{sh}}\sim -0.25\,{r}_{\mathrm{Li}}$ of Figures 2(b) and (c) (see also the magnetic energy in $\delta {B}_{x}^{2}$ and $\delta {B}_{z}^{2}$ in Figure 2(a), at the same location). The electron whistler instability provides a mechanism for electron pitch-angle scattering and thus reduces the electron temperature anisotropy, as shown in the downstream region of Figure 3(f).

Zoom In Zoom Out Reset image size

As we have already discussed, if the electron fluid were to follow the double adiabatic predictions, ${T}_{e,\perp }\propto n$ and ${T}_{e,\parallel }\,\propto$ const. The fact that the perpendicular temperature profile in Figure 3(f) (blue line) resembles the density profile (black line in Figure 3(a)) and the fact that ${T}_{e}\sim {T}_{e,\perp }$ (compare the green and blue curves in Figure 3(e)) suggest that most of the increase in electron temperature comes from adiabatic compression. However, the fact that ${T}_{e,\parallel }$ is not constant across the shock requires nonadiabatic processes. In order to quantify the degree of nonadiabatic (or "irreversible") electron heating, we compare in Figure 3(g) the mean electron temperature T_e with the adiabatic prediction

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{e,\mathrm{ad}}}{{T}_{e0}}={\left(\displaystyle \frac{{n}_{e}}{{n}_{e0}}\right)}^{2/3}.\end{eqnarray} \tag{ 25 }$$

This estimate of the adiabatic temperature assumes an isotropic gas, which is valid, given the small degree of electron anisotropy far downstream of the shock (see Figure 3(f) at $x-{x}_{\mathrm{sh}}\lesssim -1\,{r}_{\mathrm{Li}}$). Figure 3(g) shows the excess of T_e above ${T}_{e,\mathrm{ad}}$ in units of the upstream electron temperature. Most of the irreversible heating occurs at two locations: $x-{x}_{\mathrm{sh}}\sim 0$, that is, in the shock transition region; and $x-{x}_{\mathrm{sh}}\sim -2.5\,{r}_{\mathrm{Li}}$, where the density suffers another compression and strong proton-driven waves are generated (see Figure 2). These two locations are marked by the vertical dotted lines in Figures 2 and 3, and the particle and wave properties there will be further studied below. In the far downstream, the temperature excess over the adiabatic estimate saturates at ${T}_{e}-{T}_{e,\mathrm{ad}}\sim 0.3\,{T}_{e0}$ (Figure 3(g)).

Zoom In Zoom Out Reset image size

An alternative (and possibly more rigorous) estimate of the degree of irreversible electron heating is given by the specific entropy s_e (i.e., the entropy per particle), measured with the electron distribution function f_e as

$$\begin{eqnarray}&&{s}_{e}\equiv -\displaystyle \frac{\int {d}^{3}p\ {f}_{e}\mathrm{ln}{f}_{e}}{\int {d}^{3}p\ {f}_{e}}\,,\end{eqnarray} \tag{ 26 }$$

where the normalization is such that $\int {d}^{3}p\,{f}_{e}={n}_{e}/{n}_{e0}$. To construct the spatial profile of s_e(x), we first bin the particles by their x position with a width of ${\rm{\Delta }}x=100$ cells. In each spatial bin, we compute ${f}_{e}({\boldsymbol{p}})$ by constructing a three-dimensional histogram of the particle momenta. In each direction (${p}_{e,x}$, ${p}_{e,y}$, and ${p}_{e,z}$), the central bin of the histogram lies at the mean momentum, and the histogram spans four standard deviations above and below the mean. Each standard deviation is resolved with 10 momentum bins.

Figure 3(h) shows the change of electron entropy with respect to the upstream value. In analogy to Figure 3(g), the increase in electron entropy is localized around $x-{x}_{\mathrm{sh}}\sim 0$ and $x-{x}_{\mathrm{sh}}\sim -2.5\ {r}_{\mathrm{Li}}$ (indicated by the gray and pink vertical dotted lines, respectively). The increase in electron specific entropy saturates at ${\rm{\Delta }}{s}_{e}\sim 0.25$ in the far downstream.

4.2.1. Electron Whistler Waves

The physics of particle irreversible heating that we have described in Section 2 relies on two ingredients: a certain level of particle anisotropy and a mechanism to break the adiabatic invariance. As we have shown above, a large-scale magnetic field amplification (e.g., resulting from shock compression of the upstream field) will lead to electron anisotropy with ${T}_{e,\perp }\gt {T}_{e,\parallel }$. In turn, this triggers the growth of whistler waves, which scatter the electrons in pitch angle, providing a mechanism to break the adiabatic invariance and generate irreversible heat. Below, we show that the two ingredients needed for entropy increase are indeed present in the two locations where the entropy profile shows the fastest increase (vertical dotted lines in Figures 2 and 3).

At the shock (gray dotted line in Figures 2 and 3), the electron temperatures are driven to ${T}_{e,\perp }\gt {T}_{e,\parallel }$ by shock compression of the upstream field, via conservation of the first and second adiabatic invariants. In Figure 4, we show the spacetime diagram of various quantities, in the time interval $20.0\leqslant {{\rm{\Omega }}}_{{ci}}t\leqslant 27.4$ and along the y extent of the box. The x location is fixed at the shock ramp (more precisely, $x-{x}_{\mathrm{sh}}=4\,c/{\omega }_{{pe}}$). Shock compression of the upstream field (see Figure 4(a), where $| B| /{B}_{0}\sim 2.2$) leads to a temperature anisotropy ${T}_{e,\perp }/{T}_{e,\parallel }-1\simeq 0.6$ (Figure 4(d)). Both the field amplification and the degree of temperature anisotropy are nearly constant in time and uniform in y.

Zoom In Zoom Out Reset image size

As a result of the strong temperature anisotropy, magnetic waves are excited throughout the y range consistently over time. Panels (b) and (c) show the spacetime diagrams of the magnetic fluctuations $\delta {B}_{x}$ and $\delta {B}_{z}$, revealing the presence of high-frequency and short-wavelength modes (as also seen in Figures 2(b) and (c) near the shock). Figures 4(e) and (f) show the corresponding power spectra, as a function of frequency ω (horizontal axis) and wavenumber k_y (vertical axis). The power spectrum displays a pronounced peak at frequency $\omega \simeq 0.5\,{{\rm{\Omega }}}_{{ce}}$ (here ${{\rm{\Omega }}}_{{ce}}=({m}_{i}/{m}_{e}){{\rm{\Omega }}}_{{ci}}$ is the electron gyrofrequency) and wave vector ${k}_{y}\simeq 0.5\,{\omega }_{{pe}}/c$. We have compared this with linear theory of the electron whistler instability (e.g., Gary & Madland 1985; Gary & Wang 1996; Gary & Karimabadi 2006) by solving the dispersion relation

$$\begin{eqnarray}0 & = & {D}^{\pm }({k}_{y},{\rm{\Omega }})\\ & = & {{\rm{\Omega }}}^{2}-{c}^{2}{k}_{y}^{2}+{\omega }_{{pi}}^{2}\displaystyle \frac{{\rm{\Omega }}}{{k}_{y}{v}_{i}}Z({\zeta }_{i}^{\pm })\\ & & +{\omega }_{{pe}}^{2}\displaystyle \frac{{\rm{\Omega }}}{{k}_{y}{v}_{e,\parallel }^{2}}Z({\zeta }_{e}^{\pm })+{\omega }_{{pe}}^{2}\left(\displaystyle \frac{{T}_{e\perp }}{{T}_{e,\parallel }}-1\right)[1+{\zeta }_{e}^{\pm }Z({\zeta }_{e}^{\pm })]\,,\end{eqnarray} \tag{ 27 }$$

where ${\zeta }_{e}^{\pm }=({\rm{\Omega }}\pm {{\rm{\Omega }}}_{{ce}})/{k}_{y}{v}_{e,\parallel },$ ${v}_{e,\parallel }={(2{k}_{{\rm{B}}}{T}_{e,\parallel }/{m}_{e})}^{1/2},$ ${\zeta }_{i}^{\pm }\,=({\rm{\Omega }}\pm {{\rm{\Omega }}}_{{ci}})/{k}_{y}{v}_{i\parallel },$ ${v}_{i}={(2{k}_{{\rm{B}}}{T}_{i}/{m}_{i})}^{1/2},$ and $Z(\zeta )$ is the plasma dispersion function:

$$\begin{eqnarray}&&Z(\zeta )=\displaystyle \frac{1}{\sqrt{\pi }}{\int }_{-\infty }^{\infty }{dx}\displaystyle \frac{\exp (-{x}^{2})}{x-\zeta }.\end{eqnarray} \tag{ 28 }$$

The input values of ${v}_{i},{v}_{e,\parallel },{T}_{e\perp }/{T}_{e,\parallel }-1$ for the dispersion relation are taken from the time and space averages of the corresponding quantities over the same time period and spatial extent as the spacetime diagram in Figure 4. The resulting theoretical prediction for the real part of the frequency is shown with a black solid line in panels (e) and (f), and it matches extremely well the contours of the power spectrum. The imaginary part of the frequency, that is, the growth rate of the mode, is plotted with a dashed black curve. The value of k_y giving the fastest growth agrees well with the location of the peak of the power spectrum (${k}_{y}\simeq 0.5\,{\omega }_{{pe}}/c$). The excellent agreement between the simulation data and the electron whistler dispersion relation confirms that the waves in the shock ramp are produced by the electron whistler instability.

Figure 5 shows similar plots at the location indicated in Figures 2 and 3 with a vertical dotted pink line, at $x-{x}_{\mathrm{sh}}\simeq -2.5\,{r}_{\mathrm{Li}}\simeq -122\,c/{\omega }_{{pe}}$. Here, field amplification is driven by a combination of two effects: the density (and so, the frozen-in magnetic field) experiences another large-scale compression, and the proton-driven waves shown in Figure 2 further increase the local magnetic field intensity.

Zoom In Zoom Out Reset image size

As compared to Figure 4, the spacetime diagrams show now a higher degree of inhomogeneity, imprinted by the anisotropy-driven long-wavelength proton modes. These fluctuations coexist with weaker small-wavelength, high-frequency modes, which only appear in localized patches (e.g., at $x\sim 80\,c/{\omega }_{{pe}}$ and $t\sim 100\,{{\rm{\Omega }}}_{{ce}}^{-1}$ in Figures 5(b) and (c)). The high-frequency waves are generated in regions where field amplification (Figure 5(a) at $x\sim 80\,c/{\omega }_{{pe}}$ and $t\sim 100\,{{\rm{\Omega }}}_{{ce}}^{-1}$) causes the electron anisotropy (Figure 5(d)) to exceed the threshold for whistler growth (Figure 5(g)), which is given by

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{e,\perp }}{{T}_{e,\parallel }}-1\simeq \displaystyle \frac{0.21}{{\beta }_{e,\parallel }^{0.6}},\end{eqnarray} \tag{ 29 }$$

where ${\beta }_{e,\parallel }=8\pi {n}_{e}{k}_{{\rm{B}}}{T}_{e,\parallel }/{B}^{2}$ is the local value of the parallel electron beta (Gary 2005).

Figures 5(e) and (f) show the power spectra of $\delta {B}_{x}$ and $\delta {B}_{z}$. Most of the power is concentrated in low-frequency, long-wavelength modes, generated by the proton cyclotron or mirror instabilities. However, there is still an appreciable amount of power in high-frequency, short-wavelength modes peaking at $\omega \sim 0.7\,{{\rm{\Omega }}}_{{ce}}$ and ${k}_{y}\sim 0.7\,{\omega }_{{pe}}/c$. We apply a high-frequency, short-wavelength filter, in order to isolate the top right region in Figures 5(e) and (f) (the cutoff frequency and wavenumber of our filter are shown with dashed red lines). This allows us to extract (via an inverse Fourier transform) the spacetime wave patterns of high-frequency, short-wavelength modes, which are shown in panels (h) and (i). The two panels confirm that short-wavelength modes exist only in regions where the electron temperature anisotropy exceeds the electron whistler threshold (Figure 5(g) at $x\sim 80\,c/{\omega }_{{pe}}$ and $t\sim 100\,{{\rm{\Omega }}}_{{ce}}^{-1}$). We have measured the average electron and proton temperatures and densities in the region where the whistler threshold is appreciably exceeded (${T}_{e,\perp }/{T}_{e,\parallel }-1-0.21/{\beta }_{e,\parallel }^{0.6}\geqslant 0.3$), in order to obtain linear theory predictions. The real part and imaginary part of the resulting dispersion relation are plotted in Figures 5(e) and (f) with solid and dashed black lines, respectively. The good agreement with the power spectra extracted from our shock simulation confirms the presence of patches of whistler waves in the second ramp (at $x-{x}_{\mathrm{sh}}\sim -2.5\,{r}_{\mathrm{Li}}$) of the electron entropy profile.

To summarize, we have identified two major sites of electron entropy production in the shock downstream. One is at the shock ramp, and the other is at a distance of $\sim 2.5\,{r}_{\mathrm{Li}}$ behind the shock, where density compression and proton-driven waves both contribute to magnetic field amplification. Both sites show the presence of electron whistler waves triggered by electron temperature anisotropy. Whistler waves provide the pitch-angle scattering required to break electron adiabatic invariance and to generate entropy. In the following two sections, we further elucidate the physics of entropy production in these two sites, by means of periodic box simulations. The reader primarily interested in the implications for shocks can skip Sections 5 and 6 and proceed directly to Section 7, where the heating model is validated in the shock simulation.

To emulate the conditions for electrons in the shock ramp, we set up a suite of compressing box experiments, using the method introduced in Sironi & Narayan (2015) and Sironi (2015). Here, we report only its main properties. We solve Maxwell's equations and the Lorentz force in the fluid comoving frame, which is related to the laboratory frame by a Lorentz boost. In the comoving frame, we define two sets of spatial coordinates with the same time coordinate. The unprimed coordinate system has a basis of unit vectors, so it is the appropriate coordinate set to measure all physical quantities. Yet, it is convenient to redefine the unit length of the spatial axes in the comoving frame such that a particle subject only to compression stays at fixed coordinates. This will be our primed coordinate system. Then, compression with rate q is accounted for by the diagonal matrix

$$\begin{eqnarray}{\boldsymbol{L}}=\displaystyle \frac{\partial {\boldsymbol{x}}}{\partial {\boldsymbol{x}}^{\prime} }=\left(\begin{array}{ccc}{(1+qt)}^{-1} & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{array}\right),\end{eqnarray} \tag{ 30 }$$

which has been tailored for compression along the x axis, as expected in our shock.

A uniform ordered magnetic field ${{\boldsymbol{B}}}_{0}$ is initialized along the y direction (in analogy to the shock setup). We define ${{\rm{\Omega }}}_{{ci}}$ as the proton Larmor frequency in the initial field B₀. Maxwell's equations in the primed coordinate system (Sironi & Narayan 2015) prescribe that the field will grow in time as ${{\boldsymbol{B}}}_{0}(1+q\,t)$, which is consistent with flux freezing (the particle density in the box increases at the same rate). From the Lorentz force in the compressing box (Sironi & Narayan 2015), the component of particle momentum aligned with the field does not change during compression, whereas the perpendicular momentum increases as $\propto \sqrt{1+q\,t}$. This is consistent with the conservation of the first and second adiabatic invariants.

This method is implemented for 1D, 2D, and 3D computational domains, with periodic boundary conditions in all directions. In the previous section, we have shown that the whistler instability is the dominant mode in the shock ramp. Its wave vector is nearly aligned with the field direction (i.e., along $\hat{y}$). It follows that the evolution of the dominant mode can be conveniently captured by means of 1D simulations with the computational box oriented along y, which we will be employing in this section. Yet, all three components of electromagnetic fields and particle velocities are tracked. In 1D simulations, we can employ a large number of particles per cell (typically 10⁴ per cell), so we have adequate statistics for the calculation of the electron specific entropy from the phase space distribution function. In addition, in 1D simulations, we can readily extend our results up to the realistic mass ratio. Even though we only show results from 1D runs, we have checked that the main conclusions hold in 2D.

As a result of the large-scale compression encoded in Equation (30), both electrons and protons will develop a temperature anisotropy, and we should witness the development of both electron and proton anisotropy-driven modes. However, in our reference shock, no proton modes grow in the shock ramp (they only develop a few Larmor radii behind the shock). For this reason, in our compressing box runs, we artificially inhibit the update of the proton momentum (effectively, this corresponds to the case of infinitely massive protons, which only serve as a charge-neutralizing fluid).

The compression rate q is measured directly from our reference shock simulation. There, we can quantify the profile of electron density as a function of the comoving time of the electron fluid, which follows from

$$\begin{eqnarray}&&\tau \equiv \int \displaystyle \frac{{dx}^{\prime} }{{V}_{{xe}}(x^{\prime} )},\end{eqnarray} \tag{ 31 }$$

where V_xe is the electron fluid velocity in the shock frame, and the integral goes from the upstream to the downstream region. Figure 6 shows the density profile as a function of τ from our reference run. The density oscillates on a timescale comparable to the proton gyration time ${{\rm{\Omega }}}_{{ci}}^{-1}$, which is expected given that the proton dynamics controls the shock structure. At the ramp starting near $\tau =0$, the electron density increases by a factor of ∼3.5 within $\sim 1\,{{\rm{\Omega }}}_{{ci}}^{-1}$. Even though the density increase is not perfectly linear, we find that a linear approximation with $q=2.5\,{{\rm{\Omega }}}_{{ci}}$ provides a reasonable fit (see the orange dashed line in Figure 6). We remark that our electron heating model is agnostic of the exact profile of density compression, as long as the compression rate and the resulting field amplification rate are much slower than the electron gyration frequency.

Zoom In Zoom Out Reset image size

Below, we fix $q=2.5\,{{\rm{\Omega }}}_{{ci}}$. With increasing mass ratio, the separation between q and the electron cyclotron frequency ${{\rm{\Omega }}}_{{ce}}$ will increase as ${m}_{i}/{m}_{e}$. As in our reference shock run, electrons are initialized to have ${k}_{{\rm{B}}}{T}_{e0}={10}^{-2}\,{m}_{e}{c}^{2}$, and the strength of the background magnetic field is set so that ${\beta }_{p0}=16$. We resolve the electron skin depth with 10 cells, so the Debye length is marginally resolved. The box extent along the y direction is fixed at $43\ c/{\omega }_{{pe}}$, which is sufficient to capture several wavelengths of the electron whistler instability. The box size does not need to scale with the mass ratio, since we are artificially excluding the proton physics.

As in our reference shock run, we employ a reduced mass ratio ${m}_{i}/{m}_{e}=49$. In the periodic compressing box, this means that our choice of $q=2.5\,{{\rm{\Omega }}}_{{ci}}$ leads to a compression rate that is a factor of ∼20 lower than the electron gyration frequency.

To highlight the importance of the electron whistler instability in facilitating electron entropy production, in Figure 7 we compare two simulations, one with the background field ${{\boldsymbol{B}}}_{0}$ in the z direction, the other one with ${{\boldsymbol{B}}}_{0}$ along the y direction. Since our simulation box is oriented along y and the dominant wave vector of the electron whistler instability is parallel to the background field, if the field lies along z (which we shall call the "out-of-plane" case and indicate with dotted lines), we artificially suppress the growth of electron whistlers. By comparing it with the "in-plane" simulation with the field along y (solid lines), which does allow for whistler wave growth, we can demonstrate the importance of the electron whistler instability for entropy production.

Zoom In Zoom Out Reset image size

In the absence of electron-scale instabilities that would break the adiabatic invariance, the out-of-plane simulation is expected to follow adiabatic scalings. In fact, in the out-of-plane simulation, we see that ${T}_{e,\perp }\propto B\propto (1+{qt})$ (blue dotted line in Figure 7(b)), while ${T}_{e,\parallel }\propto {(n/B)}^{2}\propto$ const. (orange dotted line in Figure 7(b)), as expected from the double adiabatic theory. Since no whistler waves grow (notice that the fields stay at the noise level; see the dotted lines of Figure 7(a)), no mechanism exists that can transfer heat from the perpendicular to the parallel temperature, and the electron entropy remains constant.

The in-plane simulation shows a different behavior. Initially, ${T}_{e,\perp }$ and ${T}_{e,\parallel }$ follow the double adiabatic trends (solid lines in Figure 7(b)). At ${{\rm{\Omega }}}_{{ci}}t\sim 0.3$, the increasing temperature anisotropy (blue solid line in panel (c)) leads to the exponential growth of electron whistler waves (solid lines in Figure 7(a)). At ${{\rm{\Omega }}}_{{ci}}t\sim 0.4$, the wave energy is strong enough to pitch-angle scatter the electrons. As a result, heat is transferred from the perpendicular to the parallel direction. Both ${T}_{e,\perp }$ and ${T}_{e,\parallel }$ deviate from the adiabatic scalings, and the temperature anisotropy is reduced.

At $t\sim 0.4\,{{\rm{\Omega }}}_{{ci}}^{-1}$, with the onset of pitch-angle scattering and the consequent breaking of adiabatic invariance, the electron specific entropy starts to rise (solid blue line in Figure 7(d)). The most rapid entropy increase happens near the end of the exponential whistler growth, at $t\sim 0.4\mbox{--}0.5\,{{\rm{\Omega }}}_{{ci}}^{-1}$. Here, the electron anisotropy is still large, and at the same time whistler waves are sufficiently powerful to provide effective pitch-angle scattering. In other words, both terms in the square brackets of either Equation (13) or (14) are large. After the exponential growth, the electron whistler waves enter a secular phase where the wave energy (normalized to the compressed background field energy) stays almost constant (solid green line in Figure 7(a)). In this phase, whistler waves are continuously generated as the large-scale compression steadily pushes the electron anisotropy slighly above the threshold of marginal stability (indicated by the red dashed line in Figure 7(c)). Both the ingredients needed for entropy increase (i.e., nonzero electron anisotropy and efficient pitch-angle scattering mediated by whistler waves) persist during the secular phase, leading to a further increase in the electron entropy.

In Figure 7, we also explicitly validate the heating model described in Section 2. Following Equation (7), the electron energy per particle should change as

$$\begin{eqnarray}{{du}}_{e} & = & {k}_{{\rm{B}}}{T}_{e,\perp }d\mathrm{ln}B+{k}_{{\rm{B}}}{T}_{e,\parallel }d\mathrm{ln}\left(\displaystyle \frac{n}{B}\right)-{{de}}_{w,e}\\ & = & {k}_{{\rm{B}}}{T}_{e,\perp }d\mathrm{ln}B-{{de}}_{w,e},\end{eqnarray} \tag{ 32 }$$

where ${e}_{w,e}$ is the energy per particle in whistler waves (as we have discussed in Section 2, the electron energy transferred to electron modes stays entirely in the waves, so ${{de}}_{w,e}={{de}}_{w\mathrm{tot},e}$), and we have used the fact that $n/B\propto$ const. In Figure 7(e), the blue solid line shows the measured change of electron internal energy from the in-plane run, while the red dashed line is obtained by integrating Equation (32). We find excellent agreement between simulation results and our electron heating model.

The validation can also be extended to the entropy measurement, as we do in Figure 7(d). Again, the blue solid line shows the measured change in electron specific entropy (computed from the distribution function as in Equation (26)), while the red dashed line is obtained by integrating

$$\begin{eqnarray}&&{{ds}}_{e}=\left(\displaystyle \frac{1}{2}d\mathrm{ln}{T}_{e,\parallel }\right)\left(1-\displaystyle \frac{{T}_{e,\parallel }}{{T}_{e,\perp }}\right)-\displaystyle \frac{{{de}}_{w,e}}{{T}_{e,\perp }},\end{eqnarray} \tag{ 33 }$$

which follows from Equation (13) (an equivalent form can be obtained from Equation (14)). Once again, the model matches the simulation results extremely well.

We now extend our compressing box experiments up to the realistic mass ratio and show that the electron entropy increase is nearly insensitive to ${m}_{i}/{m}_{e}$ (as long as the mass ratio is larger than a few tens). Figure 8 compares the evolution of the whistler wave energy (panel (a)), the electron temperature anisotropy (panel (b)), the rate $-d\mathrm{ln}({T}_{e,\perp }/B)$ of breaking adiabatic invariance (panel (c)), and the electron entropy increase (panel (d)) when varying the mass ratio from ${m}_{i}/{m}_{e}=49$ up to ${m}_{i}/{m}_{e}=1600$ (from purple to red; see the legend in the second panel). Since we fix the compression rate to be $q=2.5\,{{\rm{\Omega }}}_{{ci}}$, a larger mass ratio corresponds to a lower compression rate in units of the electron gyration frequency ${{\rm{\Omega }}}_{{ce}}=({m}_{i}/{m}_{e}){{\rm{\Omega }}}_{{ci}}$.

Zoom In Zoom Out Reset image size

Initially, the electron temperature anisotropy grows linearly in time as ${T}_{e,\perp }/{T}_{e,\parallel }-1={qt}$, as a result of the large-scale compression. This proceeds until the energy in whistler waves reaches a fraction $\sim 3\times {10}^{-2}$ of the compressed background field energy (Figure 8(a)). At this point, whistler waves are sufficiently strong to scatter the electrons in pitch angle, breaking their adiabatic invariance and reducing the electron anisotropy by transferring energy from the perpendicular to the parallel component. In fact, notice that the peak in panel (c), that is, the time when the electron adiabatic invariance is most violently broken, always corresponds to the time when the electron anisotropy in panel (b) shows the sharpest decrease.

The onset of efficient pitch-angle scattering (and so, the peak time of electron anisotropy) occurs earlier at higher mass ratio, at a time that decreases from $t\sim 0.35\,{{\rm{\Omega }}}_{{ci}}^{-1}$ at ${m}_{i}/{m}_{e}=49$ down to $t\sim 0.1\,{{\rm{\Omega }}}_{{ci}}^{-1}$ at ${m}_{i}/{m}_{e}=1600$. This can be understood from the competition between the large-scale compression rate (which increases the electron anisotropy) and the growth rate of whistler waves (which try to reduce the anisotropy via pitch-angle scattering). The compression rate in units of the electron cyclotron frequency is $q=2.5({m}_{e}/{m}_{i}){{\rm{\Omega }}}_{{ce}}$, while the whistler growth rate (also in units of ${{\rm{\Omega }}}_{{ce}}$) depends on how much the anisotropy exceeds the whistler threshold in Equation (29). In order to balance the two rates, a higher anisotropy is needed for larger ${m}_{e}/{m}_{i}$, that is, for lower mass ratios. This has two consequences: first, the growth rate of the whistler instability (normalized to ${{\rm{\Omega }}}_{{ce}}$) will decrease at higher ${m}_{i}/{m}_{e}$, as indeed confirmed by the inset of Figure 8(a); second, lower peak anisotropies (and so, earlier onsets of efficient pitch-angle scattering) will be achieved at higher mass ratios, which explains the trend seen in Figure 8(b). In addition, since the energy of whistler waves ultimately comes from the free energy in electron anisotropy, higher mass ratios display weaker levels of whistler wave activity (panel (a)).

The electron entropy evolution in Figure 8(d) can be separated into two stages. In the first phase (which, for ${m}_{i}/{m}_{e}=49$, occurs at $t\sim 0.45\,{{\rm{\Omega }}}_{{ci}}^{-1}$), the electron entropy grows quickly. This stage corresponds to the late exponential phase of whistler wave growth (and so, we shall call it the "exponential phase"), when both the electron anisotropy (panel (b)) and the rate of breaking adiabatic invariance (panel (c))—that is, the two ingredients needed for efficient entropy production—are large. Since higher mass ratios reach lower levels of electron anisotropy, the entropy produced during this stage is a decreasing function of ${m}_{i}/{m}_{e}$, as seen in Figure 8(d) (compare the purple line growth around $t\sim 0.45\,{{\rm{\Omega }}}_{{ci}}^{-1}$ with the red line around $t\sim 0.15\,{{\rm{\Omega }}}_{{ci}}^{-1}$). After whistler waves have reached saturation, the electron entropy still increases, in a phase that we shall call "secular." Here, the electron anisotropy stays close to the threshold of marginal stability (indicated in Figure 8(c) by the dotted lines, with the same color coding as the solid curves). Continuous pitch-angle scattering (and so, persistent violation of adiabatic invariance) is needed to oppose the steadily driven compression and maintain the system close to marginal stability. It is then expected that entropy will continuously increase during the secular phase, albeit at a lower rate than in the exponential stage. For ${m}_{i}/{m}_{e}\gtrsim 400$, the electron anisotropy at late times is nearly insensitive to ${m}_{i}/{m}_{e}$ (compare yellow, orange, and red lines at ${{\rm{\Omega }}}_{{ci}}t\gtrsim 0.4$ in panel (b)), which explains why the entropy growth in the secular phase is nearly the same for all ${m}_{i}/{m}_{e}\gtrsim 400$ (Figure 8(d)).

From Figure 8(d), we can infer how the entropy increase in the shock ramp should scale with mass ratio. Since the compression in the shock ramp lasts about one proton gyration time, we compare the entropy curves at ${{\rm{\Omega }}}_{{ci}}t\sim 1$, as indicated by the vertical dotted black line in panel (d). When the mass ratio increases from ${m}_{i}/{m}_{e}=49$ to ${m}_{i}/{m}_{e}=1600$ (i.e., more than a factor of 32), the entropy produced until ${{\rm{\Omega }}}_{{ci}}t=1$ decreases from 0.065 to 0.048, only a $\sim 30 \%$ decrease. The dependence on mass ratio would be far more pronounced if we were only to consider the entropy produced during the exponential phase. However, higher mass ratio runs have earlier onset times, as we have explained above, so they spend more time (within the first ${{\rm{\Omega }}}_{{ci}}^{-1}$) in the secular phase, as compared to lower mass ratios. In summary, most of the entropy production at lower mass ratios happens during the exponential phase, whereas at higher mass ratios the secular phase lasts longer and thus compensates for the lower level of entropy generated during the exponential stage. The net effect is that the entropy increase in our compressing box with ${m}_{i}/{m}_{e}=1600$ is only slightly smaller than for ${m}_{i}/{m}_{e}=49$. The same conclusion should hold also for our reference shock.

In order to study the role of anisotropy-driven proton modes in producing electron irreversible heating, we set up a periodic simulation box with anisotropic protons. The simulation is initialized to approximate the conditions right after the shock transition. The protons are initialized as a bi-Maxwellian distribution with ${T}_{i0,\perp }/{T}_{i0,\parallel }\sim 7$, as observed just behind the shock in Figure 1(g). The value of ${T}_{i0,\parallel }$ is the same as in the shock upstream (in fact, the parallel proton temperature is nearly uniform across the shock; see the orange line in Figure 1(f)). The electron temperature increases roughly by a factor of two across the shock (Figure 3(e)), so the electrons in the tests here are initialized with ${T}_{e0}\sim 2\times {10}^{-2}\,{m}_{e}{c}^{2}/{k}_{{\rm{B}}}$ (note that in the shock upstream the electron temperature was ${10}^{-2}\,{m}_{e}{c}^{2}/{k}_{{\rm{B}}}$). We take electrons to be isotropic, since the fast growth of whistler waves in the shock ramp ensures that the degree of downstream electron anisotropy is low (see the postshock region in Figure 3(f)). We also take into account that both density and magnetic field strength have increased by a factor of ∼2.5 as compared to the shock upstream (Figure 1(a)).

We resolve the electron skin depth with seven cells in order to (marginally) capture the electron Debye length. Since both the proton cyclotron instability, which dominates over the mirror mode in the downstream of our reference run, and the electron whistler instability have the fastest growing wave vector aligned with the background field, we employ 1D simulation domains with the box aligned with the y direction of the background magnetic field. Thanks to the reduced dimensionality of our computational domain, we can employ a large number of particles per cell (10⁴). Therefore, we have adequate statistics for the calculation of the electron specific entropy, and we can properly control the effect of numerical heating. In addition, in 1D simulations, we can readily extend our results up to the realistic mass ratio. The length of the computational box is 1512 cells for ${m}_{i}/{m}_{e}=49$. Since the fastest growing mode of the proton cyclotron instability has a wave vector $\sim {\omega }_{{pi}}/c$, we increase the number of cells in our computational domain proportional to $\propto \sqrt{{m}_{i}/{m}_{e}}$, to include the same number of proton skin depths (and so, approximately the same number of proton cyclotron wavelengths).⁷

In order to compare the results obtained from the periodic box simulations with the reference shock run, in Figures 9 and 10 we show the evolution of our periodic system for ${m}_{i}/{m}_{e}=49$. As a result of the initial proton temperature anisotropy, the proton cyclotron instability develops, generating exponentially growing waves with $\delta {B}_{x}$ and $\delta {B}_{z}$ components (Figure 9(a)). As shown in Figures 10(b) and (c), the growing waves are dominated by modes with wavelength at the proton inertial scale (for ${m}_{i}/{m}_{e}=49$, the proton skin depth is $c/{\omega }_{{pi}}=7\,c/{\omega }_{{pe}}$) and frequency comparable to the proton gyration frequency, as expected for the proton cyclotron instability.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At $t\sim 4\,{{\rm{\Omega }}}_{{ci}}^{-1}$, when the energy in proton cyclotron waves reaches a fraction $\sim {10}^{-1}$ of the background field energy, efficient pitch-angle scattering quickly reduces the proton temperature anisotropy (Figure 9(c)). During the isotropization process, the proton specific entropy increases (Figure 9(e) at $t\sim 5\,{{\rm{\Omega }}}_{{ci}}^{-1}$).

The heating model described in Section 2 can be applied to protons, keeping in mind that in the current setup no perturbations in density or magnetic field exist on scales larger than the proton scales (so, no work is being done on the protons). It follows that the perpendicular and parallel energy per proton change as

$$\begin{eqnarray}&&{{du}}_{i,\perp }=-{{dq}}_{i,\perp \to \parallel }-{{de}}_{w,i\perp },\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{{du}}_{i,\parallel }={{dq}}_{i,\perp \to \parallel }-{{de}}_{w,i\parallel },\end{eqnarray} \tag{ 35 }$$

so the total change in proton energy per particle is

$$\begin{eqnarray}&&{{du}}_{i}=-{{de}}_{w,i\perp }-{{de}}_{w,i\parallel }\equiv -{{de}}_{w,i\mathrm{tot}},\end{eqnarray} \tag{ 36 }$$

which simply states that the energy lost by protons is transferred to proton waves. Following Section 2, the change in specific proton entropy is

$$\begin{eqnarray}&&{{ds}}_{i}=\left(\displaystyle \frac{1}{2}d\mathrm{ln}{T}_{i,\parallel }\right)\left(1-\displaystyle \frac{{T}_{i,\parallel }}{{T}_{i,\perp }}\right)-\displaystyle \frac{{{de}}_{w,i\mathrm{tot}}}{{T}_{i,\perp }},\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{{ds}}_{i}=(-d\mathrm{ln}{T}_{i,\perp })\left(\displaystyle \frac{{T}_{i,\perp }}{{T}_{i,\parallel }}-1\right)-\displaystyle \frac{{{de}}_{w,i\mathrm{tot}}}{{T}_{i,\parallel }}\,,\end{eqnarray} \tag{ 38 }$$

where the two expressions are equivalent, as with Equations (13) and (14). We now need to specify ${{de}}_{w,i\mathrm{tot}}$, that is, the energy per proton transferred to proton modes. As we anticipated in Section 2, this is not equal to the energy residing in proton waves, since some fraction of that is being used to perform work on the electrons. In the remainder of this section, we denote as n and B the density and magnetic field fluctuations induced by proton waves. Since the scale of the perturbations is much larger than the electron gyroradius, the fluctuations perform work on the electrons, so Equation (4) for electrons becomes

$$\begin{eqnarray}&&{{dw}}_{e}={T}_{e,\perp }d\mathrm{ln}B+{T}_{e,\parallel }d\mathrm{ln}\left(\displaystyle \frac{n}{B}\right)\equiv {{dw}}_{e,\perp }+{{dw}}_{e,\parallel }.\end{eqnarray} \tag{ 39 }$$

This energy increase in the electrons is at the expense of the energy in proton waves, so the residual energy per particle residing in proton waves will be

$$\begin{eqnarray}&&{{de}}_{w,i}={{de}}_{w,i\mathrm{tot}}-{{dw}}_{e,\perp }-{{dw}}_{e,\parallel },\end{eqnarray} \tag{ 40 }$$

and the energy equation for protons reads

$$\begin{eqnarray}&&{{du}}_{i}=-{{dw}}_{e,\perp }-{{dw}}_{e,\parallel }-{{de}}_{w,i},\end{eqnarray} \tag{ 41 }$$

where the three terms on the right-hand side can be explicitly measured in our simulations.

Figures 9(e) and (g) demonstrate that our heating model works remarkably well for protons (later on, we will show that it also works for electrons). In Figure 9(e), the blue solid line is the proton entropy change measured directly from the simulation, using the distribution function as we did in Equation (26). It matches extremely well the prediction obtained by integrating the right-hand side of the proton entropy equation, Equation (37) or (38) (see the orange dotted line in Figure 9(e)). The agreement is also remarkable regarding the proton energy equation, Equation (41). In Figure 9(g), the proton energy loss $-{\rm{\Delta }}{u}_{i}=-\int {{du}}_{i}$ is indicated as a green line. From Equation (41), this should be equal to ${\rm{\Delta }}{e}_{w,i}+{w}_{\perp }+{w}_{\parallel }$, where we have defined ${\rm{\Delta }}{e}_{w,i}=\int {{de}}_{w,i}$, ${w}_{\perp }=\int {{dw}}_{e,\perp }$, and ${w}_{\parallel }=\int {{dw}}_{e,\parallel }$. In fact, the green line nearly overlaps the red curve.

The growth of proton cyclotron waves provides a source of field amplification and density perturbations that can perform work on the electrons. Indeed, Figure 9(h) shows that during the exponential phase of the proton cyclotron instability ($4\lesssim {{\rm{\Omega }}}_{{ci}}t\lesssim 7$), the proton waves increase the electron perpendicular energy (i.e., ${{dw}}_{e,\perp }\gt 0$; see the blue line in Figure 9(h)) and decrease the parallel component (i.e., ${{dw}}_{e,\parallel }\lt 0$; see the orange line in Figure 9(h)). This leads to a temperature anisotropy ${T}_{e,\perp }\gt {T}_{e,\parallel }$ (compare the blue and orange lines in Figure 9(d) at ${{\rm{\Omega }}}_{{ci}}t\sim 5$), which can be equivalently explained as a result of the conservation of the first and second adiabatic invariants in the growing fields of the proton cyclotron waves. The resulting electron anisotropy is sufficiently strong to trigger the growth of whistler waves.

While the presence of whistler waves is hard to identify by eye in the spacetime diagrams of Figures 10(b) and (c), due to the dominance of proton cyclotron modes, we can extract their signature by applying a filter in frequency and wavenumber, as done in Section 4.2.1. Figures 10(e) and (f) show the power spectra of $\delta {B}_{x}$ and $\delta {B}_{z}$. Most of the power is concentrated near the origin at low frequencies and long wavelengths, associated with the proton cyclotron mode. However, we can still identify a significant peak around $\omega \simeq 13\,{{\rm{\Omega }}}_{{ci}}\simeq 0.3\,{{\rm{\Omega }}}_{{ce}}$ and ${k}_{y}\simeq 0.3\,{\omega }_{{pe}}/c$. In analogy with the discussion in Section 4.2.1, we associate this peak with electron whistler waves. When applying a high-pass filter whose frequency and wavelength cuts are shown as red dashed lines in Figures 10(e) and (f), we recover in the spacetime diagrams of Figures 10(h) and (i) the typical spatial and temporal patterns of electron whistler waves. As expected, most of the electron whistler activity takes place near the end of the exponential growth of proton waves, at $5\lesssim {{\rm{\Omega }}}_{{ci}}t\lesssim 7$ (see also the temporal evolution of the energy in whistler waves in Figure 9(b)). In this time interval, the electron anisotropy exceeds the threshold of whistler instability in the whole simulation domain (Figure 10(g)).

This period also corresponds to a rapid increase of the electron specific entropy, as measured directly from the electron distribution function (blue solid line in Figure 9(f)). This is expected, since whistler waves provide the pitch-angle scattering required to break adiabatic invariance, which (together with the sustained electron anisotropy; see Figures 9(d) and 10(d) at $5\lesssim {{\rm{\Omega }}}_{{ci}}t\lesssim 7$) drives efficient entropy generation. Based on our model in Section 2, the electron specific entropy should increase as

$$\begin{eqnarray}&&{{ds}}_{e}=\left[\displaystyle \frac{1}{2}d\mathrm{ln}\left(\displaystyle \frac{{T}_{e,\parallel }}{{(n/B)}^{2}}\right)\right]\left(1-\displaystyle \frac{{T}_{e,\parallel }}{{T}_{e,\perp }}\right)-\displaystyle \frac{{{de}}_{w,e}}{{T}_{e,\perp }},\end{eqnarray} \tag{ 42 }$$

which follows from Equation (13) (an equivalent form can be obtained from Equation (14)). Here, we have used the condition ${{de}}_{w,e\mathrm{tot}}={{de}}_{w,e}$ for electrons. The comparison of the measured entropy increase (blue solid line in Figure 9(f)) with the predicted entropy change (orange dotted line in Figure 9(f)) provides another validation of our heating model.

While most of the electron entropy production happens near the end of the exponential growth of proton waves, a moderate increase of the electron entropy also occurs during the secular stage (i.e., at ${{\rm{\Omega }}}_{{ci}}t\gtrsim 10$). Here, the oscillating cyclotron fluctuations are sloshing electrons around and can occasionally excite local patches of electron anisotropy that exceed the whistler threshold (e.g., see Figure 10(g) at ${{\rm{\Omega }}}_{{ci}}t\simeq 11$ and $y\simeq 175\,\,c/{\omega }_{{pe}}$). In the same region, we can identify a short episode of electron whistler activity, particularly in $\delta {B}_{z}$ in Figure 10(i). These sporadic episodes of anisotropy-driven whistler waves further increase the electron entropy.

It is worth noticing, though, that at late times the box-averaged electron anisotropy switches sign, with ${T}_{e,\parallel }\gtrsim {T}_{e,\perp }$ (in Figure 9(d), at ${{\rm{\Omega }}}_{{ci}}t\gtrsim 9$), so the opportunities for whistler growth are fewer. This behavior is consistent with the conservation of the first and second adiabatic invariants in the decaying field of the proton cyclotron waves, leading to an increase in ${T}_{e,\parallel }$ and a decrease in ${T}_{e,\perp }$ (as indeed seen in Figure 9(d) at late times). The same "inverted" anisotropy with ${T}_{e,\parallel }\gtrsim {T}_{e,\perp }$ is seen in the far downstream of our reference shock run (Figure 3(f)), accompanying the decay of the proton modes. We remark that electron entropy production is still possible when ${T}_{e,\parallel }\gtrsim {T}_{e,\perp }$, as long as the anisotropy is large enough to exceed the threshold of the firehose instability, which would then provide the mechanism for breaking adiabatic invariance. We have verified with expanding box simulations similar to the ones reported in Section 5 (not shown) that once the system exceeds the firehose threshold, the electron entropy rapidly increases.

Finally, by measuring directly the energy in whistler waves, we can also validate the electron energy equation

$$\begin{eqnarray}&&{{du}}_{e}={{dw}}_{e,\perp }+{{dw}}_{e,\parallel }-{{de}}_{w,e}.\end{eqnarray} \tag{ 43 }$$

Once again, the time integral of the left-hand side matches very well the time integral of the right-hand side (compare green and red curves in Figure 9(h)); that is, the change of electron internal energy can be accounted for by the total work done by the proton waves and the energy lost to generate electron whistler waves.

In summary, we have demonstrated that efficient electron entropy production can occur even in the absence of a large-scale compression. Magnetic and density fluctuations sourced by anisotropic protons drive electrons to become anisotropic, with ${T}_{e,\perp }\gt {T}_{e,\parallel }$. The electron anisotropy is relaxed by the growth of whistler waves, which break the electron adiabatic invariance and mediate the production of electron entropy. In the next subsection, we show that the resulting electron entropy increase is nearly independent of the proton-to-electron mass ratio.

In this subsection, we explore with periodic boxes initialized with anisotropic protons how the development of the proton cyclotron instability can lead to electron irreversible heating. We vary the mass ratio from 49 up to 1600, as indicated in the legend of Figure 11(a). Since the fastest growing mode of the proton cyclotron instability has wave vector $\sim {\omega }_{{pi}}/c$, we increase the number of cells in our computational domain as $\propto \sqrt{{m}_{i}/{m}_{e}}$, to include the same number of proton skin depths for all values of ${m}_{i}/{m}_{e}$ (and so, approximately the same number of proton cyclotron wavelengths).

Figure 11 compares the runs. Panel (a) shows the time evolution of the wave magnetic energy, which is dominated by proton cyclotron modes. Panel (b) shows the evolution of the proton temperature anisotropy, which reduces strongly at ${{\rm{\Omega }}}_{{ci}}t\gtrsim 4$ by pitch-angle scattering off the strong cyclotron waves. Unsurprisingly, since these two quantities are related to protons, their evolution is almost identical for different mass ratios. As long as the mass ratio is sufficiently large to adequately separate electron and proton scales, the proton cyclotron instability—whose polarization is resonant with protons, but nonresonant with electrons—is not affected by electron physics.

Zoom In Zoom Out Reset image size

As in the case ${m}_{i}/{m}_{e}=49$ discussed above, the growth of the proton cyclotron instability induces an electron temperature anisotropy with ${T}_{e,\perp }\gt {T}_{e,\parallel }$ (Figure 11(d)), which excites electron whistler waves (Figure 11(c))⁸ and facilitates electron entropy increase (Figure 11(f)) by violating the electron adiabatic invariance (Figure 11(e)). The dependence of the peak electron anisotropy on mass ratio for ${m}_{i}/{m}_{e}\lesssim 200$ can be understood from the same argument we have presented in Section 5: in units of the electron gyration period, the growth of proton waves (or the compressed field, for Section 5) is faster at lower ${m}_{i}/{m}_{e}$, which leads to an overshoot in electron anisotropy beyond the threshold of whistler marginal stability. The overshoot is more pronounced for lower ${m}_{i}/{m}_{e}$. Due to the higher electron anisotropy, more free energy is available for the growth of whistler waves at lower ${m}_{i}/{m}_{e}$ (see the trend from purple to green in Figure 11(c) at ${{\rm{\Omega }}}_{{ci}}t\sim 6$). Because of the higher electron anisotropy and stronger whistler waves, the electron entropy increases slightly more at lower mass ratios (in particular, see the purple line in Figure 11(f) for ${m}_{i}/{m}_{e}=49$).

On the other hand, the electron physics shows no appreciable dependence on mass ratio for ${m}_{i}/{m}_{e}\gtrsim 400$. The peak electron anisotropy at $5\lesssim {{\rm{\Omega }}}_{{ci}}t\sim 7$ saturates at the threshold of whistler marginal stability (indicated in Figure 11(c) by the dotted lines, with the same color coding as the solid lines). As a consequence, the peak strength of whistler waves is nearly independent of mass ratio (see Figure 11(b) in the same time interval), and the resulting entropy increase is the same for all mass ratios ${m}_{i}/{m}_{e}\gtrsim 400$ (Figure 11(f)). Even for ${m}_{i}/{m}_{e}=49$, the electron entropy increase at the final time is only $\sim 30 \%$ higher than for ${m}_{i}/{m}_{e}=1600$.

In the periodic box runs of Sections 5 and 6, we have extended our study to realistic mass ratios, showing that the entropy increase at ${m}_{i}/{m}_{e}=1600$ is only $\sim 30 \%$ smaller than for the choice ${m}_{i}/{m}_{e}=49$ of our reference shock simulation. In Figure 13, we investigate the dependence of the electron physics in our full shock simulations on the mass ratio, from ${m}_{i}/{m}_{e}=25$ up to ${m}_{i}/{m}_{e}=200$ (as indicated in the legend of panel (d)). We typically employ 32 computational particles per cell, with the exception of ${m}_{i}/{m}_{e}=200$, where we use 64 particles per cell to keep numerical heating under control. We keep the upstream electron temperature fixed at ${k}_{{\rm{B}}}{T}_{e0}={10}^{-2}\,{m}_{e}{c}^{2}$ so that electrons stay safely nonrelativistic. This implies that the plasma inflow velocity is slower with increasing mass ratio, as $\propto \sqrt{{m}_{e}/{m}_{i}}$.

Zoom In Zoom Out Reset image size

The proton physics is expected to be the same regardless of mass ratio, and in fact the profiles of density (panel (a)), proton temperature (panel (c)), and proton anisotropy (panel (d)) are nearly the same for all mass ratios. The same holds for the wave magnetic energy at $(x-{x}_{\mathrm{sh}})/{r}_{\mathrm{Li}}\lesssim -0.5$ (panel (b)), where proton-driven modes dominate (see Section 4 for details).

On the other hand, the peak electron anisotropy at the shock (panel (e)) is systematically lower for higher mass ratios, in perfect agreement with the trend observed in the periodic box experiments of Figure 8. Despite the pronounced difference in peak anisotropy, Figure 8(d) showed that the entropy increase until $t\sim {{\rm{\Omega }}}_{{ci}}^{-1}$ was only marginally lower at higher ${m}_{i}/{m}_{e}$. This trend (and the weak mass ratio dependence) is confirmed by the profiles of electron entropy in the shock ramp shown in Figure 13(h). Overall, Figure 13(h) confirms the results of our periodic box experiments, namely, the electron entropy increase is nearly independent of mass ratio (with the exception of the lowest mass ratio ${m}_{i}/{m}_{e}=25$ presented in Figure 13). Even though our shock simulations only extend up to ${m}_{i}/{m}_{e}=200$, the results of our periodic runs in Sections 5 and 6 suggest that the same conclusion should hold up to the realistic mass ratio.