Electron Heating in Perpendicular Low-beta Shocks

Electron heating in collisionless shocks—stated as post-shock electron/ion temperature ratio ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$—is not constrained by magnetohydrodynamic (MHD) shock jump conditions. How much do electrons heat, and how do they heat? A prediction for ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ can constrain models for gas accretion onto galaxy clusters (Avestruz et al. 2015) and cosmic-ray acceleration in supernova remnants (Helder et al. 2010; Yamaguchi et al. 2014; Hovey et al. 2018). Detailed study of the electron heating physics can also help us interpret new high-resolution data from the Magnetospheric Multiscale Mission (Chen et al. 2018; Goodrich et al. 2018; Cohen et al. 2019).

In the heliosphere, shocks of magnetosonic Mach number ${{ \mathcal M }}_{\mathrm{ms}}\gtrsim 2\mbox{--}3$ heat electrons beyond adiabatic compression via a two-step process: electrons accelerate in bulk along ${\boldsymbol{B}}$ toward the shock downstream, then relax into "flat-top" distributions in ${\boldsymbol{B}}$-parallel velocity (Feldman et al. 1982, 1983; Chen et al. 2018). Two mechanisms—quasi-static direct current (DC) fields and plasma waves—may drive ${\boldsymbol{B}}$-parallel acceleration. In the DC mechanism, an electric potential jump in the shock layer (i.e., a quasi-static electric field that points along shock normal) accelerates electrons in bulk (Feldman et al. 1983; Goodrich & Scudder 1984; Scudder et al. 1986; Scudder 1996; Hull et al. 2001; Lefebvre et al. 2007; Schwartz 2014). The DC electron energy gain scales with ${\cos }^{2}\theta$, where θ is the angle between ${\boldsymbol{B}}$ and shock normal (Goodrich & Scudder 1984). We expect no heating in exactly planar perpendicular shocks, but shock rippling from ion-scale waves (Lowe & Burgess 2003; Johlander et al. 2016; Hanson et al. 2019) can bend ${\boldsymbol{B}}$, alter θ, and enable DC heating. Plasma waves with nonzero ${E}_{\parallel }\,={\boldsymbol{E}}\cdot \hat{{\boldsymbol{b}}}$, such as oblique whistlers, can also provide electron bulk acceleration and thus heating (Wilson et al. 2014a, 2014b). Such plasma waves are intrinsic to shock structure (Krasnoselskikh et al. 2002; Wilson et al. 2009, 2012; Dimmock et al. 2019) and may be sustained by free energy from, e.g., shock-reflected ions (Wu et al. 1984; Matsukiyo & Scholer 2006; Muschietti & Lembège 2017).

In this Letter, we study thermal electron heating in multi-dimensional particle-in-cell (PIC) simulations of perpendicular shocks with realistic structure (requiring high ion/electron mass ratio ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}$; Krauss-Varban et al. 1995; Umeda et al. 2012a, 2014) and high grid resolution to resolve electron scattering and relaxation after ${\boldsymbol{B}}$-parallel bulk acceleration.

We simulate collisionless 2D (x–y) ion–electron shocks using the relativistic PIC code TRISTAN-MP (Buneman 1993; Spitkovsky 2005). We inject plasma with velocity $-{u}_{0}\hat{{\boldsymbol{x}}}$ and magnetic field ${B}_{0}\hat{{\boldsymbol{y}}}$ from the simulation domain's right-side (upstream) boundary. Injected plasma reflects from a conducting wall at x = 0, forming a shock that travels toward $+\hat{{\boldsymbol{x}}}$. The shocked downstream plasma has zero bulk velocity, and the upstream ${\boldsymbol{B}}$ is perpendicular to the shock normal, so θ = 90°. The simulation domain expands along $+\hat{{\boldsymbol{x}}}$ to keep the right-side boundary $\gtrsim 1.5\,{r}_{\mathrm{Li}}$ ahead of the shock front (Sironi & Spitkovsky 2009, Section 2), where ${r}_{\mathrm{Li}}={u}_{0}/{{\rm{\Omega }}}_{{\rm{i}}}$ is a characteristic ion Larmor radius; we checked that shock heating physics is not artificially affected by the right-side boundary. Upstream ions and electrons have equal density n₀ and temperature T₀. The plasma frequencies ${\omega }_{{\rm{p}}\{{\rm{i}},{\rm{e}}\}}=\sqrt{4\pi {n}_{0}{e}^{2}/{m}_{\{{\rm{i}},{\rm{e}}\}}}$ and cyclotron frequencies ${{\rm{\Omega }}}_{\{{\rm{i}},{\rm{e}}\}}={{eB}}_{0}/({m}_{\{{\rm{i}},{\rm{e}}\}}c)$, where subscripts ${\rm{i}}$ and ${\rm{e}}$ denote ions and electrons. We use Gaussian CGS units throughout.

