Electron Heating in Low Mach Number Perpendicular shocks. I. Heating Mechanism

Electron Heating in Low Mach Number Perpendicular shocks. I. Heating Mechanism

Xinyi Guo, Lorenzo Sironi, and Ramesh Narayan Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA
Department of Astronomy, Columbia University, 550 W 120th St, New York, NY 10027, USA
xinyi.guo@cfa.harvard.edu lsironi@astro.columbia.edu rnarayan@cfa.harvard.edu

Recent X-ray observations of merger shocks in galaxy clusters have shown that the post-shock plasma is two-temperature, with the protons hotter than the electrons. By means of two-dimensional particle-in-cell simulations, we study the physics of electron irreversible heating in perpendicular low Mach number shocks, for a representative case with sonic Mach number of 3 and plasma beta of 16. We find that two basic ingredients are needed for electron entropy production: (i) an electron temperature anisotropy, induced by field amplification coupled to adiabatic invariance; and (ii) a mechanism to break the electron adiabatic invariance itself. In shocks, field amplification occurs at two major sites: at the shock ramp, where density compression leads to an increase of the frozen-in field; and farther downstream, where the shock-driven proton temperature anisotropy generates strong proton cyclotron and mirror modes. The electron temperature anisotropy induced by field amplification exceeds the threshold of the electron whistler instability. The growth of whistler waves breaks the electron adiabatic invariance, and allows for efficient entropy production. We find that the electron heating efficiency displays only a weak dependence on mass ratio (less than drop, as we increase the mass ratio from up to ). We develop an analytical model of electron irreversible heating and show that it is in excellent agreement with our simulation results.

Subject headings:
galaxies: clusters: general — instabilities — radiation mechanisms: thermal — shock waves

1. Introduction

Galaxy clusters grow via mergers of subclusters. A large fraction of the kinetic energy of the infalling subclusters is dissipated at low Mach number shocks (, where is the ratio of shock speed and pre-shock sound speed), that heat the intracluster medium (ICM) and sometimes accelerate particles to relativistic energies (Sarazin, 2002; Ryu et al., 2003; Brüggen et al., 2012). Merger shocks in clusters are collisionless. Due to the high ICM temperatures (K) and low densities (), the collisional mean free path ( cm) is as large as the size of the cluster.

Galaxy cluster shocks are routinely observed in the radio and X-ray bands. X-ray measurements can quantify the density and temperature jumps between the unshocked (upstream) and the shocked (downstream) plasma (e.g., Markevitch et al., 2002; Finoguenov et al., 2010; Russell et al., 2010; Ogrean et al., 2013; Eckert et al., 2016; Akamatsu et al., 2017). The existence of shock-accelerated electrons is revealed by radio observations of synchrotron radiation (e.g., van Weeren et al., 2010; Lindner et al., 2014; Trasatti et al., 2015; Kale et al., 2017). Recently, the pressure jump associated with a merger shock at relatively high redshift has been measured through radio observations of the thermal Sunyaev-Zel’dovich (SZ) effect (Basu et al., 2016).

Since all of our observational diagnostics are based on radiation emitted by electrons, the proton properties (in particular, their temperature) are basically unconstrained. One usually makes the simplifying assumption that the electron temperature equals the mean gas temperature (and so, the proton temperature). This assumption is unlikely to hold in the vicinity of merger shocks. Ahead of the shock, the bulk kinetic energy of protons is a factor of larger than for electrons (here, and are the proton and electron masses, respectively). In the absence of a channel for efficient proton-to-electron energy transfer, a comparable ratio should persist between the post-shock temperatures of the two species.

While Coulomb collisions will eventually drive electrons and protons to equal temperatures, the collisional equilibration timescale (Spitzer, 1962) for typical ICM conditions is as long as yrs. In fact, X-ray observations by Russell et al. (2012) have shown that the electron temperature just behind a merger shock in Abell 2146 is lower than the mean gas temperature expected from the Rankine-Hugoniot jump conditions, and thus lower than the proton temperature. As a separate evidence, Akamatsu et al. (2017) has compiled a list of merger shocks, estimating their Mach number from both X-ray () and radio observations (), and noticed a slight bias of . Here, is derived by measuring the power-law slope of the synchrotron emission, which is related — via the theory of diffusive shock acceleration — to the density compression at the shock (and so, to the Mach number). On the other hand, is obtained from the electron free-free emission by measuring the jumps in density and temperature across the shock. It follows that, if electrons have a lower temperature than protons behind the shock, would have been underestimated.

In fact, it has long been thought that collisionless shocks can lead to a two-temperature structure at the outskirts of galaxy clusters (Fox & Loeb, 1997; Ettori & Fabian, 1998; Takizawa, 1999). Detailed cosmological hydrodynamic simulations have shown that this can significantly bias the X-ray and thermal SZ signatures (Wong & Sarazin, 2009; Rudd & Nagai, 2009). In the absence of a physical model for electron heating in low Mach number shocks, these studies usually employ an ad-hoc subgrid approach to prescribe the electron heating efficiency in shocks. Either electrons are heated adiabatically, or the non-adiabatic (or “irreversible”) heating efficiency is quantified by a phenomenological (and often, arbitrary) parameter. While observations from heliospheric low Mach number shocks have shown that electrons do not get heated much beyond adiabatic compression (Bame et al., 1979; Ghavamian et al., 2013), there has also been direct evidence of electron entropy production (i.e., non-adiabatic heating) at low Mach number shock fronts (Parks et al., 2012).

