# On the width of a collisionless shock and the index of the cosmic rays it
accelerates

Antoine Bret  ETSI Industriales, Universidad de Castilla-La Mancha, 13071,
Ciudad Real, Spain  Instituto de Investigaciones Energéticas y Aplicaciones
Industriales, Campus Universitario de Ciudad Real, 13071, Ciudad Real, Spain
Asaf Pe’er  Department of Physics, Bar-Ilan University, Ramat-Gan, 52900,
Israel

###### Abstract

Despite being studied for many years, the structure of collisionless shocks is
still not fully determined. Such shocks are known to be accelerators of cosmic
rays, which, in turn, modify the shock structure. The shock width $\lambda$ is
known to be connected to the cosmic rays (CRs) spectral index, $a$. Here, we
use an instability analysis to derive the shock width in the presence of CRs.
We obtain an analytical expression connecting the shock width to the CRs index
and to the fraction of upstream particles that are accelerated. We find that
when this fraction becomes larger than $\sim$ 30%, a new instability becomes
dominant. The shock undergoes a transition where its width increases by a
factor $\sim 8-10$, and the CRs acceleration effectively ends. Our analysis is
valid for strong, non-relativistic and unmagnetized shocks. We discuss the
implication of these results to the expected range of CRs spectra and flux
observed, and on the structure of non-relativistic collisionless shocks.

##  1 Introduction

Shock waves are a common phenomenon in many environments. Their study dates
back to the 19th century, when it was already found that a large-amplitude
sound wave will necessarily steepen during its propagation, giving rise to a
density discontinuity, i.e. a shock wave (Salas, 2007).

In a fluid with a small Knudsen number (ratio of the mean free path to the
length scale of the system), binary collisions between molecules constitute
the microscopic mechanism by which the density jump is accomplished. As a
consequence, the shock front is not a discontinuity, but a thin transition
layer instead, a few mean free paths thick (Zel’dovich & Raizer, 2002).
Towards the end of the 20th century, Sagdeev’s work demonstrated that in the
other extreme, namely in a plasma with a mean free path well above the
dimensions of the system, a shock wave can also form (Sagdeev, 1966). In this
case, the shock wave is not maintained by binary collisions, but rather by
collective electromagnetic effects. Such shocks have been dubbed
“collisionless shocks”.

Besides their interest as fundamental processes, shock waves have in recent
decades been the subject of much attention from the astrophysical community,
thanks to their ability to accelerate cosmic rays to high energies, e.g., via
the Fermi process (Blandford & Ostriker, 1978; Spitkovsky, 2008). As a result,
they are a key ingredient of CRs production scenarios, as well as other
phenomena that show non-thermal emission such as gamma-ray bursts or fast
radio bursts (Piran, 2004; Zhang, 2020). The spectra and flux of the CRs
produced in shock waves in various astronomical objects are the subject of
extensive research, be it to understand CRs origin or that of high energy
neutrinos resulting from CRs interactions.

Early works establishing the acceleration of particles by collisionless shocks
considered the shock front as a discontinuity (Blandford & Ostriker, 1978). In
this approximation, the CRs spectrum obeys a power law $N(p)dp\propto p^{-a}$,
the index of which depends on the density jump ratio $r$ with $a=3r/(r-1)$.
Here, $p$ is the particle’s momentum and $a$ is the index of the power-law
accelerated CRs. For a strong shock with $r=4$, it was found that $a=4$.

Shortly afterwards, a number of authors studied the change in this index when
considering a spatial extension of the front (Drury et al., 1982a; Schneider &
Kirk, 1987, 1989). It was found that the CRs spectrum is steeper when an
extended front is considered, so that a real collisionless shock is a worse
accelerator than a shock idealized by a discontinuity. An overall conclusion
of the aforementioned works is that the CRs index $a$ depends not only on the
density ratio $r$, but also on the width of the shock front $\lambda$:
$a=a(r,\lambda)$.

A natural scale for the shock width $\lambda$ is the ion inertial length,
$\lambda_{i}=u_{i}/\omega_{p,i}$, where $u_{i}$ is the ion velocity and
$\omega_{p,i}$ the ion plasma frequency (e.g., Treumann, 2009). However, the
shock width $\lambda$ itself is affected by the existence of CRs. Therefore,
one can expect $\lambda$ to be a function of both the CRs index $a$ and the
fraction of CRs in the plasma. As a result, the relation between $a$ and
$\lambda$ is not simple.

Several authors studied, mostly numerically, the CRs spectrum produced by a
finite width shock (Kruells & Achterberg, 1994; Marcowith & Kirk, 1999; Li et
al., 2013; Xu et al., 2022). In these works, the velocity field from the
upstream to the downstream is usually modelled by an hyperbolic tangent (Drury
et al., 1982a). Xu et al. (2022) proposed an analysis of the effect of the
Péclet number introduced below by Equation (2). Also, Li et al. (2013) applied
these concepts in the space plasma framework. Here, we introduce a new
ingredient to the topic: the determination of the shock width from an
instability analysis. Furthermore, our semi-analytic approach is easily
tractable.

In this respect, Tidman (1967) looked for the fastest growing unstable mode at
the middle of the shock front, deducing that the shock width must be
proportional to the inverse of the most unstable wavevector $k_{m}$. This
approach, besides providing a better approximation to the shock width than
simple scaling arguments, provides important physical insight into the
structure of collisionless shock fronts.

Although Tidman did not include CRs in his analysis, his reasoning suggests a
straightforward way to do it, namely solving the dispersion relation of the
distribution function modified by the CRs. As we shall see in the sequel, the
presence of CRs indeed modifies Tidman’s analysis, leading to a different
shock width. As we will show below, this width depends on both the CRs
spectral index $a$, and even more importantly, on the fraction of CRs, denoted
hereafter $\epsilon$.

The goal of this article is to examine how the existence of CRs affects the
shock structure, and using this input, to examine the previously derived
relation between the CRs index $a$ and the shock width $\lambda$. In fact, as
we shall see below, both the shock width and the CRs index depend on the
fraction of CRs in the plasma in a non-linear way. The results derived are
expected to be of importance to numerical experiments and CRs observations.
Indeed, in a typical numerical experiment where a plasma is sent to a
reflecting wall, the upstream parameters are the only ones determining the
properties of the subsequent shock. The width of the shock front and the index
of the CRs it accelerates, should therefore be functions of these upstream
parameters. This work is a first step toward determining them.
Observationally, we may have demonstrated the existence of an effective cutoff
in the ability of non-relativistic strong shock to accelerate CRs beyond a
certain fraction, thereby implying an upper limit on the CRs flux expected
from a given object. In addition, our results connect the CRs spectral index
to the shock properties.

Although the shock density ratio $r$ and the sonic Mach number $\mathcal{M}$
appear in some of our equations, the end results are established considering
$\mathcal{M}\gg 1$ and $r\sim 4$. Our conclusions are therefore restricted to
strong shocks. In addition, we do not account for an external magnetic field
and implement a 1D model.

