<figure><img src="image_1.png" /></figure>

# Origin of intense electron heating in relativistic blast waves

Arno Vanthieghem, Martin Lemoine, Laurent Gremillet

► To cite this version:

Arno Vanthieghem, Martin Lemoine, Laurent Gremillet. Origin of intense electron heating in relativistic blast waves. The Astrophysical journal letters, 2022, 930 (1), pp.L8. 10.3847/2041-8213/ac634f hal-03644417

HAL Id: hal-03644417

https://hal.science/hal-03644417

Submitted on 20 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Origin of Intense Electron Heating in Relativistic Blast Waves

Arno Vanthieghem<sup>1,2,3,4,5</sup> (✉), Martin Lemoine<sup>5</sup> (✉), and Laurent Gremillet<sup>6,7</sup> (✉)

Received 2021 December 19; revised 2022 March 16; accepted 2022 April 1; published 2022 May 3

## Abstract

The modeling of gamma-ray burst afterglow emission bears witness to strong electron heating in the precursor of Weibel-mediated, relativistic collisionless shock waves propagating in unmagnetized electron-ion plasmas. In this Letter, we propose a theoretical model, which describes electron heating via a Joule-like process caused by pitch-angle scattering in the decelerating, self-induced microturbulence and the coherent charge-separation field induced by the difference in inertia between electrons and ions. The emergence of this electric field across the precursor of electron-ion shocks is confirmed by large-scale particle-in-cell (PIC) simulations. Integrating the model using a Monte Carlo-Poisson method, we compare the main observables to the PIC simulations to conclude that the above mechanism can indeed account for the bulk of electron heating.

Unified Astronomy Thesaurus concepts: Shocks (2086); Plasma astrophysics (1261); Gamma-ray bursts (629); Relativistic jets (1390); High energy astrophysics (739); Particle astrophysics (96)

## 1. Introduction

The dissipation of the bulk energy of astrophysical outflows into nonthermal distributions of accelerated particles appears both generic and multifarious in the high-energy universe. Collisionless shock waves account for such dissipation in varied environments, from our own solar system to extreme relativistic events such as gamma-ray bursts (GRBs). At the external boundary of the GRB jet, a shock front indeed sweeps the external medium at a velocity close to that of light and eventually radiates part of the blast energy into multiwavelength spectra through synchrotron self-Compton radiation of shock-accelerated electrons (Kumar & Zhang 2015 and references therein). These emissions have by now been detected up to TeV energies (Abdalla et al. 2019; Acciari et al. 2019), and in one outstanding case, they have formed the electromagnetic counterpart of a gravitational-wave event (Abbott et al. 2017), thereby laying a foundation stone of multimessenger astrophysics.

Accordingly, this calls for an improved understanding of the physics of relativistic collisionless shock waves, the intricacy of which lies in the subtle intertwining of the self-generated electromagnetic microturbulence that mediates the dissipative dynamics of the dilute plasma with the nonthermal beam of accelerated particles that drives the microturbulence in which particles scatter and possibly radiate (Marcowith et al. 2016, for a review). One striking feature of weakly magnetized, relativistic collisionless shocks is their ability to convert a fraction as large as 10%–30% of the shock-dissipated energy into nonthermal high-energy electrons. This finding, confirmed in particle-in-cell (PIC) numerical experiments (Spitkovsky 2008; Martins et al. 2009; Haugbølle 2011; Sironi et al. 2013), nicely accounts for the large inferred radiative efficiency of GRB afterglows (Freedman &

Waxman 2001) and is therefore of prime interest for astrophysical phenomenology. For comparison, this fraction falls by one to two orders of magnitude in subrelativistic shock waves, such as those formed in supernova remnants (Völk et al. 2002) and, should electrons verify the shock-crossing conditions independently of ions, it would equal the electron-to-ion mass ratio $m_e/m_i$.

As viewed from the reference frame of the shock front, the incoming plasma energy is mostly carried by the ions, but the electrons can sap energy from this reservoir through their interaction with the electric fields generated in the shock precursor and transition layer. Yet the underlying mechanism remains widely debated: The scenarios proposed so far rely on the inductive electric field associated with the growth of the primary current filamentation instability (CFI) (Gedalin et al. 2012; Kumar et al. 2015) or of a secondary kink-type instability (Milosavljevic & Nakar 2006), as well as on the longitudinal (Gedalin et al. 2008; Plotnikov et al. 2013; Kumar et al. 2015) and/or transverse (Plotnikov et al. 2013) electro-static components that accompany the CFI and its oblique variants (Bret et al. 2008).

