Kinetic Simulations of Strongly-Magnetized Parallel Shocks: Deviations from MHD Jump Conditions

In a collisional fluid, dissipation at the front of a shockwave is provided by binary collisions. As a consequence, the shock front is a few mean free path thick (Zel’dovich & Raizer, 2002). Yet, in plasma, shock waves can be sustained by collective effects, with a front smaller than the mean free path by several orders of magnitude (Petschek, 1958; Buneman, 1964; Sagdeev, 1966).

A good example is the bow shock formed at the boundary of the Earth’s magnetosphere and the solar wind, in which the shock front thickness is on the order of 100 kilometers, while the proton mean free path is comparable to an AU (Bale et al., 2003; Schwartz et al., 2011).

Such shocks are omnipresent in astrophysical systems and are good potential candidates to explain the origin of cosmic rays, gamma-ray bursts and fast radio bursts; as such, they have been under intense scrutiny for decades (e.g., Blandford & Eichler, 1987; Piran, 2005; Sironi et al., 2021). Because of their collisionless characteristics the fluid formalism cannot be straightforwardly applied to study them. Consequently, Particle-In-Cell (PIC) simulations, based on a model that evolves the equations of collisionless, kinetic plasmas, has become an important tool employed in the investigation of shocks in recent years (e.g., Hoshino et al., 1992; Dieckmann et al., 2000; Spitkovsky, 2008a; Sironi & Spitkovsky, 2011a; Caprioli & Spitkovsky, 2014; Park et al., 2016). The fluid properties of these kinetic shock simulations, such as the compression ratio and downstream heating, are often measured by their corresponding MHD predictions, even though the relevance of MHD for collisionless shocks is not straightforward (Bret, 2020; Haggerty & Caprioli, 2020).

In this respect, parallel shocks are particular interesting when comparing the kinetic and MHD picture. In the MHD picture, the strength of a parallel magnetic field does not affect the hydrodynamics of shock; in this configuration, the field and the fluid are completely disconnected (Lichnerowicz (1976), or Kulsrud (2005), p. 141). Therefore, any detectable effect of a parallel field on a collisionless shock represents a departure from MHD and must be due to kinetic effects.

In a series of recent theoretical works collisionless (electron-positron) plasma shock hydrodynamics were reconsidered from a kinetic framework (Bret et al., 2017; Bret & Narayan, 2018, 2020), in which the density jump of a collisionless shocks was predicted deviate from the standard MHD, Rankine Hugoniot jump conditions; namely that an increasing the magnetic field strength in a parallel shock should decrease the compression ratio. In this paper, we constrain this prediction with state-of-the-art PIC simulations, and we examine the effects of increased field strength on the thermal and nonthermal properties of the shock. The paper is outlined as follows: In Sec. 2 we review the theory and predictions for the departure from MHD, in Sec. 3 we discuss the PIC model and the setup for the shock simulations, in Sec. 4 the fluid properties of the are examined, namely the density, the anisotropic temperature and the entropy, and are shown to be in good agreement with theory, in Sec. 5 we detail how the magnetic field strength modified the shock front’s width , and finally we conclude in Sec. 6.

For completeness we briefly summarize the model developed for the density jump in Bret & Narayan (2018), hereafter BN18. We considered a non-relativistic shock in an electron-positron (pair) plasma. The upstream is assumed isotropic. The external magnetic field $\mn@boldsymbol{B}_{0}$ is parallel to the shock normal, as shown in Fig. 1. The upstream density, pressure and temperature are $n_{1},P_{1}$ and $T_{1}$ respectively. Downstream quantities are labeled with the subscript “2”, and the downstream are assumed to be anisotropic with $T_{2x}\neq T_{2y}$. We define the downstream anisotropy parameter $A_{2}$ as,

The field strength is measured by the $\sigma$ parameter (magnetization),

where $v_{A0}$ is the Alfvén speed and $v_{1}$ the upstream velocity in the shock frame.

While the standard MHD jump conditions assume isotropy, several authors have adapted them to account for a downstream pressure anisotropy $A_{2}$ (Karimabadi et al., 1995; Vogl et al., 2001; Erkaev et al., 2000; Gerbig & Schlickeiser, 2011). Yet, in these studies, $A_{2}$ is a free parameter leaving the problem under-determined and an additional assumption is required to pinpoint $A_{2}$ in terms of the upstream quantities (Gary, 1993; Bale et al., 2009; Maruca et al., 2011).

The additional assumption made in BN18 is that the transit through the shock front is adiabatic in the perpendicular direction while the entropy increase goes into the parallel temperature. Regarding the perpendicular temperatures, the relations derived in (Chew et al., 1956) imply their conservation. Together with the macroscopic conservation equations, this assumption allows to fully determine the characteristics of the downstream, including its temperature anisotropy $A_{2}$.