The article is structured as follows. In Section 2 we recall the dependence of
the CRs index on the width of the shock, and derive a simple analytical
expression relating them. In Section 3 we recall Tidman’s analysis of the
front width in terms on instabilities and show how the result depends on the
CRs fraction and their spectral index, when they are accounted for. In
particular, we discover a new mode that becomes dominant for CRs fraction
exceeding $\sim 30\%$. We conclude in Section 4 discussing the implications of
our results on the shock structure and its observational consequences.

##  2 Relation between the CRs index and the shock width

The dependence of the CRs index on the shock width was considered by several
authors (e.g., Drury et al., 1982a; Schneider & Kirk, 1987, 1989). In the
limit of a discontinuity ($\lambda\rightarrow 0$), the relation derived in
Blandford & Ostriker (1978) for the power index $a$ reads,

$a=\frac{3r}{r-1}.$ (1)

For a strong shock with $r=4$, this gives $a=4$. Later on, the effect of a
finite front width was considered in Drury et al. (1982a); Schneider & Kirk
(1987, 1989). In particular, Schneider & Kirk (1987) considered various
velocity profiles for the transition between the upstream and the downstream.
Considering the simplest transition, namely linear, we prove in Appendix A
that for a strong shock with $r=4$, the power index $a$ reads,

$a=4+\frac{1}{6}\frac{\lambda}{D/u_{1}},$ (2)

where $u_{1}$ is the upstream velocity and $D$ the diffusion coefficient of
the CRs. The dimensionless quantify $\lambda u_{1}/D$ is known as the “Péclet
number” of the shock (Zel’dovich & Raizer, 2002).

Although derived through a Taylor expansion for $\lambda\ll D/u_{1}$, the
expression is indeed quite accurate up to $\lambda\sim 15D/u_{1}$ (see Figure
6 in Appendix A). Yet, this expression depends on the velocity profile in the
shock front. It is obtained for the simplest, linear profile, with constant
diffusivity. Other assumptions (such as constant diffusion length, $D/u_{1}$)
lead to a similar expression, though with a slightly different pre-factor in
front of the Péclet number (the second term); see Drury et al. (1982a);
Achterberg & Schure (2011); Xu et al. (2022).

##  3 Relation between the shock width and the CRs index from an instability
analysis

For the present article to be self-contained, we briefly remind the result
obtained in Tidman (1967). Tidman (1967) assessed the physics of a
collisionless shock using the so-called “Mott-Smith ansatz” (Mott-Smith, 1951;
Bret & Pe’er, 2021), which consists in writing the distribution function of
the particles along the shock profile ($z$ axis here), as a linear combination
of the upstream and downstream Maxwellians. One then writes,

$f(z,\mathbf{p})=n_{1}(z)f_{1}(\mathbf{p})+n_{2}(z)f_{2}(\mathbf{p}),$ (3)

where $n_{1,2}(z)$ are the weights of the upstream and downstream Maxwellians
$f_{1,2}(\mathbf{p})$. Note that such a distribution constitutes a counter-
streaming system. As such, it may be unstable and indeed, Tidman (1967)
analysed the shock structure in terms of these instabilities. In particular,
he concluded that the width of the front is proportional to the most unstable
_wavelength_ found for $z$ in the middle of the front, where the two
Maxwellians have about equal weights. Note that while instability analysis is
usually conducted in terms of unstable frequencies, this one considers
unstable wavelengths.

Neglecting temperature effects for simplicity, Tidman (1967) argued that the
shock width is

$\lambda=A\frac{u_{1}}{\omega_{p1}},$ (4)

where $u_{1}$ and $\omega_{p1}$ are the upstream velocity and ion plasma
frequency respectively, and $A$, the “Tidman’s constant”, a constant of order
10 at most. We shall now see how Eq. (4) is retrieved from an instability
analysis as in Tidman (1967), before modifying it accounting for CRs.

Figure 1: Imaginary (left) and real (right) parts of $Z$ in terms of $x$ for
various $r$ in the absence of CRs.

###  3.1 Instability analysis without cosmic rays (CRs)

Tidman’s result (4) can be derived from a one-dimensional model. Given a
distribution $f_{0}(p)$, the dispersion equation for small perturbations
$\propto\exp(ikz-i\omega t)$ reads (see for example Landau & Lifshitz (1981),
§29),

$\frac{4\pi q^{2}}{k}\int\frac{df_{0}/dp}{kv-\omega}dp=1,$ (5)

where $q$ is the ions’ charge and $v=p/m$. Let us consider a distribution of
the form,

$f_{0}(p)=n_{1}\delta(p-m_{i}u_{1})+n_{2}\delta(p-m_{i}u_{2}),$ (6)

where $m_{i}$ is the mass of the ions, $n_{1,2}=N_{1}/2$ with $N_{1}$ the far
upstream density, and $u_{1},u_{2}$ the upstream and downstream velocities
respectively. We also adopted another assumption of Mott-Smith (1951), namely
that the densities of the upstream and downstream incoming particles are the
same in the middle of the shock front, where the instability analysis in
conducted.

Since $\int\delta^{\prime}g=-\int\delta g^{\prime}$, inserting Equation (6)
into Equation (5) yields,

$\frac{4\pi n_{1}q^{2}}{m_{i}}\frac{1}{(\omega-ku_{1})^{2}}+\frac{4\pi
n_{2}q^{2}}{m_{i}}\frac{1}{(\omega-ku_{2})^{2}}=1.$ (7)

Denoting by $r$ the shock density ratio, $u_{2}=u_{1}/r$. With this and
$n_{1,2}=N_{1}/2$, Equation (7) becomes,

$\frac{4\pi\frac{N_{1}}{2}q^{2}}{m_{i}}\left(\frac{1}{(\omega-
ku_{1})^{2}}+\frac{1}{(\omega-ku_{1}/r)^{2}}\right)=1.$ (8)

Introducing the plasma frequency,

$\omega_{p1}^{2}=\frac{4\pi N_{1}q^{2}}{m_{i}},$ (9)

and the normalized frequency and wavevector,

$\displaystyle x$ $\displaystyle=$ $\displaystyle\frac{\omega}{\omega_{p1}},$
$\displaystyle Z$ $\displaystyle=$ $\displaystyle\frac{ku_{1}}{\omega_{p1}},$
(10)

Equation (8) reads

$\frac{1/2}{(x-Z)^{2}}+\frac{1/2}{(x-Z/r)^{2}}=1.$ (11)

The imaginary and real parts of $Z$ are plotted on Figure 1 in terms of $x$
for various shock compression ratio. For $r=4$, it peaks at ${\rm
Im}(Z_{max})\sim 1$ near $x=1$, with ${\rm Re}(Z_{max})\sim 1$. Denoting by
$k_{m}$ the real part of the most unstable wavevector $k$, we therefore find,