What is then the mechanism responsible for electron heating at collisionless shocks? This is a fundamental question of plasma physics, as the fluid-type Rankine-Hugoniot relations only predict the jump in the mean plasma temperature across the shock, without specifying how the shock-generated heat is distributed between the two species. To understand electron heating in collisionless shocks, fully-kinetic simulations with the particle-in-cell (PIC) method (Birdsall & Langdon, 1991; Hockney & Eastwood, 1981) are essential to self-consistently capture the non-linear structure of the shock and the role of electron and proton plasma instabilities in particle heating.

So far, PIC studies of electron heating in shocks have focused on the regime of high sonic Mach number () and low plasma beta () appropriate for supernova remnants. At very high Mach numbers, the Buneman instability can trap electrons in the shock transition region and heat them (Dieckmann et al., 2012). For lower Mach numbers, resonant wave-particle scattering induced by the modified two-stream instability (MTSI) can lead to significant electron heating at the shock front (Matsukiyo & Scholer, 2003; Matsukiyo, 2010).

The regime of low sonic Mach number and high beta most relevant for cluster merger shocks is still unexplored. In this paper, we study electron heating in low Mach number perpendicular shocks by means of two-dimensional (2D) PIC simulations. We focus on the results from a reference shock simulation with and . In a forthcoming paper (X. Guo et al., in preparation) we will explore the dependence of our conclusions on sonic Mach number and plasma beta. The choice of a perpendicular magnetic field geometry is meant to minimize the role of non-thermal electrons that are self-consistently accelerated in oblique configurations, as we have shown in Guo et al. (2014a, b). In the absence of shock-accelerated electrons returning upstream, the shock can settle to a steady state on a shorter time, thus allowing to focus on the steady-state electron heating physics. However, we emphasize that we have verified with a suite of PIC simulations of quasi-perpendicular shocks (not shown here) that the physics of electron heating presented in this paper also applies to quasi-perpendicular configurations, as long as the non-thermal electrons are energetically sub-dominant.

The rest of the paper is organized as follows. In Section 2 we lay out the theoretical framework for electron heating. Section 3 describes the setup of the reference shock simulation. Section 4 shows the shock structure of the reference simulation. We emphasize that efficient electron irreversible heating occurs at two main locations. With periodic box experiments meant to reproduce the shock conditions at the two major sites of entropy production, Section 5 and Section 6 investigate in detail the electron heating physics in these two locations, and validate the heating theory presented in Section 2. In Section 7, the heating model is then validated in the shock simulation. We conclude with a summary in Section 8.

2. The Physics of Electron Heating

As they pass through the shock, electrons will experience a density compression, which results in adiabatic heating. In addition, irreversible processes might operate, which will further increase the electron temperature. The purpose of this section is to present a general formalism for the physics of irreversible heating. Even though we will be primarily interested in electron heating, the model can be applied to any particle species. It relies on the presence of two basic ingredients: (i) a temperature anisotropy; and (ii) a mechanism to break the adiabatic invariance. We first describe the change in internal energy of an anisotropic fluid, and then consider the resulting change in entropy.111We point out that the model that we propose is reminiscent of the so-called “magnetic pumping” mechanism, where a periodically-varying external magnetic field is used in the laboratory to drive proton anisotropy and subsequent plasma heating (Spitzer & Witten, 1953; Berger et al., 1958; Borovsky, 1986).

2.1. The Change in Internal Energy

The work done on an isotropic gas with pressure and volume is simply . We shall generalize this expression to the case of an anisotropic gas having pressure perpendicular (parallel, respectively) to the magnetic field lines equal to (, respectively). Consider a magnetic flux tube with length , cross-sectional area , volume , and field strength . The magnetic flux through the area is . In response to a compression (or expansion) perpendicular to the magnetic field, the volume will change as


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


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


Defining the work done per particle as , we find that it can be separated into a “perpendicular” component and a “parallel” component as


It follows that will change the internal energy per particle associated with motions perpendicular to the field, while will affect the energy per particle associated with parallel motions.

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


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


where we denote the total energy per particle transferred to waves as , including magnetic, electric and bulk kinetic contributions (in practice, the magnetic term always dominates).

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

2.2. The Change in Entropy

For a non-relativistic bi-Maxwellian plasma with perpendicular temperature and parallel temperature , the entropy per particle (or specific entropy) is


where is the phase space distribution and is a normalization constant. By differentiating,


The temperature can be related to the internal energy per particle via the respective adiabatic index as


For a non-relativistic gas, (two degrees of freedom are available in the perpendicular direction), whereas (one degree of freedom). The equation above then becomes


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


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


As anticipated above, the two separate components and of the wave energy per particle do not explicitly enter the entropy equation.

In Equations (13) and (14), the first term on the right hand side typically dominates. This clearly demonstrates that two ingredients are required for entropy generation: (i) the presence of a temperature anisotropy; and (ii) a mechanism to break the adiabatic invariance. Note that the CGL double adiabatic theory of Chew et al. (1956) predicts that, for adiabatic perturbations, and , which follow from the conservation of the first and second adiabatic invariants. The form of Equations (13) and (14) is thus easy to understand. In most cases, it is the temperature anisotropy that provides the free energy for generating the waves responsible for breaking the adiabatic invariance.

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

3. Setup of the Shock Simulations

We perform numerical simulations using the three-dimensional (3D) electromagnetic PIC code TRISTAN-MP (Spitkovsky, 2005), which is a parallel version of the code TRISTAN (Buneman, 1993) that was optimized for studying collisionless shocks. In this section, we describe the setup of our shock simulations, which parallels closely what we used in Guo et al. (2014a, b). In Section 5 and in Section 6, we will study in more detail the physics of electron heating by employing periodic simulation domains, meant to represent two specific regions of the shock structure. The simulation setups adopted there are different, and are described in the respective sections.