The objective of this Letter is to clarify the physics of electron heating in relativistic electron-ion collisionless shocks. To do so, we have performed large-scale PIC simulations of relativistic, unmagnetized electron-ion shock waves (Section 2). In parallel, we have developed a theoretical model, which generalizes a previous model for pair shocks (Lemoine et al. 2019a, 2019b, 2019c; Pelletier et al. 2019), to describe the deceleration and heating dynamics of the background plasma electrons and ions (Section 3). The numerical integration of this model, through a Monte Carlo-Poisson (MCP) technique, gives results that compare satisfactorily well with the PIC simulations. In this picture, electron heating can be depicted as a collisionless Joule process in a large-scale, longitudinal electric field due to charge separation in the background plasma, itself entailed by the differential scattering strength of electrons and ions in the Weibel microturbulence. Our

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1"><sup>1</sup> International Research Collaboration Center, National Institutes of Natural Sciences, Tokyo 105-0001, Japan; vanthieghem@astro.princeton.edu</footnote>

<footnote n="2"><sup>2</sup> Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="3"><sup>3</sup> High Energy Density Science Division, SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA</footnote>

<footnote n="4"><sup>4</sup> Sorbonne Université, Institut Lagrange de Paris (ILP), 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="5"><sup>5</sup> Institut d'Astrophysique de Paris, CNRS—Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="6"><sup>6</sup> CEA, DAM, DIF, F-91297 Arpajon, France</footnote>

<footnote n="7"><sup>7</sup> Université Paris-Saclay, CEA, LMCE, F-91680 Bruyères-le-Châtel, France</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Close-up of the precursor of an unmagnetized electron–ion collisionless shock in the PIC simulation frame. The upstream plasma is injected from the right-hand side with a Lorentz factor $\gamma_{\infty} = 100$, proper temperatures $T_e = T_i = 0.01 m_e$, and mass ratio $r = m_i/m_e = 100$. The top (middle and bottom) panel shows the magnetic field (ion density and electron energy fraction) profile at time $t = 870 \, \omega_{\text{pi}}^{-1}$. The shock front, located at $x = 0$, separates the unshocked plasma ($x > 0$) from the shocked plasma ($x < 0$).</figcaption></figure>

results are summarized and discussed in Section 4. Throughout we use units in which $k_B = c = 1$.

## 2. Kinetic Simulations

Our kinetic simulations are performed in 2D3V (2D in physical space, 3D in momentum space) geometry using the finite-difference time-domain, relativistic PIC CALDER code (Lefebvre et al. 2003), which has been extensively used in the relativistic regime (e.g., Vanthieghem et al. 2018), where it has been shown to properly expunge the salient relativistic beam-grid instabilities by means of the Cole–Karkkainen electromagnetic solver (Kärkkäinen et al. 2006) coupled with the Godfrey–Vay filtering method (Godfrey & Vay 2014). The plasma is injected from the right-hand side of the box with a negative relativistic velocity $\beta_{\infty}$ (corresponding Lorentz factor $\gamma_{\infty}$). To further reduce computational time, the domain is progressively lengthened using a moving injector, keeping electromagnetic waves and reflected particles away from the right boundary. On the left-hand side, we impose perfect conductor conditions for the fields and specularly reflecting conditions for the particles. We present here the results of three large-scale simulations with various mass ratios $r = m_i/m_e = (25, 100, 100)$ and Lorentz factors $\gamma_{\infty} = (100, 100, 10)$; in each, the plasma is initialized with low proper temperatures $T_i = T_e = 10^{-2} m_e$. The simulation frame corresponds to the downstream rest frame (|d subscript), hence $\gamma_{\infty}$ represents the relative Lorentz factor between upstream and downstream. Correspondingly, the shock Lorentz factor, measured relative to upstream is $\sqrt{3}\gamma_{\infty}$ for a relativistic shock in 2D3V. The relativistic electron skin depth $c/\omega_{\text{pe}}$ of the unshocked plasma is resolved with a mesh size $\Delta x = 0.08c/\omega_{\text{pe}}$ and $\Delta t = 0.99\Delta x$. Our simulations use 10 particles per cell and per species and cover a duration $L_t \simeq 10^3 \, \omega_{\text{pi}}^{-1}$ (with $\omega_{\text{pi}}$ the upstream ion plasma frequency), corresponding to $\simeq 120,000$ iterations, at which time the numerical domain comprises $120,000 \times 12,000$ cells. To enable proper comparison with the model, we distinguish the background plasma particles from the nonthermal population as in

Lemoine et al. (2019a): The background plasma is defined as the set of particles with negative velocity ($\beta_x \le 0$), which never experienced any turnaround in the microturbulence.

All simulations exhibit the characteristic filamentary pattern of the CFI, which is seen to grow in amplitude from the far upstream toward the shock front, as a result of the interpenetration of the beam of accelerated particles with the background plasma (Figure 1, top and middle panels). Across the precursor, the background plasma electrons are continuously energized to the point of reaching quasi-equipartition when crossing the shock front—i.e., the energy partitions $e_{\alpha} = e_{\alpha}/(e_e + e_i) \sim \mathcal{O}(1)$, where $e_{\alpha}$ represents the energy density of species $\alpha$ (Figure 1, bottom panel).