${\rm
Re}(Z_{max})=\frac{k_{m}u_{1}}{\omega_{p1}}=1~{}~{}\Rightarrow~{}~{}k_{m}=\frac{\omega_{p1}}{u_{1}},$
(12)

and Tidman’s prescription $\lambda=A/k_{m}$ retrieves Equation (4). We point
out that this instability analysis is conducted in terms of the most unstable
wave vector $k$, whereas such analysis usually focusses on the most unstable
frequency $\omega$.

###  3.2 Instability analysis with CRs

We are now prepared to see how the presence of CRs modifies Equation (4). We
consider the same one-dimensional model, adding a CRs component to the
distribution function,

$f_{0}(p)=n_{1}\delta(p-m_{i}u_{1})+n_{2}(1-\epsilon)\delta(p-m_{i}u_{2})+n_{cr}f_{cr}(p),$
(13)

where $n_{cr}\equiv\epsilon N_{1}$ is the CRs density, $\epsilon$ being the
fraction of upstream particles accelerated. As a consequence of this
acceleration, the number of thermal particles in the second beam is depleted
by the same amount, hence the factor $(1-\epsilon)$.

Particle-in-cell (PIC) simulations of cosmic ray acceleration typically yield
$\epsilon\sim 10^{-2}$ (Sironi et al., 2013; Marcowith et al., 2016). However,
even the longest simulations to date, like for example those carried on by
Groselj et al. (2024), reproduce but an extremely short slice of the shock
life.

Furthermore, there are clear hints that magnetic field generation and particle
acceleration are inter-related, while the fraction of energy carried by the
accelerated particles grows with time (Keshet et al., 2009; Groselj et al.,
2024), following the development of magnetic fields at increasingly larger
coherence scale. Therefore, since the number of accelerated particles grows
with time, we shall explore values of $\epsilon$ up to $\epsilon=0.5$.

For $f_{cr}(p)$, we consider

$f_{cr}(p)=\kappa p^{-a},~{}~{}~{}p\in[P_{min},+\infty],$ (14)

where $a$ is the power index and $\kappa$ the normalizing constant, given by

$\int_{P_{min}}^{\infty}\kappa
p^{-a}dp=1~{}~{}\Rightarrow~{}~{}\kappa=\frac{a-1}{P_{min}^{1-a}}.$ (15)

It is common in literature to introduce $\zeta$, the fraction of the shock
kinetic energy imparted into CRs. Given the distribution (13), the energy
contained in CRs in our scenario reads $\int\epsilon N_{1}f_{cr}(p)E(p)dp$,
with $E(p)=p^{2}/2m_{i}$. In the present 1D model and for $a=4$ or higher, the
resulting integral converges. Still, since our instability analysis only
relies on $\epsilon N_{1}$, the number of CRs, rather than on their energy, we
shall not refer to $\zeta$ in the sequel.

Inserting now Equation (13) into Equation (5) gives,

$\frac{4\pi\frac{N_{1}}{2}q^{2}}{m_{i}}\left(\frac{1}{(\omega-
ku_{1})^{2}}+\frac{1-\epsilon}{(\omega-ku_{1}/r)^{2}}\right)+J_{cr}=1,$ (16)

with

$J_{cr}=\epsilon\frac{4\pi
N_{1}q^{2}}{k}\kappa\int\frac{(-a)p^{-a-1}}{kp/m_{i}-\omega}dp.$ (17)

Here, the integration domain is $[-\infty,-P_{min}]\cup[P_{min},+\infty]$.
After some manipulations, we find

$J_{cr}=-\epsilon\frac{a(a-1)}{xZ}\frac{m_{i}u_{1}}{P_{min}}~{}\varphi\left(a,\frac{Z}{x}\frac{P_{min}}{m_{i}u_{1}}\right),$
(18)

where

$\displaystyle\varphi(a,y)$ $\displaystyle=$ $\displaystyle
2\int_{1}^{+\infty}\frac{s^{-a-1}}{s^{2}y^{2}-1}ds,$ (19) $\displaystyle=$
$\displaystyle\frac{2}{(a+2)y^{2}}~{}_{2}F_{1}\left(1,\frac{a+2}{2};\frac{a+4}{2};\frac{1}{y^{2}}\right),~{}~{}y\in\mathbb{C},~{}a>4.$

where ${}_{2}F_{1}$ is the hypergeometric function ${}_{2}F_{1}(a,b;c;z)$.

Since the thermal particles are the seed particles for acceleration to high
energies, one has $P_{min}\sim 2m_{i}v_{sh}$ (Caprioli & Spitkovsky, 2014a;
Caprioli et al., 2015), where $v_{sh}$ is the speed of the shock in the
upstream frame. Therefore, in the present notations, $v_{sh}=u_{1}$ with, in
addition, $u_{1}=\mathcal{M}_{1}c_{s1}$, where $\mathcal{M}_{1}$ is the shock
Mach number and $c_{s1}$ the upstream speed of sound. One thus finds

$\frac{P_{min}}{m_{i}u_{1}}=2.$ (20)

The CRs-modified dispersion equation eventually reads,

$\frac{1/2}{(x-Z)^{2}}+\frac{(1-\epsilon)/2}{(x-Z/r)^{2}}-\epsilon\frac{a(a-1)}{xZ}\frac{1}{2}~{}\varphi\left(a,2\frac{Z}{x}\right)=1,$
(21)

with $x\in\mathbb{R}$ and $Z\in\mathbb{C}$.

Figure 2: Imaginary part of $Z$ in terms of $x$ given by Equation (21), for
various choices of $(\epsilon,a)$. The red dashed line is the solution without
CRs ($\epsilon=0$), for $r=4$. The curve that is close to it and falls to 0
for $x\sim 1.6$ show how CRs modify this mode. The other curve that extends to
height $x$’s is the new unstable mode triggered by the CRs. The maximum is
always around $x=1$.

Solution to Equation (21) gives both the imaginary part of $Z$, whose maximum
gives the fastest growing mode, and well as the real part of $Z$, which
correspond to the shock width. From here and below, we restrict our study to
the strong shock regime with $r=4$. Figure 2 displays the imaginary part of
$Z$ in terms of $x$, given by Equation (21). The red dashed line is the
solution without CRs ($\epsilon=0$).

The CRs modified dispersion equation now evidences 2 modes. The first one is
simply the unstable mode arising from the interaction of the two ion beams,
modified by the presence of the CRs. But another unstable mode arises from the
presence of the CRs, and extends to higher $x$’s. For values of $\epsilon\sim
0.3$ and above, it becomes the dominant mode. Given this result, one can
characterize analytically how the most unstable $Z$ varies with $\epsilon$ and
$a$, which we do now.