For shock simulations, we use a 2D simulation box in the plane, with periodic boundary conditions in the direction. Even though the simulations are two-dimensional in space, all three components of particle velocities and electromagnetic fields are tracked. The shock is set up by reflecting an upstream electron-proton plasma moving along the direction off a conducting wall at the leftmost boundary of the computational box (). The interplay between the reflected stream and the incoming plasma causes a shock to form, which propagates along . In the simulation frame, the downstream plasma is at rest.

The upstream electron-proton plasma is initialized following the procedure described by Zenitani (2015), as a drifting Maxwell-Jüttner distribution with electron temperature equal to the proton temperature (i.e. ), and bulk velocity . This gives a simulation-frame Mach number


where is the sound speed in the upstream, is the Boltzmann constant, is the adiabatic index for an isotropic non-relativistic gas, and is the proton mass. Below, we will adopt the usual definition of Mach number, as the ratio between the upstream flow velocity and the upstream sound speed in the shock rest frame (rather than in the downstream frame of the simulations, as in Equation (15)), where the upstream moves into the shock with speed . We will then parameterize our results with the Mach number


The shock-frame Mach number is related to the downstream-frame Mach number via


where the density jump across the shock, in the limit of weakly magnetized flows, is equal to


In writing these relations we have assumed an isotropic gas, which is valid upstream of the shock by our initial conditions, and is also valid sufficiently downstream of the shock, as we will see in the discussion that follows.

The incoming plasma carries a uniform magnetic field , and its associated motional electric field . The magnetic field strength is parametrized by the plasma beta


where is the number density of the incoming protons and electrons. Alternatively, one could quantify the magnetic field strength via the Alfvénic Mach number .

The magnetic field is initialized in the simulation plane along the direction, i.e., perpendicular to the shock normal. We find that the shock physics is properly captured by 2D simulations only if the field is lying in the simulation plane. A posteriori, this will be motivated by the fact that the plasma instabilities excited in the downstream region have wavevectors preferentially parallel or quasi-parallel to the background magnetic field. We have explicitly verified with an additional simulation having magnetic field initialized along (so, pependicular to both the shock normal and the simulation plane), that the electron heating efficiency is completely suppressed, just as in 1D simulation results (Appendix A). Our choice of an in-plane magnetic field configuration will be justified again in the following sections.

For accuracy and stability, PIC codes have to resolve the plasma oscillation frequency of the electrons


and the electron plasma skin depth , where and are the electron charge and mass, respectively. On the other hand, the shock structure is controlled by the proton Larmor radius


where the shock-frame Alfvénic Mach number is . Similarly, the evolution of the shock occurs on a time scale given by the proton Larmor gyration period . The need to resolve the electron scales, and at the same time to capture the shock evolution for many , is an enormous computational challenge for the realistic mass ratio . Thus we adopt a reduced mass ratio for our reference run, which is sufficient to properly separate the electron and proton scales. This allows us to follow the system for long times, until the shock reaches a steady state. We have explicitly verified that the electron heating physics in our shock simulations is nearly the same for higher mass ratios (see Section 7, where we test up to ). In addition, in Section 5 and Section 6 we demonstrate via analytical arguments and PIC simulations in periodic domains that the electron heating efficiency is nearly independent of over the range from up to the realistic mass ratio.

As in Guo et al. (2014a, b), the upstream plasma is initialized at a “moving injector”, which recedes from the wall in the direction at the speed of light. When the injector approaches the right boundary of the computational domain, we expand the box in the direction. This way both memory and computing time are saved, while following at all times the evolution of the upstream regions that are causally connected with the shock. Further numerical optimization can be achieved by allowing the moving injector to periodically jump backward (i.e. in the direction), resetting the fields to its right (see Sironi & Spitkovsky (2009)). For a perpendicular shock (i.e., with magnetic field perpendicular to the shock direction of propagation), no particles are expected to escape ahead of the shock, so we choose to jump the injector in the direction so as to keep a distance of a few proton Larmor radii ahead of the shock. This suffices to properly capture the heating physics of electrons and protons. We have checked, though only for relatively early times, that simulations with and without the jumping injector give identical results.

In the main body of this paper, we present the results from a reference run with and , as motivated by galaxy cluster shocks. The upstream plasma is initialized with and . We remark that even though our values for the plasma temperature and bulk speed are motivated by galaxy cluster shocks, the results can be readily applied to other systems (e.g., the solar wind), as long as the dimensionless ratios and are the same and all the velocities remain non-relativistic. We will investigate the dependence of the results on the Mach number and the plasma beta in a forthcoming paper (X. Guo et al., in preparation).

We employ a spatial resolution of computational cells per electron skin depth , which is sufficient to resolve the Debye length of the upstream electrons for our chosen temperature of . We have tested that a spatial resolution of 7 cells per electron skin depth can still capture the electron heating physics. We use a time resolution of . Each cell is initialized with computational particles ( per species), but we have tested that a larger number of particles per cell (up to 64 per species) does not change our results (Appendix B). For the reference run, the transverse size of the computational box is , corresponding to , but we have tested that simulations with a transverse box size up to show essentially the same results.