Depending on the value of $A_{2}$ thus obtained, the downstream can be firehose unstable if the temperature anisotropy is large relative to the magnetic field strength. For weak fields, i.e., where $\sigma$ smaller than some critical value $\sigma_{c}$, the downstream is firehose unstable, while for strong enough fields, with $\sigma>\sigma_{c}$, the downstream remains stable and constitutes the end state of the shock. The resulting density jump ($r$) can be derived,

where,

The model explained in BN18 made use of the parameter $\chi_{1}$, which is a pseudo Mach number. The reason why the model is not adapted to the use of a proper Mach number can be understood by examining Eq. (3), which shows that for $\sigma=0$ and $\chi_{1}\gg 1$, the density jump reads $r=4$, while for $\sigma>\sigma_{c}$ and $\chi_{1}\gg 1$, we find $r=2$. Therefore, the model presents a transition from a 3D strong shock with $\gamma=5/3$ at $\sigma=0$, to a 1D strong shock with $\gamma=3$ at $\sigma>\sigma_{c}$ (Bret, 2021). As a result, within the model of BN18, it is impossible to define a Mach number using a fixed adiabatic index since the effective adiabatic index of the model varies between 5/3 and 3.

In the rest of this work, the parameter $\chi_{1}$ is related to the real Mach number $\mathcal{M}$ through,

where $\gamma$ is the adiabatic index of the plasma. Accordingly, the downstream anisotropy $A_{2}$ reads,

Notably, in the strong shock limit $\chi_{1}=\infty$, Eqs. (3,6) reduce to,

and

with, as expected, $r_{\infty}(\sigma=0)=4$ and $A_{2\infty}(\sigma=0)=1$.

In order to test these predictions we perform a survey of non-relativistic 3D particle-in-cell (PIC) simulations of electron-positron (pair-plasma) shocks with TRISTAN-MP (Buneman, 1993; Spitkovsky, 2005). The choice to model a pair-plasma was made to more directly compare with the theory developed in BN18 which was derived for an electron-positron plasma. This choice was motivated by two reasons: the temperature of the two species would have the same value downstream of the shock and the instabilities, relevant for this problem, are unchanged in the electron-positron limit (Gary & Karimabadi, 2009; Schlickeiser, 2010). In the ion-electron limit, the temperature of the two species is not necessarily equal (Guo et al., 2017, 2018; Tran & Sironi, 2020), however, examining this problem in this framework would require a description of the ion/electron energy partition. While this is an important step in understanding the physics of non-relativistic collisionless shocks, it is beyond the scope of this work

The setup of the shock simulations are similar to those described in Spitkovsky (2008b); Sironi & Spitkovsky (2011b); Sironi et al. (2013), where the shock is formed by initializing the plasma with a supersonic flow towards a reflecting wall (in the negative $x$ direction) on the left hand side of the simulation. The plasma is reflected and forms a shock that is mediated by some counter streaming instabilities. The plasma is continuously injected by an injector receding away from the shock which extends the simulation domain. The boundary conditions along $y$ and $z$ are periodic.

Values in the simulations are measured in normalized code units: speeds are normalized to the speed of light, times to the inverse plasma frequency $\omega_{pe}^{-1}=\sqrt{m/4\pi ne^{2}}$, where $e$ is the fundamental charge and $n$ is the number density of the inflowing plasma, and lengths, consistently, to the electron inertial length $d_{e}=c/\omega_{pe}$. The maximum allowed simulation length (i.e., where the receding injector stops) is $L_{x}=1200c/\omega_{pe}$, and the transverse width and depth are $L_{y}=L_{z}=16c/\omega_{pe}$ for all simulations. The simulations are performed with $2$ electron/positron pairs per cell, $4$ grid cells per $c/\omega_{pe}$, and with the time step chosen such that $\Delta t=.45\Delta x/c$. A convergence simulation was run for the $M_{A}=M_{s}$ case, with double the transverse box size, twice the spatial resolution and $\times 4$ particles per cell, and no significant difference in the results found.

The inflowing electrons and positrons are initialized as a drifting Maxwellian with equal temperature, $\Delta\gamma\equiv k_{B}T/mc^{2}=10^{-4}$. Because of the reflecting wall, the downstream plasma is at rest in the simulation frame. In this frame, the upstream plasma is in-flowing at $v_{sh}/c=0.1$, which corresponds to a Mach number of $\mathcal{M}=v_{sh}/\sqrt{\gamma P/n}\approx 7.75$. However, the shock predictions and theory outlined above are determined in the rest frame of the shock. The two frames are ultimately related by the density jump $r$, i.e., $v_{1}=v_{sh}/(1-1/r)$, as discussed in many previous simulations of shocks formed by reflecting walls (Caprioli & Spitkovsky, 2014; Haggerty & Caprioli, 2019, 2020).