Our fiducial simulations have ion/electron mass ratio ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}\,=625$ and total plasma beta ${\beta }_{{\rm{p}}}=16\pi {n}_{0}{k}_{{\rm{B}}}{T}_{0}/{{B}_{0}}^{2}=0.25$. The fast magnetosonic, sonic, and Alfvén Mach numbers are ${{ \mathcal M }}_{\mathrm{ms}}\,={u}_{\mathrm{sh}}/\sqrt{{{c}_{{\rm{s}}}}^{2}+{{v}_{{\rm{A}}}}^{2}}=1\mbox{--}10$, ${{ \mathcal M }}_{{\rm{s}}}={u}_{\mathrm{sh}}/{c}_{{\rm{s}}}=3\mbox{--}20$, and ${{ \mathcal M }}_{{\rm{A}}}\,={u}_{\mathrm{sh}}/{v}_{{\rm{A}}}=1.5\mbox{--}10$. The sound speed ${c}_{{\rm{s}}}=\sqrt{2{\rm{\Gamma }}{k}_{{\rm{B}}}{T}_{0}/({m}_{{\rm{i}}}+{m}_{{\rm{e}}})}$, Alfvén speed ${v}_{{\rm{A}}}={B}_{0}/\sqrt{4\pi {n}_{0}({m}_{{\rm{i}}}+{m}_{{\rm{e}}})}$, and ${u}_{\mathrm{sh}}$ is the speed of upstream plasma in the shock's rest frame; for nonrelativistic speeds, ${u}_{\mathrm{sh}}={u}_{0}/(1-1/r)$, where $r\leqslant 4$ is the MHD shock-compression ratio. The one-fluid adiabatic index Γ is not known a priori, but it is set self-consistently by the degree of ion and electron isotropization. We report Mach numbers assuming Γ = 2, which overestimates ${{ \mathcal M }}_{\mathrm{ms}}$ by $\sim 1\mbox{--}10 \%$ for stronger shocks that isotropize ions/electrons and have ${\rm{\Gamma }}\approx 5/3$.

The grid cell size ${\rm{\Delta }}x={\rm{\Delta }}y=0.1c/{\omega }_{\mathrm{pe}}$ and the time step ${\rm{\Delta }}t=0.045{\omega }_{\mathrm{pe}}^{-1}$ so that $c=0.45{\rm{\Delta }}x/{\rm{\Delta }}t$. Upstream plasma has 16 particles per cell per species. We smooth the electric current with 32 sweeps of a three-point binomial ("1–2–1") filter at each time step (Birdsall & Langdon 1991, Appendix C). The N = 32 sweeps approximate a Gaussian filter with standard deviation $\sqrt{N/2}=4$ cells or $0.4c/{\omega }_{\mathrm{pe}}$. The filter's half-power cutoff is at wavenumber $k\approx \sqrt{2/N}{({\rm{\Delta }}x)}^{-1}=2.5{\left(c/{\omega }_{\mathrm{pe}}\right)}^{-1}$, which implies 50% damping at wavelength $\lambda \approx \pi \sqrt{2N}{\rm{\Delta }}{\rm{x}}=2.5\left(c/{\omega }_{{\rm{pe}}}\right)$. Electron-scale waves may be damped, but we will later show that electron-scale waves mainly scatter rather than heat. We simulated 2D ${{ \mathcal M }}_{\mathrm{ms}}=3.1$, ${\beta }_{{\rm{p}}}=0.25$ shocks with 4× larger or smaller sweep number N; the ratio ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ did not change much. We adjust T₀, u₀, and B₀ to control ${{ \mathcal M }}_{\mathrm{ms}}$ and ${\beta }_{{\rm{p}}}$ while keeping shocked electrons nonrelativistic, i.e., post-shock ${k}_{{\rm{B}}}{T}_{{\rm{e}}}\,\lesssim 0.05{m}_{{\rm{e}}}{c}^{2}$. The ratio $\tau ={\omega }_{\mathrm{pe}}/{{\rm{\Omega }}}_{{\rm{e}}}=2.5\mbox{--}11$ (τ ≫ 1 for solar wind and astrophysical settings). The transverse (y) width is 2.9–$5.8\,c/{\omega }_{\mathrm{pi}}\,=720\mbox{--}1440$ cells. Simulation durations are 10–$20\,{{\rm{\Omega }}}_{{\rm{i}}}^{-1}\,=932$–$2736\,{\omega }_{\mathrm{pi}}^{-1}$ so that post-shock ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ reaches steady state. The temperatures ${T}_{\{{\rm{i}},{\rm{e}}\}}$, ${T}_{\{{\rm{i}},{\rm{e}}\}\parallel }$, and ${T}_{\{{\rm{i}},{\rm{e}}\}\perp }$ are moments of the particle distribution in a 5^N cell region, where $N\in \{1,2,3\}$ is the domain dimensionality. The comoving frame boost for moment calculation uses a fluid velocity also averaged over 5^N cells. All ∥- and ⊥-subscripted quantities are taken with respect to local ${\boldsymbol{B}}$.