Figure 1.— Shock structure and proton dynamics at . The coordinate is measured relative to the shock location , and it is normalized to the proton Larmor radius . From top to bottom, we plot: (a) the -averaged 1D profiles of proton density (black, in units of the upstream value), magnetic field (green, in units of the upstream field ) and total magnetic field strength (red, in units of the upstream field ); (b) the cross-shock electrostatic potential energy , in units of the proton upstream bulk energy ; (c)-(e) the proton phase spaces , , and , where the proton momentum is in units of and the proton thermal velocity is defined as ; (f) the proton temperature perpendicular (, blue line) and parallel (, orange line) to the magnetic field, and the mean proton temperature (green line); (g) the proton anisotropy (blue line) and the anisotropy upper bound in Equation (23) (red dashed line).
Figure 2.— 1D and 2D structure of magnetic fluctuations in our reference shock run at . In panel (a), we plot the energy of magnetic field fluctuations in the , and directions (blue, orange and green lines, respectively) normalized to the magnetic energy of the frozen-in field, which is defined as . Panels (b)-(d) show the 2D structure of the field fluctuations , and , respectively. The coordinate is measured relative to the shock location , and it is normalized to the proton Larmor radius . In panels (b)-(d), the coordinate is in units of the proton Larmor radius .

4. Shock Structure

In this section, we describe the structure of our reference shock run, with , and . We first discuss the proton dynamics and the generation of magnetic fluctuations sourced by the proton temperature anisotropy. Then, we present the electron dynamics and focus on the profile of electron irreversible heating. We will identify two main locations where the electron entropy increases: the shock ramp and the site where proton-driven waves grow in the downstream. The electron heating physics in these two regions will be investigated in Section 5 and Section 6, respectively.

4.1. Proton Dynamics and Proton-Driven Instabilities

In this subsection, we describe the proton dynamics, with a focus on proton isotropization and thermalization downstream of the shock. Figure 1 shows the profile of various quantities in the shock at time , as a function of the coordinate relative to the shock location , in units of the proton Larmor radius defined in Equation (21).

Panel (a) shows the -averaged profile of the proton number density in units of the proton density in the upstream (black line). The density compression at the shock reaches over a distance of , consistent with the expectation that the thickness of a perpendicular shock should be of the order of the proton Larmor radius (Bale et al., 2003; Scholer & Burgess, 2006). The density oscillates on a typical length scale of after the overshoot and then relaxes to the Rankine-Hugoniot value of beyond a distance of behind the shock.

The density pile-up at the shock is related to the electrostatic potential that develops in the shock transition region. This phenomenon has been well studied via hybrid simulations of collisionless shocks (e.g. Leroy et al., 1981, 1982; Leroy, 1983). As shown in

Figure 3.— b

), the potential energy reaches of the incoming proton energy . As a result, a significant fraction of the incoming protons are reflected back toward the upstream, leading to a pile-up of particles just in front of the shock. The reflected protons can be identified as the ones with positive and ahead of the shock in the phase spaces of

Figure 4.— c

) and (d), respectively. As the reflected protons gyrate in the shock-compressed magnetic field, they gain energy from the upstream motional electric field. Upon their second encounter with the cross-shock potential, the reflected protons now have sufficient energy to penetrate the shock. In the downstream region just behind the shock, the protons keep gyrating in the plane perpendicular to the shock-compressed magnetic field (compare the phase spaces in

Figure 5.— c

) and (d), at ). The peaks in density seen in

Figure 6.— a

) are then correlated with the locations where the proton gyro-phase is such that most protons have small (e.g., at , and ). The amplitude of the density oscillations gets smaller as the gyrating reflected protons become more and more phase-mixed with the directly transmitted protons, at .

Since the post-shock protons gyrate in the plane perpendicular to the field, the momentum dispersion along the direction of the field is expected to be nearly the same on the two sides of the shock (see the phase space in

Figure 7.— e

) near the shock). Further behind the shock, the dispersion in increases. This can be also quantified with the -averaged profiles of the proton temperature perpendicular () and parallel () to the background magnetic field, as in

Figure 8.— f

). Here, the component of the temperature tensor is defined as , where are the particle velocities in the fluid comoving frame, is the comoving particle Lorentz factor, and the average is performed over the particle distribution at a given spatial location. As

Figure 9.— f

) shows, the mean proton temperature , defined as222The factor of two that multiplies in the definition of comes from the fact that the perpendicular motion has two degrees of freedom.


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

Figure 10.— c


The decrease in perpendicular temperature, and the resulting increase in parallel temperature, suggests that protons are being scattered in pitch angle. In fact, in the region where the variation in and is most pronounced, strong magnetic waves are observed in

Figure 11.— T

heir wavelength is comparable to the proton skin depth, indicating that they are driven by protons (as opposed to electrons). In

Figure 12.— a

), we compare the 1D profiles (averaged over the direction) of the magnetic fluctuations , and , normalized to , where is defined as the magnitude of the flux-frozen magnetic field (i.e., , where is the -averaged particle density).333The frozen-in magnetic field is also used in the definition of . The energy of proton-driven waves peaks at . In

Figure 13.— a

), they are responsible for the excess of magnetic field strength (red curve) above the flux-freezing prediction (which would correspond to the density profile, in black).

The dominant mode at in the and direction has a wavevector nearly parallel to the background field (

Figure 14.— b

) and (c)), consistent with the proton cyclotron instability (Kennel, 1966; Davidson & Ogden, 1975). The waves in are slightly weaker (compare the green line with the blue and orange curves in

Figure 15.— a

)) and have oblique wavevectors (

Figure 16.— d

)), as expected for the mirror mode (Chandrasekhar et al., 1958; Barnes, 1966; Hasegawa, 1975; McKean et al., 1993). The presence of mirror modes breaks the flux freezing condition, as shown by the fact that in