The simulations are initialized with a uniform magnetic field anti-parallel to the upstream flow (i.e., $+x$). The initial magnetic field strength varies between simulations and is characterized by the shock energy density in the downstream frame $\sigma^{\prime}=B_{1}^{2}/4\pi n_{1}mv_{sh}^{2}$ and is related to the value described by Eq. 2 by the density jump, $\sigma=\sigma^{\prime}(1-1/r)^{2}$. The survey consists of 9 simulations performed over a range of $\sigma^{\prime}$ values from $10^{-3}$ to $25$, allowing for the consideration of both super-Alfvénic and sub-Alfvénic shocks. The simulations are evolved for 1000’s of $\omega_{pe}^{-1}$ but are stopped before particles can escape the simulation.

Various hydrodynamic shock quantities can be determined from the survey of simulations, and in the following section they are presented and compared with the theory of BN18.

From the density profiles shown in Fig. 2(a), it is clear that the width of the shock has some dependence on $\sigma$. To further examine this effect, we determined the width of the shock fronts for each simulation at $6000\,\omega_{pe}^{-1}$. This was achieved by fitting the shifted, average density profiles with a ${\rm tanh}$ function of the form,

The value of $\lambda$ for each simulation is determined with a non-linear least squares fit and the results are shown in Fig. 6.

Although the width of the shock varies non-monotonically with $\sigma$, two distinct patterns appear to emerge when considering simulations with $\sigma\ll 0.1$ or $\sigma\gg 0.1$. In the small $\sigma$ case, the shock width seems to reach a constant size for decreasing $\sigma$. This is consistent with the shock formation being mediated by the two stream instability for the present parameters (Bret et al., 2008), which has a growth rate independent of $\sigma$.

For $\sigma\gg 0.1$, the shock width increases with $\sigma$. This increase is potentially due to the increasing number of pitch angle scattering required to return particles traveling up-streaming back towards the downstream. Each deflection is expected to alter a particle’s pitch angle by $\Delta mu\sim\delta B/B_{0}\sim 1/\sqrt{\sigma}$. For a larger value of $\sigma$, more deflections will be required to send particles back potentially leading to a larger shock width. This process is shown in an animation of self-consistent particle tracing for the $\sigma^{\prime}=2$ simulation (orange line in this papers figures) in the supplementary material.

In this work, using self-consistent kinetic PIC simulations, we verify the departure of parallel collisionless shocks from fluid MHD predictions; by varying the magnetization on collisionless pair plasma shocks we show that a parallel propagating shock wave is suppressed for sufficiently strong magnetic field strengths, consistent with kinetic theoretical predictions of BN18. For a fixed (parallel) shock inclination and sonic Mach number ($v_{1}/v_{th}=10$), we test the predictions of BN18 of how the shock properties depend on the strength of the upstream magnetic field relative to the upstream shock speed, a quantity parametrized by $\sigma\equiv B_{0}^{2}/4\pi n_{1}mv_{1}^{2}$. The simulations are found to be in good agreement with the theoretical predictions, namely that as $\sigma$ is increased to $\gtrsim 1$, the shocked plasma remains well magnetized, suppressing the firehose instability and leading to a large field aligned temperature anisotropy downstream. The entropy jump between the downstream and upstream is considered and found to be approximate agreement with predictions.

Additionally we examine how the kinetic properties of the shock are modified by the magnetization. The effective shock width is determined for a fixed time for each of the simulations, and the width is found to be a minimum for an upstream plasma beta of unity ($\beta\sim 1$). For decreasing and increasing values of $\sigma$ the shock width broadens, but saturating for small values of $\sigma$ as the instability mediating the shock reformation becomes two stream like and insensitive to the magnetic field strength.

In the theory and simulations of this work, only a pair plasma was considered; this choice simplified the problem by ensuring that the downstream temperatures would be the same for both species. In shocks comprised of electron-proton plasmas, it is likely that the temperature of each species will be different (Tran & Sironi, 2020). This added complexity should be further explored to insure the applicability of this work to space and astrophysical shocks and will be considered in future investigations.

These results contribute to the building consensus collisionless shocks fundamentally require a kinetic description to accurately capture both the thermal and non-thermal plasma physics.

The datasets used in this manuscript were derived using the open source particle-in-cell software Tristan-MP v2 (https://ntoles.github.io/tristan-wiki/). The the average output for each simulation requires $\sim 14$ GB of space ($\gtrsim 200$ GB for all of the simulations and convergence studies) and as such either the data underlying this article or the simulation initialization files will be shared on reasonable request to the corresponding author.

We thank the anonymous referee for thorough comments which helped improve the manuscript. This work has been achieved under projects ENE2016-75703-R from the Spanish Ministerio de Economía y Competitividad and SBPLY/17/180501/000264 from the Junta de Comunidades de Castilla-La Mancha. As well as the FDSS NSF AGS-1936393 grant to the Institute for Astronomy of the University of Hawaii. Simulations were performed on computational resources provided by the University of Chicago Research Computing Center and XSEDE TACC (TG-AST180008).