Note that for high values of $\epsilon$, the growth-rate of the mode triggered
by the CRs becomes quite flat around the maximum. In this case, the most
unstable mode does not grow so much faster than its neighbors. Could this
invalidate Tidman’s analysis according to which the fastest growing mode
governs the linear phase of the instability? We show in Appendix B that it
does not.

Figure 3: Largest Im$(Z)$ with $Z$ root of (21), with the corresponding
Re$(Z_{max})$ and $x_{max}$ for $(\epsilon,a)\in[10^{-2},0.5]\times[4,5]$.

###  3.3 Fastest growing mode in the presence of CRs: analytical expression

For each couple $(\epsilon,a)$ we searched for the most unstable mode in
Figure 2 (highest Im$(Z)$) and systematically computed its corresponding $x$
and Re$(Z)$. As stated above, we found that the value of $x$ of the fastest
growing mode is always near $x\equiv x_{max}=1$. By definition, we set
$Z(x_{max})\equiv Z_{max}$. The most important quantity for our purpose is
Re$(Z_{max})$, since it defines the width of the shock through Equation (12).
It is plotted in Figure 3, alongside $x_{max}$ and Im$(Z_{max})$ as a function
of the cosmic rays power law index, $a$ and their fraction $\epsilon$.

At low values of $\epsilon\lesssim 0.3$ (the threshold slightly varies with
$a$), the ion-ion modified branch leads, implying that the shock width is only
slightly modified by the existence of CRs. However, as higher values of
$\epsilon$, the CRs triggered instability becomes dominant. At any rate,
Re$(Z_{max})$ hardly varies for $a\in[4,5]$. Although we present results
restricted to $a\in[4,5]$, we numerically checked that this weak dependence on
$a$ extends all the way up to $a=15$.

When the CRs triggered instability leads (blue region of Figure 3-center),
Re$(Z_{max})\sim 0.2$. When the two-stream modified branch leads (yellow
region of the same plot), Re$(Z_{max})\sim 1.65$. This is written as