Figure 17.— a

) the -averaged transverse magnetic field profile deviates at from the density profile (in black, which tracks the flux freezing prediction).

Both the proton cyclotron instability and the mirror instability are sourced by proton temperature anisotropy. In fact, since the motion of downstream protons right behind the shock is mostly confined in the plane, a large temperature anisotropy arises, with (

Figure 18.— g

)). The anisotropy provides free energy for the growth of proton cyclotron waves and mirror modes, which scatter the protons in pitch angle and reduce their anisotropy back to the upper bound corresponding to marginal stability Gary et al. (1997) (see the red dashed line in

Figure 19.— g

)), which is


Here, is the local value of the proton plasma beta, computed with the parallel proton temperature.

Figure 20.— Shock structure and electron dynamics at . From top to bottom, we plot: (a) the -averaged profiles of electron density (black) and total magnetic field strength (red); (b)-(d) the electron phase spaces , , and , where the electron momentum is in units of and the electron thermal velocity is ; (e) the electron temperature perpendicular (, blue) and parallel (, orange) to the magnetic field, and the mean electron temperature (green); (f) the electron anisotropy ; (g) the excess electron temperature over the adiabatic expectation for an isotropic gas; (h) the electron entropy profile, measured as in Equation (26).

4.2. Electron Dynamics and Heating

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

Figure 20(a) shows the electron density profile (black line), which strongly resembles that of the protons (black line in Figure 1(a)), and thus ensures approximate charge neutrality. While a small degree of charge separation at the shock is responsible for establishing the electric potential shown in

Figure 21.— b

), the fact that is nearly uniform at suggests that charge neutrality is satisfied very well in the far downstream.

Figure 20(b)-(d) shows the electron phase space. Since electrons have opposite charge than protons, they are not reflected back upstream by the cross-shock potential. In fact, unlike for protons, there is no reflected electron population with just ahead of the shock (compare

Figure 22.— b

) with

Figure 23.— c


Figure 20(e) shows the temperature profile of electrons, for the perpendicular component (blue), the parallel component (orange) and the mean temperature (green), which is defined as


The profile of perpendicular temperature (blue line) follows closely the density compression (compare with the black line in

Figure 24.— a

)) and starts to rise just ahead of the shock at . This is consistent with the double adiabatic theory Chew et al. (1956), which predicts (and in flux freezing, ).444We remark, as we have already pointed out at the end of Section 2, that the field strength should include all the magnetic contributions at scales much larger than the electron Larmor radius and at frequencies much lower than the electron gyration frequency. The double adiabatic theory applies to electrons, since the density and magnetic field compression occurs on scales much larger than the electron Larmor radius. This is not true for protons, since the shock thickness and the scale length of the downstream oscillations seen in

Figure 25.— a

) are set by the proton Larmor radius .

The parallel electron temperature (orange line in

Figure 26.— e

)) initially remains unchanged, as the CGL theory predicts and as a result of flux freezing (compare the green and black lines in Figure 1(a) in the vicinity of the shock). The increase in perpendicular temperature at the shock, while the parallel temperature stays the same as in the upstream, leads to a strong electron anisotropy, up to (

Figure 27.— f

)). This excites the electron whistler instability, which creates the small-wavelength transverse magnetic waves in and seen in the region of Figure 2(b) and (c) (see also the magnetic energy in and in

Figure 28.— a

), at the same location). The electron whistler instability provides a mechanism for electron pitch angle scattering and thus reduces the electron temperature anisotropy, as shown in the downstream region of

Figure 29.— f


As we have already discussed, if the electron fluid were to follow the double adiabatic predictions, and const. The fact that the perpendicular temperature profile in

Figure 30.— f

) (blue line) resembles the density profile (black line in

Figure 31.— a

)), and the fact that (compare green and blue curves in

Figure 32.— e

)), suggests that most of the increase in electron temperature comes from adiabatic compression. However, the fact that is not constant across the shock requires non-adiabatic processes. In order to quantify the degree of non-adiabatic (or, “irreversible”) electron heating, we compare in

Figure 33.— g

) the mean electron temperature with the adiabatic prediction