The transversely averaged drift Lorentz factor and temperature profiles for the background plasma species (Figure 2) show that, in each case, the electrons are abruptly slowed down upon entering the shock precursor, and progressively heated up to $T_e \simeq 0.5 T_i$ downstream, where the ion temperature matches the (2D relativistic) hydrodynamic jump conditions: $T_i/m_i = \gamma_{\infty}/2$. For $\gamma_{\infty} = 100$ and $r = (25, 100)$, the ion Lorentz factor shows a regular decrease until reaching the (subshock) transition layer ($x \simeq 50c/\omega_{\text{pi}}$), where the $\gamma_i$ and $\gamma_e$ curves converge with each other and rapidly drop to $\sim 1$. For $\gamma_{\infty} = 10$ and $r = 100$, the ions undergo only weak deceleration over most of the precursor and are hence mainly stopped in the transition layer.

A key observation, specific to the unmagnetized regime, is the existence of a particular reference frame—the Weibel frame (R<sub>w</sub>)—in which the microturbulence is essentially magnetostatic because it is generated by the CFI (Pelletier et al. 2019). In the simulation frame, this Weibel frame strongly decelerates from a large four-velocity $|u_w| \sim \gamma_{\infty}$ at the tip of the precursor to low subrelativistic $|u_w| \sim 0.5$ values at the shock transition. The dynamics of that Weibel frame can be extracted from the PIC simulations, in particular its
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Drift Lorentz factors (top) and temperatures (bottom) for ions (blue) and electrons (red) of the background plasma, extracted from our three reference PIC simulations with parameters: (a) $\gamma_{\infty} = 100$ and $r = 25$, (b) $\gamma_{\infty} = 100$ and $r = 100$, and (c) $\gamma_{\infty} = 10$ and $r = 100$. In the top panels, the solid gray lines plot the Lorentz factor of the Weibel frame as extracted from each simulation (see text).</figcaption></figure>

three-velocity $\beta_w(x) = -\langle\delta E_y^2\rangle_y^{1/2}/\langle\delta B_z^2\rangle_y^{1/2}$ (sign fixed by our conventions).

The x-profile of the corresponding Lorentz factor $\gamma_w$ is plotted (in gray) in the top panels of Figure 2. In each simulation case, the Weibel and electron Lorentz factors exhibit similar dynamics: This corroborates that the background plasma electrons are strongly coupled to the turbulence throughout the precursor, unlike the ions whose larger inertia lets them stream well ahead. This effect will be shown to be mainly responsible for the generation of a coherent, longitudinal electric field, which will eventually lead to intense electron heating.

Our PIC simulations confirm the generation of a coherent longitudinal electric field $\langle E_x \rangle_y$ across the precursor. Figure 3(a) plots the spatial profile of the associated electric potential, $\langle \phi \rangle_y = -\int \langle E_x \rangle_y dx$, in the three simulation cases.<fnref n="8" /> The potential is shown in solid lines, and the shaded areas indicate the $\pm 1\sigma$ variations of $-\int E_x dx$ (again, weighted by the electron density, but without the average over y) in the transverse direction. Importantly, this potential exhibits a net mean value, whose magnitude reaches $\langle \phi \rangle_y \approx 0.2-0.5 \gamma_{\infty} m_i / e$. It is sufficient, in sign and in order of magnitude, to explain electron heating up to near equipartition. In the following, we develop a theoretical model to interpret the main features of those PIC simulations, i.e., the deceleration and heating of electrons and ions, the origin of this electric field, and its role in electron heating.

We interpret this coherent electric field as resulting from two causes (see also below): (i) to a large extent, the charge separation between electrons and ions in the background plasma owing to their different degree of coupling to the turbulence, and (ii) to a lesser extent, at larger distances from the shock, the (apparent) net charge density, $\rho_b$, carried by the beam of suprathermal particles. The latter charge imbalance is weak in the simulation with $r = 25$ (red line in Figure 3(b)) but significant for $r = 100$ (black and blue lines). In the

simulations, it originates primarily from the initial stage of shock formation but, in actual situations, it can be sustained by differential injection of electrons and ions at the shock front. We thus take it into account and weigh its relevance against charge separation between the background electrons and ions by comparing the three simulations.

One can expect other contributions to $\langle E_x \rangle_y$, in particular an inductive electric field associated with the CFI growth or an electrostatic field linked to the broadband nature of the CFI. Qualitatively, those are expected to give contributions of opposite signs in current filaments of opposite polarity, thus translating into a small net value when directly averaged over the transverse direction. The field strength measured without weighting by the electron density supports the idea that it instead originates from charge separation, as discussed above. The successful comparison of our model to the PIC simulations will further corroborate this idea.

The above two sources of charge imbalance act differently on the background electrons and on the electric potential: The positively charged suprathermal beam tends to repel the electrons as it moves away from the shock, while the background plasma ions moving in the other direction tend to drag them toward the shock. Their zones of influence can be read off the bottom panel of Figure 3. The charge density profile impacts the bulk dynamics of the background electrons because of approximate charge and current neutrality at every point (as observed in the simulations). To see this, consider three populations: the electrons and ions of the thermal background and a beam of suprathermal particles carrying charge density $\rho_b$. To a good approximation, the particle current density of the background ions is conserved along the precursor (verified in the simulations), $n_i u_i = n_{\infty} u_{\infty}$. Then, imposing overall charge and current (quasi-)neutrality leads to

