I N T RO D U C T I O N

In a collisional fluid, dissipation at the front of a shock wave 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 of the order of 100 km, while the proton mean free path is comparable to an AU (Bale, Mozer & Horbury 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 E-mail: colbyh@hawaii.edu these kinetic shock simulations, such as the compression ratio and downstream heating, are often measured by their corresponding magnetohydrodynamic (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 (electronpositron) 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 shock was predicted to 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 non-thermal properties of the shock. The paper is outlined as follows: In Section 2, we review the theory and predictions for the departure from MHD, in Section 3 we discuss the PIC model and the set-up for the shock simulations, in Section 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 Section 5 we detail how the magnetic field strength modified the shock front's width, and finally we conclude in Section 6.

T H E O RY

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 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 labelled with the subscript '2', and the downstream are assumed to be anisotropic with T 2x = T 2y . We define the downstream anisotropy parameter A 2 as,

The field strength is measured by the σ 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, Krauss-Varban & Omidi 1995;Erkaev, Vogl & Biernat 2000;Vogl et al. 2001;Gerbig & Schlickeiser 2011). Yet, in these studies, A 2 is a free parameter leaving the problem underdetermined and an additional assumption is required to pinpoint A 2 in terms of the upstream quantities (Gary 1993;Bale et al. 2009;Maruca, Kasper & Bale 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, Goldberger & Low (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 σ smaller than some critical value σ c , the downstream is firehose unstable, while for strong enough fields, with σ > σ 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 χ 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 equation (3), which shows that for σ = 0 and χ 1 1, the density jump reads r = 4, while for σ > σ c and χ 1 1, we find r = 2. Therefore, the model presents a transition from a 3D strong shock with γ = 5/3 at σ = 0, to a 1D strong shock with γ = 3 at σ > σ 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 χ 1 is related to the real Mach number M through,

where γ is the adiabatic index of the plasma. Accordingly, the downstream anisotropy A 2 reads,

Notably, in the strong shock limit

and

S I M U L AT I O N S

In order to test these predictions we perform a survey of nonrelativistic 3D 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, Sironi & Narayan 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 set-up of the shock simulations are similar to those described in Spitkovsky (2008b), Sironi & Spitkovsky (2011b), Sironi, Spitkovsky & Arons (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 ω -1 pe = m/4π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/ω pe . The maximum allowed simulation length (i.e. where the receding injector stops) is L x = 1200c/ω pe , and the transverse width and depth are L y = L z = 16c/ω pe for all simulations. The simulations are performed with 2 electron/positron pairs per cell, 4 grid cells per c/ω pe , and with the time-step chosen such that t = .45 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 × 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, γ ≡ 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 M = v sh / √ γ P /n ≈ 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 antiparallel 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 σ = B 2 1 /4πn 1 mv 2 sh and is related to the value described by equation ( 2) by the density jump, σ = σ (1 -1/r) 2 . The survey consists of nine simulations performed over a range of σ 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 ω -1 pe but are stopped before particles can escape the simulation.

H Y D RO DY NA M I C S

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.

Density jump

Fig. 2 shows the shock density profile integrated in the transverse directions, for σ up to 25. Data are plotted at time 6000ω -1 pe to insure that the shock is fully formed. The jump obtained for values of σ 1 is 4 consistent with the MHD prediction for strong shocks. Although the critical magnetization is close to unity (σ c ∼ 1), a clear departure from the standard MHD prediction is observed for σ small as 0.1. In this regime, where 0.1 ࣠ σ ࣠ 1, the Mach number based on the Alfvén speed would range between 1 and 10, and so the decreasing compression ratios seems consistent with the decreasing Mach number.

The lower plot in Fig. 2 displays the density jump for a range of σ values, contrasted against equation (7) represented by the black dashed line. Simulations show the downstream density enhancement decreasing from 4 to 2, consistent with prediction presented above. The theoretical predictions attribute the departure from the MHD prediction to the stabilization of downstream temperature anisotropy, which can be directly examined from the summations and is discussed in Section 4.4.

Entropy

For simulations with relatively small values of σ σ c , it is clear that a shock in the typical sense is forming, i.e. the upstream and downstream regions are separated by a sharp discontinuity, and the downstream plasma is compressed consistent with the MHD jump conditions; however, for simulations with σ σ c it is less clear if a bona fide shock has developed, or we are witnessing interpenetrating fluids. Ultimately, shocks must increase entropy, and as such the issue can be addressed by the entropy difference between the upstream and the downstream. While there remain unsettled issues in the interpretation of entropy and reversibility in a collisionless system, or even effective collisions introduced by poor counting statistics of the PIC model (e.g. Montgomery & Nielson 1970;Kadanoff 2017;Shalaby et al. 2017;Yang et al. 2017;Liang et al. 2019;Du et al. 2020, and references therein), we limit considerations to the definition of kinetic entropy density per particle:

where S is the kinetic entropy density, f is the distribution function. This calculation was explicitly performed in section 4 of BN18 and considered again in Bret ( 2021) where the proof is made that for large values of χ 1 , the entropy jump for both the large and small σ regimes is definitely positive. Since the marginal firehose stability is observed in the present simulations at small σ , while the perpendicular temperature conservation is also observed far larger σ 's (see Section 4.4), the calculations of BN18 and Bret (2021) should apply, ensuring an entropy increase at the transition. The measured entropy jump ( s) for all of the simulations as a function of σ is shown in panel (a) of Fig. 3, with the black dashed line showing the theoretical prediction. The value is calculated by taking the difference of integrated discretized distribution functions averaged over a region of 100 d e in the downstream and upstream. There is good quantitative agreement for simulations with σ < σ c , and as σ is increased and the downstream remains firehose stable across the shock, there is qualitative with theory, i.e. the difference in entropy reduces. It is worth noting that the entropy increase for the σ ≈ 6 simulation agrees well with theory; however, in the simulation with an even larger magnetization, σ ≈ 14, the entropy goes to nearly zero ( s ≈ 0). This difference can be understood by considering the parallel downstream distribution functions for these simulations, shown in panel (b) of Fig. 3. In the entropy derivation of BN18, the distribution was assumed to have the form of a bi-maxwellian, while it is clear from the dark red line in the lower panel, that the distribution has two clear peaks associated with counter-streaming beams. However, it is likely that -if the simulation were run for longer times in a larger box -the double-peaked distribution function would tend towards the bi-maxwellian state, as the distribution function has already gaussianized relative to the initial distribution. From this it seems reasonable that the entropy will increase across the discontinuity, thus constituting a collisionless shock even in large-σ systems.

Firehose stability

The density jump measured in the simulations is in good agreement with the theoretical predictions from BN18 and to further verify the consistency between simulations and BN18, we also test the assumptions that form the basis of the model. The main ingredient 10) (left side, A, green line; right side, firehose, maroon) superimposed on the density normalized to the downstream value (black line) for σ = 0.1 (σ ≈ 0.056). The system is firehose unstable in regions where 1 -1/β ࢱ > A 2 , corresponding to the yellow shaded region in the upstream, in contrast to the marginal stability in the downstream. Panel (c):

The averaged downstream temperature anisotropy (right axis, red line with squares) and firehose parameter (left axis, black line with triangles) shown as a function of σ . The red dashed line corresponds to prediction for A 2 in the strong shock limit, as defined in equation ( 6).

of the BN18 model is that for weak magnetic field strengths, the plasma turns firehose unstable at the shock front crossing. As a result, it migrates to the stability threshold so that the downstream is marginally firehose stable. For strong fields, we expect the downstream to never become unstable, the condition for which is (Gary & Karimabadi 2009;Schlickeiser 2010),

with

Fig. 4(a) shows the anisotropy parameter A(x) = T y /T x with the x-axis re-centred on the shock front for various values of σ (colour coded). The anisotropy far upstream for all simulations, namely x 700c/ω pe , is 1 but it decreases in the shock foot. 1 In contrast with what was assumed in BN18, the structure of the shock is not a discontinuity, which will be discussed further in Section 5. The anisotropy parameter downstream behaves as expected; at low σ , the field cannot sustain a stable anisotropy and the downstream is nearly isotropic and A 2 ∼ 1. As the field strength increases and σ > σ c , A 2 progressively goes to 0. This is further demonstrated in Fig. 4(c), in which the red line (corresponding to the right y-axis) shows the average downstream value of A 2 as a function of σ . In Fig. 4(c), the red dashed line shows the prediction for A 2 as a function of σ (equation 8), and the two are found to be in good agreement.

Another ingredient of BN18 was that for weak fields, the plasma turns firehose unstable when crossing the shock front before it migrates to marginal stability. In this respect, Fig. 4(b) shows both sides of equation ( 10) superimposed on the shock density profile for σ = 0.056 (σ = 0.1). While the downstream shows marginal stability with A 2 ∼ 1 -1/β ࢱ2 , we find A 2 < 1 -1/β ࢱ2 for 100c/ω pe < x < 250c/ω pe . As expected, the plasma turns firehose unstable as it passes to the downstream, before it settles to marginal stability.

Finally, we can determine the downstream firehose condition for the different σ simulations; Fig. 4(c) shows the averaged downstream values from the two sides of equation ( 10) (with A 2 denoted by the red line with squares measured by the right axis and the black line with triangles measured by the left axis). For small values of σ Both terms are close to 1. As σ increases to σ 0.1, the LHS of equation ( 10) drops and even becomes negative implying that the downstream is firehose stable, while simultaneously the temperature ratio (A 2 ), also decreases from 1 to 0 showing that a stable field aligned temperature anisotropy has developed.

Structure of the downstream temperature anisotropy

The predictions laid out in BN18 ultimately derived from the temperature anisotropy generated downstream of the shock. The anisotropy is predicted to grow for strong magnetic fields, as the inflowing plasma is expected to remain magnetized, such that the perpendicular temperature should not change appreciably across the shock. This assumption is examined in Fig. 5, which shows the parallel (T xx , panel a) and out-of-plane (T zz , panel b) temperature profiles for the different shock simulations. Note that the temperatures in this figure are normalized to twice the kinetic energy of the inflowing plasma, rather than to the rest mass energy, and the x-axis is centred around the shock. For the simulations with smaller values of σ (࣠ 0.3, blue to green) both the parallel and perpendicular temperatures increase to roughly equal values downstream, resulting in an isotropic downstream temperature, consistent with Fig. 4(a). For increasing magnetization (σ 0.3, yellow to red), the parallel temperature increases dramatically, while the perpendicular temperature becomes nearly constant across the shock, verifying the assumption of BN18.

Further details about the downstream anisotropy and be found from Fig. 5(c), which analyses how the thermal energy is split among the degrees of freedom. In the plot, the average downstream temperatures are shown as a function of σ , with the parallel temperature represented by the pink line with squares, perpendicular with a green line with triangles (y)/purple line with circles (z) and the average with a grey line with 'x's (denoted by Tr[T] = T i ). 2 For σ < 0.1, 2 The choice of normalization makes it clear that the total downstream temperature is equal to the second moment of two counter streaming beam populations moving at ±v sh . This fact is likely a consequence of the simulation design, as the downstream population must be at rest in simulation frame, and so all of the energy in this frame must be thermal. It seems natural that the downstream temperature must be equal to the composite energy of the counter streaming cold beams. the different components of the temperature are equal. As the only stable configuration has A 2 ∼ 1, the downstream is isotropic with its thermal energy evenly distributed along the 3 degrees of freedom. As a consequence, each one carries 1/3 of the total, and all together they sum up to mv 2 sh which is the total amount of energy initially available in the 2 cold counter streaming plasmas. For σ > 0.1, the downstream becomes increasingly anisotropic with increasing σ according to the mechanism described in BN18. The parallel temperature T x rises while T y and T z decrease in such a way that Tr[T] is kept constant and equal to mv 2 sh (momentum conservation).

S H O C K F RO N T W I D T H

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

The value of λ 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 σ , two distinct patterns appear to emerge when considering simulations with σ 0.1 or σ 0.1. In the small σ case, the shock width seems to reach a constant size for decreasing σ . 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 σ .

For σ 0.1, the shock width increases with σ . This increase is potentially due to the increasing number of pitch angle scattering required to return particles travelling up-streaming back towards the downstream. Each deflection is expected to alter a particle's pitch angle by mu ∼ δB/B 0 ∼ 1/ √ σ . For a larger value of σ , 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 σ = 2 simulation (orange line in this paper's figures) in the supplementary material.

C O N C L U S I O N

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 σ ≡ B 2 0 /4πn 1 mv 2 1 . The simulations are found to be in good agreement with the theoretical predictions, namely that as σ is increased to 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 (β ∼ 1). For decreasing and increasing values of σ the shock width broadens, but saturating for small values of σ 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.