Introduction
Shock waves are fundamental phenomena in fluids and plasmas. And collisionless shock waves represent a special kind of shock wave. Whereas in a shock wave in a fluid, the mean free path is very small compared to the dimensions of the system, in a collisionless shock wave, the mean free path is very large. For example, in situ measurements showed the width of the Earth bow shock in the solar wind is about 100 km, while the proton mean free path at this location is about the Sun–Earth distance (Refs Reference Bale, Mozer and Horbury3, Reference Schwartz, Henley, Mitchell and Krasnoselskikh28). This type of shock waves can only exist in a plasma, since they are mediated by collective electromagnetic effects instead of binary collisions (Refs Reference Balogh and Treumann4, Reference Sagdeev26).
In a collisional medium, any excess energy given to one particle is quickly shared with the others through collisions, on a time scale of the collision frequency. In a collisionless plasma, on the other hand, it is possible to give energy to one particle without it being immediately shared with the others.
It was thus realized in the late 70s that collisionless shocks can accelerate particles with a non-thermal spectrum of the kind p −a, where p is the momentum of the accelerated particles and a is the power index (Refs Reference Axford, Leer and Skadron2, Reference Bell5–Reference Blandford and Ostriker7).
Since the spectrum of cosmic rays detected on Earth from space obeys a power law (Ref. Reference Cronin14), or more precisely a succession of power laws, collisionless shocks in astrophysical media, such as those found in supernova remnant (Ref. Reference Drury15), are excellent candidates for explaining the origin of cosmic rays. As such, they have been the subject of theoretical, numerical and experimental studies for decades (see (Ref. Reference Marcowith, Bret, Bykov, Dieckman, O’CDrury, Lembège, Lemoine, Morlino, Murphy, Pelletier, Plotnikov, Reville, Riquelme, Sironi and Stockem Novo24) and references therein).
Producing this kind of shocks in the laboratory, and observing accelerated particles, currently requires installations like the National Ignition Facility (Ref. Reference Fiuza, Swadling, Grassi, Rinderknecht, Higginson, Ryutov, Bruulsema, Drake, Funk, Glenzer, Gregori, Li, Pollock, Remington, Ross, Rozmus, Sakawa, Spitkovsky, Wilks and Park19). Numerical studies, on the other hand, can only probe a short part of the shock’s life after its formation (Refs Reference Grošelj, Sironi and Spitkovsky21, Reference Keshet, Katz, Spitkovsky and Waxman22). It is therefore important to pursue theoretical studies which, while they can only address simplified models of the real process, allow to explore scenarios that are currently beyond the reach of experiments or simulations.
While cosmic rays acceleration in collisionless shocks has been explored since Fermi’s times (Ref. Reference Fermi18), the total amount of energy that goes into cosmic rays is still unclear. An open question of interest to both theoreticians and observers who search for cosmic rays from various objects.
The aim of this paper is to propose, on a theoretical basis, a mechanism that could stop the acceleration process within a collisionless shock and set a limit to the maximum fraction of particles promoted to cosmic rays.
Sections 2 and 3 introduce the mechanism under scrutiny, which relies on an instability analysis. Section 4 then derives the dispersion equation of the instability. It is solved in Section 5, with the consequences discussed in Section 6 and following.
The ingredients of a collisionless shock and their connections
In order to introduce the mechanism we are proposing, we now present the various ingredients of a collisionless shock. They are schematically pictured on Figure 1.
1. The density jump r represents the ratio between the upstream and downstream densities. It is 4 in a strong sonic shock.
2. The velocity profile represents the way the plasma velocity evolves spatially from the upstream to the downstream.
3. The width of the shock front λ represents the distance over which this transition takes place.
4. Finally, cosmic rays represent the particles accelerated by the shock, mainly characterized by their power index a.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_fig1.png?pub-status=live)
Figure 1. The various ingredients of a collisionless shock and their connections. The red arrows represent the new connections proposed here. The snowflakes picture the connections we froze in the present work, namely, that we did not consider (see ‘Dispersion equation’ section).
As it happens, these four ingredients are interconnected. The power index a depends on the density jump r (Ref. Reference Blandford and Ostriker7). But it also depends on the width of the shock front λ (Ref. Reference Drury, Axford and Summers16). For example, the greater the width of the shock front, the greater the power index. Furthermore, as we have just seen, a change in the power index a is equivalent to a change in r, but the very presence of cosmic rays can also change r, when the energy they carry becomes substantial (Refs Reference Bret10, Reference Eichler17). Finally, there is a bidirectional relationship between a and the velocity profile, since the latter can change the value of a (Ref. Reference Schneider and Kirk27), while cosmic rays, when their pressure becomes significant, in turn affect the velocity profile (Ref. Reference Blasi8).
We propose here some as yet unexplored links between these four ingredients, which could provide a mechanism to terminate the acceleration process. It is a causal relationship between the density jump r, the cosmic ray power index a and the width of the shock front λ. These new connections are pictured by the red arrows on Figure 1.
The origin of the new link proposed
The theoretical study of the density profile of a shock is a difficult subject. In 1951, Harold Mott-Smith proposed an ansatz for studying such a profile within the framework of kinetic theory (Ref. Reference Mott-Smith25). He hypothesized that the particle distribution function along the shock profile is a linear combination of the upstream and downstream Maxwellians. Indeed this hypothesis has been validated by particle-in-cell (PIC) simulations (Ref. Reference Spitkovsky32).
In the Mott-Smith picture, see Figure 2, the parameters of the upstream Maxwellian (density, drift speed, temperature) constitute the problem inputs; those of the downstream Maxwellian are given by the Rankine–Hugoniot relations, and only the respective weights of these two Maxwellians depend on the position z along the shock profile. The weight of the upstream Maxwellian vanishes as one progresses into the shock, while the weight of the downstream one grows.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_fig2.png?pub-status=live)
Figure 2. The Mott-Smith ansatz. The particle distribution function along the shock profile is a linear combination of the upstream and downstream Maxwellians. The weight of the former vanishes as one progresses into the shock, while the weight of the latter grows. The far upstream and downstream flows have velocities and density
$U_{1,2}$ and
$N_{1,2}$ respectively.
In 1967, Derek Tidman used Mott-Smith’s model to study a collisionless shock (Ref. Reference Tidman34). Among other results, he concluded that the thickness of the shock front can be given by an instability analysis at a location where the two Maxwellians have approximately the same weight. Since the upstream Maxwellian is centred around the upstream flow velocity U 1, while the downstream one is centred around the downstream flow velocity
$U_2 = U_1/r$, this velocity shift results in an unstable system. Tidman analysed the instability of this system in terms of the most unstable wavelength, not the most unstable frequency, and related this wavelength to the thickness of the front, since it corresponds to the length over which the upstream flow is disrupted.
Now, Tidman’s analysis did not account for the presence of cosmic rays, that is, the particles accelerated by the shock. It only accounted for the unstable interaction between the two Maxwellians.
Here, we shall revisit Tidman’s analysis including the population of cosmic rays accelerated by the shock. As we shall see, beyond a critical amount of upstream particles promoted to cosmic rays, the instability analysis provides a radically different answer to the question of the width of the front. In turn, such a change in the front width can radically change the power index of the cosmic rays.
Since Tidman’s analysis was conducted for an unmagnetized shock, we only consider here the same kind of shocks, propagating without an external magnetic field.
Dispersion equation
Figure 1 makes it clear that it is not possible to elaborate a theoretical model accounting for all the connections involved in the problem. We therefore chose to ‘freeze’ those signalled by the snowflakes on Figure 1. Hence, we set the density ratio r to 4 in the sequel, ignoring the back reaction of the cosmic rays on the same ratio. We also considered a linear velocity profile between the upstream and the downstream, ignoring the back reaction of the cosmic rays on the same profile.
We therefore analyse the unstable system composed of the upstream and downstream Maxwellians, plus a population of cosmic rays. We implement a 1D model where location along the shock profile is identified by the z coordinate. In order to simplify the phase space of parameters, we neglect temperature effects. The distribution function under scrutiny is, therefore,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn1.png?pub-status=live)
where mi is the ion mass, δ the Dirac delta function and
$U_{1,2}$ the upstream and downstream flow speeds respectively. The first term refers to the upstream flow, the second one to the downstream flow and the third one to the cosmic rays. They obey a power law spectrum of the form,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn2.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn3.png?pub-status=live)
is a normalization constant and Pmin the injection momentum.
We shall conduct the instability analysis at a location z where
$n_1 \sim n_2 = N_1/2$. Yet, we shall consider a fraction ϵ of the upstream flow has been ‘promoted’ to cosmic rays so that we set
$n_2 =$
$N_1/2 - \epsilon$ and
$n_{cr}=\epsilon N_1$.
The derivation of the dispersion equation is standard (Ref. Reference Landau and Lifshitz23) and yields
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn4.png?pub-status=live)
where q is the elementary charge. We now introduce,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn5.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn6.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn7.png?pub-status=live)
Considering
$P_{min} \sim 2 m_i U_1$ (Refs Reference Caprioli, Pop and Spitkovsky12, Reference Caprioli and Spitkovsky13) and noting that
$U_2 = U_1/r$, where r is the compression ratio of the shockFootnote 1, Eq. (4) eventually reads,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn8.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn9.png?pub-status=live)
We now solve Eq. (8) for
$(Z,x) \in \mathbb{C} \times \mathbb{R}$, since we are interested in the most unstable wavelength.
Resolution of the dispersion equation
Setting r = 4, the numerical resolution of the dispersion equation reveals the presence of one unstable mode when cosmic rays are turned off, and two when they are turned on, as explained on Figure 3.
In the absence of cosmic rays, namely for ϵ = 0 (red curves on Figure 3), there is only one unstable mode, arising from the interaction of the two flows. Its maximum spatial growth rate has Im
$(Z) \sim 1$ and also Re
$(Z) \sim 1$ (not shown). The corresponding most unstable wavelength has, therefore,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn10.png?pub-status=live)
Within Tidman’s analysis, the width of the shock λ is proportional to this quantity, up to a factor A, the ‘Tidman’s constant’, of order 10. We thus recover Tidman’s result, established without cosmic rays, namely,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn11.png?pub-status=live)
Suppose that instead of
$\mathrm{Re}(Z) \sim 1$, we have
$\mathrm{Re}(Z) \sim l$. Then,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn12.png?pub-status=live)
which will be relevant shortly.
When cosmic rays are turned on, this two-flows unstable mode gets modified. It is pictured by the black circles on Figure 3. But now, another unstable mode appears (green circles), fully triggered by the presence of the cosmic rays. For small values of ϵ, this new unstable mode grows slower than the two-flows one, and the width of the shock is not modified. But for higher values of ϵ, like ϵ = 0.5 on Figure 3-right, this new unstable mode grows faster. It is, therefore, the real part of this most unstable Z which now defines the width of the shock front.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_fig3.png?pub-status=live)
Figure 3. Imaginary part of Z solution of Eq. (8). Dashed red curve: two-flows unstable mode, that is, solution without cosmic rays, namely for ϵ = 0. Black circles: two-flows unstable mode, modified by the cosmic rays. Green circles: new unstable mode, triggered by the presence of the cosmic rays.
We scanned the parameters phase space
$(\epsilon, a) \in [0,1]\times [4,15]$. For each couple
$(\epsilon, a)$, we computed the most unstable Z and plotted its real part on Figure 4.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_fig4.png?pub-status=live)
Figure 4. Real part of the most unstable Z solution of the dispersion equation (8). For small values of ϵ, the most unstable mode remains the one triggered by the interaction of the upstream and downstream flows. But for ϵ > 0.3, the most unstable modes arises from the presence of cosmic rays, with a Re(Z) significantly smaller than before. The red arrows picture the temporal evolution of the system (see discussion in Section 7).
One can navigate this plot following the time evolution of a strong collisionless shock. As it starts to propagate,
$\epsilon \ll 1$ and
$a \sim 4$, that is, the power index of a strong shock. The amount of accelerated particles then increases with time, but Re(Z) does not, as evidenced by the plateau at low ϵ.
Yet, once ϵ reaches
$\sim 0.3$, an abrupt transition occurs. The most unstable mode switches from the two-flows one to the cosmic rays triggered one, and Re(Z) suddenly falls from 1.65 to 0.2, almost 10 times smaller. As a result, and according to Eqs. (10,11,12), the front width λ increases almost 10 folds.
Consequences on the power index a
Back in Section 2, we referred to the connection between the front width and the power index. In (Ref. Reference Blandford and Ostriker7), the front was considered a discontinuity. Within this approximation, it was found that a shock could accelerate particles with a power index a given by
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn13.png?pub-status=live)
where r is the density ratio. For a strong sonic shock with r = 4, this gives a = 4. It was then found in (Refs Reference Drury, Axford and Summers16, Reference Schneider and Kirk27) that considering a finite width λ ≠ 0 for the front implies a larger value of a, because it is then more difficult for particles to go back and forth around the front and undergo Fermi cycles. Considering a simple linear transition over a distance λ for the velocity field and r = 4, the numerical result derived in (Ref. Reference Schneider and Kirk27) can be expressed as (Ref. Reference Bret and Pe’er11),
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20250207131310847-0886:S0263034624000053:S0263034624000053_eqn14.png?pub-status=live)
where D is the diffusion coefficient of the cosmic rays. Strictly speaking, this quantity depends on the location along the shock profile and of the momentum of the particles. Here, like in (Refs Reference Achterberg and Schure1, Reference Schneider and Kirk27), we consider an average value of D.
The dimensionless quantity
$\lambda U_1/D$ is the shock Péclet number. Equation (14) is accurate up to
$\lambda U_1/D \sim 15$, beyond which the exact value keeps growing, though slower than the linear trend (Ref. Reference Bret and Pe’er11).
Discussion
We can now draw the conclusion of these calculations following the temporal evolution of the shock on Figure 4.
At the beginning of its history, the density jump is r = 4, with a power index a = 4. Then, as it propagates, ϵ grows, as the shock accelerates more and more particles. It then follows the trajectory indicated by the red arrow on the upper plateau.
When the fraction of accelerated particles reaches about 30%, the most unstable mode becomes the one triggered by the cosmic rays, and the system falls onto the lower plateau, following the vertical red arrow downwards.
At this point, and according to Eqs. (10,11,12), the width of the front increases. In turn, according to Eq. (14), the power index increases. Suppose the front width λ was such that before the transition we had
$\lambda U_1/6D = 1$ in Eq. (14), hence a = 5. After the transition, with λ increased 10 fold,
$\lambda U_1/6D$ jumps from 1 to 10, and a from 5 to 14. Such a high power index means acceleration stops, as the extension of the shock front becomes so large that it forbids Fermi cycles.
To our knowledge, the maximum fraction of upstream particles promoted to cosmic rays is an open question, mainly studied through PIC simulations. Some long ran simulations show values of ϵ reaching about 5% (Refs Reference Caprioli, Pop and Spitkovsky12, Reference Gargaté and Spitkovsky20, Reference Sironi and Spitkovsky30). Yet, the longest simulations to date only capture a short fraction of the shock lifetime (Refs Reference Grošelj, Sironi and Spitkovsky21, Reference Keshet, Katz, Spitkovsky and Waxman22). Since the fraction ϵ grows with time (Refs Reference Caprioli and Spitkovsky13, Reference Keshet, Katz, Spitkovsky and Waxman22, Reference Sironi, Spitkovsky and Arons31), it could reach the threshold commented here in real settings.
Strictly speaking, Eq. (14), which relates the power index a to the front width λ, assumes
$\epsilon \ll 1$ since it is derived ignoring the back-reaction of the cosmic rays on the velocity profileFootnote 2. Values of
$\epsilon \sim 0.3$ invalidate this assumption. Therefore, as specified from Figure 1 and in Section 4, we ignored, ‘froze’, such back-reaction in order to obtain an analytically tractable model. Yet, we think the mechanism presented here could be robust enough to withstand a more realistic scenario.
Noteworthily, Tidman’s initial analysis was performed for an electrostatic shock whereas electromagnetic shocks mediated by the Weibel instability have a different structure (Ref. Reference Stockem, Fiuza, Bret, Fonseca and Silva33). However, as far as the shock width is concerned, the key ingredient here is the instability of the two, namely upstream and downstream, shifted Maxwellians. This coexistence of the two Maxwellians has also been observed in PIC simulations of electromagnetic shocks (Ref. Reference Spitkovsky32). We, therefore, think we can apply it for this kind of shocks.
Besides, the acceleration mechanism considered in our work is the Fermi process contemplated in (Refs Reference Blandford and Ostriker7, Reference Drury, Axford and Summers16, Reference Schneider and Kirk27) and for which Eq. (14) stands. It differs from laser shock acceleration where the upstream ions gain energy from only one interaction with the shock (Refs Reference Boella, Fiúza, Novo, Fonseca and Silva9, Reference Silva, Marti, Davies, Fonseca, Ren, Tsung and Mori29, Reference Stockem, Fiuza, Bret, Fonseca and Silva33).
In summary, we have revisited Tidman’s analysis of the width of a collisionless shock. While Tidman did not account for the shock accelerated cosmic rays, we did include them in our instability analysis. The result is that the most unstable wavelength is dramatically increased when the faction of accelerated particles reaches
$\epsilon \sim 30\%$. Hence, we have uncovered an unexplored mechanism that could put an end to particle acceleration in a collisionless shock. In this picture, the width of the shock front is given by an instability analysis between the upstream and downstream flows. When this analysis takes into account the cosmic rays accelerated by the shock, the most unstable wavelength is radically altered when the fraction of accelerated particles exceeds around 30%. This result is a potentially abrupt increase of the index of accelerated particles, which could simply mean the end of the acceleration process. Very long-term PIC simulations could help validate our scheme.
Data availability statement
The data used to support the findings of this study are available from the corresponding author upon request.
Acknowledgements
A.B. acknowledges the support from the Ministerio de Economía y Competitividad of Spain (Grant No. PID2021-125550OB-I00). A.P. acknowledges the support from the European Research Council via ERC Consolidating Grant No. 773062 (acronym O.M.J.). Thanks are due to Luis Silva for helpful inputs.
Conflicts of interest
The authors declare that they have no conflicts of interest.