Figure 1 shows a representative simulation. Ions transmit or reflect at the shock ramp (Figures 1(a)–(b)). The shock front is rippled (Figure 1(c)). A net E_x potential exists across the shock, and E_x is also modulated by the ${\boldsymbol{B}}$ rippling wavelength (Figures 1(b)–(c)). Reflected ions accelerate in the motional field ${E}_{z}={u}_{0}{B}_{0}/c$ before reentering the shock (Leroy et al. 1982); this lowers ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ in the shock foot (Figure 1(d)). Electrons heat above adiabatic expectation in the shock ramp and settle to ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}\approx 0.4;$ no appreciable heating occurs after the shock ramp (Figure 1(d)). The fluid adiabatic prediction in Figure 1(d) is ${T}_{{\rm{e}},\mathrm{ad}}/({T}_{{\rm{i}}}+{T}_{{\rm{e}}}-{T}_{{\rm{e}},\mathrm{ad}})$, using measured ${T}_{{\rm{i}}}$ and ${T}_{{\rm{e}}}$ and assuming ${T}_{{\rm{e}},\mathrm{ad}}={T}_{0}[1+2{(n/{n}_{0})}^{{\rm{\Gamma }}-1}]/3$ with Γ = 2.

We measure post-shock ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ (Figure 1(c)) as a function of ${{ \mathcal M }}_{\mathrm{ms}}$ for many simulations with varying dimensionality, magnetic field orientation θ, ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}$, and ${\beta }_{{\rm{p}}}$. We also adjust domain width, particle resolution, and current smoothing to control noise and computing cost. In simulations with $\theta \lt 90^\circ$, the right-side boundary expands at $\max \left({u}_{\mathrm{sh}},0.5c\cos \theta \right)$ to retain shock-reflected electrons streaming along ${\boldsymbol{B}}$.

We show the post-shock ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ for our fiducial 2D ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}=625$ shocks with in-plane upstream magnetic field ${B}_{0}\hat{{\boldsymbol{y}}}$ in Figure 2(a). These fiducial simulations are converged in ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ with respect to transverse (y) width. For perpendicular shocks, we find that electron heating beyond adiabatic compression requires 2D geometry with in-plane ${\boldsymbol{B}}$. Corresponding 2D simulations with out-of-plane ${\boldsymbol{B}}$ (along $\hat{{\boldsymbol{z}}}$) and 1D simulations heat electrons by compression alone (Figure 2(a)). At ${{ \mathcal M }}_{\mathrm{ms}}\sim 5\mbox{--}10$, the 2D simulations with out-of-plane ${\boldsymbol{B}}$ and 1D simulations show weak super-adiabatic heating in the shock layer, but the ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ measurement is also less precise due to numerical heating. Shimada & Hoshino (2000, 2005) saw strong electron heating in 1D perpendicular ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}=20$ shocks due to Buneman instability between shock-reflected ions and incoming electrons. The higher ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}=625$ suppresses the Buneman instability in our 1D shocks.

Zoom In Zoom Out Reset image size

Can DC heating in a 2D rippled shock—i.e., varying local magnetic field angle due to self-generated waves—explain the super-adiabatic electron heating seen in our fiducial 2D simulations? To estimate the DC heating from varying θ, we perform 1D oblique shock simulations with varying $\theta \lt 90^\circ$ (Figure 2(b)); recall that θ is the angle between ${\boldsymbol{B}}$ and shock normal. The 1D setup keeps quasi-static shock structure (averaged over shock reformation cycles) and should retain DC heating while excluding waves oblique to the shock normal. We do find super-adiabatic heating in 1D oblique shocks. Electrons heat more for lower θ, which is qualitatively consistent with DC field heating (Goodrich & Scudder 1984). For our representative ${{ \mathcal M }}_{\mathrm{ms}}=3.1$ shock, which has local ripple $\theta \gtrsim 80^\circ$ (Figure 4(f)), the DC heating inferred from 1D oblique shock simulations appears too low to explain the full amount of super-adiabatic heating (Figure 2(b), box).

Our fiducial 2D perpendicular shocks appear converged in mass ratio at ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}\,\sim$ 200−625 (Figure 2(c)), consistent with prior simulations (Umeda et al. 2012a, 2014) and theory (Krauss-Varban et al. 1995). For ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}=20\mbox{--}625$, 2D shocks agree on ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ to within a factor of 2–3. A set of 3D ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}=49$ simulations with narrower transverse width, $2.7c/{\omega }_{\mathrm{pi}}$, shows good agreement too. Agreement between 2D and 3D for ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}=49$ suggests that 2D simulations with in-plane ${\boldsymbol{B}}$ include the essential physics for electron heating.

To see how heating depends on ${\beta }_{{\rm{p}}}$, we reduce ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}$ to 49 and sweep ${\beta }_{{\rm{p}}}$ over 0.125–4 (Figure 2(d)). Electron heating increases above adiabatic at ${{ \mathcal M }}_{\mathrm{ms}}\sim 2\mbox{--}3$ for all ${\beta }_{{\rm{p}}}$. At ${{ \mathcal M }}_{\mathrm{ms}}\sim 3\mbox{--}5$ and ${\beta }_{{\rm{p}}}\leqslant 1$, two-step ${\boldsymbol{B}}$-parallel electron heating (which we describe below) operates for all ${\beta }_{{\rm{p}}}\lesssim 1$. At ${{ \mathcal M }}_{\mathrm{ms}}\sim 3\mbox{--}5$ and ${\beta }_{{\rm{p}}}\gtrsim 2$, a distinct electron cyclotron whistler instability is expected to heat electrons instead (Guo et al. 2017, 2018). At ${{ \mathcal M }}_{\mathrm{ms}}\gtrsim 5$, shock structure is more complex, which we do not explore here. The relationship between ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ and ${{ \mathcal M }}_{\mathrm{ms}}$ does not appear to depend on ${\beta }_{{\rm{p}}}$ for ${{ \mathcal M }}_{\mathrm{ms}}\gtrsim 4$.