$$ \gamma_e^2 = \frac{\gamma_i^2}{\left[ 1 + 2\gamma_i^2 (1 + \beta_b) \frac{\rho_b}{\gamma_{\infty} n_{\infty} e} \right]}, \quad (1) $$

hence $\gamma_e^{-2} \propto \rho_b$ in the far precursor where $\rho_b/\gamma_{\infty} n_{\infty} e \gtrsim \gamma_{\infty}^{-2}$. Because $\gamma_e \sim \gamma_w$ (see Figure 2), the above relation implies $\rho_b \propto \gamma_w^{-2}$. Figure 3(b) demonstrates that this scaling law is nicely verified in our simulations.

<footnote n="8">⁸ To extract the electric potential as seen by the electrons, we weight the electric contribution at each point by the electron density of the background plasma. This compensates for the slight transverse heating inhomogeneities visible in Figure 1. The potentials calculated without this weighting techniques take lower values at the shock, but the coherent longitudinal electric field remains.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 3. (a) Estimate of the coherent mean electric potential $-\int \langle E_x \rangle_y dx$ weighted by the electron density (solid curves) and in the range of $\pm 1\sigma$ variation (shaded area). The different colors respectively correspond to the three PIC simulations with $(r, \gamma_\infty) = (100, 10)$ (blue), $(25, 100)$ (black), and $(100, 100)$ (red). The dotted-dashed lines correspond to the electrostatic potential reconstructed from the Monte Carlo-Poisson integration of our theoretical model. (b) Charge density of the beam component as measured in the PIC simulations with the same color code as in panel (a). The suprathreshold beam is defined as the ensemble of particles moving with longitudinal velocity $\beta_x > 0$. The dotted lines show the relationship between the (apparent) beam charge density $\rho_b$ and $\gamma_w$, as extracted from Equation (1), namely, $\rho_b \propto \gamma_\infty n_\infty e/\gamma_w^2$, which matches well the observed profile.</figcaption></figure>

The drift speed of $\mathcal{R}_w$ in which the background electrons relax through scattering depends nontrivially on the physical characteristics of the beam and the background plasma (Pelletier et al. 2019); in the presence of a net charge, this relation becomes even more complex because of charge compensation, as discussed above. The drift velocity of $\mathcal{R}_w$ is predicted to be smaller (in magnitude) than that of the background plasma; this notably explains the sharp deceleration of background electrons at the tip of the precursor,<fnref n="9" /> where they first penetrate this microturbulence.

### 3. Physics of Electron Heating

We model the physics of electron heating in the microturbulence of the shock precursor as follows. In the decelerating Weibel frame $\mathcal{R}_w$, particles are subject to an effective gravity, to the coherent electric field $\langle E_x \rangle_y$, and to angular scattering off the microturbulence. This combination gives rise to efficient heating through a collisionless Joule process, in which the gravity and the electric field serve as driving forces along $x$, while scattering redistributes the energy gained or lost in the transverse directions. The heating induced by the effective gravity can be understood, in the simulation frame, as due to the perpendicular motional electric fields carried by the filamentary structures (Lemoine et al. 2019c). We show below

that this contribution accounts for the bulk of ion heating, while the longitudinal coherent electric field will be responsible for most of the electron heating.<fnref n="10" />

Of course, the electric field does not only heat the electrons, it also slows down the ions and accelerates the electrons toward the shock. Therefore, not all of the magnitude of $e\phi$ at the shock is converted into electron heating. To understand the detailed contribution of the electric field to electron heating, we extend the Fokker-Planck description of Lemoine et al. (2019c) by incorporating this longitudinal component. To first order in the relative velocity $|\beta_{e|w}|$ between the background electrons and $\mathcal{R}_w$, we find that this Joule heating can be characterized by the momentum diffusion coefficient (for ultrarelativistic electrons)

$$D_{pp} = \frac{\beta_{w|s}^2 P_{|w}^2}{3\nu_{e|w}} \left[ \frac{du_{w|s}}{dx|_s} \left( 1 - \frac{1}{3\beta_{w|s}^2} \right) - \frac{qE_x}{\beta_{w|s}P_{|w}} \right]^2 \quad (2)$$

with the following notations: quantities indexed with $|s$ ($|w$) are understood to be defined in the rest frame of the shock ($\mathcal{R}_w$); $P_{|w}$ represents the particle momentum in $\mathcal{R}_w$. The first term in the brackets is the inertial term associated with the deceleration of $\mathcal{R}_w$ (effective gravity), while the second represents the contribution of the coherent longitudinal electric field to the driving force. Such a Fokker-Planck analysis is well suited to the electrons here because their drift velocity in the Weibel frame $|\beta_{e|w}| \ll 1$; see Section 2.

