Electron Heating in 2D Particle-in-Cell Simulations of Quasi-Perpendicular Low-Beta Shocks

We measure the thermal electron energization in 1D and 2D particle-in-cell (PIC) simulations of quasi-perpendicular, low-beta ($\beta_p=0.25$) collisionless ion-electron shocks with mass ratio $m_i/m_e=200$, fast Mach number $\mathcal{M}_{ms}=1$-$4$, and upstream magnetic field angle $\theta_{Bn} = 55$-$85^\circ$ from shock normal $\hat{\boldsymbol{n}}$. It is known that shock electron heating is described by an ambipolar, $\boldsymbol{B}$-parallel electric potential jump, $\Delta\phi_\parallel$, that scales roughly linearly with the electron temperature jump. Our simulations have $\Delta\phi_\parallel/(0.5 m_i {u_\mathrm{sh}}^2) \sim 0.1$-$0.2$ in units of the pre-shock ions' bulk kinetic energy, in agreement with prior measurements and simulations. Different ways to measure $\phi_\parallel$, including the use of de Hoffmann-Teller frame fields, agree to tens-of-percent accuracy. Neglecting off-diagonal electron pressure tensor terms can lead to a systematic underestimate of $\phi_\parallel$ in our low-$\beta_p$ shocks. We further focus on two $\theta_{Bn}=65^\circ$ shocks: a $\mathcal{M}_s=4$ ($\mathcal{M}_A=1.8$) case with a long, $30 d_i$ precursor of whistler waves along $\hat{\boldsymbol{n}}$, and a $\mathcal{M}_s=7$ ($\mathcal{M}_A=3.2$) case with a shorter, $5d_i$ precursor of whistlers oblique to both $\hat{\boldsymbol{n}}$ and $\boldsymbol{B}$; $d_i$ is the ion skin depth. Within the precursors, $\phi_\parallel$ has a secular rise towards the shock along multiple whistler wavelengths and also has localized spikes within magnetic troughs. In a 1D simulation of the $\mathcal{M}_s=4$, $\theta_{Bn}=65^\circ$ case, $\phi_\parallel$ shows a weak dependence on the electron plasma-to-cyclotron frequency ratio $\omega_{pe}/\Omega_{ce}$, and $\phi_\parallel$ decreases by a factor of 2 as $m_i/m_e$ is raised to the true proton-electron value of 1836.