Our fiducial ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}-{{ \mathcal M }}_{\mathrm{ms}}$ data are order-of-magnitude consistent with measurements from solar wind bow shocks (Figure 2(e)), replotted from Ghavamian et al. (2013). The Saturn data are uncertain in both ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ and ${{ \mathcal M }}_{\mathrm{ms}}$ due to a lack of ion temperature measurements from Cassini (Masters et al. 2011), so ${{ \mathcal M }}_{\mathrm{ms}}=0.671{{ \mathcal M }}_{{\rm{A}}}$ (equivalent to ${\beta }_{{\rm{p}}}\sim 1.5$) is assumed following Ghavamian et al. (2013); we note ${\beta }_{{\rm{p}}}\,\sim 1.5$ is a typical value (Richardson 2002). The Earth data have ${\beta }_{{\rm{p}}}\sim 0.1\mbox{--}1$ and use directly measured ion and electron temperatures from the ISEE spacecraft (Schwartz et al. 1988). Both data sets are mostly quasi-perpendicular, with a majority of shocks having $50^\circ \lt \theta \lt 90^\circ$ (Schwartz et al. 1988; Masters et al. 2011).

For further study, we choose the weakest 2D ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}=625$, ${\beta }_{{\rm{p}}}=0.25$ shock with significant super-adiabatic heating: our representative ${{ \mathcal M }}_{\mathrm{ms}}=3.1$ simulation (Figures 1 and 2(a)). We redo this simulation with higher resolution: ${\rm{\Delta }}x={\rm{\Delta }}y\,=0.05\,c/{\omega }_{\mathrm{pe}}$ (keeping $c=0.45{\rm{\Delta }}x/{\rm{\Delta }}t$), 64 particles per cell per species, and 64 current filter passes per time step. The current filter approximates a Gaussian with standard deviation $\sim 5.7$ cells or $0.28c/{\omega }_{\mathrm{pe}};$ the filter's half-power cutoff is at wavenumber $k\approx 3.5{\left(c/{\omega }_{\mathrm{pe}}\right)}^{-1}$, which means 50% damping at wavelength $\lambda \approx 1.8c/{\omega }_{\mathrm{pe}}$. We then select 15,898 electron particles between $x=8.00\mbox{--}8.02c/{\omega }_{\mathrm{pi}}$ at $t^{\prime} \equiv t-324\,{\omega }_{\mathrm{pi}}^{-1}=0$, located $4\,c/{\omega }_{\mathrm{pi}}\,=2\,{r}_{\mathrm{Li}}$ ahead of the shock ramp, and monitor their phase-space evolution (Figure 3) and energy gain (Figure 4) through the shock. The perpendicular upstream ${\boldsymbol{B}}$ confines particles within a narrow magnetic flux tube and prevents particle drift from downstream to upstream.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Elongated, ion-scale E_∥ waves accelerate electrons along ${\boldsymbol{B}}$ in the shock foot and ramp. These waves have $| {E}_{\parallel }| /({u}_{0}{B}_{0}/c)\,\sim 0.2\mbox{--}0.6$ and wavelength ${\lambda }_{y}\sim 2c/{\omega }_{\mathrm{pi}}\sim {r}_{\mathrm{Li}}$ (Figure 3(a), arrow); we attribute this E_∥ to very oblique whistler waves (i.e., magnetosonic/lower hybrid branch) with fluctuating ${E}_{x}\gg {E}_{y},{E}_{z}$ and ${B}_{y},{B}_{z}\gt {B}_{x}$, as identified by prior PIC studies (Matsukiyo & Scholer 2003, 2006; Hellinger et al. 2007; Umeda et al. 2012b). A stronger bipolar ion-scale $| {E}_{\parallel }| /({u}_{0}{B}_{0}/c)\gtrsim 0.5$ (Figure 3(b), arrow) straddles clumps of shock-reflected ions and also accelerates electrons. Accelerated electrons appear as coherent deflections in $\gamma {\beta }_{\parallel }$–y phase space (Figures 3(g)–(h), arrows) that disrupt and relax via two-stream-like instability. Local y-regions develop asymmetric and transiently unstable $\gamma {\beta }_{\parallel }$ distributions (Figures 3(q), (v), (x)). Electron relaxation generates strong and rapid electron-scale E_∥ waves and phase-space holes with ${\lambda }_{y}\sim c/{\omega }_{\mathrm{pe}}$ (Figures 3(b), (g), (i), boxes; compare An et al. 2019) Landau damping is evidenced by flattened distributions at $\gamma {\beta }_{\parallel }\sim 0.2$ (Figures 3(k)–(l)). Electrons relax to near isotropy by $t^{\prime} \sim 140\,{\omega }_{\mathrm{pi}}^{-1}$ (Figures 3(j), (o), (t), (y)). Prior 2D PIC simulations have shown similar two-step ${\boldsymbol{B}}$-parallel heating in a shock foot setup (periodic interpenetrating beams; Matsukiyo & Scholer 2006) and in full shocks (Umeda et al. 2011, 2012b).