This estimate of the adiabatic temperature assumes an isotropic gas, which is valid, given the small degree of electron anisotropy far downstream of the shock (see

Figure 34.— f

) at ). Figure 20(g) shows the excess of above in units of the upstream electron temperature. Most of the irreversible heating occurs at two locations: , i.e., in the shock transition region; and , where the density suffers another compression, and strong proton-driven waves are generated (see

Figure 35.— .

These two locations are marked by the vertical dotted lines in

Figure 36.— n


Figure 37.— a

nd the particle and wave properties there will be further studied below. In the far downstream, the temperature excess over the adiabatic estimate saturates at (

Figure 38.— g


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


where the normalization is such that . To construct the spatial profile of , we first bin the particles by their position with a width of cells. In each spatial bin, we compute by constructing a three-dimensional histogram of the particle momenta. In each direction (i.e., , and ), the central bin of the histogram lies at the mean momentum, and the histogram spans four standard deviations above and below the mean. Each standard deviation is resolved with 10 momentum bins.

Figure 20(h) shows the change of electron entropy with respect to the upstream value. In analogy to Figure 20(g), the increase in electron entropy is localized around and (indicated by the gray and pink vertical dotted lines, respectively). The increase in electron specific entropy saturates at in the far downstream.

Figure 41.— Space-time diagrams and power spectra at a distance of ahead of the shock (as indicated by the vertical dotted grey lines in
Figure 42.— n
Figure 43.— ,
during the time interval . For this plot, the unit of time is the electron cyclotron time (the corresponding unit of frequency is ) and the unit of distance along is the electron skin depth (the corresponding unit for the wavevector is ). Panels (a)-(d) are the space-time diagrams of: (a) total magnetic field strength , in units of the upstream field ; (b)-(c) transverse magnetic field fluctuations and ; (d) electron anisotropy . Panels (e) and (f) show the power spectra of the field fluctuations presented in panels (b) and (c), respectively. In panels (e) and (f), the solid black line is the predicted real part of the frequency of electron whistler modes, whereas the dashed black line is the predicted imaginary part (i.e., the growth rate). The agreement between the prediction and our measurement confirms that the fluctuations in panels (a)-(c) are whistler waves.
Figure 48.— Space-time diagrams and power spectra at behind the shock (as indicated by the vertical dotted pink lines in
Figure 49.— n
Figure 50.— ,
during the same time interval as in
Figure 51.— F
or panels (a)-(f), see
Figure 52.— t
he only difference being that the predictions in panels (e) and (f) (solid black line for the whistler wave frequency, dashed black line for the growth rate) are computed considering the plasma properties only in regions where the electron anisotropy is well above the whistler threshold, more specifically . In panel (g), where we indeed plot , this would correspond to the dark green areas. Since panels (a)-(c) are dominated by long-wavelength slowly-propagating proton modes, we isolate electron waves via a high-pass filter in the power spectra of panels (e) and (f), keeping only the high- high- region delimited by the red dashed lines. The resulting space-time wave pattern is shown in panels (h) and (i), which reveal the presence of electron whistler waves.

4.2.1 Electron Whistler Waves

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

Figure 53.— n


Figure 54.— .

At the shock (grey dotted line in

Figure 55.— n


Figure 56.— ,

the electron temperatures are driven to by shock-compression of the upstream field, via conservation of the first and second adiabatic invariants. In

Figure 57.— w

e show the space-time diagram of various quantities, in the time interval and along the extent of the box. The location is fixed at the shock ramp (more precisely, ). Shock-compression of the upstream field (see Figure 43(a), where ) leads to a temperature anisotropy (Figure 43(d)). Both the field amplification and the degree of temperature anisotropy are nearly constant in time and uniform in .

As a result of the strong temperature anisotropy, magnetic waves are excited throughout the range consistently over time. Panels (b) and (c) show the space-time diagrams of the magnetic fluctuations and , revealing the presence of high-frequency and short-wavelength modes (as also seen in

Figure 58.— b

) and (c) near the shock). Figure 43(e) and (f) show the corresponding power spectra, as a function of frequency (horizontal axis) and wavenumber (vertical axis). The power spectrum displays a pronounced peak at frequency (here is the electron gyrofrequency) and wavevector . We have compared this with linear theory of the electron whistler instability (e.g. Gary & Madland, 1985; Gary & Wang, 1996; Gary & Karimabadi, 2006) by solving the dispersion relation


where and is the plasma dispersion function


The input values of for the dispersion relation are taken from the time- and space-averages of the corresponding quantities over the same time period and spatial extent as the space-time diagram in Figure 43. The resulting theoretical prediction for the real part of the frequency is shown with a black solid line in panels (e) and (f), and it matches extremely well the contours of the power spectrum. The imaginary part of the frequency, i.e., the growth rate of the mode, is plotted with a dashed black curve. The value of giving the fastest growth agrees well with the location of the peak of the power spectrum (). The excellent agreement between the simulation data and the electron whistler dispersion relation confirms that the waves in the shock ramp are produced by the electron whistler instability.

Figure 52 shows similar plots at the location indicated in

Figure 59.— n


Figure 60.— i

th a vertical dotted pink line, at . Here, field amplification is driven by a combination of two effects: the density (and so, the frozen-in magnetic field) experiences another large-scale compression; in addition, the proton-driven waves shown in

Figure 61.— u

rther increase the local magnetic field intensity.

As compared to

Figure 62.— t

he space-time diagrams show now a higher degree of inhomogeneity, imprinted by the anisotropy-driven long-wavelength proton modes. These fluctuations co-exist with weaker small-wavelength high-frequency modes, which only appear in localized patches (e.g., at and in

Figure 63.— b

) and (c)). The high-frequency waves are generated in regions where field amplification (

Figure 64.— a

) at and ) causes the electron anisotropy (

Figure 65.— d

)) to exceed the threshold for whistler growth (

Figure 66.— g

)), which is given by


where is the local value of the parallel electron beta Gary (2005).

Figure 52(e) and (f) show the power spectra of and . Most of the power is concentrated in low-frequency long-wavelength modes, generated by the proton cyclotron or mirror instabilities. However, there is still an appreciable amount of power in high-frequency short-wavelength modes peaking at and . We apply a high-frequency short-wavelength filter, in order to isolate the top right region in

Figure 67.— e

) and (f) (the cutoff frequency and wavenumber of our filter are shown with dashed red lines). This allows us to extract (via an inverse Fourier transform) the space-time wave patterns of high-frequency short-wavelength modes, which are shown in panels (h) and (i). The two panels confirm that short-wavelength modes exist only in regions where the electron temperature anisotropy exceeds the electron whistler threshold (Figure 52(g) at and ). We have measured the average electron and proton temperatures and densities in the region where the whistler threshold is appreciably exceeded (), in order to obtain linear theory predictions. The real part and imaginary part of the resulting dispersion relation are plotted in

Figure 68.— e

) and (f) with solid and dashed black lines, respectively. The good agreement with the power spectra extracted from our shock simulation confirms the presence of patches of whistler waves in the second ramp (at ) of the electron entropy profile.