1. INTRODUCTION Shocks heat and compress their host medium.In an ordinary fluid, collisions mediate heating; different particle species in the fluid share the same temperature after passing through a shock.In a low-density collisionless plasma, like the solar wind or interstellar medium, shocks partition their energy between various plasma waves and sub-populations of ions and electrons in a more complex way.We wish to know how much energy goes to electrons, and how and where within the shock the energy is gained, in order to (1) provide sub-grid electron heating prescriptions that may be used in two-fluid or hybrid fluid/kinetic simulations, and (2) help interpret observations of heliospheric and astrophysical shocks.
In this study, we take Earth's quasi-perpendicular, low-beta magnetospheric bow shock as an exemplar system.Typical parameters are upstream plasma beta (thermal/magnetic pressure ratio) β p ∼ 0.1-10, upstream magnetic field angle Corresponding author: Aaron Tran, Lorenzo Sironi aaron.tran@columbia.edu,lsironi@astro.columbia.edu* Now at Department of Physics, University of Wisconsin-Madison, USA from shock normal θ Bn ∼ 55-90 • , and fast magnetosonic Mach number M ms ∼ 1-10 (Farris et al. 1993, Figure 1).At such shocks, electrons may reflect from the shock front or leak from downstream to upstream, whereas thermal ions are generally confined to the downstream plasma after crossing the shock.
Electron heating can be described using a well-established cross-shock potential model (Goodrich & Scudder 1984;Scudder 1995;Hull et al. 1998).In a shock's de Hoffmann-Teller (HT) frame, electron and ion bulk flows are parallel to B both upstream and downstream.Entering the shock, a B-parallel electric potential ϕ ∥ (x) varying along the shocknormal coordinate x boosts electrons' parallel velocities.Small-scale scattering within the shock is expected to cool and smooth the post-shock electron distribution as compared to the adiabatic Liouville-mapping limit (i.e., test-particle evolution in static, macroscopic (ion-scale) B(x) and ϕ ∥ (x)) (Scudder et al. 1986a;Schwartz 2014).
The cross-shock potential model works and recovers trends in observed post-shock temperature anisotropy (Scudder et al. 1986a;Hull et al. 1998;Lefebvre et al. 2007).But, the model is silent on the origin of small-scale scattering waves required to create stable post-shock electron distributions (Schwartz 2014), e.g., flat-tops in Earth's magne-tosheath (Feldman et al. 1983b).Scattering waves diffuse electron velocities and may alter electron pressure; the electron pressure tensor divergence in Ohm's law sets ϕ ∥ , which in turn sets the free particle energy available to drive scattering waves (Veltri et al. 1990;Veltri & Zimbardo 1993a,b).
Wave-particle interactions also mediate electron heating, in a complementary view.Ion-scale whistler waves oblique to B are a prime suspect.Low-Mach shocks host largeamplitude whistler precursors with a variety of propagation angles (Ramírez Vélez et al. 2012;Wilson et al. 2017;Davis et al. 2021).Whistlers may radiate from the shock ramp (Krasnoselskikh et al. 2002;Dimmock et al. 2019) or be generated by instabilities in the shock foot (Wu et al. 1984;Winske et al. 1985;Muschietti & Lembège 2017).These whistlers have low-frequency B-parallel electric fields that may trap and heat electrons, and strong electron trapping (parallel electric field energy comparable to the electrons' thermal energy) induces phase mixing and short-wavelength electrostatic structures (McBride et al. 1972;Matsukiyo & Scholer 2003, 2006;Matsukiyo 2010).
The connection between scattering waves, electron heating, and the cross-shock potential (or, ion-scale fields and waves more generally) is not yet fully understood.Recent Magnetospheric Multiscale (MMS) mission studies have given new insights.For example, Chen et al. (2018) suggest that the cross-shock potential could arise from the cumulative effect of oblique, non-linear whistler wave E-fields; Sun et al. (2022) suggest that strong, Debye-scale double layers may help form the cross-shock potential.Cohen et al. (2019) have measured the electrostatic potential jump across an interplanetary shock, along with collocated, strong, highfrequency electrostatic fields.Hull et al. (2020) showed the onset of electron heating as a transition from adiabatic to non-adiabatic electron distributions over multiple wavelengths of a shock's oblique whistler precursor.Much recent work has been done on electron scattering by coherent lowand high-frequency whistlers at shocks (Hull et al. 2012;Wilson et al. 2016;Oka et al. 2017Oka et al. , 2019;;Page et al. 2021;Artemyev et al. 2022;Shi et al. 2023), including the construction and use of general frameworks such as stochastic shock drift acceleration (Katou & Amano 2019;Amano et al. 2020) and magnetic pumping (Lichko & Egedal 2020), with particular success in explaining non-thermal distributions.
We seek to address a few general questions.First, how much do electrons heat (∆T e , T e /T i ) across quasiperpendicular shocks?Does heating occur at the ramp (steepest density rise), foot (region in which shock-reflected ions execute a single Larmor gyration), or farther upstream within a precursor?How does the heating amount and location depend on shock parameters (M ms , θ Bn ) and simulation parameters (mass ratio, domain dimensionality)?Second, what is the origin and nature of the cross-shock potential in PIC simulations?What structure does ϕ ∥ (x) possess along the shock-normal coordinate x, beyond approximations such as ϕ ∥ (x) ∝ B(x) (Hull et al. 2000;Hull & Scudder 2000;Artemyev et al. 2022)?How wide is the potential jump in x? Third, because ϕ ∥ is not easy to directly measure in the heliosphere, how faithfully do various observational proxies (e.g., the HT-frame electric field) reproduce ϕ ∥ ? 2. SHOCK SETUP 2.1.Simulation Methods We simulate 1D and 2D ion-electron shocks using the relativistic particle-in-cell (PIC) code TRISTAN-MP (Buneman 1993;Spitkovsky 2005).Our setup follows Sironi & Spitkovsky (2009); Guo et al. (2014Guo et al. ( , 2017)); Tran & Sironi (2020).We use Gaussian CGS units throughout.
To form a shock, we inject upstream plasma with velocity −u 0 x that reflects off a conducting wall at x = 0 in a uniform Cartesian (x, y, z) grid.Coordinate unit vectors are x, ŷ, ẑ.The reflected plasma interacts with upstream plasma to drive a shock traveling towards + x.The domain's left-side (x = 0) boundary specularly reflects all particles in the simulation frame (v x → −v x ).The domain's rightside x boundary expands ahead of the shock, continuously adding new plasma to the upstream; the expansion speed is manually chosen to out-run waves, thermal electrons leaking from downstream to upstream, and mirror-reflected (shockdrift accelerated) particles that may all travel ahead of the shock.Particles and waves exiting the right-side x boundary are deleted.The y and z boundaries are periodic, and we simulate 2D shocks in the x-y plane.All three Cartesian components of particle velocities and electromagnetic fields are tracked in both 1D and 2D shocks.
The injected upstream plasma comprises Maxwellian ions and electrons with equal density n 0 and temperature T 0 , carrying a magnetic field B with magnitude B 0 and angle θ Bn with respect to shock normal ( n = x).The upstream B lies in the x-y plane.The ion and electron plasma frequencies ω p{i,e} = 4πn 0 e 2 /m {i,e} , inertial lengths d {i,e} = c/ω p{i,e} , and cyclotron frequencies Ω c{i,e} = eB 0 /(m {i,e} c) are defined with respect to upstream plasma properties.Subscripts i and e denote ions and electrons respectively.Subscripts ⊥ and ∥ are generally defined with respect to local B. A shock is specified by several dimensionless parameters: M s , (or M A , M ms ), β p , θ Bn , Γ, and c s /c.The total plasma beta β p = 16πn 0 k B T 0 /B 0 2 .The fast magnetosonic, sonic, and Alfvén Mach numbers are M ms = u sh / √ c s 2 + v A 2 , M s = u sh /c s , and M A = u sh /v A respectively.1 Here, u sh is upstream plasma speed in the shock's rest frame.The sound speed c s = 2Γk B T 0 /(m i + m e ) and the Alfvén speed v A = B 0 / 4πn 0 (m i + m e ), with Γ the one-fluid adiabatic index.We scale all velocities with respect to c by choosing any of the ratios How do we set u 0 to obtain a desired Mach number?In the simulation frame, the downstream (post-shock) plasma has bulk u x = 0, and the shock speed is u sh /r, where r is the post-shock compression ratio.We relate u 0 and u sh by an explicit boost: (Tidman & Krall 1971), we must guess the effective one-fluid adiabatic index Γ, which is not known a priori for our kinetic plasma.We adopt Γ = 5/3 to set u 0 for all simulations in this manuscript.We expect Γ ≈ 5/3 for stronger or more oblique shocks with good coupling between ion velocities parallel and perpendicular to B. But Γ ≈ 2, corresponding to a 2D non-relativistic gas, may be more appropriate for weak, nearly perpendicular shocks with little ion isotropization.Our choice of Γ = 5/3 to initialize shocks that are best described by Γ ≈ 2 can cause our stated Mach numbers to be inaccurate by ∼10% or less.
We define the shock position x sh to be located where the magnetic fluctuation B/B 0 − 1 declines to one-half its maximum value, measured rightwards (increasing x) from the maximum value of B/B 0 just after the shock ramp.
In oblique MHD shocks (i.e., θ Bn ̸ = 90 • ) with B in the x-y plane, upstream plasma with bulk velocity u y = 0 acquire a transverse drift u y ̸ = 0 after transiting the shock.Our left-side (x = 0) boundary accommodates this drift by imposing E z = (u y /c)B 0 cos θ Bn in the 10 left-most x-cells of the Yee mesh.Particles reflect 5 x-cells rightwards of the imposed-field region.The boundary u y is computed from the MHD R-H jump conditions, and it is typically ∼10% of u sh (Tidman & Krall 1971, Figures 1.5-1.6).The boundary mimics a conducting wall that slides in the y direction.
Our choice of particle reflection procedure at the left boundary can alter post-shock electron distributions.Specular reflection, v x → −v x , mixes field-parallel and perpendicular particle momenta in the downstream rest frame and thus may alter particle distributions' anisotropy.Other choices are possible, e.g., Krauss-Varban et al. (1995) let particles escape the domain at both upstream and downstream domain boundaries, with new particles continuously injected from Maxwellian distributions.In nature, the post-shock boundary condition may be set by global system effects; e.g., Mitchell et al. (2012); Mitchell & Schwartz (2013, 2014); Schwartz et al. (2019); Horaites et al. (2021) construct collisionless magnetosheath models wherein electron distributions are regulated by the two points at which magnetosheath B field lines penetrate the bow shock.Because we do not consider global transport effects upon shock structure, our simulations can only inform us about electron heating local to a shock, under the assumption that post-shock electrons catching up to the shock (i.e., traveling from downstream to upstream) have a similar distribution as those streaming away from the shock.

Shock Parameters
We fix m i /m e = 200 and β p = 0.25; we vary the magnetic obliquity θ Bn = 55-85 • and sonic Mach number M s = 3-10.The transverse shock width is 5.4 d i for 2D runs.We choose the upstream electron temperature Θ e = k B T 0 /(m e c 2 ) to keep post-shock electrons marginally non-relativistic; lower Θ e is more realistic and costlier.For M s = 3-4 we take Θ e = 0.01; for M s = 5, Θ e = 0.007; for M s = 7, Θ e = 0.005; for M s = 10, Θ e = 0.003.
The grid spacing ∆x = ∆y = 0.1d e and the timestep ∆t = 0.045ω pe −1 , so c = 0.45∆x/∆t.We use 128 particles per cell for 2D runs, and 2048 for 1D runs.The upstream thermal electron gyroradius r Le = (Γβ p /4) 1/2 d e is resolved with ∼3 cells at β p = 0.25.The upstream Debye length e d e is marginally resolved with one cell for M s = 3-5 and not resolved for M s = 7-10.At each step in the PIC algorithm, the electric current is smoothed with N = 32 sweeps of a 3-point ("1-2-1") digital filter in each axis (Birdsall & Langdon 1991, Appendix C); the filter has a halfpower cut-off at k ≈ 2/N (∆x) −1 = 2.5d e −1 .
3. REGIME MAP OF SHOCK PARAMETERS Our simulations inhabit various regimes of collisionless shock behavior.Figure 2 shows approximately where precursors of leaking thermal electrons (green shaded), waves (blue shaded), and shock-reflected electrons (convex, black dotted contours) may appear upstream of a shock as a function of M s and θ Bn .The top panel is for our simulated shock parameters.The bottom panel is for more realistic, solar windlike parameters; adopting the true proton/electron mass ratio m i /m e = 1836, in particular, permits thermal electron and whistler wave precursors for a large region of (M s , θ Bn ) parameter space.Figure 2  First, some notation.We rescale the parallel potential by the incoming ion bulk kinetic energy: (1) The tilde (∼) will mean the same normalization for other electric potentials throughout this manuscript.And, the prefix ∆ (e.g., ∆ φ∥ ) will generally mean a quantity's crossshock jump value, except for the grid step sizes ∆x, ∆y, and ∆t.
In the blue-shaded region of Figure 2, whistlers propagating along shock normal have phase and group speeds faster than the shock speed, so whistler wave trains may "phase stand" ahead of the shock (Tidman & Northrop 1968;Kennel et al. 1985;Krasnoselskikh et al. 2002;Oka et al. 2006).with solar wind-like T0 = 10 eV (bottom panel).Black dots in top panel mark our simulations; circled black dots are studied in more depth in Section 7. Black shaded region is Mms < 1, calculating MHD fast speed with θBn dependence (unlike definition and use of Mms elsewhere in manuscript; definitions agree at θBn = 90 • ).Blue shaded region is whistler sub-critical, MA < Mw. Green shaded region shows which shocks may permit thermal electrons to escape from downstream to upstream, for varying ∆ φ∥ = 0, 0.05, and 0.1 (Equation ( 1)).Red dotted line marks where Mms equals the critical Mach number (Marshall 1955;Kennel et al. 1985).Orange dashed line is sub/super-luminal boundary for Θe = 0.01; note that the Ms = 7, 10 and θBn = 85 • simulations are marginally sub-luminal, not super-luminal, because of their lower Θe values.Black dotted contours bound the regions where the SDA reflection efficiency is > 1% for varying ∆ φ∥ = 0, 0.05, and 0.1 (larger to smaller regions).
Such a precursor requires M A < M w , where the whistler Mach number and v w is the maximum phase speed of an oblique whistler based on an approximate cold-plasma dispersion relation (Krasnoselskikh et al. 2002).The orange dashed line of Figure 2 denotes the sub/superluminal boundary.Shocks are subluminal for u sh tan θ Bn < c, where θ Bn is measured in the shock frame (also called normal incidence frame, NIF); we use θ Bn measured in the simulation (downstream rest) frame, since B is nearly frame invariant when the relative velocity between NIF and simulation frame is non-relativistic.Shocks right of the boundary are sub-luminal; left, super-luminal.All sub-luminal shocks may be boosted into the de-Hoffmann Teller (HT) frame.All super-luminal shocks may be boosted into a perpendicular shock frame wherein the shock is stationary and both upstream and downstream B fields are oriented 90 • from shock normal n (Drury 1983;Kirk & Heavens 1989;Begelman & Kirk 1990).For smaller Θ e , ceteris paribus, the super-/sub-luminal boundary shifts leftward so that super-luminal shocks occupy a smaller region of the plot, while the other curves and regions in the regime plot are unchanged.Our M s = 5, M s = 7, and M s = 10 simulations use smaller values of Θ e than shown in Figure 2 (Table 1).The two simulations with M s = 7, 10 and θ Bn = 85 • are both marginally sub-luminal.
In the green-shaded regions of Figure 2, downstream thermal electrons may escape from the shock into the upstream.
For a given (θ Bn , M s , ∆ φ∥ ), we define a thermal escape criterion as follows.Consider a thermal electron with downstream rest frame four-velocity γβ ∥ equal to the mean value for a Maxwell-Jüttner distribution of temperature T e2 , and γβ ⊥ = 0 (i.e., zero pitch angle).Here γ is the Lorentz factor; β is three-velocity in units of c.We compute an estimated post-shock electron temperature T e2 from the non-relativistic oblique MHD R-H conditions (Tidman & Krall 1971) assuming post-shock equipartition, T e2 = T i2 .We choose the sign of γβ ∥ > 0 to indicate travel towards + x in the downstream rest frame.Now, boost into the HT frame.If the electron has HT-frame β ∥ > 0 and HT-frame kinetic energy exceeding the cross-shock potential jump, i.e., (γ − 1)m e c 2 > ∆ϕ ∥ , then we conclude that "typical" thermal electrons may escape from downstream to upstream.Equivalently, for nonrelativistic HT-frame boosts, the escape criterion is approximately equal to where c se2 is the downstream electron thermal speed interpreted here as a B-parallel velocity, and θ Bn2 is the postshock magnetic field angle with respect to shock normal.0, 0.05, 0.1, with successively darker green shading corresponding to larger ∆ φ∥ .Our assumptions somewhat overestimate the possibility of post-shock thermal electron escape; the quasi-perpendicular shocks in this manuscript typically have post-shock T e2 < T i2 and ∆ φ∥ ∼ 0.1-0.2.Our escape criterion also neglects mirroring as electrons travel from high to low B magnitude, which would aid the escape of electrons with non-zero pitch angle as perpendicular energy transfers to parallel energy.
The red dotted line of Figure 2 denotes the critical Mach number (Marshall 1955;Edmiston & Kennel 1984;Kennel et al. 1985); for M ms exceeding this Mach number, ion reflection is expected to become important for shock dissipation.
The black dotted contours in Figure 2 denote where shock drift acceleration (SDA) may reflect > 1% of incoming thermal electrons.Successive (smaller) contours correspond to ∆ φ∥ = 0, 0.05, and 0.1.We calculate the SDA reflection efficiency by applying Guo et al. (2014, Equations ( 13)-( 14)) to a Monte-Carlo sampling of upstream Maxwell-Jüttner distributions for a grid of shock parameters (M s , θ Bn ).The reflection efficiency is defined as the ratio between the number of particles with upstream HT-frame v ∥ < 0 (approaching the shock) that would mirror and reflect upon reaching the shock front, compared to the total number of particles with upstream HT-frame v ∥ < 0.

OVERVIEW OF SHOCK SIMULATIONS
Figure 3 shows the ion density structure for our (M s , θ Bn ) parameter sweep of 2D shocks.Shocks change from laminar to turbulent as the Mach number rises.All M s = 3 shocks appear laminar.As θ Bn decreases from 85 • for the weaker shocks (M s = 3-5), shock precursors appear as coherent waves and filaments at x > x sh (Figure 3), and also spatiallycoincident rises in T e (Figure 4).Towards higher M s , ion reflection becomes important; we observe discrete clumps of ions and distorted B-field lines within 2d i of the shock ramp; this reflection is seen for all θ Bn .To limit the scope of this work, we cannot firmly identify all the instabilities or linear modes in Figure 3, but let us speculate briefly.Ion/electron beam drifts coupling to the whistler mode (various forms of modified two-stream instability, aka MTSI) (Wu et al. 1984;Hellinger & Mangeney 1999;Matsukiyo & Scholer 2003, 2006;Muschietti & Lembège 2017) could be responsible for fine compressive filamentation with wavevector nearly perpendicular to B seen in the M s = 5, θ Bn = 65, 55 • cases.At M s = 10 and θ Bn = 55 • , the localized field-line bending and less coherent fluctuations (i.e., not clearly wave-like) may be a weaker form of the turbulence and reconnection in 2D and 3D PIC quasi-parallel shocks simulated by Bessho et al. (2020Bessho et al. ( , 2022)); Ng et al. (2022).At M s = 7-10 and θ Bn ∼ 85-75 • , the reflected ion clumps and shock ramp fluctuations are likely similar to those reported by Tran & Sironi (2020).
Figure 4 shows the 1D y-averaged electron temperature T e for both 1D and 2D shocks.The temperature in units of m e c2 is the trace of γβ i β j f e (r, p)d 3 p for electron distribution function f e and tensor indices i,j.In the near-perpendicular θ Bn = 85 • and 75 • shocks, electron heating is much greater in 2D than in 1D; a 2D domain is required for non-adiabatic heating to occur within the shock ramp.Ions' effective Γ and hence the shock propagation speed both change when comparing 1D and 2D, especially at M s ≳ 5. Towards lower θ Bn , the trend changes.The weakest, oblique shocks generally have similar 1D and 2D heating (e.g.M s = 3, θ Bn = 55 • ).Some stronger oblique 2D shocks heat electrons less than their 1D counterparts (M s = 5-7, θ Bn = 65-55 • ).
Figure 5 plots the downstream electron/ion temperature ratio T e /T i for the 2D runs.We measure T e /T i in the downstream plasma within x − x sh = −18 to −3d i ; details are given in Appendix A. For comparison, we also show exactlyperpendicular shock data at m i /m e = 200 and 625 that were previously presented in Tran & Sironi (2020).The dotted black line corresponds to an adiabatic Γ = 2 compression for electrons based on the MHD R-H shock jump prediction at θ Bn = 90 • .We observe that the electron-ion temperature ratio at θ Bn = 85 • is similar to that at θ Bn = 90 • .And, T e /T i increases towards unity as θ Bn decreases.

PARALLEL POTENTIAL PRAXIS
In quasi-perpendicular shocks where electrons are wellmagnetized, and the de Hoffmann-Teller (HT) frame is accessible (i.e., the shock is subluminal), the space physics community has an established description of electron heating.Electron motion, at 0th order, is assumed to obey quasistatic, macroscopic fields B(x), E(x) varying only along x with no time dependence (Feldman et al. 1983b;Goodrich & Scudder 1984;Scudder et al. 1986a;Hull et al. 1998Hull et al. , 2000Hull et al. , 2001)).By "macroscopic", we mean that B and E are assumed to vary on ion length scales.Individual electron trajectories within the HT frame may be modeled as conserving two invariants, energy ε and magnetic moment µ (Hull et al. 1998): where ϕ ∥ is an electrostatic potential that forms within the shock transition due to the B-parallel electric field component, E ∥ .The parallel potential reads: where b is the magnetic field unit vector, and E ∥ is set by the electron momentum equation (i.e., generalized Ohm's law): Here, V e is the electron bulk velocity, P e is the electron pressure tensor, n e is the electron density, and d/dt = ∂/∂t + V e • ∇ is the total derivative.Collisions are neglected.The inertial term dV e /dt is small, 2 so the ambipolar term ∇ • P e sets the parallel electric field in all frames: In the HT frame, V e is parallel to B, up to usually-small correction terms localized within the shock layer (Scudder 1987).So, E HT ≈ E amb as well.In the following discussion, we treat simulation-frame measurements of E ∥ , ϕ ∥ , and B as equal to their HT-frame values, because all three quantities are frame invariant to order O(u/c) in boost velocity u.This quasi-static, macroscopic description is not complete, for shocks host high-frequency and short-wavelength fields that will scatter electrons and break the invariants µ and ε (Veltri et al. 1990;Hull et al. 1998;See et al. 2013;Schwartz 2014).Within the context of this macroscopic-field model, scattering is required to obtain physical (stable) electron distributions, and scattering will alter P e and therefore ϕ ∥ as well.Nevertheless, the macroscopic fields do successfully describe the electron dynamics observed by spacecraft (Scudder et al. 1986a;Lefebvre et al. 2007).
Let us now compare various ways to measure ϕ ∥ , which we summarize in Figures 6 and 7 for one example shock simulation.See also Schwartz et al. (2021, Section 3) for a recent, similar discussion of ϕ ∥ measurement methods applied to MMS data.

Liouville Mapping Measurement of the Parallel Potential
We can estimate ∆ϕ ∥ using Liouville mapping fits of the electron distribution function, as is commonly done to interpret satellite data (Scudder et al. 1986a;Schwartz et al. 1988;Hull et al. 1998;Lefebvre et al. 2007;Johlander et al. 2023).If electrons evolve adiabatically across a shock, their distribution at small pitch angle will inflate in energy following Liouville's theorem; by comparing upstream and downstream electron distributions, the amount of energy gain can be fitted to infer ∆ϕ ∥ .
We extract both upstream (unshocked) and downstream (shocked) electron distribution function cuts along p ∥ by selecting particles with pitch angle < 15 • , all in the HT frame.Downstream distributions are sampled from x − x sh = −18 to −3d i .Upstream distributions are sampled from x − x sh = 65d i to 105d i for θ Bn = 55 • shocks, 45d i to 85d i for θ Bn = 65 • shocks, and 25d i to 65d i for θ Bn = 75 • shocks.
Next, we map the sampled upstream particles to a notional downstream state, specified by a magnetic compression ratio B 2 /B 0 and potential jump ∆ϕ ∥ .We fix B 2 /B 0 to its volume-averaged value in the region x−x sh = −18 to −3d i , while ∆ϕ ∥ is a free parameter.Let u ⊥ = γβ ⊥ and u ∥ = γβ ∥ be the dimensionless four-velocity components of each particle.The Liouville-mapped particle momenta u ⊥2 and u ∥2 , written in terms of the upstream momenta u ⊥0 and u ∥0 , are: The mapping applies to both incoming (u ∥0 < 0) and outgoing (u ∥0 > 0) electron populations, so the outgoing electrons are mapped backwards in time.We then select particles with pitch angle < 15 • and bin in p ∥ to form a Liouvillemapped distribution (Figure 6, orange curves), which can be compared to the true downstream distribution (Figure 6, blue curve).The Liouville-mapped distributions are normalized to both momentum-and real-space volumes.The momentum-space volume is a cone fixed by the pitch-angle selection.The real-space volume changes across the shock due to Bperpendicular compression and B-parallel velocity change across the shock; the ratio of 3D real-space volume elements In the non-relativistic limit, If two particles are separated along a B field line as they cross the potential jump ∆ϕ ∥ , their time-staggered accelerations will tend to increase their real-space separation.This ∆ϕ ∥ -induced dilation is most important for small pitch angles and small initial parallel velocities, i.e., β ∥0 ∼ 0.
In contrast to our particle-based procedure, some studies (e.g., Chen et al. 2018;Cohen et al. 2019) perform Liouville mapping by shifting the distribution f in energy coordinates, f (ε) → f (ε + e∆ϕ ∥ ).We find general agreement between the two procedures in one-off experiments with our simulation data.
The Liouville-mapped distributions are fitted to the downstream distributions using a least-squares fit that minimizes a cost function where N is the number of momentum bins.We separately fit the p ∥ < 0 and p ∥ > 0 parts of f within multiple bands of f values, manually chosen to lie below the downstream flat-top f value and above counting noise at high energies.
Figure 6 shows the Liouville mapping procedure for the 2D, M s = 4, θ Bn = 65 • shock.The mean best-fit ∆ φ∥ = 0.12 for incoming electrons (HT-frame p ∥ < 0, upstream to downstream).The mean best-fit ∆ φ∥ = 0.061 for outgoing electrons (HT-frame p ∥ > 0, downstream to upstream).Fit values have a ±∼20% range for incoming electrons and a larger range for outgoing electrons.Appendix B shows the Liouville-mapped distributions for all the 2D shocks in our parameter sweep.

Direct Measurement of the Parallel Potential
In PIC simulations, we can measure ϕ ∥ from global particle and field information that is not always available in satellite measurements.We consider only 1D, y-averaged measures of ϕ ∥ and neglect the possibility of different potential jumps along different B field lines in 2D.For this section of the manuscript, angle brackets ⟨• • •⟩ denote transverse (yaxis) averaging.Integrals have the same limits as in Equation (4).
First, we may measure Second, we may measure E ∥ from the electron pressure tensor divergence (Figure 7 Liouville mapping: by comparing upstream (black curve) and downstream (blue curve) electron distribution cuts along the p ∥ axis, we fit Liouville-mapped upstream distributions (orange curves) to the downstream distributions in order to estimate ∆ φ∥ for different electron populations (incoming p ∥ < 0 versus outgoing p ∥ > 0) and for different bands in f .For incoming p ∥ < 0, the smallest and largest best-fit values of ∆ φ∥ = 0.100 and 0.150 are highlighted in orange alongside their corresponding Liouvillemapped distributions.For outgoing p ∥ > 0, the smallest and largest best-fit values of ∆ φ∥ = 0.045 and 0.073 are similarly highlighted.Momentum is scaled using the upstream electron thermal speed cse = ΓkBT0/me.The shock is the same 2D, Ms = 4, θBn = 65 • case in Figure 7.
We compute P e on a grid by depositing particles with a 5 cell flat-top kernel (5 × 5 in 2D), and we evaluate ∇ using a 2nd-order centered finite difference.
Third, we can simplify Equation ( 11).Assume a gyrotropic pressure tensor P e = P e∥ bb + P e⊥ (I − bb ), following Goodrich & Scudder (1984), to obtain (Figure 7 (12) All d/dy and d/dz terms are neglected.Note that the spatial average is taken over the entire integrand, unlike Equations ( 10) and ( 11).
We have chosen a particular order for the transverseaverage operations, ⟨ b • E⟩/⟨ x • b⟩ for Equation ( 10) and the same with E amb replacing E for Equation (11).We could have instead averaged prior to all the vector projections, e.g., integrate ⟨ b⟩ • ⟨E⟩/⟨ x • b⟩ in Equation (10) ("early" averaging).Or, we could have deferred averaging until the entire integrand is formed, e.g., integrate ⟨( b • E)/( x • b)⟩ in Equation (10) ("late" averaging).We compared the three averaging procedures across the shock parameter range.The "late" averaging gives poor results; it may be too sensitive to local regions where x • b approaches zero, whereas ⟨ x • b⟩ is constrained to be non-zero in the other procedures due to B x conservation across the shock.The "early" averaging agrees well with our adopted procedure in some cases, but not all.
There is not a clear reason to favor or disfavor our adopted procedure as compared to "early" averaging, due to the inherent approximation of describing a 2D shock with a 1D yaveraged profile; disagreement between the two procedures is due solely to 2D effects.

Indirect Measurement of the Parallel Potential
The HT-frame cross-shock potential (Figure 7(e)), will approximately equal the parallel potential ϕ ∥ (x) if E amb points towards + x in the HT frame (Goodrich & Scudder 1984).Why should E amb point towards + x? Suppose that P e is isotropic (scalar) and that the shock structure varies solely along the shock-normal coordinate x.Then, the ambipolar field and the parallel potential Since E HT ≈ E amb , then ϕ HT ≈ ϕ ∥ .
Next, let us consider a situation in which E amb may not point towards + x.Suppose that P e is gyrotropic and not isotropic.Then, P e has off-diagonal terms proportional to (P e∥ − P e⊥ ) because b does not coincide with a coordinate axis in the HT frame.Let us still assume a shock varying solely along x; i.e., neglect d/dy and d/dz terms.Then, the ambipolar field We see ϕ ∥ ̸ = ϕ HT due to the non-x components of E amb in this situation. 3This electron anisotropy effect is just one of several higher-order corrections to the quasi-static field model as discussed by Scudder (1987).
We test the importance of the off-diagonal P e terms in the simulation coordinate system by computing shown as a blue dotted line in Figures 7(c) and 8.
If E HT = E amb and the shock is 1D-like (d/dy and d/dz terms negligible), then Equation ( 15) resolves to Equation ( 14), which should equal ϕ HT .If in addition E amb points towards + x, then Equation ( 15) should equal both ϕ HT and ϕ ∥,amb .We shall see that these assumptions do not fully hold in our simulations.

Parallel Potential Praxis -Results
Figure 7 compares different proxies for ϕ ∥ in our example M s = 4, θ Bn = 65 • shock (same as Figure 1), All potentials are scaled to the shock-frame upstream ion bulk kinetic energy (Equation ( 1)).We use 50 snapshots with close time spacing to smooth out short-timescale variation.Individual snapshots are plotted as faint colored lines, while their average is plotted as a thicker solid line.The time spacing between snapshots is ∆x/ k B T 0 /m e , where k B T 0 /m e is an upstream electron thermal velocity; the snapshots span t = 40.026 to 40.123Ω ci −1 .Our choice of time spacing allows thermal electrons to translate by ≳ 1 grid cell between snapshots, in order to decorrelate thermal fluctuations in consecutive snapshots while keeping the shock stationary on ion scales.We set ϕ = 0 at x − x sh = 45d i for each ϕ ∥ proxy and for each snapshot before averaging the data together.
What do we learn from Figure 7? We observe an extended rise in the potential ahead of the shock, spanning tens of d i , corresponding to a precursor region of increased T e as seen in Figure 4. Local spikes in ϕ ∥ (potential wells for electrons) appear at the position of magnetic troughs in the precursor wave train; a strong potential spike appears in the magnetic trough immediately adjacent to the shock ramp at x = x sh (Figure 7(a)).The E ∥ measured from the PIC mesh is noisy; multiple time snapshots must be averaged to obtain a good measurement of ϕ ∥,grid .Both PIC field-based estimators ϕ ∥,grid and ϕ HT show more fluctuation than the particle-based estimators ϕ ∥,amb and ϕ ∥,gyr .The particlebased estimators can give a good estimate of the potential ϕ ∥ with a single snapshot.Both particle-based estimators are in reasonable agreement with each other, suggesting that gyrotropy is a good assumption.The particle-based estimators are close to, albeit ∼5% lower than, ϕ ∥,grid .
The net cross-shock jump in ϕ HT (Figure 7(e)) agrees with the other proxies we consider, but the potential shape 3 Analogous to the disagreement between normal-incidence and HT frame cross-shock potentials explained by Goodrich & Scudder (1984).along x shows pronounced d i -scale fluctuation that does not appear in the other proxies; e.g., in the post-shock region x − x sh = −10 to 0d i and the foreshock region x − x sh = 15 to 30d i .Why is this?If the off-diagonal terms in P e were the sole cause for disagreement between ϕ ∥ and ϕ HT , as discussed in Section 5.3, then we would expect Equation 15(blue dotted) to coincide with ϕ HT (orange).But that is not so.The disagreement between ϕ ∥,amb and ϕ HT must come from other effects, e.g., E HT ̸ = E amb due to lowfrequency motional fluctuations outside the shock ramp that are not fully cancelled by the global HT frame boost.Our calculation of Equation ( 15) also shows that the nondiagonal terms in P e contribute appreciably to ϕ ∥,amb in our shocks.The non-diagonal P e gradients are responsi- ble for (i) potential spikes within the shock precursor region x − x sh = 0 to 30d i , and (ii) a ∼20% increase in ϕ ∥ across the shock ramp at x = x sh , compared to integrating only the x • (∇ • P e ) piece (Figure 7(c)).We expect that our shocks, having relatively low β p , will show larger pressure anisotropies than higher-β p shocks wherein pressure anisotropy may be bounded to lower magnitude by various microinstabilities.
Figure 8 proceeds to a wider shock parameter range.We show the same ϕ ∥ proxies for θ Bn = 75 • to 55 • and varying M s , following the same procedure as in Figure 7.All runs use 50 evenly-spaced snapshots within a short time interval.For M s = 3-5, t ≈ 40.02-40.12Ωci −1 (with ±0.01Ω ci −1 variation on the exact timing for individual runs).For M s = 7, t = 30.04-30.14Ωci −1 .For M s = 10, t = 20.23-20.33Ωci −1 .The potential is pinned to ϕ = 0 at x − x sh = 45d i for all shocks and time snapshots.
What do we learn from Figure 8? The HT-frame potential is noisy, but it does a reasonable job of replicating the magnitude of the ϕ ∥ jump across the shock ramp, as measured by other proxies.The HT-frame potential deviates from other ϕ ∥ estimators within the precursor wave trains of low M s and low θ Bn shocks.For example, in the M s = 3, θ Bn = 55 • case, the HT potential shows a gradual rise with less evident fluctuations than ϕ ∥ .In the M s = 4, θ Bn = 55 • case, the HT potential shows more fluctuation than ϕ ∥ .Disagreement between ϕ HT and ϕ ∥ is reasonable in such precursors because the HT-frame boost's cancellation of motional electric fields may be imperfect within a shock transition of finite width.
The particle-based ϕ ∥ proxies agree well with each other for θ Bn = 65-55 • (Figure 8); more disagreement is seen in θ Bn = 75 • shocks, which are closer to perpendicular and which we expect to deviate from the quasi-static electron heating model description (Goodrich & Scudder 1984).The direct grid measure ϕ ∥,grid can be noisy for weak shocks, but our averaging procedure reduces the uncertainty across most of the shock parameter range considered.We suggest that agreement or disagreement between our proxy measurements in Figures 7 and 8 can be broadly attributed to systematic error: 2D effects, neglected terms in the generalized Ohm's law (i.e., effective frictional force from collisionless scattering), or finite-Larmor radius drift contributions to electron heating that involve B-perpendicular electric fields (Northrop 1963).
The calculation of Equation ( 15), which tests whether E amb may be approximated as lying along x, returns a crossshock potential jump that is lower than the other ϕ ∥ proxies in all but the weakest (M s = 3) shocks.
Lastly, the Liouville-mapping inferred values for ∆ φ∥ (red triangles) show multiple trends across the (M s , θ Bn ) parameter range.The reader may inspect the detailed Liouville mapping fits of f (r, p) in Appendix B. Weaker, lower-θ Bn shocks (upper right of Figure 8) return a Liouville-mapping ∆ φ∥ value for incoming electrons that is in general agreement with Equations ( 10)-( 13) and has a range of tens of percent, when fitted to different bands in f (see Figure 6).The inferred potential for outgoing (p ∥ > 0) electrons is consistently smaller than that for incoming (p ∥ < 0) electrons.Towards stronger and higher-θ Bn shocks, the asymmetry between outgoing-and incoming-inferred potentials becomes larger, and the Liouville-inferred potential jump may disagree with Equations ( 10)-( 13) by a factor of 2× or more.In cases where the outgoing (p ∥ > 0) electrons return a value of ∆ϕ ∥ significantly smaller than other estimators such as ϕ ∥,amb , the Liouville mapping may not have a straightforward physical interpretation.Downstream electrons with p ∥ > 0 could be well confined by ϕ ∥ , and other mechanisms in the shock ramp may be needed to explain the outgoing upstream electron distributions.In the M s = 10, θ Bn = 75 • case (Figure 8, bottom left), no Liouville mapping is performed for outgoing electrons because too few upstream electrons have HT-frame p ∥ > 0.
In Figure 8, we have omitted the near-luminal θ Bn = 85 • shocks.In these shocks, it is less meaningful to describe electron heating in terms of a global, 1D parallel potential.Conduction along x becomes comparable to, or slower than, bulk advection of flux tubes into the shock.And, we have not accounted for finite-Larmor radius drifts in the electron heating model (Goodrich & Scudder 1984).Electrons will still gain energy by parallel electric fields in sufficiently strong shocks, but this parallel energization becomes a local process that is not well described by a 1D, y-averaged potential (Tran & Sironi 2020).
Based on Figures 7 and 8, we adopt ϕ ∥,amb as our preferred estimator for ϕ ∥ in the rest of this manuscript, as it can provide a reasonable estimate of ϕ ∥ from a single time snapshot.
6. CROSS-SHOCK POTENTIAL SCALING WITH SHOCK PARAMETERS

Parallel Potential Scaling
What is the relationship between ϕ ∥ and T e as a function of shock parameters?If we assume isotropic electrons with a polytropic equation of state P e ∝ n Γ e , where Γ is an effective adiabatic index, Equations ( 11) and ( 12) simplify to The effective adiabatic index Γ should reflect the underlying kinetic physics and simulation setup (e.g., a local planar shock versus a global magnetosphere model), so what is its value?
We measure the cross-shock potential and electron temperature jumps in the post-shock region x − x sh = −18 to −3d i as space-and density-weighted averages respectively (Appendix A).Besides ϕ ∥ , we also measure the normal incidence frame (NIF) potential ϕ NIF = E x dx as a commensal "add-on"; ϕ NIF is not directly related to our discussion of electron heating, but it is convenient to also present here.
Figure 9 suggests that both 1D and 2D simulations show an effective Γ ≈ 5/3 relation between ∆ϕ ∥ and ∆T e , although the stronger shocks may deviate towards a higher effective Γ.The 1D simulations can attain higher T e and ∆ϕ ∥ , yet the effective Γ = 5/3 remains similar to our 2D simulations.We caution that the effective Γ for post-shock electrons could be sensitive to our choice of domain left-side boundary condition, which may modify the post-shock electron isotropization.

Normal Incidence Frame Potential Scaling
We further consider how the normal incidence frame (NIF) cross-shock potential, ϕ NIF = E x dx, scales with M ms and θ Bn in our simulations.Although less relevant to electron heating, ϕ NIF is still important to a kinetic shock's internal structure (Burgess & Scholer 2015).The NIF cross-shock potential's dependence on Mach number has also been studied in prior hybrid (Leroy et al. 1982;Quest 1986) and PIC (Shimada & Hoshino 2005) shock simulations.
Our results are shown in the right column of Figure 10, which compares ∆ϕ NIF to a low-M A model from Gedalin (1996, Equation ( 40)) (and Gedalin (1997, Equation ( 28)); Jones & Ellison (1991, Sec. 5.1)), Equation ( 16) comes from integrating E x of generalized Ohm's Law across the shock, assuming scalar P e ∝ n 2 e and B ∝ n e .The upstream β e = β p /2 in our simulations.To evaluate the compression ratio r in Equation ( 16), we use the oblique MHD shock jump conditions, similar to Bale et al. (2008).Our calculation of r introduces a weak θ Bn dependence shown by the closely-spaced curves of varying color in Figure 10; an oblique θ Bn = 55 • decreases the model prediction by ≲ 10% compared to the nearly-perpendicular θ Bn = 85 • .
Figure 10 shows that ∆ϕ NIF is of order one-half the incident ion bulk kinetic energy, 0.5m i u 2 sh , as is well understood (Leroy et al. 1982;Leroy 1983;Burgess & Scholer 2015).Both 2D and 1D shocks show that ∆ϕ NIF varies with θ Bn , with more oblique θ Bn ∼ 55 • yielding higher ∆ϕ NIF than more perpendicular θ Bn .The 1D shocks also show stronger variation in ∆ϕ NIF with θ Bn .We observe that going from M s = 3 to 4, the potential increases for all θ Bn considered.Towards higher M s , the potential appears to level off or decrease slightly; the Gedalin (1997) model does not apply for these stronger shocks.

CASE STUDIES OF ELECTRON DYNAMICS
WITHIN THE PARALLEL POTENTIAL

Electron Phase Space
What electron distributions form within and generate the parallel potential?We focus upon two shock cases that are weak (M A < M w , sub-critical) or strong (M A > M w , super-critical) and therefore have quite different structure.The weaker case is the M s = 4, θ Bn = 65 • shock previously shown in Figures 1 and 7.The stronger case is the M s = 7, θ Bn = 65 • shock.Within Section 7, we refer to each case by its sonic Mach number M s alone.
Figure 11 surveys electron phase space evolution along the shock-normal coordinate for our 2D M s = 4 case, in the manner of Scudder et al. (1986b) and Chen et al. (2018).We show electron momenta γβ in units of m e c, where γ is the relativistic Lorentz factor and β is three-velocity in units of c.The momentum coordinates ( n∥ , n⊥1 , n⊥2 ) are an orthogonal, right-handed triad defined by n∥ = b, n⊥1 ∝ ( b × V e ) × b, and n⊥2 ∝ b × V e ; here, V e is local electron bulk three-velocity, and b is computed locally at each particle position.This momentum coordinate system follows Chen et al. (2018).By construction, γβ ⊥1 includes both bulk motion and Larmor gyration, and γβ ⊥2 includes only Larmor gyration.
Figure 11(a) shows that ion density oscillates in sync with magnetic field strength, as expected for the compressible, fast (whistler-branch) mode.In panel (b), we adopt ϕ ∥,amb as our proxy for the parallel potential.Panels (c)-(r) show how electron momentum distributions evolve through the shock, with momenta measured in the simulation frame.The distributions are thermal Maxwellians far upstream of the shock, at x − x sh = 40-45d i (panels (j),(r)).Within the shock precursor, the distributions become asymmetric in γβ ∥ ; we see a beam-like component with γβ ∥ < 0 on top of a broader,  "flat-top" component; the asymmetric form is most prominent in strong magnetic troughs with δB/B 0 ∼0.1 (panels (m),(o)).Distributions within magnetic peaks are closer to a single component in form, but with a skew towards negative γβ ∥ (panels (n),(p)).Downstream of the shock, the beam-like component disappears and the distribution shows a smooth, flat region in γβ ∥ centered near zero (panels (c),(k)).
In some 2D momentum distributions (panels (e),(g)), particles with γβ ∥ > 0 show a hint of anisotropy, with some particles at γβ ∥ ∼ 0.25 having large γβ ⊥2 ≳ γβ ∥ ; the right side of the 2D distributions is stretched along the γβ ⊥2 axis.We expect that these back-streaming electrons were magnetic mirror-reflected within the shock ramp or precursor.Backstreaming electrons are also visible in the far-upstream 1D distributions (panels (p),(q),(r)) when plotted with a logarithmic scale.
The evolution of parallel distributions through the shock, showing a transient beam and flat-top structure, is broadly consistent with satellite observations (Feldman et al. 1983b,a) and prior 1D and 2D PIC simulations (Savoini & Lembege 1994).It can be noted that Feldman et al. (1983a) showed that weak, interplanetary shocks with M s ∼ 1-4 exhibit less beam/flat-top structure as compared to the stronger bow shock at Earth's magnetosphere.Although our M s = 4 case has M A < M w and is sub-critical, it shows parallel electron behavior similar to space measurements of stronger shocks.This is partly an artifact of the reduced mass ratio adopted for our simulations; we return to this point in Section 7.2.
Electron distributions in the M s = 7 case (Figure 12) show similar distortions as the M s = 4 case: namely, a onesided beam with γβ ∥ < 0 accelerating towards the shock and eroding through the shock ramp to leave a broad, flattened post-shock distribution in γβ ∥ .The population of backstreaming electrons is more prominent than in the M s = 4 case.
Figure 13 shows how the HT-frame γβ ∥ electron distributions correlate with the potential ϕ ∥ for the M s = 4 and M s = 7 case studies in 1D and 2D.In Figure 13(e-h), the blue and black solid curves show ϕ ∥ recast as a parallel velocity ± 2eϕ ∥ /m e , and black dashed lines show the approximate HT-frame downstream bulk velocity, In Figure 13(i-l), trapped electrons are defined to have total energy ε < 0 (Equation ( 2)), while backstreaming electrons have ε > 0 and v ∥ > 0. The electron fractions in Figure 13(il) are computed within 0.1d i bins along x using the 1D yaveraged ϕ ∥ (x) = ϕ ∥,amb (x).Some downstream electrons with v ∥ < 0 are not counted as trapped in Figure 13, but most such electrons are de facto trapped because reflection at the left boundary or scattering within the downstream rest frame generally will not give them ε > 0 and v ∥ > 0 as needed to escape back upstream.We may understand an electron's parallel velocity evolution by casting Equations ( 2) and (3) into dimensionless form: with v ∥0 and v ⊥0 being the electron's initial HT-frame velocity components far upstream, where ϕ ∥ (x) = 0 and B(x) = B 0 .The upstream electron thermal speed c se = Γk B T 0 /m e , and the normalized φ∥ = eϕ ∥ /(0.5m i u 2 sh ) is typically 0.1-0.2 in our simulated shocks (Figure 10).Equation (18) shows that electrons gain parallel energy from ϕ ∥ and lose energy from magnetic mirroring, and M 2 s moderates the relative importance of ϕ ∥ versus mirroring.
Let us trace how electrons evolve within the potential ϕ ∥ and the corresponding velocity floor/barrier ± 2eϕ ∥ /m e , neglecting scattering.Suppose an electron enters the shock with v ∥0 < 0 and v ⊥0 = 0 (no magnetic mirroring) in the HT frame; further suppose that |v ∥0 | ≪ 2eϕ ∥ /m e .That electron will follow the negative potential floor v ∥ (x) = − 2eϕ ∥ (x)/m e into the shock.It reflects at the left boundary by reversing the sign of v x in the simulation frame, i.e. the downstream rest frame.The post-shock v ∥ distribution is thus centered on the HT bulk flow velocity v ∥ = −u HT2 (Equation ( 17)), shown by dashed horizontal lines in Figure 13(e-h).Post-shock electrons may escape upstream if they have v ∥ ≳ 2eϕ ∥ /m e (more precisely, ε > 0 for non-zero pitch angle).The fraction of downstream electrons that may escape upstream is ≲ 10% of the total downstream population in all case-study shocks (Figure 13(i-l)); most shocked electrons are confined to x − x sh < 0, consistent with the prediction of Figure 2.
If an electron entering the shock instead has non-zero v ⊥ and a large-enough pitch angle, it may mirror reflect specularly within the HT frame (i.e., v ∥ → −v ∥ ).If the macroscopic B(x) and ϕ ∥ (x) are time-stationary and there is no scattering, the reflected electron will freely escape back upstream with ε > 0 and v ∥ > 0. In practice, reflection combined with scattering may lead to local electron trapping and energization within the precursor wave train and shock ramp (Katou & Amano 2019).
In each panel of Figure 13, the potential floor (solid blue and black lines with v ∥ < 0) tracks the deformation of the upstream thermal electron beam entering the shock from right to left.In the post-shock region, the potential floor roughly corresponds to the anti-parallel (v ∥ < 0) edge of the electron distribution.The agreement of the potential floor and the electron distribution edge in Figure 13 corresponds to good agreement at p ∥ < 0 in the Liouville mapping procedure of Figure 6.
Let us now compare 1D versus 2D for the M s = 4 case in Figure 13.The overall structure of ϕ ∥ is similar in both 1D and 2D, with a gradual rise in the shock precursor region and jump at the shock ramp to a final post-shock value of ∆ φ∥ ≈ 0.1.The 1D shock shows sharper d i -scale potential spikes embedded within the precursor wave train; in both 1D and 2D the spikes occur within magnetic troughs.Electrons clump at local magnetic maxima, especially in 1D, which we interpret as due to magnetic mirroring.In 1D, ϕ ∥ spikes also occur at magnetic maxima within x − x sh = 0-5d i (Figure 13, top left panel) just ahead of the shock ramp.
The 1D and 2D cases differ in the post-shock region x < x sh for both the M s = 4 and 7 cases.The 2D-shock distributions are diffuse and smooth in v ∥ , whereas the 1D-shock distributions show localized (in x) beams or clumps that may correspond to transient phase space holes.The 1D case also shows larger precursor fluctuations in the trapped and untrapped electron fractions (and, larger magnetic and ϕ ∥ fluctuations) than the 2D case.It is possible that similarly strong fluctuations occur in 2D but are hidden by y-averaging; however, inspection of B/B 0 images suggests that the 2D precursor fluctuations are coherent and so of lower amplitude than in 1D.
In the M s = 7 case, the 1D jump in ϕ ∥ is ∼2× larger than the equivalent 2D shock.The 2D shock shows more v ∥ diffusion in the shock foot as electrons approach the ramp (x − x sh ∼ 1d i ).We suggest that strong, non-adiabatic scattering embedded within the shock ramp drives so-called "infilling" of the parallel distribution (Hull et al. 1998) and hence a net cooling, as compared to the 1D case.

Effects of Electron Plasma-Cyclotron Frequency Ratio and Ion-Electron Mass Ratio
Let us now explore the effects of the electron plasma-tocyclotron frequency ratio ω pe /Ω ce and the mass ratio m i /m e upon the parallel potential, for fixed Mach number and β p .In this Subsection, we focus solely on the 1D M s = 4 case for two reasons.First, the extended ϕ ∥ structure over many d i in the M s = 4 shock precursor has not been studied before in a quasi-perpendicular shock, to our knowledge.Second, computing cost and numerical noise both rise as either m i /m e or ω pe /Ω ce are raised towards realistic values; the use of 1D simulations helps us limit both cost and noise.In this Subsection (Figures 14, 15) all simulations use 8192 upstream particles per cell, four times larger than our standard runs (Section 2.2).
Figure 14 shows that ω pe /Ω ce does not have a strong effect on φ∥ or the electron phase-space behavior, which we interpret as follows.Suppose that ϕ ∥ scales with some electrostatic fluctuation strength, δE ∥ , and is integrated over an ion skin depth: When normalized to ion kinetic energy, Let us now suppose that δE ∥ scales like a whistler wave's electromagnetic fluctuation δE ⊥ , which will be motional for low frequencies: δE ⊥ ∼ δuδB/c (where δu and δB are the wave's velocity and magnetic fluctuations).Then The fluctuations δu/v A and δB/B 0 should not depend on u sh /c for the non-relativistic solar wind (Verscharen et al. 2020).Then, In this scaling, we find that φ∥ shows no explicit dependence on u sh /c or ω pe /Ω ce .Our heuristic argument has caveats.The whistler wave strength δB/B 0 > 0.1 is outside the regime of linear fluctuations, especially at the shock ramp (δB/B 0 ∼ 1).And, the supposed proportionality δE ∥ ∼ δE ⊥ is suspect; the strength of δE ∥ could be regulated by secondary electrostatic modes that could introduce m i /m e or ω pe /Ω ce dependence into the scaling argument.
Figure 15 shows with increasing mass ratio m i /m e : ϕ ∥ decreases, the electron v ∥ phase space is less distorted, and magnetic mirroring from µ conservation becomes more important.Both downstream and upstream distributions are more centered on v ∥ = 0 with increasing m i /m e .Is the changed electron response due solely to the decrease in ϕ ∥ ?Or, would the electron response change even if ϕ ∥ (x) were identical for both m i /m e = 200 and 1836?
Suppose that we have two shocks with the same M s , φ∥ (x), and B(x) profiles, and the shocks differ only in their mass ratio m i /m e .The right-hand side of Equation ( 18) is not directly affected by m i /m e , since (v ⊥0 /c se ) 2 ∼ 1 in the HT frame, so the relative importance of ϕ ∥ and mag-netic mirroring is unchanged.In other words, the Liouville mapping of a single trajectory starting at any (v ∥0 , v ⊥0 ) does not change with the mass ratio.Instead, m i /m e influences the electron dynamics via v .When the HT-frame bulk upstream velocity u HT = u sh / cos θ Bn is much larger than the electron thermal speed c se , then In the opposite limit u HT ≪ c se , v ∥0 /c se ∼ 1.The mass ratio m i /m e shifts the distribution along v ∥0 /c se in the HT frame.Although m i /m e does not change electron trajectories and separatrices in phase space, which are fully determined by φ∥ (x) and B(x) (Equation ( 18)), m i /m e does change the initial sampling of said trajectories, and hence the relative fraction of electrons that pass into the shock downstream as opposed to being reflected by magnetic mirroring.
Figure 16 summarizes the just-preceding discussion by showing the test-particle electron Liouville mapping for two idealized shocks having identical M s , φ∥ (x), and B(x) and differing only in mass ratio m i /m e ; see, e.g., Yuan et al. (2008) for a more sophisticated example of such test-particle mapping.The Liouville map is performed following Equations ( 5) to ( 8).For a fixed starting point (v ⊥0 , v ∥0 ), the Liouville mapping (phase space flow) from pre-to post-shock is identical and has no dependence upon mass ratio, when specified in terms of upstream electron thermal velocity c se .But, the starting electron distribution is offset farther from v ∥0 = 0 in the m i /m e = 200 case; this leads to a stronger v ∥ < 0 beam and a more asymmetric post-shock distribution in the HT frame, as compared to the m i /m e = 1836 case.
It is not a good assumption that φ∥ is the same for two shocks of different m i /m e .Indeed, φ∥ decreases with m i /m e in our 1D M s = 4 example by a factor of 2 between m i /m e = 200 to 400, and then appears unchanged for higher m i /m e .Nevertheless, as we have just shown, it is helpful to separate the effects of u HT /c se and φ∥ upon electron energization.

Precursor Wave Properties
In this Subsection, we verify that the precursor waves in our 2D M s = 4 and M s = 7 shocks are upstreampropagating whistler waves, oblique to both B and n, with phase and group speeds exceeding the shock speed.It is already well established that weak, low-beta, quasiperpendicular solar wind shocks are formed of, regulated by, and radiate low frequency (ω ≪ Ω ce ) whistler-branch modes; evidence is given by theory (Bickerton et al. 1971;Tidman & Krall 1971;Krasnoselskikh et al. 2002), experiments (Robson & Sheffield 1969), observations (Fairfield & Feldman 1975;Greenstadt et al. 1975;Mellott & Greenstadt 1984;Oka et al. 2006;Hull et al. 2012;Wilson et al. 2012;Wilson 2016;Wilson et al. 2017;Oka et al. 2019;Hull et al. 2020;Lalti et al. 2022), and fully-kinetic simulations (Liewer et al. 1991;Savoini & Lembege 1994, 2010; Riquelme &  18).Red line bounds the region wherein electrons will mirror reflect upon encountering the shock.Bottom row: downstream electron distribution in (v ∥ , v ⊥ ).The colormap shows phase-space density in arbitrary units, but matched in all panels to allow quantitative comparison.
Spitkovsky 2011).The goal is to frame our case-study shocks within a broader context.
Figure 17 shows the 2D structure of our case study shocks.Both M s = 4 and 7 cases show electromagnetic precursors oblique to B (Figure 17(a-f)).The electron temperature map shows hot filaments at x − x sh = 0 in the M s = 4 case (Figure 17(g-h)), some of which appear connected to strong bipolar electrostatic structures at y = 2d i and 5d i within the shock ramp (Figure 17(i-l)).Small-scale (< d i ) electrostatic structures appear in E ∥ , or equivalently in b • E amb , and they inhabit wave troughs of local B minima within the shock precursor region x > x sh (Figure 17(i-l)); many such structures are bipolar electron holes (positive electric potential).
Figure 18 presents Fourier power spectra of B z within the shock precursor regions in order to identify the ion-scale precursor wave modes and to measure their propagation angles and phase speeds.We choose B z to include all electromagnetic modes in the 2D simulation domain.We sample B z at 0.25Ω −1 ci intervals for a duration 5Ω −1 ci towards each simulation's end.The spatial grid of B z used for analysis is 6× downsampled in all directions.The spatial windows are x − x sh = 2-45d i and 2-15d i in the M s = 4 and M s = 7 cases respectively, at end of each simulation.The spatial window is stationary in the simulation frame; at earlier time snapshots, the window is farther from the shock, and the precursor gradually advances into the spatial window over 5Ω −1 ci .Before computing the Fourier transform, we apply a Blackman-Harris window to both x and t coordinates to reduce spectral leakage.The frequencies and wavenumbers in Figure 18 are presented in the simulation frame, and θ Bk is the angle between wavevector k and the upstream magnetic field at θ Bn = 65 • .The plotted frequency range is mildly asymmetric because we use an even number of time snapshots, so Nyquist-frequency power is assigned to the ω = +12.6Ω−1 ci bin.In Figure 18(c-d), we plot an approximate cold whistler dispersion relation from Krasnoselskikh et al. (1985Krasnoselskikh et al. ( , 2002)): with ω in the plasma rest frame.We plot Equation ( 20) with a Doppler shift into the simulation frame, ω → ω − k x u 0 ; the simulation frame's upstream flow is also shown as ω = −k x u 0 (dotted white line).Due to our coarse time sampling, wave power can alias in frequency ω, and we account for this by also aliasing one of the dispersion curves (Figure 18(d), cyan dashed).
Wave phase velocities may be computed from Figure 18  as: with upstream v A /c = 0.028 and 0.020 in our M s = 4 and M s = 7 shocks respectively; θ kn is the angle between between k and shock normal n = x.Equation ( 21) is only a unit conversion and vector projection and does not use Equation ( 20).
The M s = 4 shock shows two coherent precursor wave trains (Figure 18(a),(c)).The strongest precursor travels along the shock normal n; a second precursor travels at an oblique angle ∼35 • with respect to both B and n.Both wave trains have projected phase speed ω/k x near or above the shock speed ω/k x = u sh /r = 0.028c (Figure 18(c), dashed white line).The shock normal-aligned train has ω = 2.5Ω ci , k x = 0.21d −1 e , and phase speed ω/k = ω/k x = 0.024c.The oblique wave train has ω = 3.7Ω ci , k x = 0.14d −1 e , and phase speed ω/k = (ω/k x ) cos θ kn = 0.046c.We take ω, k x , and k y to be the location of 1D power spectrum maxima within slices of the 2D (ω, k x ) and (k x , k y ) spectra.Measurement uncertainty comes from the coarsest Fourier transform bins, ∆k y = 0.08d −1 e and ∆ω = 1.3Ω −1 ci , so our phase speed estimates are only good to tens of percent.Nevertheless, we conclude that the shock normalaligned train's ω/k ≈ u sh /r permits it to phase stand in the shock frame.The oblique train's phase-velocity vector has shock-normal ( n = x) projected component equal to (ω/k) cos θ kn = (ω/k x ) cos 2 θ kn = 0.040c ≳ u sh /r, permitting it to co-move with or out-run the shock.
The M s = 7 shock shows multiple modes that we attribute to forward-propagating whistlers with propagation angles θ Bk ∼ 0-40 • .In the k x -k y spectrum (Figure 18 power within k x ∼ 0.1-0.4d−1 e , at all frequencies within the Nyquist-limited band −11.3 to 12.6Ω ci .We interpret these blobs as frequency-aliased whistler modes at near-parallel θ Bk ≲ 15 • propagation angles (Figure 18(d)).The k xk y power spectrum in conjunction with Equation (20) helps break the frequency-aliasing degeneracy and supports our identification of the waves as whistlers.As an example, the aliased dispersion for θ Bk = 15 • (dashed cyan) crosses two blobs of negative frequency (gray arrows) in Figure 18(d), which we attribute to k x and k y of the same θ Bk in the k xk y spectrum (Figure 18(b), gray arrows).The angle θ Bk is uncertain; the same power could be explained by, e.g., doubly-aliased θ Bk = 0 • whistlers.We cannot cleanly identify frequencies and wavevector angles due to both aliasing and the coarse frequency-and wavenumber-space resolution.But, we can conclude that the supercritical M s = 7 shock hosts oblique (θ Bk ≲ 40 • ) whistler modes with simulationframe frequencies ω ≳ 5Ω ci around and above the lower hybrid range.For the θ Bk = 25 • and 40 • modes, we estimate that their shock normal-projected phase velocities (ω/k) cos θ kn = (ω/k x ) cos 2 θ kn ≈ 0.026c and 0.036c respectively, at or above the shock speed u sh /r = 0.023c in the simulation frame.
Figure 19 shows that the precursor waves are right-hand polarized, for both shock-normal aligned and oblique modes, in both the M s = 4 and M s = 7 2D shocks.Following Stix (1992, Ch. 1), polarization is defined by the rotation sense of a wave's electric field about its background (upstream) magnetic field B 0 at a fixed point in space; counterclockwise rotation about B 0 is right-handed polarization.We project E along the right-handed coordinate triad of unit vectors ( n⊥1 , n⊥2 , b0 ), where b0 is the upstream magnetic field for M s = 7 (Fig. 18(b)); these modes have either one or two standing wavelengths along y.In the M s = 4 case, we see that the n-parallel precursor wave has higher amplitude and longer x extent than the n-oblique wave; in the M s = 7 case, the opposite holds.All four wavetrains show E ⊥1 one quarter cycle ahead of E ⊥2 and are therefore right-hand polarized, consistent with a fast mode having k at acute angle to x (i.e., forward propagating with respect to x).

Precursor Wave Interpretation
What is the origin of the precursor whistlers?A full answer to this question is beyond the scope of our work, but we shall give some comments.In the M s = 4 case, the shock-normal whistlers are consistent with a phase-standing wave train expected to form for M A < M w (Kennel et al. 1985;Krasnoselskikh et al. 2002).For the whistlers oblique to n, we consider two plausible mechanisms: nonlinear wave steepening within the shock ramp (Krasnoselskikh et al. 2002), or beam resonance with reflected ions gyrating within the shock foot, which may be called "modified two stream instability" (MTSI) or "kinetic cross-field streaming instabil- ity" (Wu et al. 1984;Hellinger et al. 1996;Hellinger & Mangeney 1997;Muschietti & Lembège 2017).There also exists a less-studied wave-wave decay instability that may generate daughter whistlers having k mis-aligned with respect to their parent whistlers (Galeev & Karpman 1963;Decker & Robson 1972), noted by Mellott & Greenstadt (1984, pg. 2158), but we do not consider it for now.Lower hybrid drift instabilities have been invoked to explain shock structure and electron heating (Krall & Liewer 1971;Wu et al. 1984;Stasiewicz & Eliasson 2020); but, drift waves propagate orthogonal to both B and ∇P e (i.e., n) and are thus suppressed in our 2D simulations.We also do not consider "slow" MTSI caused by a relative drift between incoming ions and electrons (Muschietti & Lembège 2017).Both lower hybrid drifts and "slow" MTSI are expected to drive waves nearly perpendicular to B that are mainly electrostatic and so would not explain the oblique, electromagnetic power that we are now considering.
Both M s = 4 and 7 shocks have M A below the non-linear whistler Mach number √ 2M w as defined by Krasnoselskikh et al. (2002), and Figure 18 shows that whistlers in range of θ Bk can outrun both shocks.Of note, our M s = 7 shock having M w < M A < √ 2M w is similar to the Cluster shock crossing presented by Dimmock et al. (2019), wherein d escale shock ramp fluctuation was presented as evidence for nonlinear steepening as the origin of oblique whistlers.More broadly speaking, it seems plausible that 2D shock ramp rippling or filamentation, coupled with wave steepening in the ramp, may cause whistler wave emission at various angles (cf.Sundkvist et al. 2012).
Can reflected ions resonate with the oblique precursor waves?We assess this by comparing the reflected ions' xy plane velocity β xy = v xy /c against the precursor whistler phase speeds (Equation ( 21)), all measured in the simulation frame.We separate reflected ions from "core" ions using an x-z momentum-space threshold γβ z > 0.02 − 0.5γβ x , which selects ions with high v x and v z being accelerated by the upstream motional field near the shock ramp.The same momentum-space cut is applied for both M s = 4 and M s = 7.We compute the mean and standard deviation of β xy within 0.5d i bins along x, noting that the reflected ion distributions are not Gaussian or particularly symmetric in phase space.
Figures 20 and 21 show the ion distributions and measurements of β xy .Reflected ions appear within ±1d i of x sh for both the M s = 4 and 7 shocks, with a shorter shock-foot gyration (≲ 0.5d i along x) for the M s = 4 shock.In the last simulation snapshot, the reflected ions have a number density contrast ∼ 3-8% (M s = 4) and ∼20% (M s = 7) with respect to the total ion population.
The M s = 4 shock's reflected ion beam has mean v xy ∼0.02-0.026c(standard deviation 0.01c).It appears difficult for the reflected ions to resonate with the θ Bk = 35 • whistler mode at ω/k ∼ 0.046c, but the ions could (in principle) lie in beam resonance with the shock-normal θ Bk = 65 • whistler at ω/k ∼ 0.024c.
The M s = 7 shock's reflected ion beam has mean v xy ∼0.03-0.04c(standard deviation 0.02c).The reflected ions might resonate with the lower-frequency θ Bk = 40 • and 25 • whistler modes having ω/k ∼ 0.033 and 0.04c respectively.The reflected ions may be too slow to resonate with the higher-frequency θ Bk ≲ 15 • whistlers, which have ω/k ≳ 0.05c.
Our analysis is limited because we do not consider any of (1) time or spatial (2D) variation in the reflected ion properties, (2) angle between reflected ions' velocity vector and the precursor wavevectors, (3) energy or power balance, (4) more detailed kinetic stability analysis.A proper treatment of (1) may loosen our constraints if time-intermittent or localized, gyrophase-bunched ion reflection can attain a larger range of beam speeds and higher density with respect to the incoming flow.A proper accounting of (2-4) should help tighten our constraints by culling observed wave modes that could not plausibly be driven by a resonant ion beam.Now, (1-4) are not unusually difficult to assess, but we emphasize that the purpose of Sections 7.3 and 7.4 is to provide broader context for our primary focus of electron heating physics.We can at least conclude that reflected ions' beam resonance likely does not explain all of the oblique wave precursor power observed in our simulations, and perhaps nonlinear shock ramp processes are needed.
We warn that the propagation angle θ Bk of our whistlers is subject to both a mode selection effect from the simulation domain's periodic boundary in y, as well as angle-dependent Landau and transit-time damping that may differ between our simulated m i /m e = 200 and the true m i /m e = 1836.We also warn that our angles θ Bk are by-eye estimates, rather than being obtained by a formal fitting procedure; we deem this acceptable since neither precise wave angle determination nor wave driving mechanism determination is the primary purpose of this manuscript, and there is considerable systematic uncertainty due to the frequency aliasing in our M s = 7 shock analysis anyways.
Our most robust conclusion is that forward-propagating precursor whistlers appear at a variety of propagation angles between θ Bk = 0 • to θ Bk = θ Bn in our case study shocks.These whistlers are important in setting the ion-scale structure of the electron's parallel potential.We expect such precursor whistlers to occur in a wide range of simulated parameter space as fully-kinetic simulations catch up to decades of spacecraft observations.8. CONCLUSIONS 8.1.Summary We have measured electrons' ambipolar B-parallel potential ϕ ∥ in a survey of 2D quasi-perpendicular PIC shocks spanning M s = 3-10 (M A ∼ 1-5) and θ Bn = 85-55 • , with upstream total plasma beta β p = 0.25.Different particle-and field-based measurements of ϕ ∥ agree in the overall ϕ ∥ jump at the ∼10% level, and we find that ϕ ∥ is most robustly estimated in our simulations from the electron pressure tensor divergence.Off-diagonal terms of the electron pressure tensor P e are not negligible for our lowβ p shock simulations; the integral of the ambipolar electric field E amb projected along shock normal can underestimate ϕ ∥ by ∼10-30% as compared to other measurement methods.Both the magnitude of ϕ ∥ , and the correlation between ϕ ∥ and ∆T e , are similar to prior reports based on observational data (Schwartz et al. 1988;Hull et al. 1998Hull et al. , 2000)).We also measure the normal incidence frame (NIF) potential ϕ NIF and show its variation with θ Bn and Mach number.
We show that a quasi-perpendicular shock with Alfvén Mach number M A below a critical whistler Mach number, M A < M w , can host a ϕ ∥ (x) profile extending over many tens of d i within a whistler wave precursor.The potential shows transient spikes within magnetic troughs, as well as a secular increase towards the shock over many wave cycles.We speculate that the potential spikes (bipolar electric fields) are due to non-linear steepening of the large-amplitude whistler wave's electrostatic field (Vasko et al. 2018a;An et al. 2019), but more work is needed to properly explore and test this hypothesis.
In two case-study shocks (m i /m e = 200, M s = 4 and 7, θ Bn = 65 • ) without backstreaming ions, we find that shock precursor whistlers span various propagation angles between θ Bk = 0 • to θ Bk = θ Bn , including a phase-standing wave train along shock-normal in the M s = 4 case (Tidman & Northrop 1968;Kennel et al. 1985).The whistler precursors also host small-scale electrostatic structures, including electron holes, with wavevector k parallel to B. We tentatively find that shock-reflected (gyrating, not backstreaming) ions may lie in beam resonance with some, but not all, of the precursor wave power.Other mechanisms, such as nonlinear steepening within the shock ramp (Omidi & Winske 1988, 1990;Krasnoselskikh et al. 2002) or wave-wave decay instability (Galeev & Karpman 1963;Decker & Robson 1972), may help explain some precursor wave power that cannot be attributed to beam ions.

Future Directions
Our description of electron energization has been macroscopic, focusing upon electric fields at ion scales; we have not quantified the microscopic scattering and dissipation that is needed to provide true irreversible heating and to regulate the strength of ϕ ∥ via the electron pressure tensor.Moreover, the case-study simulations of Section 7 have electron holes as the main electrostatic structure, whereas MMS observations have revealed electron and ion holes, double layers, ion acoustic waves, and electron cyclotron harmonics (Goodrich et al. 2018;Wang et al. 2021).Towards higher β p , small-scale parallel whistlers will likely become important (Hull et al. 2012;Page et al. 2021).Much recent work has been done with MMS and some tailored simulations (Goodrich et al. 2018;Vasko et al. 2018b;Wang et al. 2021;Shen et al. 2021;Kamaletdinov et al. 2021;Sun et al. 2022).Fewer fully-kinetic shock simulations have captured the macroscopic, ion-scale coupling to microscopic electrostatic scales (Wilson et al. 2021).In such shock simulations, we would like to know: what types of electrostatic structures appear, and where in the shock do said structures appear, as a function of shock parameters M s , θ Bn , and β p ?What is their contribution to B-parallel scattering and pitch-angle scattering of electrons?As the latter two questions are answered, can we find any qualitative or quantitative trends in shock parameter space to test against MMS data?
A few caveats must be noted for further study of these secondary electrostatic structures.First, our simulations are limited by spatial resolution and numerical noise.The upstream electron Debye length is marginally resolved with one cell, and not resolved when PIC current filtering is accounted for; small-scale structures may be biased towards lower k.Discrete PIC macroparticles create numerical scattering (Birdsall & Langdon 1991) whose contribution needs to be separated from that of more-physical collisionless scattering.High mass ratio m i /m e , high particle sampling, and high space and time resolution may be needed to unveil physics at the electron Debye scale in future simulations; Vlasov solvers may be another way to make progress (e.g.Juno et al. 2018).Second, we have only used equal-temperature ionelectron Maxwellian distributions to model upstream plasma in our shocks.Solar wind shocks have multiple populations: core, halo, and strahl electrons, pick-up ions, et cetera.The electron halo and strahl, in particular, may alter Landau damping rates.Third, artificially low ω pe /Ω ce can modify the types and Fourier spectra of electrostatic waves, e.g., for the electron cyclotron drift instability (ECDI) (Muschietti & Lembège 2013;TenBarge et al. 2021), and hence also modify the energy exchange between ions and electrons.Lastly, the use of 2D versus 3D simulations will matter -just as the electron trapping is weaker in 2D versus 1D, so we also expect that 3D simulations may reduce electron trapping, admit a broader spectrum of precursor wave modes, and permit electrons to scatter and gain energy in different and possibly new ways (Trotta & Burgess 2019).
We are grateful to anonymous referees whose comments have helped to stimulate much of this study.We acknowledge helpful conversations, brief or extended, with Takanobu Amano, Artem Bohdan, Collin Brown, Li-Jen Chen, Luca Comisso, Greg Howes, Arthur Hull, Jimmy Juno, Matthew Kunz, Steve Schwartz, Navin Sridhar, Jason TenBarge, Ivan Vasko, Lynn B. Wilson III, participants and organizers of the Geospace Environment Modeling (GEM) Focus Group "Particle Heating and Thermalization in Collisionless Shocks in the MMS Era", and many other members of the scientific community.We especially thank Lynn for alerting us to a Liouville mapping normalization error in an earlier version of this manuscript.We are grateful for the use of analysis and simulation code developed by Xinyi Guo.
AT Figure 22 shows the 1D, y-averaged profiles of φNIF , φ∥ , ∆T e , and ∆T i for our full set of 2D and 1D shocks.All potentials and temperature profiles are computed using 50 finely-spaced snapshots following the alignment and averaging procedure as described in Section 5.4.We compute the volume-averaged mean of ϕ ∥ and ϕ NIF to obtain the singlepoint jumps ∆ϕ ∥ and ∆ϕ NIF plotted in Figures 9-10.The mean electron and ion temperatures are weighted by their corresponding species density.We compute the mean ϕ ∥ , ϕ NIF , and ∆T e in the region x − x sh = −18 to −3d i , highlighted by light green in Figure 22.The normal-incidence frame (NIF) potential ϕ NIF is measured in the simulation frame, since the boost to the NIF frame does not alter E x .
B. BEST-FIT LIOUVILLE-MAPPED DISTRIBUTIONS In Figure 23, we show Liouville-mapped distributions for many shocks in our parameter sweep.Figure 23 is structured similarly to Figure 6, except that it shows only one Liouvillemapped distribution for the mean best-fit ∆ φ∥ value on each side of p ∥ = 0. We observe that weaker, more oblique shocks (upper right of Figure 23) show more adiabatic electron behavior, while stronger and more perpendicular shocks are not well described by Liouville mapping for incoming (p ∥ < 0) electrons.

C. ELECTRON KINETIC ENERGY IN THE HT FRAME
Figure 24 shows electron kinetic energy, multiplied by sgn(v ∥ ), as an alternative to x-v ∥ phase space in Figure 13.The electron kinetic energy so plotted has a direct mapping to the "Trapped" and "Backstreaming" fractions in

D. NUMERICAL CONVERGENCE OF SHOCK STRUCTURE
We use our case study shocks to assess numerical convergence with respect to the upstream number of particles per cell (PPC), both ions and electrons.Recall that we adopted 2048 and 128 PPC for 1D and 2D simulations respectively.
The overall shock structure is well converged for all of our case studies (Figure 25(a)-(d)).In 2D, the shock ramp and precursor waves are not matched in phase for simulations with different PPC; the 1D simulations, constrained to a single wave train along x, are more coherent.The electron energy gain, as measured by both T e and ϕ ∥,amb , appears converged to within 10% (Figure 25(e)-(l)).For the 2D M s = 7 shock, the upstream to downstream jumps in T e and ϕ ∥,amb do not vary monotonically with PPC, which we interpret as statistical fluctuation rather than a lack of convergence.
The trapped and backstreaming electron fractions (Figure 25(m)-(t)) also appear converged in all of our case study shocks.The trapped fraction, in particular, serves as a coarse proxy for local phase mixing and scattering of both v ∥ and pitch angle within the whistler precursor.The electrostatic power in numerical noise is not small, being within an order of magnitude of the power in electron holes and smallscale structures (Figure 17(i)-(l)).It is not clear if the relative contributions of different particle scattering mechanisms are converged.But, at least, Figure 17(i)-(l) suggests that the total phase-space flow in our shocks is insensitive to PPC sampling.
E. SIMULATION PARAMETERS Table 1 provides simulation input parameters and some derived parameters for the main parameter sweep of 1D and 2D shocks, as well as a set of 1D simulations with varying ω pe /Ω ce and m i /m e presented in Section 7.2.The parameters σ, Θ e , and u 0 are actual code inputs, from which the shock parameters M s , M A , M ms , β p , u sh , u sh /r can be derived assuming a single-fluid adiabatic index Γ = 5/3.The columns of Table 1 not already defined in the main text are as follows.
• tΩ ci is the simulation duration.
• u y is the post-shock transverse drift along y predicted from the R-H conditions and used to set the electromagnetic fields at the simulation's left boundary.
• u right is the manually-chosen expansion speed of the domain's right-side boundary.
(a) shows the structure of an example 2D, M s = 4, θ Bn = 65 • shock.Incoming ions with ⟨v x ⟩ = −u 0 = −0.0238cwere reflected at x = 0 to form a shock now at x ∼ 40 d i , traveling from left to right.The ion v x distribution oscillates ahead of the shock (x ≳ 40 d i ) as part of a precursor whistler wave train.Figure 1(b) shows the bulk ion velocities ⟨v y ⟩ and ⟨v z ⟩ for our example shock; the black dotted line shows the transverse u y deflection predicted by the MHD R-H conditions and also used to set the left-side x = 0 wall's transverse velocity.The transverse ⟨v y ⟩ drift is nearly constant throughout the post-shock region (0 < x < 40d i ). Figure 1(c) shows the shock structure in ion density and ion/electron temperatures.The electron temperature rises ahead of the shock within the same region (x ∼ 40 to 70 d i ) as the ion bulk-velocity oscillations.
(top panel) serves as an interpretive key to the presence or absence of shock precursors in Figures 3 and 4. Let us explore each piece of Figure 2 in turn.

MFigure 2 .
Figure2.Regimes of βp = 0.25 shock behavior in Ms-θBn space for our simulations (top panel) and true mass ratio mi/me = 1836 with solar wind-like T0 = 10 eV (bottom panel).Black dots in top panel mark our simulations; circled black dots are studied in more depth in Section 7. Black shaded region is Mms < 1, calculating MHD fast speed with θBn dependence (unlike definition and use of Mms elsewhere in manuscript; definitions agree at θBn = 90 • ).Blue shaded region is whistler sub-critical, MA < Mw. Green shaded region shows which shocks may permit thermal electrons to escape from downstream to upstream, for varying ∆ φ∥ = 0, 0.05, and 0.1 (Equation (1)).Red dotted line marks where Mms equals the critical Mach number(Marshall 1955;Kennel et al. 1985).Orange dashed line is sub/super-luminal boundary for Θe = 0.01; note that the Ms = 7, 10 and θBn = 85 • simulations are marginally sub-luminal, not super-luminal, because of their lower Θe values.Black dotted contours bound the regions where the SDA reflection efficiency is > 1% for varying ∆ φ∥ = 0, 0.05, and 0.1 (larger to smaller regions).

Figure 3 .
Figure 3. Ion density at end of each 2D shock simulation, with ending time t = 20 to 40 Ωci −1 labeled in each panel.Columns vary magnetic obliquity from θBn = 85 • (left) to 55 • (right); rows vary shock strength from Ms = 3 (top) to 10 (bottom).Faint black lines trace magnetic field lines.See online journal for animated time evolution from t = 0 to end of all simulations.

Figure 5 .
Figure 5. Post-shock electron-ion temperature ratio Te/Ti from our 2D runs as a function of sonic Mach number Ms (left) and magnetosonic Mach number Mms (right).The upstream plasma beta βp = 0.25 for all points.Dotted black line is baseline from adiabatic Γ = 2 compression of electrons, without other ion-electron energy exchange, based on MHD R-H shock jump prediction at θBn = 90 • .
Figure 6.Liouville mapping: by comparing upstream (black curve) and downstream (blue curve) electron distribution cuts along the p ∥ axis, we fit Liouville-mapped upstream distributions (orange curves) to the downstream distributions in order to estimate ∆ φ∥ for different electron populations (incoming p ∥ < 0 versus outgoing p ∥ > 0) and for different bands in f .For incoming p ∥ < 0, the smallest and largest best-fit values of ∆ φ∥ = 0.100 and 0.150 are highlighted in orange alongside their corresponding Liouvillemapped distributions.For outgoing p ∥ > 0, the smallest and largest best-fit values of ∆ φ∥ = 0.045 and 0.073 are similarly highlighted.Momentum is scaled using the upstream electron thermal speed cse = ΓkBT0/me.The shock is the same 2D, Ms = 4, θBn = 65 • case in Figure7.
Here, b x , b y , b z are components of b; notice b x = x• b.Since E HT ≈ E amb , ϕ HT = − 1 en e dP e,xx dx dx .

Figure 7 .
Figure 7. Parallel electrostatic potential measured four different ways in an example 2D, Ms = 4, θBn = 65 • shock, using 50 evenly-spaced snapshots within t = 40.026 to 40.123Ωci −1 .(a) Magnetic field magnitude scaled to its upstream value, B0.(b) Potential ϕ ∥,grid measured from PIC grid E ∥ .(c) Potential ϕ ∥,amb measured from electron pressure tensor divergence, ∇ • Pe.Dotted line is the contribution from x • (∇ • Pe) (Equation (15)), which samples only the shock-normal component of E amb , to show that the non-x components of E amb contribute measurably to ϕ ∥,amb .(d) Potential ϕ ∥,gyr measured from electron pressure assuming gyrotropy (Equation (12)).(e) Potential ϕHT measured from the HT-frame electric field by integrating the shock-normal component EHT,x.In all panels, faint colored lines are single time snapshots; thick lines are averages.Multiple time snapshots appear in panels (c)-(d), but they agree well enough so as to be indistinguishable.

Figure 8 .
Figure 8. Parallel electrostatic potential measured in different ways for 2D shocks.Black solid line is ϕ ∥,grid , Equation (10).Blue solid line is ϕ ∥,amb , Equation (11).Dotted blue line is same as Equation (11), except the integrand is replaced with ⟨ x • (∇ • Pe)/(ene)⟩.Green solid line is ϕ ∥,gyr , Equation (12).Orange dashed line is ϕHT, Equation (13).Shaded regions (gray, orange) represent uncertainty in the grid-based proxies ϕ ∥,grid and ϕHT; region widths are the expected standard error of the mean, ±σ/ √ N with N = 50 for the number of time snapshots and σ the standard deviation, computed at each x position.Red triangles are mean of best-fit ∆ φ∥ from Liouville mapping; left-facing triangle is from incoming (HT-frame v ∥ < 0) electrons, and right-facing triangle is from outgoing (HT-frame v ∥ > 0) electrons.Red vertical bars show the range of ∆ φ∥ values inferred from different f -band fits of the same electron populations; see Figure 6 and Section 5.1.

Figure 11 .Figure 12 .
Figure 11.Electron momentum distribution functions sampled in intervals along x for a 2D, Ms = 4, θBn = 65 • shock at t = 40Ω −1 ci .(a): Ion density ni and total magnetic field strength B normalized to their upstream values, n0 and B0.(b): Parallel potential φ∥,amb along x computed as in Equation (11) and normalized to 0.5miu 2 sh .(c)-(j): Reduced 2D distributions in parallel and perpendicular momenta, γβ ∥ and γβ ⊥2 , normalized so that the histogram integral is unity.The normalization thus does not show variation in electron density ne along the shock.Momentum coordinate system is defined in manuscript text.(k)-(r): Reduced 1D distribution in parallel momentum γβ ∥ .Like in panels (c)-(j), the 1D histogram bins are normalized so that their integral equals one.The electron distributions of (c)-(r) are sampled from gray-shaded x intervals in (a)-(b), shown by light gray lines connecting panels (c)-(j) to regions in (b).

Figure 15 .
Figure 15.Effect of mi/me, ranging from 200 to 1836, upon ϕ ∥ (x) and electron x-v ∥ phase space behavior in the HT frame for 1D Ms = 4, θBn = 65 • shock.The upstream electron temperature Θe = 0.01 for all runs.Top panel: φ∥ with colors indicating mi/me value.Fiducial mi/me = 200 black curve.Bottom panels: electron x-v ∥ phase space, same organization as Figure 13(e), with mi/me increasing as rows descend.Second row from top is fiducial case, like Figure 13(e) but simulated using 4× more particles.

Figure 16 .
Figure 16.Liouville mapping for B/B0 = 1.8 and φ∥ = 0.10, holding Ms = 4, θBn = 65 • fixed.We take 10 5 samples of a Maxwellian, boost into the HT frame, and map electron trajectories from upstream to downstream; all velocities and boosts are assumed non-relativistic.Left column shows reduced mi/me = 200, right column shows true proton-electron mi/me = 1836.Top row: initial upstream distribution in (v ∥ , v ⊥ ), normalized to cse; this corresponds to v ∥0 and v ⊥0 in Equation (18).Red line bounds the region wherein electrons will mirror reflect upon encountering the shock.Bottom row: downstream electron distribution in (v ∥ , v ⊥ ).The colormap shows phase-space density in arbitrary units, but matched in all panels to allow quantitative comparison.
direction, n⊥1 = ẑ, and n⊥2 ∝ ẑ × b0 . 4All precursor wavevectors lie at an acute angle with respect to x (Figure 18(a),(c)), so wave fronts advance toward + x.Therefore, at a fixed time, right-handedness is shown by E ⊥1 offset one quarter cycle ahead of E ⊥2 when plotted as a function of x in Figure 19.We separate the n-parallel and n-oblique modes in order to check their polarizations.The n-parallel modes are shown in Figure 19(a-b) by averaging along y to capture only the k y = 0 mode.The n-oblique modes are shown in Figure 19(c-d) by Fourier transforming E, setting its k y = 0 Fourier coefficient to zero, undoing the transform, and then averaging either the bottom one-half or one-fourth of the simulation domain along y.The choice of one-half and one-fourth allows us to capture wave power at, respectively, k y = 0.082d −1 e for M s = 4 (Fig. 18(a)) and k y = 0.164d −1 e

Figure 19 .
Figure 19.Shock precursor electric fields show right-handed polarization for both 2D Ms = 4 (left column) and 2D Ms = 7 (right column), when boosted into the upstream rest frame and projected into B0-perpendicular components.Right-handedness is shown by E ⊥1 component offset in x one quarter cycle ahead of E ⊥2 component.(a-b): Electric fluctuations E ⊥1 and E ⊥2 , y-averaged (i.e., only ky = 0 mode) and scaled to the upstream motional field strength u0B0/c in the simulation frame.(c-d): Like (a-b), but show modes oblique to shock normal by computing y average for one-half (c) or one-fourth (d) of the simulation domain, using Fourier-filtered fields with ky = 0 component forced to zero.
The HT potential has been used to interpret electron heating in satellite measurements, and it has been shown to agree with other proxies for electron heating (most notably, Liouville mapping), within the large systematic uncertainties of measuring low-frequency E fields in space(Cohen et al. (Comis ¸el et al. 2015;Marghitu et al. 2017by an explicit boost of simulation-frame fields with a constant global boost velocity; this may be contrasted with adaptive methods that account for local velocity and B variations(Comis ¸el et al. 2015;Marghitu et al. 2017).
was supported by NASA FINESST 80NSSC21K1383 and a NASA Cooperative Agreement awarded to the New York Space Grant Consortium; AT acknowledges travel support from Columbia's Arts & Sciences Graduate Council and Department of Astronomy.AT and LS acknowledge partial support by NASA ATP 80NSSC20K0565 and NSF AST-1716567.LS acknowledges support from the Cottrell Scholars Award, and DoE Early Career Award DE-SC0023015.