Interestingly, Equation (2) indicates that the sign of the electric field does not matter much as it relates to heating through the Joule effect. This agrees, at least qualitatively, with the observation that electrons are systematically heated over the precursor, even in regions in which $\langle E_x \rangle_y$ is negative because of the influence of the nonneutral particle beam (Figure 2).

To validate our theoretical model against PIC simulations, we numerically integrate it using an MCP method, which solves the transport equation, including the requisite physical ingredients of our theoretical model while discarding other kinetic effects. We then extract the main observables, namely the profiles of the four-velocity and temperature of each species along the shock normal, and compare them to those seen in the PIC simulation.

Specifically, the model describes the kinematics of the background plasma particles (in $\mathcal{R}_w$) through the (discretized) stochastic equations

$$\Delta\mu_{\alpha|w} = \sqrt{2\nu_{\alpha|w}\Delta t_{|w}} \mathcal{X}, \quad (3)$$

$$\Delta p_{\alpha|w}^x = p_{\alpha|w} \Delta\mu_{\alpha|w} + q_{\alpha} E_x \Delta t_w - (\beta_{w|s} \epsilon_{\alpha|w} + p_{\alpha|w}^x) \frac{du_{w|s}}{dx|_s} \Delta t_w, \quad (4)$$

where $\mu_{\alpha|w}$ denotes the pitch-angle cosine relative to the shock normal, $p_{\alpha|w}^x \equiv \mu_{\alpha|w} P_{\alpha|w}$ is the longitudinal momentum, and $\mathcal{X} \sim \mathcal{N}(0, 1)$ represents white noise, whose role is to simulate pitch-angle scattering in the microturbulence. This numerical scheme is thus similar to that of a 1D electrostatic PIC code. At

<footnote n="9">⁹ Note that Figure 2 displays the drift Lorentz factors and temperatures. When plotting the average Lorentz factor, one recovers the behavior shown in Spitkovsky (2008) and Haugbølle (2011).</footnote>

<footnote n="10">¹⁰ The transverse electromagnetic fields, as extracted from our simulations to define $\mathcal{R}_w$, encompass perpendicular electric fields of any nature. While their contribution is subdominant relative to Weibel modes, the perpendicular electric fields associated with oblique instabilities and their effect are thus, to some extent, captured in our analysis.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** In shaded bands, distributions of the electron (red) and ion (blue) scattering frequencies across the shock precursor as extracted from the MCP solution to the transport equation using the continuous formula (5) with ($r, \gamma_{\infty}$) = (100, 100). For comparison, we overlay the scattering frequencies in the trapped (dotted lines) and untrapped (solid lines) limits as extracted from the corresponding PIC simulation using Equations (34) and (35) of Lemoine et al. (2019c).</figcaption></figure>

every time step, the momentum of each particle is advanced in $\mathcal{R}_w$ through Equations (3) and (4), where the discretized electric field is interpolated at the particle position. Moving to the shock frame, the particles' positions are evolved using their updated momenta, and the electrostatic field is computed by a standard Poisson solver from the grid-projected particle charges. A complete description of the numerical method will be provided elsewhere (A. Vanthieghem 2022, in preparation).

The above MCP model depends on the following ($x$-dependent) parameters: the scattering frequencies $\nu_{\alpha|w}(x, p)$, the law of deceleration of the Weibel frame $u_w(x)$, and the external charge $\rho_b(x)$ imposed by the beam. Once these functional forms are fixed, the numerical integration gives the law of deceleration $u_{e,i}(x)$ and heating $T_{e,i}(x)$ of the electrons and ions of the background plasma, which can be compared with the outputs of the PIC simulations.

The model equations are expressed in the shock frame, assuming stationarity, while our PIC simulations are time dependent and run in the downstream rest frame. Given their long timescales, however, the precursor has reached, at least close to the shock, a near self-similar profile that depends only on $x - \beta_{sh}t \propto x_i$. A fully consistent Lorentz transform from the downstream frame to the shock-front frame would require the PIC plasma profiles to be recorded at multiple time steps. Here, we simply approximate the Lorentz boost by multiplying the space coordinate by $\gamma_{sh} = \sqrt{4/3}$.

The scattering frequency $\nu_{|w}$ in the Weibel frame is predicted to take on different forms depending on whether the particle is trapped or not in the magnetic filaments; known expressions are summarized in Equations (34) and (35) of Lemoine et al. (2019c). For a unified description, our MCP model uses the following modified continuous formulation that retains the momentum dependencies in the trapped and untrapped limits:

$$ \nu_{|w}(p_{|w}) = \begin{cases} \nu_{0|w} \beta_{|w} (p_{|w}/m_i)^{-2}, & \text{if } p_{|w} \ge p_{0|w}, \\ \nu_{|w} (p_{0|w}) (\gamma_{|w}/\gamma_{0|w})^{-1}, & \text{else} \end{cases}, \quad (5) $$

