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. 1. The density jump r represents the ratio between the upstream and downstream densities. It is 4 in a strong sonic shock.

2. 2. The velocity profile represents the way the plasma velocity evolves spatially from the upstream to the downstream.

3. 3. The width of the shock front λ represents the distance over which this transition takes place.

4. 4. Finally, cosmic rays represent the particles accelerated by the shock, mainly characterized by their power index a.

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.

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 ₁, 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,

(1)\begin{equation} f(p) = n_1(z)\delta(p-m_iU_1) + n_2(z)\delta(p-m_iU_2) + n_{cr}f_{cr}(p), \end{equation}

where m_i 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,

(2)\begin{equation} f_{cr}(p) = \kappa p^{-a}, \quad p\in [P_{min} , +\infty], \end{equation}

where

(3)\begin{equation} \kappa = \frac{a-1}{P_{min}^{1-a}} \end{equation}

is a normalization constant and P_min 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

(4)\begin{align} &\frac{4\pi \frac{N_1}{2}q^2 }{m_i}\frac{1}{(\omega - kU_1)^2} + \frac{4\pi (1-\epsilon) \frac{N_1}{2}q^2 }{m_i}\frac{1}{(\omega - kU_2)^2} \nonumber\\& \quad +\frac{4\pi \epsilon N_1 q^2 }{m_i} \kappa \int \frac{a p^{-1-a}}{\omega - kp/m_i}dp = 1, \end{align}

where q is the elementary charge. We now introduce,

(5)\begin{equation} \omega_{p1} = \frac{4\pi N_1 q^2}{m_i}, \end{equation}

(6)\begin{equation} x = \frac{\omega}{\omega_{p1}}, \end{equation}

(7)\begin{equation} Z = \frac{kU_1}{\omega_{p1}}. \end{equation}

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 ¹, Eq. (4) eventually reads,

(8)\begin{equation} \frac{1/2}{(x-Z)^2} + \frac{(1-\epsilon)/2}{(x-Z/r)^2} - \epsilon\frac{a(a-1)}{2xZ}\varphi\left(a,\frac{2Z}{x} \right) = 1, \end{equation}

where

(9)\begin{equation} \varphi(\alpha,\beta) = 2\int_1^{+\infty}\frac{t^{-\alpha-1}}{t^2\beta^2-1}dt. \end{equation}

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,

(10)\begin{equation} \mathrm{Re}(Z) \sim 1 \quad \Rightarrow \quad \frac{1}{\mathrm{Re}(k)} \sim \frac{U_1}{\omega_{p1}}. \end{equation}

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,

(11)\begin{equation} \lambda \sim A \frac{U_1}{\omega_{p1}}, \quad A = \mathcal{O}(10). \end{equation}

Suppose that instead of $\mathrm{Re}(Z) \sim 1$, we have $\mathrm{Re}(Z) \sim l$. Then,

(12)\begin{equation} \lambda \sim A \frac{U_1}{\omega_{p1}}\frac{1}{l}, \end{equation}

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.

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.

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

(13)\begin{equation} a = \frac{3r}{r-1}, \end{equation}

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),

(14)\begin{equation} a \sim 4 + \frac{1}{6}\frac{\lambda}{D/U_1}, \end{equation}

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).