To summarize, we have identified two major sites of electron entropy production in the shock downstream. One is at the shock ramp, and the other is at a distance of behind the shock, where density compression and proton-driven waves both contribute to magnetic field amplification. Both sites show the presence of electron whistler waves triggered by electron temperature anisotropy. Whistler waves provide the pitch-angle scattering required to break electron adiabatic invariance and to generate entropy. In the following two sections, we further elucidate the physics of entropy production in these two sites, by means of periodic box simulations.

5. Electron Heating in the Shock Ramp

The first increase in electron entropy happens in the shock ramp. As a result of the shock-compression of the upstream field ( by flux freezing), electrons become anisotropic and they trigger whistler waves. Below, we model the shock compression in a periodic box using a novel form of the PIC equations introduced in Sironi & Narayan (2015); Sironi (2015), which incorporates the effect of a large-scale compression of the system. We briefly describe the simulation setup in Section 5.1, we discuss periodic box simulations applicable to our reference shock run in Section 5.2, and we describe the dependence on mass ratio in Section 5.3.

Figure 69.— As a function of the comoving time of the electron fluid defined in Equation (31), we present the density profile experienced by electrons as they propagate from upstream to downstream (solid blue line). The time axis is shifted such that just ahead of the shock. The shock-compression felt by incoming electrons can be approximated as with compression rate (orange dashed line).

5.1. Simulation Setup

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


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

A uniform ordered magnetic field is initialized along the direction (in analogy to the shock setup). We define as the proton Larmor frequency in the initial field . Maxwell’s equations in the primed coordinate system (Sironi & Narayan, 2015) prescribe that the field will grow in time as , which is consistent with flux freezing (the particle density in the box increases at the same rate). From the Lorentz force in the compressing box (Sironi & Narayan, 2015), the component of particle momentum aligned with the field does not change during compression, whereas the perpendicular momentum increases as . This is consistent with the conservation of the first and second adiabatic invariants.

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

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

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


where is the electron fluid velocity in the shock frame, and the integral goes from the upstream to the downstream region. Figure 69 shows the density profile as a function of from our reference run. The density oscillates on a timescale comparable to the proton gyration time , which is expected given that the proton dynamics controls the shock structure. At the ramp starting near , the electron density increases by a factor of within . Even though the density increase is not perfectly linear, we find that a linear approximation with provides a reasonable fit (see the orange dashed line in

Figure 70.— .

We remind that our electron heating model is agnostic of the exact profile of density compression, as long as the compression rate and the resulting field amplification rate are much slower than the electron gyration frequency.

Below, we fix . With increasing mass ratio, the separation between and the electron cyclotron frequency will increase as . As in our reference shock run, electrons are initialized to have , and the strength of the background magnetic field is set so that . We resolve the electron skin depth with 10 cells, so the Debye length is marginally resolved. The box extent along the direction is fixed at , which is sufficient to capture several wavelengths of the electron whistler instability. The box size does not need to scale with the mass ratio, since we are artificially excluding the proton physics.

Figure 71.— Time evolution of various space-averaged quantities in a 1D periodic box whose compression rate is chosen to mimic the effect of the shock ramp. We compare two field geometries, with background field lying either along the axis of the simulation box (“in-plane” configuration, solid lines) or along the direction perpendicular to the box (“out-of-plane” configuration, dotted lines): (a) energy in magnetic field fluctuations, normalized to the energy of the compressed magnetic field (the legend is appropriate for the in-plane configuration, whereas for the out-of-plane case the orange line refers to ); (b) electron temperature perpendicular (, blue lines) and parallel (, orange lines) to the background field; (c) electron temperature anisotropy (blue lines), and comparison with the threshold of the electron whistler instability, as in Equation (29) (dashed red line); (d) electron entropy change, measured from the electron distribution function as in Equation (26) (blue solid) or predicted from Equation (33) (red dashed); (e) electron energy increase in units of , measured directly (blue solid) or predicted using Equation (32) (red dashed).

5.2. Application to the Reference Shock

As in our reference shock run, we employ a reduced mass ratio . In the periodic compressing box, this means that our choice of leads to a compression rate that is a factor of lower than the electron gyration frequency.

To highlight the importance of the electron whistler instability in facilitating electron entropy production, in Figure 71 we compare two simulations, one with the background field in the direction, the other one with along the direction. Since our simulation box is oriented along and the dominant wavevector of the electron whistler instability is parallel to the background field, if the field lies along (which we shall call “out-of-plane” case, and indicate with dotted lines) we artificially suppress the growth of electron whistlers. By comparing it with the “in-plane” simulation with the field along (solid lines), which does allow for whistler wave growth, we can demonstrate the importance of the electron whistler instability for entropy production.

In the absence of electron-scale instabilities that would break the adiabatic invariance, the out-of-plane simulation is expected to follow adiabatic scalings. In fact, in the out-of-plane simulation, we see that (blue dotted line in

Figure 72.— b

)), while const (orange dotted line in

Figure 73.— b

)), as expected from the double adiabatic theory. Since no whistler waves grow (notice that the fields stay at the noise level, see the dotted lines of

Figure 74.— a

)), no mechanism exists that can transfer heat from the perpendicular to the parallel temperature, and the electron entropy remains constant.