with $p_{|w}$ the particle momentum in the Weibel frame ($\gamma_{|w}$ its Lorentz factor); $\nu_{0|w} \propto \epsilon_B/k_{\perp}$ a reference scattering frequency expressed in terms of $\epsilon_B = \langle B^2 \rangle_y / 4\pi m_i n_i \gamma_{\infty}^2$, the magnetic

energy density fraction; and $k_{\perp}$ the dominant wavenumber in the transverse direction, i.e., the inverse of the filament size. Particles transit from trapped to untrapped populations when their Larmor radius becomes comparable with the typical size of a filament, i.e., $p_{0|w}/m_i \approx \epsilon_B^{1/2} \omega_{pi}/k_{\perp}$. As further detailed in the Appendix, to solve this transport model, we fix the constant of proportionality between $\nu_{0|w}$ and $\epsilon_B/k_{\perp}$; we also extract the profile of $u_w$ and $\rho_b$, which we use directly in the MCP solver.

We now discuss the results of this procedure, first comparing the scattering frequencies, as reconstructed along the history of the background plasma in the MCP model, with those extracted from the PIC simulation. For the latter, we apply known formulae (Equations (34) and (35) of Lemoine et al. 2019c). The comparison is shown in Figure 4, which confirms that the reconstruction provides a satisfactory match to the scattering frequencies estimated from the PIC simulations. As one moves closer to the shock front, $\nu_{e|w}$ decreases while $\nu_{i|w}$ increases, by an order of magnitude or more, as a result of the evolution of the CFI spectrum from electron kinetic scales at the tip of the precursor to ion kinetic scales near the shock front, in conjunction with strong electron heating; $\nu_{e|w}$ and $\nu_{i|w}$ indeed meet near the shock front where the electrons reach near equipartition.

Our model indeed reproduces the large amount of electron heating observed at the shock for our three different simulations; this is illustrated in Figure 5. More precisely, for a scaling $\nu_{0|w} \sim \epsilon_B/k_{\perp}$ (meaning, with a prefactor not far from unity), it captures the spatial profiles of the temperature and four-velocity for both ions and electrons of the background plasma. This satisfactory reconstruction holds across the whole precursor, including the abrupt electron slowdown at its tip, down to the fast dynamics inside the subshock layer. This finding represents the main result of our work.

This numerical integration also reproduces fairly well the profile of the electric potential across the precursor; see the dashed lines in Figure 3(a). Near the shock, where most of the electron heating occurs, the reconstructed values lie within a factor $\sim 2$ of the observed ones, within the shaded bands of variation.

That the “best-fit” values tying $\nu_{0|w}$ to $\epsilon_B/k_{\perp}$ in our numerical MCP reconstructions lie within an order of magnitude of unity provides additional support to our model because $\nu_{0|w} \sim \epsilon_B/k_{\perp}$ is expected on theoretical grounds. We note some discrepancy between the prefactors used for the $\gamma_{\infty} = 100$ simulations (6 for both) and those used for $\gamma_{\infty} = 10$ (0.2). We do not expect this prefactor to depend on $\gamma_{\infty}$. However, we note that this “best-fit” value is subject to some uncertainty; it depends in particular on the extracted profiles of $k_{\perp}$, $\epsilon_B$, $u_w$, and $\rho_b$ from the PIC simulations (as discussed in the Appendix). We also note that the notion of “best fit” itself is somewhat vague: Acceptable, albeit less satisfactory reconstructions of $u_{i,e}$ and $T_{i,e}$ can be obtained for values of the prefactor differing by up to an order of magnitude, as illustrated in Figure 8 of the Appendix.

We find that ions and electrons are heated through different mechanisms. For ions, the dominant contribution comes from collisionless Joule heating driven by the effective gravity in the decelerating turbulence frame, or equivalently, by stochastic interactions with the perpendicular motional electric fields if seen in the simulation frame: The light blue line in Figure 5 (middle panel) indeed shows that $T_i$ and $u_i$ can be reproduced
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of the drift Lorentz factor (top) and temperature (bottom) between the PIC simulation with the respective ($r, \gamma_{\infty}$) = (25, 100), (100, 100), (100, 10) (gray) and the MCP solution to the transport equation for the ions (blue) and electrons (red). The transport equation is integrated using the scattering frequency given by Equation (5) with $\nu_{0|w} \simeq 6, 6, 0.2 \epsilon_B/k_{\perp}$. The light red and blue curves in the middle panel correspond to the solution to the transport equation ignoring the contribution of the electrostatic field; while they globally reproduce the observed amount of ion heating, they indicate insufficient electronic heating. We thus find that the shock dynamics is correctly reproduced for $\nu_{0|w} \sim \epsilon_B/k_{\perp}$ when accounting for the electrostatic field contribution.</figcaption></figure>