${\rm Re}(Z_{max})\equiv
G(\epsilon)=\left\\{\begin{array}[]{l}0.2,~{}~{}~{}\epsilon\gtrsim\epsilon_{0},\\\
1.65,~{}~{}~{}\epsilon\lesssim\epsilon_{0},\end{array}\right.$ (22)

with,

$\epsilon_{0}\sim 0.36-\frac{a-4}{16}.$ (23)

Finally, following Tidman’s analysis, the width of the front is still obtained
from Equation (12),

${\rm
Re}(Z_{max})=\frac{k_{m}u_{1}}{\omega_{p1}}~{}~{}\Rightarrow~{}~{}k_{m}=\frac{\omega_{p1}}{u_{1}}{\rm
Re}(Z_{max}),$ (24)

yielding now, using $\lambda=A/k_{m}$, a modified Eq. (4), namely,

$\lambda=A\frac{u_{1}}{\omega_{p1}}~{}\frac{1}{G(\epsilon)}.$ (25)

We are finally in a position to contrast the two analysis conducted above.
Replacing $\lambda$ in Equation (2) by its value from Equation (25) gives

$\displaystyle a$ $\displaystyle=$ $\displaystyle
4+\frac{1}{6}\frac{u_{1}}{D}~{}A\frac{u_{1}}{\omega_{p1}}~{}\frac{1}{G(\epsilon)},$
(26) $\displaystyle=$ $\displaystyle 4+\frac{\theta}{G(\epsilon)},$

with

$\displaystyle\theta$ $\displaystyle\equiv$
$\displaystyle\frac{A}{6}\frac{u_{1}^{2}}{D\omega_{p1}}=\frac{A}{6}\frac{u_{1}}{\omega_{p1}}\frac{u_{1}}{D}=\frac{A}{6}\frac{\mathcal{M}_{1}c_{s1}}{\omega_{p1}}\frac{u_{1}}{D}$
(27) $\displaystyle\sim$
$\displaystyle\frac{A}{6}\frac{\mathcal{M}_{1}v_{th1}}{\omega_{p1}}\frac{u_{1}}{D}=\frac{A}{6}\mathcal{M}_{1}\frac{\lambda_{D1}u_{1}}{D},$

where $\lambda_{D1}u_{1}/D$ is eventually a modified Péclet number, where the
width of the shock $\lambda$ is replaced by the upstream Debye length
$\lambda_{D1}$.

##  4 Conclusions

In this work, we derived expressions relating the shock width $\lambda$ to the
CRs spectral index $a$ and the fraction $\epsilon$ of CRs in the upstream
region. Thereby, our results also connect $\lambda$ with $a$, which we found
have only a weak connection on each other. A key result is a very sharp
transition that occurs when $\epsilon\simeq 30\%$, as both the shock width and
the CR index abruptly change, due to a new instability mode becoming dominant.

We implemented two different approaches, which enabled us to derive analytical
expressions for these two quantities in terms of the upstream shock
properties. The first approach dates back to Drury et al. (1982a); Schneider &
Kirk (1987, 1989), where the index $a$ was computed for a given front width
$\lambda$. The second approach relies on Tidman (1967) who estimated $\lambda$
through an instability analysis. In his original work, though, CRs were
ignored. As we showed here, it turns out that when accounting for the CRs, a
similar instability analysis shows that the shock width is sensitive to the
fraction $\epsilon$ of upstream particles accelerated to become CRs, and to a
lesser extent, to their spectral index $a$. Combined together, we derived
closed analytical expressions for $a$ and $\lambda$ as functions of
$\epsilon$, namely our Equations (25, 26). Together with Figure 2, these are
the main results of this work.

The expression for the diffusion coefficient $D$ used here represents an
average diffusion over the shock front. While clearly one expects the
diffusion coefficient to vary along the shock, this approach enables
analytical calculations. Indeed, this is the same approach that was used by
Schneider & Kirk (1987); Achterberg & Schure (2011). It can be physically
justified if the variation of the diffusion coefficient along the shock front
is not large.

A key result of this work is a jump in the value of the CRs power law index
$a$ that occurs when $\epsilon$ exceeds $\sim 0.3$: from $4+\theta/1.65$ to
$4+\theta/0.2$. This is further accompanied by an abrupt change of the CRs
spectral index, $a$. If this approximation holds for $\theta$ as large as
$\theta\sim 15/6=2.5$, as suggested by the results in Appendix A, it implies a
sharp increase in the value of $a$, from $a\simeq 5.5$ to $a\simeq 16.5$. This
high value of $a$ effectively implies a cutoff in the CRs spectrum - the shock
will reach a structure that prevents it from further accelerating particles to
high energies. This result therefore puts strong limits on the possible values
of both $\epsilon$, as well as the CRs spectral index, $a$.

The determination of the fraction of upstream particles “promoted” to CRs is
still an open question. Due to the complexity of the problem and to the lack
of a complete theory, this problem is currently mainly studied through PIC
simulations. Such simulations typically yield $\epsilon\sim$ a few % (Gargaté
& Spitkovsky, 2012; Sironi et al., 2013; Marcowith et al., 2016). It is
established by now that while this fraction does not exceed a few % for
relativistic shocks, it is higher for non-relativistic ones, and was
numerically observed to be $\sim 5$ % (Gargaté & Spitkovsky, 2012). However,
existing simulations are computationally limited and reproduce but a small
fraction of the shock lifetime in astronomical objects. Furthermore, numerical
evidence points towards a continuous increase of CRs number (Keshet et al.,
2009; Groselj et al., 2024), following developments of magnetic fields at
increasingly larger coherence scale. Therefore, it is likely that values of
$\epsilon$ of tens of % could be achieved in nature.

The result of our work suggests an upper limit to the fraction of CRs in the
plasma, of $\epsilon\simeq 0.3$. If confirmed, this could be used to constrain
astronomical sources of CRs, and their dynamical effects in objects where non-
relativistic collisionless shocks may be present, such as stellar winds or
accretion disks. It could further put strong constraints on the CRs spectral
index $a$, which is expected to be limited to a relatively narrow range of
$4\leq a\leq 5.5$. This is easy to validate observationally.

As particles are continuously accelerated, two scenarios can be envisioned.
The first is a steady state, in which $\epsilon$ saturates close to, but less
than, $0.3$, while the CRs spectral index is fixed. A second scenario where
the acceleration continues until $\epsilon$ reaches the limiting value of
$0.3$. Once this happens, the shock structure is so heavily modified that it
is effectively destroyed, and CRs acceleration stops. CRs already present
likely escape, and the shock starts to re-build itself. Currently, we are not
able to determine which of these scenarios is more plausible, and we believe
future PIC simulations could shed light on this question.

Figure 4: Schematic representation of the causal relations between the four
main ingredient of a collisionless shock: front width, cosmic rays, velocity
profile and upstream/downstream density ratio. The arrows indicate causal
effects, with the references cited here. “BO78” stands for Blandford &
Ostriker (1978) and “SK87” for Schneider & Kirk (1987). “Drury82” and “SK87”
stand in bold because they are considered here. The present work adds 2 causal
effects, pictured in red, as we explore how cosmic rays could regulate the
front width via plasma instabilities. In this work, the density ratio
regulates the ratio $u_{1}/u_{2}$ in Eq. (13), hence the result of the
instability analysis.

Several important ingredients are omitted in the current analysis.

  1. 1.

The effect of a magnetic field on the shock structure should be accounted for.
Magnetic fields can be of external or internal origin. External fields can
originate from a magnetized central object (e.g., a star) accreting or
ejecting material, which gives rise to collisionless shocks. Internal fields
are expected to be self-generated, alongside the acceleration of CRs. PIC
simulations indicate a gradual increase in the characteristic coherence scale
of these fields, accompanying the acceleration of CRs to increasingly higher
energies (Keshet et al., 2009; Sironi & Spitkovsky, 2009, 2011; Groselj et
al., 2024). However, whether they affect the structure of a shock wave remains
to be studied.

  2. 2.

Our model is 1D. Yet, and still in connection with magnetization, what may
happen in 3D is that once injected, CRs start to efficiently generate a
turbulence that corrugates the shock front, modifying the local magnetic
obliquity and the local injection rate (Caprioli & Spitkovsky, 2014b). In such
circumstances, the front width becomes position dependent, which could
profoundly affect our instability analysis.

  3. 3.

We chose the simplest velocity profile for analytical tractability, as we want
to highlight the importance of the instability analysis carried here as a tool
in determining the shock width. However, it is known that the velocity profile
depends on the CRs when their fraction becomes important (Eichler, 1984;
Blasi, 2002). Likewise, the shock density ratio of 4 is likely to increase
with the CRs fraction (Ellison & Eichler, 1985; Bret, 2020).

A velocity change upstream could be accounted for within the framework of our
calculation, by changing the value of $u_{1}$ in Equation (13). We numerically
checked the validity of our results for a large range of CRs index, namely
$a\in[4,15]$, which can result from different velocity profiles and
compression ratios (see Equations 1 & 2). Therefore, we expect a more
realistic upstream velocity profile to lead to a result similar to the one
presented here. As demonstrated by the cartoon in figure 4, a more exact
evaluation requires a holistic approach of the problem, as the four main
ingredients of a shock are interrelated. Such an assessment is beyond the
scope of our work.

  4. 4.

Finally, our analysis assumes a constant Péclet number, whereas in reality it
could vary in space and with the energy of the CRs through the diffusion
coefficient $D$.

The present work eventually introduces a new ingredient to the physics of
collisionless shocks, namely the regulation of the front width through
counter-streaming instabilities when the CRs fraction rises (see red arrow on
figure 4). It will be interesting to see if this new ingredient dissolves in
the others, or ends up modifying the whole recipe for collisionless shocks.

##  5 Acknowledgments

We thank Anatoly Spitkovsky, Siddhartha Gupta, Damiano Caprioli and Lorenzo
Sironi for enriching inputs. A.B. acknowledges support from the Ministerio de
Economía y Competitividad of Spain (Grant No. PID2021-125550OB-I00). A.P.
acknowledges support from the European Research Council via ERC consolidating
grant No. 773062 (acronym O.M.J.).

##  Appendix A Proof of Equation (2)

Figure 5: Thick front model, with a linear transition in velocity space
between the upstream and the downstream.

We consider the velocity profile pictured in Figure 5 with a shock front
extending from $z=-\lambda/2$ to $z=+\lambda/2$ and,

$\displaystyle z<-\lambda/2~{}~{}$ $\displaystyle\Rightarrow$
$\displaystyle~{}~{}u(z)=u_{2},$ $\displaystyle-\lambda/2<z<\lambda/2~{}~{}$
$\displaystyle\Rightarrow$
$\displaystyle~{}~{}u(z)=\frac{u_{1}+u_{2}}{2}+\frac{u_{1}-u_{2}}{\lambda}z,$
$\displaystyle z>\lambda/2~{}~{}$ $\displaystyle\Rightarrow$
$\displaystyle~{}~{}u(z)=u_{1}.$ (A1)

The flow goes from right to left, so $u_{1,2}<0$. The equation for $f(p,z)$
for isotropic CRs reads (Skilling, 1971; Drury et al., 1982a; Schneider &
Kirk, 1987),

$u\frac{\partial f}{\partial z}=\frac{\partial}{\partial
z}\left(D\frac{\partial f}{\partial
z}\right)+\frac{1}{3}\frac{du}{dz}p\frac{\partial f}{\partial p}.$ (A2)

We now assume that the distribution $f$ is of the form $f(p,z)=\alpha
g(z)p^{-a}$, with $(\alpha,a)\in\mathbb{R}^{2}$. We also consider $D>0$ is a
constant. Eq. (A2) simplifies to,

$u(z)\frac{\partial g}{\partial z}=D\frac{\partial^{2}g}{\partial
z^{2}}-\frac{a}{3}\frac{du}{dz}g(z).$ (A3)

Note that $\alpha$ simplifies everywhere: the solution is defined up to a
numerical factor $\alpha$.

We now write this equation in the 3 domains for $u(z)$. For $z<-\lambda/2$ and
$z>+\lambda/2$, the algebra can be found for example in Kulsrud (2005), from
page 381. It does not depend on $\lambda$ and is therefore identical to the
case $\lambda=0$ (sharp front).

###  A.1 Case $z<-\lambda/2$

There, $u(z)=u_{2}<0$, so,

$\frac{\partial}{\partial z}\left(D\frac{\partial g}{\partial
z}-u_{2}g\right)=0.$ (A4)

The term between parenthesis is therefore a constant that we name $c$. The
equation $Dg^{\prime}-u_{2}g=c$ has solution,

$g=c_{1}\exp\left(\frac{u_{2}}{D}z\right)-\frac{c}{u_{2}},~{}~{}~{}c_{1}\in\mathbb{R}.$
(A5)

Since $u_{2}z\rightarrow+\infty$ when $z\rightarrow-\infty$, we must have
$c_{1}=0$. Therefore,

$\displaystyle g(z<-\lambda/2)$ $\displaystyle=$
$\displaystyle-\frac{c}{u_{2}},$ (A6) $\displaystyle\left.D\frac{\partial
g}{\partial z}\right|_{z=-(\lambda/2)^{-}}$ $\displaystyle=$ $\displaystyle
0.$

###  A.2 Case $z>+\lambda/2$

There, $u(z)=u_{1}<0$, so,

$\frac{\partial}{\partial z}\left(D\frac{\partial g}{\partial
z}-u_{1}g\right)=0.$ (A7)

Again the term between parenthesis is constant. It is 0 far upstream since
there are no CRs there111We assume no pre-existing CRs are present far
upstream., so it is 0 for any $z>+\lambda/2$. The solution is,

$\displaystyle g(z>+\lambda/2)$ $\displaystyle=$ $\displaystyle
c_{2}\times\exp\left(\frac{u_{1}}{D}z\right),~{}~{}~{}c_{2}\in\mathbb{R},$
(A8) $\displaystyle\left.D\frac{\partial g}{\partial
z}\right|_{z=+(\lambda/2)^{+}}$ $\displaystyle=$ $\displaystyle
u_{1}g(+\lambda/2).$

Since the full solution is defined up to a multiplying constant, we set here
$c_{2}=1$.

###  A.3 Determination of the power index $a$

We here provide a direct expression for the power index $a$. We start
integrating (A2) from $z=-\lambda/2$ to $z=+\lambda/2$,

$\displaystyle\int_{-\lambda/2}^{+\lambda/2}u(z)\frac{\partial g}{\partial
z}dz$ $\displaystyle=$ $\displaystyle\underbrace{\left.D\frac{\partial
g}{\partial
z}\right|_{z=+\lambda/2}}_{=u_{1}g(+\lambda/2)}-\underbrace{\left.D\frac{\partial
g}{\partial
z}\right|_{z=-\lambda/2}}_{=0}-\frac{a}{3}\int_{-\lambda/2}^{+\lambda/2}\frac{du}{dz}g(z)dz$
(A9) $\displaystyle=$ $\displaystyle
u_{1}g(+\lambda/2)-\frac{a}{3}\int_{-\lambda/2}^{+\lambda/2}\frac{du}{dz}g(z)dz$

We then integrate by part the left-hand-side and replace $du/dz$ by its value
in the integrand on the right-hand-side,

$u_{1}g(+\lambda/2)-u_{2}g(-\lambda/2)-\frac{u_{1}-u_{2}}{\lambda}\int_{-\lambda/2}^{+\lambda/2}g(z)dz=u_{1}g(+\lambda/2)-\frac{a}{3}\frac{u_{1}-u_{2}}{\lambda}\int_{-\lambda/2}^{+\lambda/2}g(z)dz.$
(A10)

This gives, introducing the density ratio $r$ with $u_{1}=ru_{2}$,

$\displaystyle a$ $\displaystyle=$ $\displaystyle
3\left(\frac{u_{2}g(-\lambda/2)}{\frac{u_{1}-u_{2}}{\lambda}\int_{-\lambda/2}^{+\lambda/2}g(z)dz}+1\right),$
(A11) $\displaystyle=$ $\displaystyle 3\frac{\chi-1+r}{r-1},$

with,

$\chi\equiv\frac{\lambda
g(-\lambda/2)}{\int_{-\lambda/2}^{+\lambda/2}g(z)dz}.$ (A12)

Clearly $\lim_{\lambda=0}\chi=1$ so the sharp front result (1) is retrieved.

###  A.4 Resolution for $-\lambda/2<z<+\lambda/2$ and expression of the power
index $a$

With profile (A), Eq. (A3) now reads,

$\displaystyle\left(\frac{u_{1}+u_{2}}{2}+\frac{u_{1}-u_{2}}{\lambda}z\right)\frac{\partial
g}{\partial z}$ $\displaystyle=$ $\displaystyle D\frac{\partial^{2}g}{\partial
z^{2}}-\frac{a}{3}\frac{u_{1}-u_{2}}{\lambda}g(z),$
$\displaystyle\partial_{z}g(-\lambda/2^{+})$ $\displaystyle=$ $\displaystyle
0,$ $\displaystyle g(+\lambda/2^{-})$ $\displaystyle=$
$\displaystyle\exp\left(\frac{u_{1}}{D}\frac{\lambda}{2}\right),$ (A13)

with boundary conditions coming from (A6,A8).

Now we use the notation $\partial_{z}g=g^{\prime}$. We also set $z=\lambda\xi$
so that $\partial_{z}=\lambda^{-1}\partial_{\xi}$. Introducing $u_{1}=ru_{2}$,
we get for $g(\xi)$ (derivatives are now with respect to $\xi$),

$\displaystyle\frac{D}{u_{2}\lambda}g(\xi)^{\prime\prime}-\left(\frac{r+1}{2}+(r-1)\xi\right)g(\xi)^{\prime}-a\frac{r-1}{3}g(\xi)$
$\displaystyle=$ $\displaystyle 0,$ (A14) $\displaystyle g^{\prime}(-1/2)$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle g(+1/2)$ $\displaystyle=$
$\displaystyle\exp\left(\frac{r}{2}\frac{u_{2}\lambda}{D}\right).$

For a strong shock with $r=4$, this gives,

$\displaystyle-\frac{1}{\Lambda}g(\xi)^{\prime\prime}-\left(\frac{5}{2}+3\xi\right)g(\xi)^{\prime}-a~{}g(\xi)$
$\displaystyle=$ $\displaystyle 0,$ (A15) $\displaystyle g^{\prime}(-1/2)$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle g(+1/2)$ $\displaystyle=$
$\displaystyle\exp\left(-2\Lambda\right),$

where222The “-” sign in the exponential comes from the definition of
$\Lambda$, where $|u_{2}|$ is considered instead of $u_{2}$.,

$\Lambda\equiv\frac{\lambda|u_{2}|}{D}.$ (A16)

The solution is found with _Mathematica_ , in terms of $a$ and $\Lambda$. It
reads

$\displaystyle g(\xi)$ $\displaystyle=$
$\displaystyle\frac{e^{-2\Lambda}(\mathcal{A}+\mathcal{B})}{\mathcal{C}+\mathcal{D}},~{}~{}~{}\mathrm{where},$
(A17) $\displaystyle\mathcal{A}$ $\displaystyle=$
$\displaystyle\sqrt{-\Lambda}\,_{1}F_{1}\left(1+\frac{a}{6};\frac{3}{2};\frac{-\Lambda}{6}\right)H_{\frac{a}{3}}\left(\frac{\sqrt{-\Lambda}(6\xi+5)}{2\sqrt{6}}\right),$
$\displaystyle\mathcal{B}$ $\displaystyle=$
$\displaystyle\sqrt{6}H_{-\frac{a}{3}-1}\left(\frac{\sqrt{-\Lambda}}{\sqrt{6}}\right)\,_{1}F_{1}\left(\frac{a}{6};\frac{1}{2};\frac{1}{24}-\Lambda(6\xi+5)^{2}\right),$
$\displaystyle\mathcal{C}$ $\displaystyle=$
$\displaystyle\sqrt{-\Lambda}H_{-\frac{a}{3}}\left(2\sqrt{\frac{2}{3}}\sqrt{-\Lambda}\right)\,_{1}F_{1}\left(1+\frac{a}{6};\frac{3}{2};\frac{-\Lambda}{6}\right),$
$\displaystyle\mathcal{D}$ $\displaystyle=$
$\displaystyle\sqrt{6}H_{-\frac{a}{3}-1}\left(\frac{\sqrt{-\Lambda}}{\sqrt{6}}\right)\,_{1}F_{1}\left(\frac{a}{6};\frac{1}{2};\frac{8-\Lambda}{3}\right),$

where $H_{n}$ and ${}_{1}F_{1}$ are the Hermite polynomial of order $n$ and
the Kummer confluent hypergeometric function, respectively. The following
method then allows to derive $a$,

  * •

Compute $\int_{-1/2}^{1/2}g(\xi)d\xi$ in terms of $a,\Lambda$.

  * •

Then compute $\chi$ from Eq. (A12).

  * •

From the value of $\chi$, compute the corresponding value of $a$ from Eq.
(A11). Call it $a_{1}$, function of $a,\Lambda$.

  * •

Solve $a=a_{1}(a,\Lambda)$ to obtain the final $a(\Lambda)$.

Figure 6: Numerical determination of $a(\Lambda/4)$ with
$\Lambda=\lambda/(D/|u_{2}|)$ where $\lambda$ is the front width. The
horizontal axis considers $\Lambda/4$ instead of $\Lambda$ for a direct
comparison with Figure 2a of Schneider & Kirk (1987), where the $\eta$
parameter in Figure 2a is $\lambda u_{1}/D$ instead of $\lambda u_{2}/D$.

The power index $a$ can be computed numerically. A Taylor expansion can be
obtained for $\Lambda\ll 1$. It reads,

$a=4+\frac{2}{3}\Lambda+\mathcal{O}(\Lambda^{2}),$ (A19)

which, accounting for $\Lambda=\lambda u_{2}/D$ and $u_{2}=u_{1}/4$, is the
result of Eq. (2). Figure 6 compares the numerical value of $a$ with
expression (A19). The horizontal axis considers $\Lambda/4$ instead of
$\Lambda$ for a direct comparison with Figure 2a of Schneider & Kirk (1987),
where the $\eta$ parameter in Figure 2a is $\lambda u_{1}/D$ instead of
$\lambda u_{2}/D$.

The agreement between the numerical calculation and Equation (A19) is indeed
very good up to $\Lambda/4\sim 15$. For the numerical calculation, Figure
2a/profile (a) of Schneider & Kirk (1987) is retrieved.

##  Appendix B Instability analysis accounting for the vicinage of the most
unstable mode

Let us model the flatness of the function Im$(Z)(x)$ around its maximum
Im$(Z_{max})$, reached for $x_{m}$, through

$\mathrm{Im}(Z)(x)\sim\mathrm{Im}(Z_{max})-b(x-x_{m})^{2},$ (B1)

where

$b=\frac{1}{2}\left.\frac{d^{2}\mathrm{Im}(Z)}{dx^{2}}\right|_{x=x_{m}}.$ (B2)

$b\ll 1$ gives a flat spectrum, and $b\gg 1$ an extremely peaked one around
the maximum. Here, $b$ is of order unity. For example, for $a=4$ and
$\epsilon=0.5$, we find numerically $b\sim 0.36$.

Translated to $(\omega,k)$ variables, Eq. (B1) gives,

$k_{i}(\omega)\sim k_{i,m}-\frac{b}{u_{1}\omega_{p1}}(\omega-\omega_{m})^{2},$
(B3)

where $k_{i}$ is the imaginary part of $k$ and $k_{i,m}$ its maximum value,
reached for $\omega=\omega_{m}$. Assume now all modes from
$\omega_{m}-\Delta\omega$ to $\omega_{m}+\Delta\omega$ are excited with
similar amplitude. The wave packet evolves on average like

$\frac{1}{2\Delta\omega}\int_{\omega_{m}-\Delta\omega}^{\omega_{m}+\Delta\omega}\exp\left[\left(k_{i,m}-b^{\prime}(\omega-\omega_{m})^{2}\right)z\right]d\omega=e^{k_{i,m}z}~{}\frac{\sqrt{\pi}}{2}\frac{\mathrm{erf}\left(\sqrt{(\Delta\omega)^{2}b^{\prime}z}\right)}{\sqrt{(\Delta\omega)^{2}b^{\prime}z}},$
(B4)

where erf$(x)$ is the error function and $b^{\prime}=b/u_{1}\omega_{p1}$.
Therefore, the wave packet grows like $e^{k_{i,m}z}$, times a factor function
of $(\Delta\omega)^{2}b^{\prime}z$. With $\Delta\omega\sim\omega_{p1}$, we
obtain $(\Delta\omega)^{2}b^{\prime}z\sim b\omega_{p1}z/u_{1}$ so that when
the most unstable mode will have propagated 1 $e$-folding length
$z=k_{i,m}^{-1}\sim u_{1}/\omega_{p1}$, the quantity
$(\Delta\omega)^{2}b^{\prime}z$ will be of order unity, together with the
correction to the exponential growth $e^{k_{i,m}z}$.

The situation is eventually very similar to the “vanilla” two-stream
instability, found between 2 symmetric cold electron beams over a background
of fixed protons. With the present dimensionless variables $(x,Z)$, the
dispersion equation reads $(x-Z)^{-2}+(x+Z)^{-2}=1$ (Boyd & Sanderson (2003),
§6.5.1). It can be solved exactly, giving a second derivative at the maximum
growth rate such that the constant $b$ defined by Eq. (B1) is $b=3/4$.

Therefore, in analogue to the two-stream instability where the fastest growing
mode governs the linear phase, while the nearby modes only slightly modify the
sin structure of the growing wave, we expect here a similar pattern. As in
Tidman’s analysis, the shock structure will be mainly determined by the
fastest growing mode, with the nearby modes only slightly modifying its
structure.

## References

  * Achterberg & Schure (2011) Achterberg, A., & Schure, K. M. 2011, MNRAS, 411, 2628, doi: 10.1111/j.1365-2966.2010.17868.x
  * Blandford & Ostriker (1978) Blandford, R., & Ostriker, J. 1978, Astrophysical Journal, 221, L29
  * Blandford & Ostriker (1978) Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29, doi: 10.1086/182658
  * Blasi (2002) Blasi, P. 2002, Astroparticle Physics, 16, 429, doi: 10.1016/S0927-6505(01)00127-X
  * Boyd & Sanderson (2003) Boyd, T. J. M., & Sanderson, J. J. 2003, The Physics of Plasmas (Cambridge University Press)
  * Bret (2020) Bret, A. 2020, ApJ, 900, 111, doi: 10.3847/1538-4357/aba68e
  * Bret & Pe’er (2021) Bret, A., & Pe’er, A. 2021, Journal of Plasma Physics, 87, 905870204, doi: 10.1017/S002237782100012X
  * Caprioli et al. (2015) Caprioli, D., Pop, A.-R., & Spitkovsky, A. 2015, ApJ, 798, L28, doi: 10.1088/2041-8205/798/2/L28
  * Caprioli & Spitkovsky (2014a) Caprioli, D., & Spitkovsky, A. 2014a, ApJ, 783, 91, doi: 10.1088/0004-637X/783/2/91
  * Caprioli & Spitkovsky (2014b) —. 2014b, ApJ, 794, 47, doi: 10.1088/0004-637X/794/1/47
  * Drury et al. (1982a) Drury, L. O., Axford, W. I., & Summers, D. 1982a, MNRAS, 198, 833, doi: 10.1093/mnras/198.3.833
  * Drury et al. (1982b) —. 1982b, MNRAS, 198, 833, doi: 10.1093/mnras/198.3.833
  * Eichler (1984) Eichler, D. 1984, ApJ, 277, 429, doi: 10.1086/161711
  * Ellison & Eichler (1985) Ellison, D. C., & Eichler, D. 1985, Phys. Rev. Lett., 55, 2735, doi: 10.1103/PhysRevLett.55.2735
  * Gargaté & Spitkovsky (2012) Gargaté, L., & Spitkovsky, A. 2012, ApJ, 744, 67, doi: 10.1088/0004-637X/744/1/67
  * Groselj et al. (2024) Groselj, D., Sironi, L., & Spitkovsky, A. 2024, arXiv e-prints, arXiv:2401.02392, doi: 10.48550/arXiv.2401.02392
  * Keshet et al. (2009) Keshet, U., Katz, B., Spitkovsky, A., & Waxman, E. 2009, ApJ, 693, L127, doi: 10.1088/0004-637X/693/2/L127
  * Kruells & Achterberg (1994) Kruells, W. M., & Achterberg, A. 1994, A&A, 286, 314
  * Kulsrud (2005) Kulsrud, R. M. 2005, Plasma physics for astrophysics (Princeton, NJ: Princeton Univ. Press)
  * Landau & Lifshitz (1981) Landau, L., & Lifshitz, E. 1981, Course of Theoretical Physics, Physical Kinetics, Vol. 10 (Elsevier, Oxford)
  * Li et al. (2013) Li, G., Kong, X., Zank, G., & Chen, Y. 2013, ApJ, 769, 22, doi: 10.1088/0004-637X/769/1/22
  * Marcowith & Kirk (1999) Marcowith, A., & Kirk, J. G. 1999, A&A, 347, 391, doi: 10.48550/arXiv.astro-ph/9905176
  * Marcowith et al. (2016) Marcowith, A., Bret, A., Bykov, A., et al. 2016, Reports on Progress in Physics, 79, 046901, doi: 10.1088/0034-4885/79/4/046901
  * Mott-Smith (1951) Mott-Smith, H. M. 1951, Physical Review, 82, 885, doi: 10.1103/PhysRev.82.885
  * Piran (2004) Piran, T. 2004, Rev. Mod. Phys., 76, 1143
  * Sagdeev (1966) Sagdeev, R. Z. 1966, Reviews of Plasma Physics, 4, 23
  * Salas (2007) Salas, M. D. 2007, Shock Waves, 16, 477, doi: 10.1007/s00193-007-0084-z
  * Schneider & Kirk (1987) Schneider, P., & Kirk, J. G. 1987, ApJ, 323, L87, doi: 10.1086/185063
  * Schneider & Kirk (1989) —. 1989, A&A, 217, 344
  * Sironi & Spitkovsky (2009) Sironi, L., & Spitkovsky, A. 2009, ApJ, 698, 1523, doi: 10.1088/0004-637X/698/2/1523
  * Sironi & Spitkovsky (2011) —. 2011, ApJ, 726, 75, doi: 10.1088/0004-637X/726/2/75
  * Sironi et al. (2013) Sironi, L., Spitkovsky, A., & Arons, J. 2013, ApJ, 771, 54, doi: 10.1088/0004-637X/771/1/54
  * Skilling (1971) Skilling, J. 1971, ApJ, 170, 265, doi: 10.1086/151210
  * Spitkovsky (2008) Spitkovsky, A. 2008, ApJ, 682, L5, doi: 10.1086/590248
  * Tidman (1967) Tidman, D. A. 1967, Physics of Fluids, 10, 547, doi: 10.1063/1.1762148
  * Treumann (2009) Treumann, R. A. 2009, A&A Rev., 17, 409, doi: 10.1007/s00159-009-0024-2
  * Xu et al. (2022) Xu, Y. D., Li, G., & Yao, S. 2022, ApJ, 932, 87, doi: 10.3847/1538-4357/ac6ce4
  * Zel’dovich & Raizer (2002) Zel’dovich, Y. B., & Raizer, Y. P. 2002, Physics of shock waves and high-temperature hydrodynamic phenomena (Dover Publications), xxvii, 916 p.
  * Zhang (2020) Zhang, B. 2020, Nature, 587, 45, doi: 10.1038/s41586-020-2828-1