E_∥ is the main source of non-adiabatic electron heating (Figure 4(a)). We decompose the sample electrons' mean energy gain W_tot into parallel and perpendicular work, ${W}_{\parallel }\,=-e\langle \int {E}_{\parallel }{v}_{\parallel }{dt}\rangle$ and ${W}_{\perp }=-e\langle \int {E}_{\perp }{v}_{\perp }{dt}\rangle$, integrated for every particle over every code time step Δt such that ${W}_{\mathrm{tot}}={W}_{\parallel }+{W}_{\perp }$. Angle brackets are particle averages. We estimate the adiabatic heating as ${W}_{\perp ,\mathrm{ad}}\,=\langle {\sum }_{n}({\gamma }_{n\to n+1,\mathrm{ad}}-{\gamma }_{n}){m}_{{\rm{e}}}{c}^{2}\rangle$, where

$$\begin{eqnarray}&&{\gamma }_{n\to n+1,\mathrm{ad}}=\sqrt{1+{\left(\gamma {\beta }_{\parallel }\right)}_{n}^{2}+{\left(\gamma {\beta }_{\perp }\right)}_{n}^{2}\left({B}_{n+1}/{B}_{n}\right)}\end{eqnarray} \tag{ 1 }$$

captures electron heating from compression between time steps n and n+1.

In Equation (1), we assume $(\gamma {\beta }_{\perp }{)}_{n}^{2}/{B}_{n}$ and ${(\gamma {\beta }_{\parallel })}_{n}$ are constant during compression; γ, ${\beta }_{\parallel }$, and ${\beta }_{\perp }$ are evaluated in the electron fluid's rest frame. The sum in ${W}_{\perp ,\mathrm{ad}}$ uses a coarse output time step ${\rm{\Delta }}{t}_{\mathrm{out}}=400\,{\rm{\Delta }}t=9\,{\omega }_{\mathrm{pe}}^{-1}$.

The non-adiabatic perpendicular work can be explained by high-frequency scattering of electron parallel energy into perpendicular energy (Figure 4(b)). We Fourier-space filter E_x, E_y, and E_z to isolate wavenumbers ${k}_{y}\gt {\left(c/{\omega }_{\mathrm{pe}}\right)}^{-1}$ and then construct ${E}_{\perp ,\mathrm{HF}}$ and ${E}_{\parallel ,\mathrm{HF}}$ by projecting the Fourier-filtered fields onto local ${\boldsymbol{B}}$. Then, ${W}_{\perp ,\mathrm{HF}}=-e\langle \sum {E}_{\perp ,\mathrm{HF}}{v}_{\perp }{\rm{\Delta }}{t}_{\mathrm{out}}\rangle$ and ${W}_{\parallel ,\mathrm{HF}}=-e\langle \sum {E}_{\parallel ,\mathrm{HF}}{v}_{\parallel }{\rm{\Delta }}{t}_{\mathrm{out}}\rangle$. We find that ${W}_{\perp ,\mathrm{HF}}$ and ${W}_{\perp }-{W}_{\perp ,\mathrm{ad}}$ agree to $\sim 10 \%$, suggesting that non-adiabatic ${W}_{\perp }$ comes from electron-scale scattering of parallel energy. Exact agreement is not expected due to the coarse time step Δt_out and the arbitrary k_y cut.

The particle-averaged adiabatic moment ${{p}_{\perp }}^{2}/B$ grows in steps that correlate with increases in W_tot and ${W}_{\parallel }$. Bulk acceleration at $t^{\prime} =84.6\,{\omega }_{\mathrm{pi}}^{-1}$ and $t^{\prime} =107.3\,{\omega }_{\mathrm{pi}}^{-1}$ coincides with momentarily constant ${{p}_{\perp }}^{2}/B$ and increasing ${W}_{\parallel }$ prior to a scattering episode. Then, ${{p}_{\perp }}^{2}/B$ increases during strong electron scattering at $t^{\prime} =89.6\,{\omega }_{\mathrm{pi}}^{-1}$ and $112.3\,{\omega }_{\mathrm{pi}}^{-1}$ while ${W}_{\parallel }$ flattens off (Figures 3 and 4(a), (c)).

Figure 4(e) shows mean work from grad B, inertial, and polarization drifts, as well as $\partial B/\partial t$ induction work; see Northrop (1961, 1963), Goodrich & Scudder (1984), Dahlin et al. (2014), and Rowan et al. (2019). Each ${W}_{\mathrm{drift}}=-e\langle \sum {\boldsymbol{E}}\cdot {{\boldsymbol{v}}}_{\mathrm{drift}}{\rm{\Delta }}{t}_{\mathrm{out}}\rangle$, where ${{\boldsymbol{v}}}_{\mathrm{drift}}$ is one of

$$\begin{eqnarray*}&&{{\boldsymbol{v}}}_{\{{\rm{\nabla }}B,\mathrm{inert},\mathrm{pol}\}}=-\displaystyle \frac{\gamma {m}_{{\rm{e}}}c}{{eB}}\hat{{\boldsymbol{b}}}\times \left\{\displaystyle \frac{{{v}_{\perp }}^{2}}{2B}{\rm{\nabla }}B,\,{v}_{\parallel }\displaystyle \frac{d\hat{{\boldsymbol{b}}}}{{dt}},\,\displaystyle \frac{d{{\boldsymbol{v}}}_{E}}{{dt}}\right\},\end{eqnarray*}$$