fairly well even when ignoring longitudinal electric fields. This is not altogether surprising as the ions are the dominant carriers of inertia and, in the case of a pair shock, for which this coherent electric field vanishes (Pelletier et al. 2019), heating proceeds in the same manner. On the other hand, the electric field appears to provide the dominant source of heating for electrons. This is here illustrated by the light red line, which shows the corresponding evolution for electrons in the absence of $\langle E_x \rangle_y$; in that case, $T_e$ at the shock lies $\gtrsim 2$ orders of magnitude below that observed in PIC simulations.

## 4. Discussion and Conclusions

Our results thus indicate that the bulk of electron heating, up to near equipartition at relativistic, unmagnetized shock waves, is associated with the self-consistent generation of a coherent (along the transverse direction) longitudinal electric field across the shock precursor. This field mainly originates from the charge separation imposed by the differential dynamics of the background plasma ions and electrons in the Weibel microturbulence. In our theoretical description, electron heating results from a collisionless Joule process imparted by the coherent electric field in conjunction with pitch-angle scattering off the microturbulence, as expressed in the $R_w$ frame in which this microturbulence can be seen to be essentially magnetostatic. The effective gravity force, felt by background ions and electrons in the decelerating $R_w$ frame, provides an additional source of heating; while subdominant for electrons, it accounts for most of the heating and deceleration of the ions and therefore for the shock-crossing conditions.

Our conclusions rest on the successful comparison between a numerical integration of the above model and large-scale PIC simulations of relativistic electron-ion shocks, which we have conducted for values of the shock Lorentz factor up to 100 and for mass ratios up to 100. Those simulations show the progressive heating of the ions and electrons up to quasi-equipartition, and they reveal the existence of a longitudinal coherent electric field across the shock precursor. To test our model, we have extracted the velocity profile of the microturbulence frame $R_w$ and the beam charge profile. Using those estimates together with a general law for the scattering of particles in the microturbulence, we have performed a numerical

integration of our model, which describes the stochastic heating and deceleration of a plasma in a decelerating Weibel-type turbulence, in the presence of a longitudinal electric field, by means of an MCP solver. This integration compares satisfactorily to the PIC simulations: In particular, it reproduces fairly well the evolution of the electric potential as a function of distance to the shock, as well as the profiles of deceleration and heating of ions and electrons across the precursor. This suggests that this model captures the main features of electron and ion heating at weakly magnetized, relativistic shock waves.

In conclusion, let us point out some similarities and differences with previous theoretical ideas on electron heating. Our model shares some features with that of Gedalin et al. (2008), which describes electron heating in a cross-shock potential; in our scenario, however, the potential is rather found to extend over the precursor, and heating is mostly stochastic, not directly associated with a DC jump as in that reference. Moreover, the electric field originates from the difference in inertia, not from the growth of the CFI, as proposed in Gedalin et al. (2012) and Kumar et al. (2015), nor from the broadband nature of that instability and its oblique variants (Gedalin et al. 2008; Plotnikov et al. 2013; Kumar et al. 2015). We also note that the geometric configuration of the shock precursor is key to the emergence of this longitudinal electric field; a symmetric counterstreaming configuration, with two interpenetrating plasmas sharing similar physical characteristics, as simulated by Kumar et al. (2015), could not observe its presence. Finally, the origin of that electric field also differs from that envisaged in Milosavljevic & Nakar (2006) and Naseri et al. (2018), which attribute heating to the inductive electric field generated by the disruption of Weibel filaments.

A.V. thanks Frederico Fiuza and the anonymous referee for insightful suggestions and acknowledges support by the US DOE Early Career Research Program under FWP 100331. The authors acknowledge financial support from the ILP Labex (reference ANR-10-LABX-63) as part of the Idex SUPER (reference ANR-11-IDEX-0004-02); the DIWINE Emergence SU 2019 program; the ANR-14-CE33-0019 MACH project; the ANR-20-CE30-0030 UnRIP project; the NSF grants PHY-1804048 and AST-1814708. Simulations were run on the HPC
resources of TGCC/CCRT under the allocations 2018-A0030407666 and 2019-A0030407666 made by GENCI and on CORI at the National Energy Research Scientific Computing Center (NERSC) through the ALCC award.

## Appendix
### Estimating the Scattering Frequency

Here, we specify how the profiles of the physical quantities ($k_{\perp}$, $\epsilon_B$, $u_w$, $\rho_B$) used in our MCP model are extracted from the

