On the Width of a Collisionless Shock and the Index of the Cosmic Rays It Accelerates

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 (CRs), which, in turn, modify the shock structure. The shock width λ is known to be connected to the CR 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 CR index and to the fraction of upstream particles that are accelerated. We find that when this fraction becomes larger than ∼30%, a new instability becomes dominant. The shock undergoes a transition where its width increases by a factor ∼8–10, and the CR acceleration effectively ends. Our analysis is valid for strong, nonrelativistic, and unmagnetized shocks. We discuss the implication of these results on the expected range of CR spectra and flux observed and on the structure of nonrelativistic collisionless shocks.


INTRODUCTION
Shock waves are a common phenomenon in many environments.Their study dates back to the 19 th 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 20 th 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 ∝ 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 λ: a = a(r, λ).
A natural scale for the shock width λ is the ion inertial length, λ i = u i /ω p,i , where u i is the ion velocity and ω p,i the ion plasma frequency (e.g., Treumann 2009).However, the shock width λ itself is affected by the existence of CRs.
Therefore, one can expect λ 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 λ 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 ǫ.
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 λ.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 M appear in some of our equations, the end results are established considering M ≫ 1 and r ∼ 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 ∼ 30%.We conclude in Section 4 discussing the implications of our results on the shock structure and its observational consequences.

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 (λ → 0), the relation derived in Blandford & Ostriker (1978) for the power index a reads, 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, where u 1 is the upstream velocity and D the diffusion coefficient of the CRs.The dimensionless quantify λu 1 /D is known as the "Péclet number" of the shock (Zel'dovich & Raizer 2002).Although derived through a Taylor expansion for λ ≪ D/u 1 , the expression is indeed quite accurate up to λ ∼ 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).

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, where n 1,2 (z) are the weights of the upstream and downstream Maxwellians f 1,2 (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 where u 1 and ω 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.

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 ∝ exp(ikz − iωt) reads (see for example Landau & Lifshitz (1981), §29), where q is the ions' charge and v = p/m.Let us consider a distribution of the form, 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 δ ′ g = − δg ′ , inserting Equation (6) into Equation (5) yields, Denoting by r the shock density ratio, u 2 = u 1 /r.With this and n 1,2 = N 1 /2, Equation ( 7) becomes, Introducing the plasma frequency, and the normalized frequency and wavevector, 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 Im(Z max ) ∼ 1 near x = 1, with Re(Z max ) ∼ 1. Denoting by k m the real part of the most unstable wavevector k, we therefore find, and Tidman's prescription λ = 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 ω.

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, where n cr ≡ ǫN 1 is the CRs density, ǫ 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 Particle-in-cell (PIC) simulations of cosmic ray acceleration typically yield ǫ ∼ 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 ǫ up to ǫ = 0.5.
For f cr (p), we consider where a is the power index and κ the normalizing constant, given by It is common in literature to introduce ζ, the fraction of the shock kinetic energy imparted into CRs.Given the distribution (13), the energy contained in CRs in our scenario reads ǫ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 ǫN 1 , the number of CRs, rather than on their energy, we shall not refer to ζ in the sequel.
Since the thermal particles are the seed particles for acceleration to high energies, one has P min ∼ 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 = M 1 c s1 , where M 1 is the shock Mach number and c s1 the upstream speed of sound.One thus finds The CRs-modified dispersion equation eventually reads, with x ∈ R and Z ∈ C. 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 (ǫ = 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 ǫ ∼ 0.3 and above, it becomes the dominant mode.Given this result, one can characterize analytically how the most unstable Z varies with ǫ and a, which we do now.
Note that for high values of ǫ, 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.

Fastest growing mode in the presence of CRs: analytical expression
For each couple (ǫ, 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 ≡ x max = 1.By definition, we set Z(x max ) ≡ 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 ǫ.
At low values of ǫ 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 ǫ, the CRs triggered instability becomes dominant.At any rate, Re(Z max ) hardly varies for a ∈ [4, 5].Although we present results restricted to a ∈ [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 ) ∼ 0.2.When the two-stream modified branch leads (yellow region of the same plot), Re(Z max ) ∼ 1.65.This is written as with, Finally, following Tidman's analysis, the width of the front is still obtained from Equation ( 12), yielding now, using λ = A/k m , a modified Eq. ( 4), namely, We are finally in a position to contrast the two analysis conducted above.Replacing λ in Equation ( 2) by its value from Equation ( 25) gives with where λ D1 u 1 /D is eventually a modified Péclet number, where the width of the shock λ is replaced by the upstream Debye length λ D1 .

CONCLUSIONS
In this work, we derived expressions relating the shock width λ to the CRs spectral index a and the fraction ǫ of CRs in the upstream region.Thereby, our results also connect λ with a, which we found have only a weak connection on each other.A key result is a very sharp transition that occurs when ǫ ≃ 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 λ.The second approach relies on Tidman (1967) who estimated λ 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 ǫ 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 λ as functions of ǫ, 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.

Drury82
Present work . 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 u1/u2 in Eq. ( 13), hence the result of the instability analysis.
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 ǫ exceeds ∼ 0.3: from 4 + θ/1.65 to 4 + θ/0.2.This is further accompanied by an abrupt change of the CRs spectral index, a.If this approximation holds for θ as large as θ ∼ 15/6 = 2.5, as suggested by the results in Appendix A, it implies a sharp increase in the value of a, from a ≃ 5.5 to a ≃ 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 ǫ, 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 ǫ ∼ 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 ∼ 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 ǫ 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 ǫ ≃ 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 ≤ a ≤ 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 ǫ saturates close to, but less than, 0.3, while the CRs spectral index is fixed.A second scenario where the acceleration continues until ǫ 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.
Several important ingredients are omitted in the current analysis.
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. 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. 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 ∈ [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. 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.We consider the velocity profile pictured in Figure 5 with a shock front extending from z = −λ/2 to z = +λ/2 and, 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), ) We now assume that the distribution f is of the form f (p, z) = αg(z)p −a , with (α, a) ∈ R 2 .We also consider D > 0 is a constant.Eq. (A2) simplifies to, Note that α simplifies everywhere: the solution is defined up to a numerical factor α.
We now write this equation in the 3 domains for u(z).For z < −λ/2 and z > +λ/2, the algebra can be found for example in Kulsrud (2005), from page 381.It does not depend on λ and is therefore identical to the case λ = 0 (sharp front).
Since the full solution is defined up to a multiplying constant, we set here c 2 = 1.
1 We assume no pre-existing CRs are present far upstream.

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 = −λ/2 to z = +λ/2,

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

Figure 2 .
Figure 2. Imaginary part of Z in terms of x given by Equation (21), for various choices of (ǫ, a).The red dashed line is the solution without CRs (ǫ = 0), for r = 4.The curve that is close to it and falls to 0 for x ∼ 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.

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