with γ and ${v}_{\perp }$ evaluated in the electron fluid's rest frame. We take ${{\boldsymbol{v}}}_{E}=\langle c{\boldsymbol{E}}\times {\boldsymbol{B}}/{B}^{2}\rangle$ to reduce noise; otherwise, the ${{\boldsymbol{v}}}_{\mathrm{drift}}$ terms use ${\boldsymbol{E}}$ and ${\boldsymbol{B}}$ fields seen by individual particles. The $d/{dt}$ terms are one-sided finite differences, e.g., $d{{\boldsymbol{v}}}_{E}/{dt}=[{({{\boldsymbol{v}}}_{E})}_{n+1}-{({{\boldsymbol{v}}}_{E})}_{n}]/{\rm{\Delta }}{t}_{\mathrm{out}}$. And, ${W}_{\mathrm{induct}}=\gamma {m}_{{\rm{e}}}{v}_{\perp }^{2}(\partial B/\partial t)/(2B)$, with $\partial B/\partial t\,=[{B}_{n+1}({{\boldsymbol{r}}}_{n})-{B}_{n-1}({{\boldsymbol{r}}}_{n})]/(2{\rm{\Delta }}{t}_{\mathrm{out}})$ and ${{\boldsymbol{r}}}_{n}$ the particle position at time step n. We find that grad B drift and induction together give fluid-like adiabatic compression. Inertial and polarization drifts give less work, but some other electron samples have W_inert comparable to ${W}_{{\rm{\nabla }}B}$ (Appendix A). We compare the summed drifts to ${W}_{\mathrm{tot},\mathrm{out}}=-e\langle \sum {\boldsymbol{E}}\cdot {\boldsymbol{v}}{\rm{\Delta }}{t}_{\mathrm{out}}\rangle$. We conclude that W_tot agrees with both the summed drift work and W_tot,out, given uncertainty from both the guiding-center drift approximation and the coarse integration time step.

Earlier, we argued that DC heating alone may not explain all super-adiabatic heating in our fiducial 2D shock, based on downstream volume-averaged ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ (Figures 1(c) and 2(b)). So, Figure 4(f) estimates DC-like work as ${\bar{W}}_{\mathrm{drift}}=-e\sum \langle {\boldsymbol{E}}\rangle \,\cdot \langle {{\boldsymbol{v}}}_{\mathrm{drift}}\rangle {\rm{\Delta }}{t}_{\mathrm{out}}$ and ${\bar{W}}_{\parallel }=-e\sum \langle {E}_{\parallel }\rangle \langle {v}_{\parallel }\rangle {\rm{\Delta }}{t}_{\mathrm{out}}$. The ${\boldsymbol{E}}$ average removes waves along $\hat{{\boldsymbol{y}}}$ to keep only 1D-like shock fields. The ${{\boldsymbol{v}}}_{\mathrm{drift}}$ average gives a mean drift trajectory and mostly discards gyration. The DC-like parallel work ${\bar{W}}_{\parallel }$ goes to zero, and the DC-like contribution to super-adiabatic heating appears small. Fluid-like adiabatic compression is preserved in ${W}_{\mathrm{induct}}+{\bar{W}}_{{\rm{\nabla }}B}$. Figure 4(g) separates E_x, E_y, and E_z contributions to ${W}_{\parallel }$ as ${W}_{\parallel ,{Ei}}=-e\langle \sum {E}_{i}{b}_{i}{v}_{\parallel }{\rm{\Delta }}{t}_{\mathrm{out}}\rangle$, where $i=x,y,z$ and b_i is the ith component of $\hat{{\boldsymbol{b}}}$. E_y gives parallel heating, whereas E_x and E_z cause parallel cooling.

The quantities ${W}_{\perp ,\mathrm{ad}}$, ${W}_{\{\parallel ,\perp \},\mathrm{HF}}$, ${W}_{\mathrm{drift}}$, W_tot,out, ${\bar{W}}_{\mathrm{drift}}$, ${\bar{W}}_{\parallel }$, and ${W}_{\parallel ,{Ei}}$ are integrated with coarse time step Δt_out and converged at the $\sim 10 \%$ level. The error regions in Figures 4(a), (b), (e) are defined in Appendix B.

We have measured ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ in 2D PIC simulations of perpendicular shocks to inform models of astrophysical systems lacking direct ${T}_{{\rm{e}}}$ or ${T}_{{\rm{i}}}$ measurements. In a ${{ \mathcal M }}_{\mathrm{ms}}\,=3.1$, ${\beta }_{{\rm{p}}}=0.25$ rippled shock, quasi-static DC fields provide fluid-like adiabatic heating, and most super-adiabatic heating is from ion-scale E_∥ waves.