The in-plane simulation shows a different behavior. Initially, and follow the double adiabatic trends (solid lines in

Figure 75.— b

)). At , the increasing temperature anisotropy (blue solid line in panel (c)) leads to the exponential growth of electron whistler waves (solid lines in

Figure 76.— a

)). At , the wave energy is strong enough to pitch-angle scatter the electrons. As a result, heat is transferred from the perpendicular to the parallel direction. Both and deviate from the adiabatic scalings and the temperature anisotropy is reduced.

At , with the onset of pitch-angle scattering and the consequent breaking of adiabatic invariance, the electron specific entropy starts to rise (solid blue line in Figure 71(d)). The most rapid entropy increase happens near the end of the exponential whistler growth, at . Here, the electron anisotropy is still large, and at the same time whistler waves are sufficiently powerful to provide effective pitch-angle scattering. In other words, both terms in the square brackets of either Equation (13) or Equation (14) are large. After the exponential growth, the electron whistler waves enter a secular phase where the wave energy (normalized to the compressed background field energy) stays almost constant (solid green line in Figure 71(a)). In this phase, whistler waves are continuously generated as the large-scale compression steadily pushes the electron anisotropy slighly above the threshold of marginal stability (indicated by the red dashed line in

Figure 77.— c

)). Both the ingredients needed for entropy increase (i.e., nonzero electron anisotropy and efficient pitch angle scattering mediated by whistler waves) persist during the secular phase, leading to further increase in the electron entropy.


Figure 78.— w

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


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

The validation can also be extended to the entropy measurement, as we do in

Figure 79.— d

). Again, the blue solid line shows the measured change in electron specific entropy (computed from the distribution function as in Equation (26)), while the red dashed line is obtained by integrating


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

Figure 80.— Dependence on mass ratio (up to ) of various space-averaged quantities in a 1D periodic box with compression rate (the legend is in panel (b)). The background field is aligned with the box (in-plane configuration). We plot: (a) energy in magnetic field fluctuations, normalized to the energy of the compressed field; (b) electron temperature anisotropy (solid lines) and threshold condition for the electron whistler instability (dotted lines with the same color coding as the solid lines); (c) rate of violation of adiabatic invariance ; (d) electron entropy change, measured from the electron distribution function as in Equation (26). After (vertical dotted black line in panel (d)), which corresponds to the end of the compression phase in the shock ramp, the entropy change is nearly independent of the mass ratio.

5.3. Dependence on the Mass Ratio

We now extend our compressing box experiments up to the realistic mass ratio and show that the electron entropy increase is nearly insensitive to (as long as the mass ratio is larger than a few tens). Figure 80 compares the evolution of the whistler wave energy (panel (a)), the electron temperature anisotropy (panel (b)), the rate of breaking adiabatic invariance (panel (c)) and the electron entropy increase (panel (d)) when varying the mass ratio from up to (from purple to red, see the legend in the second panel). Since we fix the compression rate to be , a larger mass ratio corresponds to a lower compression rate in units of the electron gyration frequency .

Initially, the electron temperature anisotropy grows linearly in time as , as a result of the large-scale compression. This proceeds until the energy in whistler waves reaches a fraction of the compressed background field energy (

Figure 81.— a

)). At this point, whistler waves are sufficiently strong to scatter the electrons in pitch angle, breaking their adiabatic invariance and reducing the electron anisotropy by transferring energy from the perpendicular to the parallel component. In fact, notice that the peak in panel (c), i.e., the time when the electron adiabatic invariance is most violently broken, always corresponds to the time when the electron anisotropy in panel (b) shows the sharpest decrease.

The onset of efficient pitch-angle scattering (and so, the peak time of electron anisotropy) occurs earlier at higher mass ratio, at a time that decreases from at down to at . This can be understood from the competition between the large-scale compression rate (which increases the electron anisotropy) and the growth rate of whistler waves (which try to reduce the anisotropy via pitch angle scattering). The compression rate in units of the electron cyclotron frequency is , while the whistler growth rate (also in units of ) depends on how much the anisotropy exceeds the whistler threshold in Equation (29). In order to balance the two rates, a higher anisotropy is needed for larger , i.e., for lower mass ratios. This has two consequences: first, the growth rate of the whistler instability (normalized to ) will decrease at higher , as indeed confirmed by the inset of

Figure 82.— a

); second, lower peak anisotropies (and so, earlier onsets of efficient pitch angle-scattering) will be achieved at higher mass ratios, which explains the trend seen in

Figure 83.— b

). In addition, since the energy of whistler waves ultimately comes from the free energy in electron anisotropy, higher mass ratios display weaker levels of whistler wave activity (panel (a)).

The electron entropy evolution in

Figure 84.— d

) can be separated into two stages. In the first phase (which, for , occurs at ), the electron entropy grows quickly. This stage corresponds to the late exponential phase of whistler wave growth (and so, we shall call it “exponential phase”), when both the electron anisotropy (panel (b)) and the rate of breaking adiabatic invariance (panel (c)) — i.e., the two ingredients needed for efficient entropy production — are large. Since higher mass ratios reach lower levels of electron anisotropy, the entropy produced during this stage is a decreasing function of , as seen in

Figure 85.— d

) (compare the purple line growth around with the red line around ). After whistler waves have reached saturation, the electron entropy still increases, in a phase which we shall call “secular”. Here, the electron anisotropy stays close to the threshold of marginal stability (indicated in

Figure 86.— c

) by the dotted lines, with the same color coding as the solid curves). Continuous pitch-angle scattering (and so, persistent violation of adiabatic invariance) is needed to oppose the steadily-driven compression and maintain the system close to marginal stability. It is then expected that entropy will continuously increase during the secular phase, albeit at a lower rate than in the exponential stage. For , the electron anisotropy at late times is nearly insensitive to