PIC simulations. All other quantities then follow from the numerical integration. The dominant wavevector $k_{\perp}$ of the microturbulence is obtained from the Fourier space decomposition along the y-axis across the shock precursor, as illustrated in Figure 6. One observes a progressive transition from electron skin depth scales at the tip of the precursor to an ion skin depth scale close to the shock transition. The magnetization fraction $\epsilon_B$ is obtained from the y-averaged magnetic field energy density. The magnetic field profiles for the different simulations are shown on the left of Figure 7

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Longitudinal profile of the transverse Fourier modes of the out-of-plane ($B_z$) magnetic field across the shock precursor for ($r, \gamma_{\infty}$) = (25, 100), (100, 100), (100, 10) (from left to right). The two horizontal dotted lines correspond to the electron (top) and ion (bottom) plasma scales. The black line follows the dominant transverse wavenumber across the precursor, which is used to estimate the scattering frequency $\nu_{0|w}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Left: out-of-plane magnetic field component ($B_z$) in the simulation plane, for ($r, \gamma_{\infty}$) = (25, 100), (100, 100), (100, 10) (from top to bottom). Right: one-dimensional profile along the shock normal of the magnetic energy fraction $\epsilon_B$, for the three simulations. This profile is used in the estimate of $\nu_{0|w}$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison of the drift Lorentz factor (top) and temperature (bottom) between the PIC simulation with ($r, \gamma_{\infty}$) = (100, 100) (gray) and the MCP solution to the transport equation (including the electric field contribution) for the ions (blue) and the electrons (red). The transport equation is integrated using the scattering frequency given by Equation (5) with, from left to right, $\nu_{0|w}$ = 0.6, 6, 60 $\epsilon_B/k_{\perp}$.</figcaption></figure>
while the associated magnetization level is shown on the
right. Those quantities then serve to estimate the leading
dependency of the scattering frequency, namely $\nu_{0|w} \propto \epsilon_B/k_\perp$;
see Equation (5). The remaining prefactor, which we expect to
be of the order of unity, is adjusted to reproduce the
simulation results.

To illustrate the dependence of our results on this prefactor,
we focus on the simulation with the largest Lorentz factor and
largest mass ratio, ($r, \gamma_\infty$) = (100, 100). We plot in Figure 8
the reconstructed history of the plasma four-velocity and
temperature, as in Figure 5, for values of $\nu_{0|w}$ alternatively
larger and smaller by an order of magnitude. As this prefactor
departs from the value (6) used in the main text, the profiles
start to differ from those measured in the PIC simulation,
although the global amount of electron heating remains of the
correct order of magnitude. This indicates that our results are
not strongly sensitive to the choice of $\nu_{0|w}$.

ORCID iDs

Arno Vanthieghem https://orcid.org/0000-0002-3643-9205

Martin Lemoine https://orcid.org/0000-0002-2395-7812

Laurent Gremillet https://orcid.org/0000-0003-0116-5248

References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJL, 848, L12

Abdalla, H., Adam, R., Aharonian, F., et al. 2019, Natur, 575, 464

Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2019, Natur, 575, 455

Bret, A., Gremillet, L., Bénisti, D., & Lefebvre, E. 2008, PhRvL, 100, 205008

Freedman, D. L., & Waxman, E. 2001, ApJ, 547, 922

Gedalin, M., Balikhin, M. A., & Eichler, D. 2008, PhRvE, 77, 026403

Gedalin, M., Smolik, E., Spitkovsky, A., & Balikhin, M. 2012, EPL, 97, 35002

Godfrey, B. B., & Vay, J.-L. 2014, JCoPh, 267, 1

Haugbølle, T. 2011, ApJL, 739, L42

Kärkkäinen, M., Gjonaj, E., Lau, T., & Weiland, T. 2006, in Proc. Int. Computational Accelerator Physics Conf. (Online) 35, https://accelconf.web.cern.ch/icap06/PAPERS/MOM2IS03.PDF

Kumar, P., & Zhang, B. 2015, PhR, 561, 1

Kumar, R., Eichler, D., & Gedalin, M. 2015, ApJ, 806, 165

Lefebvre, E., Cochet, N., Fritzler, S., et al. 2003, NucFu, 43, 629

Lemoine, M., Gremillet, L., Pelletier, G., & Vanthieghem, A. 2019a, PhRvL, 123, 035101

Lemoine, M., Pelletier, G., Vanthieghem, A., & Gremillet, L. 2019b, PhRvE, 100, 033210

Lemoine, M., Vanthieghem, A., Pelletier, G., & Gremillet, L. 2019c, PhRvE, 100, 033209

Marcowith, A., Bret, A., Bykov, A., et al. 2016, RPPh, 79, 046901

Martins, S. F., Fonseca, R. A., Silva, L. O., & Mori, W. B. 2009, ApJL, 695, L189

Milosavljevic, M., & Nakar, E. 2006, ApJ, 641, 978

Naseri, N., Bochkarev, S. G., Ruan, P., et al. 2018, PhPl, 25, 012118

Pelletier, G., Gremillet, L., Vanthieghem, A., & Lemoine, M. 2019, PhRvE, 100, 013205

Plotnikov, I., Pelletier, G., & Lemoine, M. 2013, MNRAS, 430, 1280

Sironi, L., Spitkovsky, A., & Arons, J. 2013, ApJ, 771, 54

Spitkovsky, A. 2008, ApJL, 673, L39

Vanthieghem, A., Lemoine, M., & Gremillet, L. 2018, PhPl, 25, 072115

Völk, H. J., Berezhko, E. G., Ksenofontov, L. T., & Rowell, G. P. 2002, A&A, 396, 649