Xinyi Guo shared and provided helpful assistance for some software used to perform and analyze these simulations. Alex Bergier and colleagues provided excellent assistance with Columbia's Habanero cluster. Adam Masters and Parviz Ghavamian kindly shared Saturn bow shock data. We thank Matthew W. Abruzzo, Luca Comisso, Greg Howes, Anatoly Spitkovsky, Vassilis Tsiolis, and Lynn B. Wilson III for discussion. We thank the anonymous referees for integral comments and critiques. L.S. and A.T. were supported by the Sloan Fellowship to L.S., NASA ATP-80NSSC20K0565, and NSF AST-1716567. Some work was done at UCSB KITP, which is supported by NSF PHY-1748958. Simulations were run on Habanero (Columbia University), Edison (NERSC), and Pleiades (NASA HEC). Columbia University's Shared Research Computing Facility is supported by NIH Research Facility Improvement grant 1G20RR030893-01 and the New York State Empire State Development, Division of Science Technology and Innovation (NYSTAR) Contract C090171. NERSC is a U.S. Department of Energy Office of Science User Facility operated under Contract DE-AC02-05CH11231. The NASA HEC Program is part of the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.

Facilities: NERSC - , Pleiades -

Figures A1 and A2 present two animations of our electron sample traversing the shock front.

In Figure A3, we show the work decomposition from Figure 4 for many electron samples. At $t^{\prime} =0\,{\omega }_{\mathrm{pi}}^{-1}$, we selected all electrons in the regions $x\in [5.60,5.62]c/{\omega }_{\mathrm{pi}}$, $x\,\in [6.00,6.02]c/{\omega }_{\mathrm{pi}}$, and so on with even spacing $0.4c/{\omega }_{\mathrm{pi}}$ to get 17 electron particle samples of similar size.

Several quantities in Figure 4 are summed with a coarse time step ${\rm{\Delta }}{t}_{\mathrm{out}}=400{\rm{\Delta }}t=9\,{\omega }_{\mathrm{pe}}^{-1}$, namely: ${W}_{\perp ,\mathrm{ad}}$, ${W}_{\perp }-{W}_{\perp ,\mathrm{ad}}$, ${W}_{\perp ,\mathrm{HF}}$, ${W}_{\parallel ,\mathrm{HF}}$, ${W}_{\mathrm{induct}}$, ${W}_{{\rm{\nabla }}{\rm{B}}}$, W_inert, W_pol, W_tot,out, ${\bar{W}}_{{\rm{\nabla }}{\rm{B}}}$, ${\bar{W}}_{\mathrm{inert}}$, ${\bar{W}}_{\mathrm{pol}}$, ${\bar{W}}_{\parallel }$, ${W}_{\parallel ,{Ex}}$, ${W}_{\parallel ,{Ey}}$, and ${W}_{\parallel ,{Ez}}$. To check convergence, we downsample in time each quantity's summand by 4 and compute an error $\delta f$ at discrete time t_n as

$$\begin{eqnarray}&&\delta f({t}_{n})=\max \left\{\left|f({t}_{m})-{f}_{1/4}({t}_{m})\right|| 0\leqslant m\leqslant n\right\},\end{eqnarray} \tag{ B1 }$$

where ${f}_{1/4}$ is the 4× downsampled version of f. Thus $\delta f$ is strictly non-decreasing with t_n. The error regions defined by Equation (B1) are plotted as shaded areas in Figures 4 and A3.

Figure B1 shows curves with 2×, 4×, 8× downsampling for the 17 distinct electron samples of Figure A3, including the sample shown in Figure 4.

Figure C1 shows that our fiducial 2D simulations are converged with respect to transverse width. Most of the varying transverse width simulations are not listed in Table D1. For ${{ \mathcal M }}_{\mathrm{ms}}=9.1$ only, the 1D simulation uses a slightly higher upstream temperature than the 2D simulations, so the ratio ${{\rm{\Omega }}}_{{\rm{i}}}$/${\omega }_{\mathrm{pi}}$ differs between 1D and 2D. In this case, we matched times based on ${{\rm{\Omega }}}_{{\rm{i}}}^{-1}$ rather than ${\omega }_{\mathrm{pi}}^{-1}$.

Table D1 contains input parameters, derived shock parameters, run durations, and downstream temperature measurements for all simulations in our Letter. The first row is the high-resolution run used in Figures 3 and 4; the remaining rows are presented in Figure 2. A machine-readable version of Table D1, in comma-separated value (CSV) ASCII, is available. Below, we define all table columns.

• 1.  
  mi_me is the ion–electron mass ratio ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}$.
• 2.  
  theta and phi specify the upstream magnetic field orientation, measured in the simulation frame. θ is the angle between ${\boldsymbol{B}}$ and the x-coordinate axis. φ is the angle between the y–z plane projection of ${\boldsymbol{B}}$ and the z-coordinate axis. To visualize these angles, see Figure 1 of Guo et al. (2014), but note that their ${\varphi }_{{\rm{B}}}=\pi /2-\varphi$ is the complement of our φ. For all our simulations, θ corresponds to the angle between ${\boldsymbol{B}}$ and shock normal. The 2D simulations with in-plane ${\boldsymbol{B}}$ have θ = 90° and φ = 90°. The 2D simulations with out-of-plane ${\boldsymbol{B}}$ (i.e., ${\boldsymbol{B}}$ along $\hat{{\boldsymbol{z}}}$) have θ = 90° and φ = 0°. The 1D simulations with oblique ${\boldsymbol{B}}$ have have $\theta \lt 90^\circ$.
