1. Introduction
Shock waves are fundamental phenomena in fluids. In an electrically neutral medium, only binary collisions can mediate the transition between the upstream and the downstream. Consequently, the width of the shock front is of the order of a few mean free paths (Zel’dovich & Raizer Reference Zel’dovich and Raizer2002). In a plasma where the mean free path is greater than the dimensions of the system, shock waves can also propagate, mediated by collective electromagnetic effects (Sagdeev Reference Sagdeev1966). Such shock waves are called ‘collisionless shocks’.
When it comes to analysing such shocks, the jump equations of magnetohydrodynamics (MHD) are frequently used, even though this formalism assumes a small mean path since it ultimately relies on fluid mechanics. It is therefore important to know to what extent the use of MHD is legitimate, or not, when it comes to analysing collisionless shocks.
We previously (Bret & Narayan Reference Bret and Narayan2018) developed a model capable of predicting the density jump of a collisionless shock in the presence of a parallel magnetic field. The interest of this set-up is that according to MHD, the magnetic field should play no role in the density jump of a parallel shock (Kulsrud Reference Kulsrud2005). However, our model showed that for a strong enough field, the density jump of a strong sonic shock can go from 4 to 2. This theory has been validated by Particle-in-Cell (PIC) simulations by Haggerty, Bret & Caprioli (Reference Haggerty, Bret and Caprioli2022).
The model of Bret & Narayan (Reference Bret and Narayan2018) was non-relativistic. However, even in the relativistic case, a sharp reduction in the density jump, compared with MHD predictions, was also observed by Bret et al. (Reference Bret, Pe’er, Sironi, Sadowski and Narayan2017).
The aim of the present work is to develop the relativistic version of the theory described by Bret & Narayan (Reference Bret and Narayan2018), and compare the result with the PIC simulations of Bret et al. (Reference Bret, Pe’er, Sironi, Sadowski and Narayan2017).
In the non-relativistic regime, in line with Buckingham’s
$\pi$
theorem (Buckingham Reference Buckingham1914), the number of parameters and basic units allows to compute the density jump in terms of two dimensionless parameters only: the upstream sonic Mach number
$M_s$
and the magnetic parameters
$\sigma$
, defined later by (3.5).Footnote
1
In the present relativistic regime, the speed of light is added to the list of parameters. Conversely, the problem must be solved in terms of 2 + 1 dimensionless parameters:
$M_s$
,
$\sigma$
and
$\gamma _1$
, the upstream Lorentz factor. In the following, we shall only consider the strong sonic shock case, namely
$M_s=\infty$
, having only
$\sigma$
and
$\gamma _1$
left as free parameters.
Since the non-relativistic limit has been dealt with by Bret & Narayan (Reference Bret and Narayan2018), we shall here only explore the regime
$\gamma _1 \gg 1$
.
Sections 2 and 3 summarise the physics of the theory and the relativistic conservation equations needed. Then, the problem is solved in §§ 4–6, in the front frame. Yet, because the PIC simulations were performed in the downstream frame, our solution must be boosted to this frame. This is achieved in § 7.
Figure 1. Densities
$\rho _i$
and pressures
$P_i$
are measured in the fluids (upstream and downstream) rest frame. Lorenz factors
$\gamma _i$
are measured in the front frame. The upstream is assumed isotropic, but not the downstream. Only the case of a strong shock, namely
$P_1=0$
(sonic Mach number
$M_s=\infty$
) and
$\gamma _1 \gg 1$
, is studied in this work.
2. Method
For the present work to be self-contained, the method described by Bret & Narayan (Reference Bret and Narayan2018) is recalled.
The basic idea is that as the upstream plasma goes through the front, it undergoes an anisotropic compression: it is mainly compressed in the direction parallel to the motion. The resulting downstream is therefore anisotropic, with an increased temperature perpendicular to the motion and the field, while the perpendicular temperature has been conserved.
This stage of the downstream is called ‘Stage 1’. In a collisional fluid, binary collisions would quickly restore isotropy on a time scale given by the collision frequency. In a non-magnetised collisionless medium, the Weibel instability would also restore isotropy, on a time scale given by the growth rate (Weibel Reference Weibel1959; Silva, Afeyan & Silva Reference Silva, Afeyan and Silva2021). However, in the presence of an ambient magnetic field, both theory and solar wind observations show that an anisotropy can definitely be stable (Gary Reference Gary1993; Bale et al. Reference Bale, Kasper, Howes, Quataert, Salem and Sundkvist2009; Maruca, Kasper & Bale Reference Maruca, Kasper and Bale2011). This is the very reason why the behaviour of a collisionless shock can differ from its fluid counterpart. Note that shock accelerated particles, neglected here, can also trigger a departure from MHD (Bret Reference Bret2020; Haggerty & Caprioli Reference Haggerty and Caprioli2020).
If the field is strong enough to stabilise it, then Stage 1 is the end state of the downstream. In contrast, the plasma migrates to the instability threshold,Footnote 2 namely, by definition ‘Stage 2’.
Whether we characterise Stage 1 or 2, the density jump is computed from the anisotropic MHD equations (Hudson Reference Hudson1970). Yet, within these equations, the anisotropy degree of the downstream is left as a free parameter. The process previously described specifically aims at determining this parameter. More precisely, there are three MHD jump equations for the parallel case. However, the downstream state involves four unknowns.Footnote
3
Therefore, an additional equation is needed to solve the problem. For Stage 1, this missing equation is simply
$T_{\perp 2}=T_{\perp 1}$
. For Stage 2, the marginal stability requirement provides the missing equation.
The formalism was developed for a pair plasma. It enables consideration of only one single parallel temperature and one single perpendicular temperature instead of four different temperature parameters for ions and electrons, as species of different mass can be heated differently at the front crossing (Guo et al. Reference Guo, Sironi and Narayan2017, Reference Guo, Sironi and Narayan2018). Yet, our strategy eventually relies mainly on the macroscopic MHD formalism. This enables the thought that it could be extended unchanged to electron/ion plasmas, only re-scaling some dimensionless parameters. PIC simulations are currently performed to confirm this point (Shalaby, Bret & Fraschetti Reference Shalaby, Bret and Fraschetti2025). Yet, in the present work, we shall stick to pair plasmas.
3. Conservation equations
For a parallel shock,Footnote 4 the relativistic conservation equations in the front frame have been derived by Double et al. (Reference Double, Baring, Jones and Ellison2004) and Gerbig & Schlickeiser (Reference Gerbig and Schlickeiser2011). They read
\begin{align} \left [\rho \gamma \beta \right ] & = 0,\nonumber\\ \left [\gamma ^2\beta (e+P_\parallel )\right ] & = 0, \nonumber \\ \left [\gamma ^2\beta ^2(e+P_\parallel ) +P_\parallel \right ] & = 0, \nonumber \\ \left [ B \right ] & = 0, \end{align}
where
$e, \rho , P$
are the internal energy, the mass density, and the pressure of the upstream and downstream fluids in their rest frame. Here,
$\gamma$
is the Lorentz factor of the fluid in the front frame and
$\beta =\sqrt {1+\gamma ^{-2}}$
. The last equation shows the field is conserved at the interface, as is the case for a non-relativistic system (Majorana & Anile Reference Majorana and Anile1987; Kulsrud Reference Kulsrud2005). We shall therefore drop it when writing these equations in the following. In addition, these equations imply that for such an orientation of the field, the MHD problem is eventually equivalent to the fluid problem.
The internal energy
$e$
for an anisotropic relativistic gas reads (Double et al. Reference Double, Baring, Jones and Ellison2004; Gerbig & Schlickeiser Reference Gerbig and Schlickeiser2011)
\begin{equation} e = \frac {2P_\perp + P_\parallel }{3(\varGamma - 1)} + \rho c^2. \end{equation}
Strictly speaking, the adiabatic index
$\varGamma$
is a function of the upstream quantities since depending on them, the downstream can be relativistic in its rest frame, or not.Footnote
5
A non-relativistic downstream has
$\varGamma = 5/3$
, and
$\varGamma = 4/3$
if it is relativistic. Because we here consider
$\gamma _1 \gg 1$
, we shall consider
$\varGamma =4/3$
in the following.
The first three equations contain four unknowns:
$\rho _2$
,
$\gamma _2$
,
$P_{\parallel 2}$
and
$P_{\perp 2}$
. The Stage 1 and 2 constraints on
$P_{\perp ,\parallel 2}$
then provide the fourth equation allowing to solve the problem.
The problem will be dealt with in terms of the following dimensionless variables:
\begin{align} r &= \frac {\rho _2}{\rho _1}, \end{align}
\begin{align} \theta _{\perp ,\parallel } &= \frac {P_{\perp ,\parallel }}{\rho _1c^2}, \end{align}
\begin{align} \sigma &= \frac {B_0^2/4\pi }{\gamma _1\rho _1c^2}, \end{align}
\begin{align} A_2 &= \frac {P_{\perp 2}}{P_{\parallel 2}}, \end{align}
which stand for the density ratio, the dimensionless pressures, the magnetic field parameterFootnote 6 and the downstream anisotropy, respectively.
Accounting for
$P_1=0$
, (3.1) reads
\begin{align} \rho _1 \gamma _1\beta _1 &= \rho _2 \gamma _2\beta _2, \end{align}
\begin{align} \gamma _1^2\beta _1 &= \gamma _2^2\beta _2\left (\frac {2\theta _{\perp 2} + \theta _{\parallel 2}}{3(\varGamma - 1)} + \frac {\rho _2}{\rho _1} +\theta _{\parallel 2}\right )\!, \end{align}
\begin{align} \gamma _1^2\beta _1^2 &= \gamma _2^2\beta _2^2\left (\frac {2\theta _{\perp 2} + \theta _{\parallel 2}}{3(\varGamma - 1)} + \frac {\rho _2}{\rho _1} +\theta _{\parallel 2}\right ) +\theta _{\parallel 2}. \end{align}
Using the first one to derive
\begin{equation} r = \frac {\rho _2}{\rho _1} = \frac {\gamma _1\beta _1}{\gamma _2\beta _2}, \end{equation}
and using it to eliminate
$\rho _2$
in the second and the third, we obtain the system we shall solve for both Stage 1 and Stage 2,
\begin{eqnarray} \gamma _1^2\beta _1 &\,=\,& \gamma _2^2\beta _2\left (\frac {2\theta _{\perp 2} + \theta _{\parallel 2}}{3(\varGamma - 1)} + \frac {\gamma _1\beta _1}{\gamma _2\beta _2} +\theta _{\parallel 2}\right )\!, \end{eqnarray}
\begin{eqnarray} \gamma _1^2\beta _1^2 &\,=\,& \gamma _2^2\beta _2^2\left (\frac {2\theta _{\perp 2} + \theta _{\parallel 2}}{3(\varGamma - 1)} + \frac {\gamma _1\beta _1}{\gamma _2\beta _2} +\theta _{\parallel 2}\right ) +\theta _{\parallel 2}. \end{eqnarray}
4. Stage 1
Here, we assume
$T_{\perp }$
is conserved, so that
$T_{\perp 2}=T_{\perp 1}=0$
. We therefore set
$\theta _{\perp 2}=0$
in (3.11) and (3.12). Eliminating
$\theta _{\parallel 2}$
between (3.11) and (3.12) gives the equation for
$\beta _2$
,
\begin{equation} (3 \varGamma -2) \gamma _1 \gamma _2^2 \beta _2 (\beta _2-\beta _1)+3 (\varGamma -1) (\gamma _1- \gamma _2)=0. \end{equation}
Considering in addition
$\beta _1 \sim 1$
and
$\varGamma =4/3$
, it simplifies to
\begin{equation} \sqrt {1-\beta _2^2}-(\beta _2-1)^2 \gamma _1 = 0. \end{equation}
Figure 2. Density ratio for Stage 1, normalised to
$\gamma _1$
.
Equation (3.3) then gives the density ratio for Stage 1,
\begin{equation} rS1 \equiv \frac {\rho _2}{\rho _1 } = \frac {\gamma _1\beta _1}{\gamma _2\beta _2} \sim \frac {\gamma _1}{\gamma _2\beta _2}, \end{equation}
We plot in figure 2 the quantity
$rS1/\gamma _1$
for a direct comparison with the isotropic case (see Appendix A), namely,
\begin{equation} \left .\frac {\rho _2}{\rho _1 }\right |_{\mathrm{Iso}} = 2^{3/2}\gamma _1 \sim 2.82\gamma _1. \end{equation}
With
$\gamma _1 \gg 1$
, we find the expansion
\begin{equation} rS1 = (2\gamma _1)^{2/3} + \mathcal{O}(1). \end{equation}
Figure 3. Stability diagram of Stage 1. Since
$T_{\perp 2} = 0$
, it lies on the red line and is firehose stable within the shaded area.
4.1. Stability of Stage 1
Since
$T_{\parallel 2} \neq 0$
while
$T_{\perp 2} = 0$
, the downstream anisotropy parameter in Stage 1 is
\begin{equation} A_2 = \frac {P_{\perp 2}}{P_{\parallel 2} } = 0. \end{equation}
As is the case for the non-relativistic system, Stage 1 can therefore be firehose unstable. The threshold of this instability is the same for the relativistic case than for the non-relativistic one (Noerdlinger & Yui Reference Noerdlinger and Yui1968; Barnes & Scargle Reference Barnes and Scargle1973), that is,
\begin{equation} A_2 = 1 - \frac {1}{\beta _{p\parallel 2}}, \end{equation}
whereFootnote 7
\begin{equation} \beta _{p\parallel 2} = \frac {P_{\parallel 2}}{B_0^2/4\pi }. \end{equation}
The stability diagram so defined is pictured in figure 3. With
$T_{\perp 2} = 0$
, Stage 1 lies on the red line. It is firehose stable for
$\beta _{p\parallel 2} \leq 1$
, that is,
\begin{equation} P_{\parallel 2} \leq \frac {B_0^2}{4\pi }. \end{equation}
Using the dimensionless field parameter (3.5), criteria (4.7) reads
\begin{equation} \frac {\theta _{\perp 2}}{\theta _{\parallel 2}} = 1 - \gamma _1\frac {\sigma }{\theta _{\parallel 2}}. \end{equation}
Figure 4. Critical value
$\sigma _c$
of
$\sigma$
defined by (4.13), above which the magnetic field stabilises Stage 1.
Solving the right-hand side = 0, that is,
\begin{equation} \gamma _1\sigma _c \equiv \theta _{\parallel 2}, \end{equation}
allows to define
$\sigma _c$
, the critical
$\sigma$
below which Stage 1 is unstable. To this extent, the value of
$\theta _{\parallel 2}$
is extracted from (3.11). With
$\theta _{\perp 2} =0$
, and accounting for
$\varGamma =4/3$
and
$\beta _1 \sim 1$
, it reads
\begin{equation} \theta _{\parallel 2} = \frac {\gamma _1}{\gamma _2}\frac {\gamma _1- \gamma _2}{2 \beta _2 \gamma _2}. \end{equation}
Substituting in (4.11) then gives
\begin{equation} \sigma _c = \frac {\gamma _1(1-\beta _2^2) - \sqrt {1-\beta _2^2}}{2 \beta _2} = \gamma _1(1-\beta _2), \end{equation}
where
$\beta _2$
is function of
$\gamma _1$
via (4.2). Figure 4 shows
$\sigma _c$
plotted in terms of
$\gamma _1$
.
For
$\gamma _1 \gg 1$
, we find the expansion
\begin{equation} \sigma _c = (2\gamma _1)^{1/3} + \mathcal{O}(\gamma _1^{-1/3}). \end{equation}
As the upstream kinetic energy grows, it takes an increasing field to stabilise Stage 1.
5. Stage 2
We need now assess the end state of the downstream in case Stage 1 is unstable, namely, Stage 2. To do so, we come back to (3.11) and (3.12), and add the marginal firehose stability requirement (4.10). The resolution follows the following lines.
-
(i) Use (4.10) to express
$\theta _{\parallel 2}$
in terms of
$\theta _{\perp 2}$
and
$\sigma$
, namely,(5.1)
\begin{equation} \theta _{\parallel 2} = \theta _{\perp 2} + \gamma _1 \sigma . \end{equation}
-
(ii) Replace the resulting expression of
$\theta _{\parallel 2}$
in (3.11) and (3.12), which gives (5.2)
\begin{eqnarray} \beta _2 \gamma _2^2 \left (\frac {\gamma _1}{\beta _2 \gamma _2}+2 (\gamma _1 \sigma +\theta _{\perp 2})+2 \theta _{\perp 2}\right )-\gamma _1^2 &\,=\,& 0, \end{eqnarray}
(5.3)
\begin{eqnarray} \beta _2^2 \gamma _2^2 \left (\frac {\gamma _1}{\beta _2 \gamma _2}+2 (\gamma _1 \sigma +\theta _{\perp 2})+2 \theta _{\perp 2}\right )-\gamma _1^2+\gamma _1 \sigma +\theta _{\perp 2} &\,=\,& 0. \end{eqnarray}
-
(iii) Then, eliminate
$\theta _{\perp 2}$
between these two equations.
For
$\gamma =4/3$
and
$\gamma _1 \gg 1$
, we eventually obtain an equation for
$\beta _2$
only,
\begin{equation} \sqrt {1-\beta _2^2} + \gamma _1(4\beta _2 - 3\beta _2^2 - 1) = 2 \beta _2\sigma . \end{equation}
Figure 5. Solutions of (3.11), (3.12) and (4.10) for (a)
$\beta _2$
, (b)
$A_2$
and (c)
$r/\gamma _1$
in Stage 2. The blue branch is the physical one as its
$r/\gamma _1$
merges with the fluid result for
$\sigma =0$
(i.e.
$r/\gamma _1=2^{3/2}\sim 2.82$
). The lower branch in panel (c), the orange one, starts from
$r\lt 1$
.
The downstream anisotropy is then computed as
\begin{equation} A_2=\gamma _1 \frac {-\beta _2^3+\beta _2+\beta _2^2 \Omega (\gamma _1+\sigma )+\Omega (\sigma -\gamma _1)} {\Omega \sigma \left (\beta _2^2 (\gamma _1-3)+\gamma _1-1\right )-\left (\beta _2^2-1\right ) \gamma _1 (\beta _2-\Omega \gamma _1)}, \end{equation}
with
$\Omega = \sqrt {1-\beta _2^2}$
. It is easily checked that
$A_2=1$
for
$\sigma =0$
.
Figure 5 displays the
$\beta _2$
solutions of (5.4),
$A_2$
given by (5.5) and
$r/\gamma _1$
given by (3.10), in terms of
$(\sigma ,\gamma _1)$
. Because (5.4) gives two branches for
$\beta _2$
,
$A_2$
and
$r/\gamma _1$
display two branches as well. Figure 5(c) shows that only the upper branch, the blue one, merges with the fluid result for
$\sigma =0$
(i.e.
$r/\gamma _1=2^{3/2}\sim 2.82$
). In contrast, the lower branch, the orange one, starts from
$r\lt 1$
, a non-physical result. Therefore, the branch we need to consider is the blue one in panel (c), that is, the lower one in the
$\beta _2$
plot.
Figure 6 shows a cut of figure 5(b) for
$\gamma _1=10$
. While it obviously goes to unity for
$\sigma =0$
, it never goes to 0 on any branch. Indeed, it hardly goes down to 0.8 with the blue branch, namely, the physical one. In other words, our model does not allow for Stage 2 solutions down to
$A_2=0$
, where it would merge with Stage 1.
Noteworthily, the non-relativistic result of Bret & Narayan (Reference Bret and Narayan2018) displayed a continuous transition between Stages 1 and 2, but only for a strong sonic shock.Footnote 8
Figure 6. Cut of figure 5(b) for
$\gamma _1=10$
.
Figure 7. Function
$r(\sigma )/\gamma _1$
defined by (6.1) for three values of
$\gamma _1$
. The horizonal lines pertain to Stage 1, where the density ratio does not depend on
$\sigma$
. The dashed lines picture the Stage 1 or 2 solutions which are not relevant because Stage 1 is stable, or not.
6. Putting Stages 1 and 2 together
We can now put Stages 1 and 2 together, to offer a complete solution, in the front frame, for the density jump and any
$(\gamma _1 \gg 1, \sigma )$
.
Because in our scenario, the plasma goes first to Stage 1, Stage 1 is the end state of the downstream as long as it is stable. The function
$r(\sigma )$
is therefore defined by Stage 1 for
$\sigma \gt \sigma _c$
. If
$\sigma \lt \sigma _c$
, then the downstream settles in Stage 2, which therefore defines the density ratio. Hence,
\begin{equation} r(\sigma ) = \left \{ \begin{array}{cc} \mathrm{Stage\,2}, &\qquad \sigma \lt \sigma _c, \,\,\mathrm{Figure}\, 5, \\ \mathrm{Stage \,1}, &\qquad \sigma \gt \sigma _c,\,\,\mathrm{Figure}\, 2. \end{array} \right . \end{equation}
Figure 7 displays the result. We plotted the function
$r(\sigma )$
defined previously for three values of
$\gamma _1$
. The horizonal lines pertain to Stage 1, where the density ratio does not depend on
$\sigma$
. The dashed lines picture the Stage 1 or 2 solutions which are not relevant because Stage 1 is stable, or not. For
$\gamma _1=15$
and 30, Stage 2 picks up (in terms of
$\sigma$
), when Stage 1 becomes unstable. For
$\gamma _1=10$
, Stage 2 does not offer solutions up to
$\sigma _c$
, which is why there is a
$\sigma$
-gap without solution for
$\sigma \in [2,2.6]$
. Such has also been found to be the case for some oblique orientations of the field in the non-relativistic regime (Bret & Narayan Reference Bret and Narayan2022; Bret, Haggerty & Narayan Reference Bret, Haggerty and Narayan2024).
7. Switching to the downstream frame and comparison with 2-D3V PICs
The PIC simulations of Bret et al. (Reference Bret, Pe’er, Sironi, Sadowski and Narayan2017) were performed using the code TRISTAN-MP (Spitkovsky Reference Spitkovsky2005). It is a parallel version of the code TRISTAN (Buneman Reference Buneman, Matsumoto and Omura1993) optimised to study collisionless shocks. The spatial computational domain was two-dimensional (2-D), while the velocities and fields were rendered in three dimensions. The simulations implemented the ‘mirror’ method, sending a cold (
$k_BT \ll mc^2$
) relativistic pair plasma towards a reflecting wall. Since the shock is formed from the reflected part of the plasma interacting with the incoming part, the simulation was eventually performed in the shock downstream.
The simulation box was rectangular in the
$x,y$
plane, with periodic boundary conditions in the
$y$
direction and
$x$
direction of the flow. Each simulation cell was initialised with 16 electrons and 16 positrons. The time step of the simulations was
$\Delta t = 0.045 \omega _p^{-1}$
, with
$\omega _p^2=4\pi n_0 q^2/\gamma _{1,df} m$
, where
$n_0$
is the incoming plasma density in its own frame, and
$m,q$
the electron mass and charge, respectively. The plasma skin depth
$c/\omega _p$
was resolved with 10 simulation cells. The computational domain was
$\sim 102 c/\omega _p$
in the
$y$
direction and
$3600 c/\omega _p$
in the
$x$
direction, that is, approximately
$36\,000$
cells long. Simulations typically extended up to
$t=3600 \omega _p^{-1}$
.
In (3.1), all fluid quantities are defined in their own rest frame. In turn, the Lorentz and
$\beta$
factors are measured in the front frame. Yet, the 2-D3V PIC simulations of Bret et al. (Reference Bret, Pe’er, Sironi, Sadowski and Narayan2017) were performed in the downstream frame. Therefore, we need to boost figure 7 to the downstream frame for a comparison with simulations.
A key quantity in this respect is
$\gamma _{1,df}$
, the Lorentz factor of the upstream seen from the downstream. With
\begin{equation} \beta _{1,df} = \frac {\beta _1-\beta _2}{1-\beta _1\beta _2}, \end{equation}
it reads
\begin{equation} \gamma _{1,df} = \frac {1}{\sqrt {1-\beta _{1,df}^2}}= \gamma _1\gamma _2(1-\beta _1\beta _2), \end{equation}
where the subscript
$df$
stands for ‘downstream frame’.
Boosting figure 7 to the downstream frame implies boosting the density ratio, the
$\sigma$
parameter and the
$\gamma _1$
parameter.
-
(i) Regarding the density ratio, we have
(7.3)Now,
\begin{equation} r_{df} = \frac {\rho _{2, df}}{\rho _{1, df}}. \end{equation}
$\rho _{2, df}$
is simply
$\rho _2$
since it is the downstream density in its own rest frame. Regarding
$\rho _{1, df}$
, it reads
$\rho _{2, df}= \gamma _{1,df}\rho _1$
. Therefore, and still for
$\gamma _1 \gg 1$
,(7.4)Hence,
\begin{equation} r_{df} = \frac {\rho _{2, df}}{\rho _{1, df}} = \frac {\rho _2}{\gamma _{1,df}\rho _1} = \frac {1}{\gamma _{1,df}}\frac {\gamma _1\beta _1}{\gamma _2\beta _2} \sim 1+\frac {1}{\beta _2}. \end{equation}
(7.5)where
\begin{equation} r_{df} = \left \{ \begin{array}{cc} 1 + 1/\beta _{2,S2}, &\qquad \sigma \lt \sigma _c, \\ 1 + 1/\beta _{2,S1}, &\qquad \sigma \gt \sigma _c, \end{array} \right . \end{equation}
$\beta _{2, S1,S2}$
is
$\beta _2$
in Stage 1 and 2, given by (4.2) and (5.4), respectively.
-
(ii) The
$\sigma$
parameter also needs a boost. In (3.5), the field remains unchanged since it is parallel to the boost. Only
$\rho _1$
and
$\gamma _1$
are modified. With
$\rho _{1,df}=\gamma _{1,df}\rho _1$
, we find(7.6)In this
\begin{equation} \sigma _{df} = \frac {B_0^2/4\pi }{\gamma _{1,df}\,\rho _{1,df}c^2} = \frac {\sigma }{\gamma _1\gamma _2^2(1-\beta _1\beta _2)^2} \sim \frac {\sigma }{\gamma _1\gamma _2^2(1-\beta _2)^2} = \frac {1+\beta _2}{1-\beta _2}\frac {\sigma }{\gamma _1}. \end{equation}
$\sigma$
re-scaling,
$\beta _2$
is given by its Stage 1 expression when considering Stage 1, and by its Stage 2 expression when considering Stage 2.
With respect to the critical
$\sigma _c$
of
$\sigma$
below which Stage 1 is unstable, it reads, according to (4.13),
$\sigma _c = \gamma _1(1-\beta _2)$
. Using (7.6), we obtain its value in the downstream frame,(7.7)where
\begin{equation} \sigma _{c, df} = 1 + \beta _2, \end{equation}
$\beta _2$
must be given by its Stage 1 expression since it defines the stability threshold for Stage 1. Therefore, the transition from Stage 1 to Stage 2, seen from the downstream frame, always occurs for
$\sigma _{c, df} \in [1,2]$
.
-
(iii) Finally, the PIC simulations presented by Bret et al. (Reference Bret, Pe’er, Sironi, Sadowski and Narayan2017) considered two values for the upstream Lorentz factors, 10 and 30, defined in the downstream frame. In other words, simulations varying
$\sigma _{df}$
were performed at constant
$\gamma _{1,df}$
. It means that to match a PIC ran at
$\gamma _{1,\mathrm{PIC}}$
, we need to run the abovementioned calculations for
$\sigma$
and
$\gamma _1$
, fulfilling(7.8)where we used
\begin{equation} \gamma _{1,df} = \gamma _1\gamma _2(1-\beta _2) = \gamma _{1,\mathrm{PIC}}, \end{equation}
$\beta _1 \sim 1$
. Figure 8 shows the value of
$\gamma _{1,df}$
in terms of
$(\sigma , \gamma _1)$
, for Stages 1 and 2. Superimposed are the contours of constant
$\gamma _{1,df}=10$
and
$30$
. They are flat for Stage 1, where
$\gamma _{1,df}$
is independent of
$\sigma$
, and slightly bent upwards for Stage 2. As a result, we find the following correspondence.
In Stage 1,
(7.9)
\begin{align} \gamma _{1,df} = 10 &\Rightarrow \gamma _1 = 43, \nonumber \\[3pt] \gamma _{1,df} = 30 &\Rightarrow \gamma _1 = 230. \end{align}
In Stage 2,
(7.10)It has been checked that for Stage 2, considering a constant
\begin{align} \gamma _{1,df} = 10 &\Rightarrow \gamma _1 \sim 16, \nonumber \\[3pt]\gamma _{1,df} = 30 &\Rightarrow \gamma _1 \sim 43. \end{align}
$\gamma _1$
hardly affects the result.
With (7.2) reading
$\gamma _{1,df} = \gamma _1\gamma _2(1-\beta _1\beta _2)$
, one could expect in (7.9) for Stage 1, a value of
$\gamma _1$
lower than 230. Indeed, if, as hinted in the Appendix,
$\beta _2 \sim 1/3$
and
$\gamma _2 \sim 1$
, then with
$\gamma _1=230$
,
$\gamma _{1,df} \sim 230 \times 1 \times 2/3 = 153 \gt 30$
. However, in the present case, the value of
$\beta _2$
has to be that derived in §4, for the anisotropic downstream in Stage 1, not the isotropic case treated in the Appendix. Additionally, for
$\gamma _1 = 230$
, Stage 1 has
$\beta _2 \sim 0.967$
, which explains how the product
$\gamma _1\gamma _2(1-\beta _1\beta _2)$
is reduced to just 30.
Figure 8. Value of
$\gamma _{1,df}$
in terms of
$(\sigma , \gamma _1)$
, for Stages 1 and 2, with the contours of constant
$\gamma _{1,df}=10$
and 30.
Figure 9. Density ratio measured in the downstream frame,
$r_{df}$
, in terms of
$\sigma _{df}$
. The transition between the two stages occurs for the critical
$\sigma _{c, df}$
given by (7.7). The squares show the results of the PIC simulations performed by Bret et al. (Reference Bret, Pe’er, Sironi, Sadowski and Narayan2017).
Figure 9 is the end result of this article. It displays the density ratio measured in the downstream frame,
$r_{df}$
, in terms of
$\sigma _{df}$
. The transition between the two stages occurs for the critical
$\sigma _{c, df}$
given by (7.7). The squares shows the results of the PIC simulations performed by Bret et al. (Reference Bret, Pe’er, Sironi, Sadowski and Narayan2017).Footnote
9
Note that Stage 1 always defines the density ratio for large fields
$\sigma _{df} \gt \sigma _{c, df}$
. It is therefore instructive to derive an expansion of (7.4) in Stage 1. We find
\begin{align} r_{df, S1} & = 1 + \frac {1}{\beta _{2,S1}} = 2 + \left ( \frac {\sqrt {2}}{\gamma _1} \right )^{2/3} + \mathcal{O}(\gamma _1^{-4/3}) \nonumber \\ & = 2 + \frac {1}{\gamma _{1, \mathrm{PIC}}} + \mathcal{O}(\gamma _{1, \mathrm{PIC}}^{-2}) . \end{align}
8. Conclusion
In this work, we have implemented the relativistic version of the model presented by Bret & Narayan (Reference Bret and Narayan2018). The system considered is a collisionless shock in pair plasma, with a parallel magnetic field. Since numerical simulations have been carried out by Bret et al. (Reference Bret, Pe’er, Sironi, Sadowski and Narayan2017), we are able to compare the results of the present theory with these simulations.
After adapting to the relativistic case the theory presented by Bret & Narayan (Reference Bret and Narayan2018), and boosting the results to the downstream frame, we reach, in figure 9, the comparison of the relativistic theory with the PIC simulations.
The agreement is satisfactory. Our model reproduces the fall of the density ratio over the correct
$\sigma _{df}$
range, with an asymptotic value close to 2 in the strong field regime
$\sigma _{df} \gt \sigma _{c, df}$
.
We observe two kinds of discrepancies in figure 9, between the theory and the simulations. (i) There is a systematic underestimation of the density ratio of the theory with respect to simulations and (ii) the PIC simulations find solutions where the model does not offer any, like for
$\sigma _{df} = 1.5$
with
$\gamma _{1,df}=30$
. Let us briefly comment on each point.
-
(i) In the absence of a magnetic field, the density jump measured in the downstream is given by (see Appendix A)
(8.1)For an adiabatic index
\begin{equation} r_{df} = \frac {\varGamma }{\varGamma -1}. \end{equation}
$\varGamma = 4/3$
, this gives
$r_{df}=4$
. Yet, as seen in figure 9, PIC simulations for a weak field give a slightly higher value,
$r_{df}\sim 4.3$
. Such an overshoot could be interpreted as the PIC effectively simulating an adiabatic index
$\varGamma =1.3$
, slightly lower than the fully relativistic
$\varGamma =4/3 \sim 1.33$
considered in the theory.Footnote
10
Yet, this solution seems improbable, for if the downstream were not fully relativistic in the simulation, its adiabatic index would rather lean towards a higher value, namely
$\varGamma =5/3 \sim 1.66$
. Why, then, is the density ratio in the PIC simulations higher than (8.1) predicts? The most probable answer is related to accelerated particles. It has been known for long that collisionless shocks can accelerate particles (Blandford & Ostriker Reference Blandford and Ostriker1978). Additionally, it has been equally known for long than their effect on the shock is to raise the density ratio (Berezhko & Ellison Reference Berezhko and Ellison1999; Stockem et al. Reference Stockem, Fiúza, Fonseca and Silva2012). To date, this is the only known mechanism capable of raising the density ratio above its MHD value (see Bret (Reference Bret2020) and references therein). Therefore, strictly speaking, the PICs presented here are not a simulation of the theory, since the theory ignores particles acceleration, which simulations necessarily include. The PIC simulations eventually simulate the system discussed here, embedded in a bath of cosmic rays.
-
(ii) In addition to ignoring particle acceleration, other differences between the PIC simulations and the theory are noticeable. For example, figures 3(b), 4(b) and 6(b) of Bret et al. (Reference Bret, Pe’er, Sironi, Sadowski and Narayan2017) show that while the upstream temperature is nearly 0, the downstream is definitely anisotropic in the PIC simulations, but its perpendicular temperature is never 0, even in the strong field regime. This means that Stage 1, as described here, is not fully realised. Probably due to the influence of cosmic rays, the simulations do not stick to the present theory, even though they follow its trend. It is therefore not surprising that the simulations find downstream solutions even when our rigid, ‘cosmic rays-less’ theory does not allow any. Future work could focus on some detailed analysis of the particle distribution in the PIC. Yet, as mentioned, it could be challenging to run some PIC simulations cancelling the cosmic rays they necessarily produce.
There is a striking similarity between the non-relativistic and relativistic theories. In the non-relativistic theory, the density jump in strong field tends to
$r=2$
. This is an important departure from the MHD prediction, since the latter predicts a field-independent density jump for a parallel shock of
$r=4$
(strong sonic shock).
The same pattern is observed here. Relativistic MHD still predicts a field-independent density jump for this parallel geometry, which amounts to a density jump measured in the downstream
$r_{df} \sim 4$
. In contrast, the relativistic theory presented here, and confirmed by simulations, gives a density jump also tending towards
$r_{df}=2$
.
It is likely that, as in the non-relativistic case (Bret & Narayan Reference Bret and Narayan2022; Bret et al. Reference Bret, Haggerty and Narayan2024), the departure from MHD is maximum for this obliquity of the magnetic field. This could be confirmed in future works.
Extending this work to finite upstream pressures would require modifying the left-hand sides of (3.8) and (3.9), before implementing similar calculations following the same road map. While this could be the goal of some future work, the result would probably be a lowering of the density ratio with diminishing sonic Mach number, as was observed in the non-relativistic case (Bret & Narayan Reference Bret and Narayan2018).
Acknowledgements
Editor Luís O. Silva thanks the referees for their advice in evaluating this article.
Funding
A.B. acknowledges support by the Ministerio de Economía y Competitividad of Spain (Grant No. PID2021-125550OB-I00). R.N. was partially supported by the Black Hole Initiative at Harvard University, which is funded by the Gordon and Betty Moore Foundation grant 8273, and the John Templeton Foundation grant 62286.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Jump conditions for an isotropic fluid with
$P_1=0$
Although the simple case of a relativistic shock in an isotropic fluid has been abundantly discussed in the literature (Taub Reference Taub1948; Anile Reference Anile1990; Landau & Lifshitz Reference Landau and Lifshitz2013), we here present a simple derivation of the jump conditions, for
$P_1=0$
and
$\gamma _1 \gg 1$
, using the present notation (Narayan Reference Narayan2012).
We start from (3.7)–(3.9), where only
$P_1=0$
has been assumed. Setting
$\theta _{\perp 2}=\theta _{\parallel 2} = \theta _{2}$
to reflect isotropy gives
\begin{eqnarray} \rho _1 \gamma _1\beta _1 &\,=\,& \rho _2 \gamma _2\beta _2, \end{eqnarray}
\begin{eqnarray} \gamma _1^2\beta _1 &\,=\,& \gamma _2^2\beta _2\left (\frac {\varGamma }{\varGamma - 1}\theta _2 + \frac {\rho _2}{\rho _1} \right ), \end{eqnarray}
\begin{eqnarray} \gamma _1^2\beta _1^2 &\,=\,& \gamma _2^2\beta _2^2\left (\frac {\varGamma }{\varGamma - 1}\theta _2 + \frac {\rho _2}{\rho _1} \right ) +\theta _2. \end{eqnarray}
Computing (A.3)
$- \beta _2 \times$
(A.2) gives
\begin{equation} \gamma _1^2\beta _1(\beta _1-\beta _2)=\theta _2. \end{equation}
Using this expression for
$\theta _2$
in (A.2) gives
\begin{equation} \gamma _1 - \gamma _2 = \frac {\varGamma }{\varGamma - 1}\gamma _1\gamma _2^2\beta _2(\beta _1-\beta _2). \end{equation}
Since
$\gamma _1 \gg 1$
,
$\beta _1 \sim 1$
. In addition, because the downstream has to be subsonic while the speed of sound in a relativistic gas is
$c/\sqrt {3}$
, we know
$\beta _2 \lt 1/\sqrt {3}$
. Hence,
$\gamma _2 \sim 1$
so that
$\gamma _1 \gg \gamma _2$
. Equation (A.5) therefore simplifies to
\begin{eqnarray} \gamma _1 - \gamma _2 &\sim & \gamma _1 \nonumber \\ &\sim & \frac {\varGamma }{\varGamma - 1}\gamma _1\gamma _2^2\beta _2(1-\beta _2) \nonumber \\ & \Rightarrow & \,\, \beta _2 = \varGamma - 1, \end{eqnarray}
where
$\gamma _2^2 = (1-\beta _2^2)^{-1}$
and
$\beta _1 \sim 1$
has been used. This expression for
$\beta _2$
is identical to (16) of Kirk & Duffy (Reference Kirk and Duffy1999). Then,
\begin{equation} \gamma _2 = \frac {1}{\sqrt {1-\beta _2^2}} = \frac {1}{\sqrt {(2-\varGamma ) \varGamma }}, \end{equation}
and from (3.10),
\begin{equation} r = \frac {\rho _2}{\rho _1} \sim \frac {\gamma _1}{\beta _2\gamma _2} = \gamma _1 \frac {\sqrt {\varGamma (2-\varGamma )}}{\varGamma }. \end{equation}
Finally,
\begin{equation} \gamma _{1,df} \sim \gamma _1(1-\beta _2) = \gamma _1\sqrt {\frac {2-\varGamma }{\varGamma }} \end{equation}
and
\begin{equation} r_{df} = \frac {r}{\gamma _{1,df}} = \frac {\varGamma }{\varGamma -1}. \end{equation}
For
$\varGamma =4/3$
, (A.8) and (A.10) give
$r = 2^{3/2}\gamma _1$
and
$r_{df} = 4$
, respectively. Note that PIC simulations, where velocities are rendered in only two dimensions, give a density ratio of
$r_{df} = 3$
, corresponding to
$\varGamma =3/2$
.
Bret et al. (Reference Bret, Pe’er, Sironi, Sadowski and Narayan2017) rendered velocities in three dimensions, as explained in §7. If they had been rendered in 2-D, the density jump would have been 3. For a direct comparison of the 2-D and 3-D cases in PIC simulations, see figure 3 of Bret et al. (Reference Bret, Stockem, Narayan and Silva2014).
Open access