• 3.  
  my and mz are the numbers of grid cells along $\hat{{\boldsymbol{y}}}$ and $\hat{{\boldsymbol{z}}}$. Our 2D simulations have ${\mathtt{mz}}=1$, and 1D simulations have ${\mathtt{my}}={\mathtt{mz}}=1$.
• 4.  
  betap, Ms, Ma, and Mms are the shock plasma beta ${\beta }_{{\rm{p}}}$, sonic Mach number ${{ \mathcal M }}_{{\rm{s}}}$, Alfvén Mach number ${{ \mathcal M }}_{{\rm{A}}}$, and fast magnetosonic Mach number ${{ \mathcal M }}_{\mathrm{ms}}$. These numbers are derived from TRISTAN-MP input parameters sigma, delgam, and u0 (defined just below). First, the total plasma beta is

  $$\begin{eqnarray*}&&{\beta }_{{\rm{p}}}=\displaystyle \frac{4{\gamma }_{0}{\rm{\Delta }}{\gamma }_{{\rm{i}}}}{\sigma \left({\gamma }_{0}-1\right)\left(1+{m}_{{\rm{e}}}/{m}_{{\rm{i}}}\right)},\end{eqnarray*}$$

  where ${\gamma }_{0}=1/\sqrt{1-{\left({u}_{0}/c\right)}^{2}}$ is the Lorentz factor of the upstream flow in the simulation frame. The sonic Mach number depends on the upstream plasma speed in the shock's rest frame:

  $$\begin{eqnarray*}&&{{ \mathcal M }}_{{\rm{s}}}=\displaystyle \frac{{u}_{\mathrm{sh}}}{{c}_{{\rm{s}}}}=\displaystyle \frac{{u}_{0}}{\left(1-1/r\left({{ \mathcal M }}_{{\rm{s}}}\right)\right){c}_{{\rm{s}}}},\end{eqnarray*}$$

  and we solve this implicit expression for ${{ \mathcal M }}_{{\rm{s}}}$ (and thus also u_sh) using an input ${u}_{0}$ and assumed fluid adiabatic index Γ = 2 (note that Γ enters into both the Rankine–Hugoniot expression for MHD shock-compression ratio r and the sound speed c_s). Once ${{ \mathcal M }}_{{\rm{s}}}$ and u_sh are known, ${{ \mathcal M }}_{{\rm{A}}}$ and ${{ \mathcal M }}_{\mathrm{ms}}$ are known as well. This procedure for estimating shock parameters is taken directly from Guo et al. (2017).
• 5.  
  sigma is the magnetization, a ratio of upstream magnetic and kinetic enthalpy densities:

  $$\begin{eqnarray*}&&{\mathtt{sigma}}=\sigma \equiv \displaystyle \frac{{{B}_{0}}^{2}}{4\pi \left({\gamma }_{0}-1\right)\left({m}_{{\rm{i}}}+{m}_{{\rm{e}}}\right){n}_{0}{c}^{2}}.\end{eqnarray*}$$
• 6.  
  delgam is the upstream plasma temperature, in units of ion rest mass energy:

  $$\begin{eqnarray*}&&{\mathtt{delgam}}={\rm{\Delta }}{\gamma }_{{\rm{i}}}\equiv \displaystyle \frac{{k}_{{\rm{B}}}{T}_{0}}{{m}_{{\rm{i}}}{c}^{2}}.\end{eqnarray*}$$
• 7.  
  u0 is the upstream plasma velocity, in units of speed of light:

  $$\begin{eqnarray*}&&{\mathtt{u}}{\mathtt{0}}\equiv {u}_{0}/c.\end{eqnarray*}$$
• 8.  
  ppc0 is number of particles (both electrons and ions) per cell in the upstream plasma.
• 9.  
  c_omp is the number of grid cells per electron skin depth c/${\omega }_{\mathrm{pe}}$.
• 10.  
  ntimes is the number of current filter passes.
• 11.  
  dur is the simulation duration in units of upstream inverse ion cyclotron frequency ${{\rm{\Omega }}}_{{\rm{i}}}^{-1}$.
• 12.  
  Te_Ti is our measurement of downstream temperature ratio ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$. As described in the main text, we manually choose a downstream region that is minimally affected by the left-side reflecting wall and the right-side shock front relaxation. Our measurement of ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ uses downsampled grid output of the particle temperature tensor; however, the temperature tensor itself is calculated for each grid cell using the full particle distribution in a 5^N cell region, where N ∈ {1, 2, 3} is the domain dimensionality.
• 13.  
  Te_Ti_std is the standard deviation of ${T}_{{\rm{e}}}/{T}_{{\rm{i}}}$ within the downstream region that we consider. Like Te_Ti, downsampled grid output is used for this estimate.
• 14.  
  Te and Ti are the downstream electron and ion temperatures scaled to their respective rest mass energies, i.e., ${k}_{{\rm{B}}}{T}_{{\rm{e}}}/({m}_{{\rm{e}}}{c}^{2})$ and ${k}_{{\rm{B}}}{T}_{{\rm{i}}}/({m}_{{\rm{i}}}{c}^{2})$. We measure all of Te, Ti, and Te_Ti in the same manually chosen downstream region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (17 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (17 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size