© The Authors 2026

1. Introduction

The influence of cosmic rays (CRs) over their own transport is crucial to explain numerous phenomena, from the diffusive properties of astrophysical plasmas, to the maximum energy reached by accelerators and the subsequent spectra we measured at the Earth (see Blasi 2019 for a review of these effects). In recent decades, substantial progress has been made in understanding how CRs shape the environment in which they are produced. This is especially true for the vicinity of shockwaves, which are expected to be key sites of CR acceleration across a wide range of astrophysical settings, from supernova remnants (SNRs) (Blasi 2013; Amato 2014) to galaxy clusters (Blasi et al. 2007).

Despite the success of early theoretical approaches to diffusive shock acceleration (DSA), in which CRs were mostly treated as test particles (Krymskii 1977; Bell 1978), it soon became clear that to account for CRs up to the knee region or even higher energy CRs, efficient shock acceleration and strong magnetic fields–likely powered by CR-driven instabilities–were needed. In particular, the resonant streaming instability was invoked as a process responsible for magnetic-field amplification (MFA) upstream of nonrelativistic shocks (Skilling 1975; Bell 1978) and proven to lead to acceleration up to tens of teraelectronvolts in typical SNRs (Lagage & Cesarsky 1983). This process has the advantage of producing Alfvén waves on the right scale to scatter the same particles responsible for exciting the instability. On the other hand, the saturation of the instability to δB/B ≲ 1 implies that this process cannot account for particle acceleration to the knee. More recently, a nonresonant branch of the streaming instability (Bell 2004) received much attention: this instability is excited by the current induced by the CR particles escaping the upstream region of a shock and can lead to δB/B ≫ 1 due to its nonresonant nature. The instability is expected to saturate when particles start scattering on the self-generated perturbations. The implications of the excitation of the nonresonant modes in terms of the maximum energy of particles accelerated in SNR shocks have been discussed in a number of articles (Bell et al. 2013; Schure & Bell 2013; Cristofari et al. 2020, 2021). The consensus among these works is that although the instability appears to be very effective in amplifying the magnetic field on the relevant scales for very fast shocks in dense environments, the maximum energy at the beginning of the Sedov phase of SNRs typically falls short of a petaelectronvolt by at least one order of magnitude. This leaves open the possibility that other instabilities may lead to more effective MFA and perhaps more effective particle acceleration.

One such mechanism is associated with the pressure gradient that unavoidably forms upstream of a shock due to energy-dependent particle transport, which may induce the so-called Drury or acoustic instability (AI) in the presence of density inhomogeneities. This instability was first proposed by Drury (1984) and studied in detail by Drury & Falle (1986) (hereafter DF) using a two-fluid prescription in the absence of magnetic fields. This instability can overcome dissipation by CR diffusion (Ptuskin 1981) if the scale height of the CR pressure is smaller than the ratio between the local diffusion coefficient and sound speed, a condition naturally met at shockwave precursors. This results in an exponential growth of small inhomogeneities present in the interstellar medium (ISM) as they approach the shock surface. Numerical solutions to the equations in DF confirm the onset of the instability and the growth of density perturbations to nonlinear levels, if enough time is granted.

These findings have prompted further investigation of AI, especially in the context of CR-modified shocks (see Zank et al. 1990; Kang et al. 1992; Ryu et al. 1993; Begelman & Zweibel 1994; Webb et al. 1999, and references therein), although the nonlinear evolution of the system and the implications for MFA remain poorly understood. Beresnyak et al. (2009) proposed that the CR pressure gradient acting on large-amplitude density inhomogeneities upstream of a strong shock could induce turbulence, and that this might be reflected in higher maximum energies of the accelerated particles. In this picture, the solenoidal fluid motions act as a small-scale dynamo, amplifying magnetic fields by at least one order of magnitude via the stretch-twist-fold mechanism.

Drury & Downes (2012) (hereafter, DD) and Downes & Drury (2014) ran 2D magnetohydrodynamical (MHD) simulations of the precursor region, assuming that CRs were only responsible for the formation of an assigned pressure gradient. These simulations confirmed that MFA actually occurs. In 2D simulations, slightly more effective magnetic energy amplification is found, with a weak predominance of power on larger scales compared to the results of 3D simulations. More recently, del Valle et al. (2016) (hereafter, VLS) explored and expanded on this idea by means of 3D MHD simulations of SNR shock precursors. Their work provided quantitative insights into the development of small-scale dynamo, turbulence, and scales where equipartition between fluid and magnetic turbulence should be expected. We comment further on these results throughout this paper.

Much of this previous work is based on a few critical assumptions: the first is that a very large part (typically 60%) of the ram pressure is channeled into accelerated particles. Although early work on the nonlinear theory of DSA appeared to suggest that such an energy conversion efficiency would be feasible (Malkov & Drury 2001; Blasi 2002), later work, based on nonlinear DSA models (Morlino & Caprioli 2012) and simulations (Caprioli & Spitkovsky 2014), suggests a typical efficiency of ≲10%. The question of MFA in higher Mach numbers with more realistic shock precursors remains open. The second assumption is that the pre-existing density fluctuations in the upstream plasma are already relatively large: this leaves little room for AI, and in fact much of the previous work on the topic focused on the possibility that the CR pressure gradient may induce turbulence in an inhomogeneous plasma. Both of these assumptions were adopted in order to minimize the demands on the numerical machinery employed to simulate the instability.

In this work, we focused on three main goals: 1) providing an analytic framework to assess the importance of AI in the precursor in the presence of a magnetic field and confirming the linear evolution of the instability with 2D MHD simulations; 2) bridging the gap with previous work, by starting with small density perturbations and investigating the possibility that turbulence may be formed due to AI first and warping of the magnetic-field lines later; 3) using more realistic CR acceleration efficiencies and Mach numbers that are more appropriate to SNRs in the beginning of the Sedov phase (M∼ few hundreds). In this paper, we take special care to discuss the scales accessible to the simulations (namely with negligible numerical damping), in 2D and 3D, compared with the wavenumbers that are excited by AI.

This article is organized as follows. In Sect. 2, we present an analytic derivation of the dispersion relation and the growth rate of AI. We derived simple expressions for the maximum number of e-folds by which small upstream perturbations can grow before reaching the shock. In Sect. 3, we describe how we used a setup similar to that of DD to perform MHD simulations of the instability, comparing the observed growth rates with the expected ones obtained analytically. We also studied the nonlinear stage of the system, discuss MFA, and calculate the power spectrum of magnetic perturbations. Finally, in Sect. 4 we comment on the comparison between MFA due to the CR pressure gradient and the Bell instability, speculating about a possible cooperation between the two processes. We summarize our findings and present our conclusions in Sect. 5.

2. Acoustic instability

We considered a plasma with density ρ, velocity u, pressure P, adiabatic index γ, and magnetic field B. Since we aim to describe the plasma dynamics upstream of a nonrelativistic shock, where a CR pressure gradient is expected, the MHD equations should include a force term due to the presence of a CR pressure gradient, ∇P_CR:

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathbf{u})=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}·\mathrm{\nabla })\mathbf{u}=-\frac{\mathrm{\nabla }P}{\rho }-\frac{\mathrm{\nabla }{P}_{\mathrm{CR}}}{\rho }+\frac{(\mathrm{\nabla }\times \mathbf{B})\times \mathbf{B}}{4\pi \rho },\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial }{\partial t}\left(\frac{P}{{\rho }^{\gamma }}\right)+(\mathbf{u}·\mathrm{\nabla })\left(\frac{P}{{\rho }^{\gamma }}\right)=0,\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathbf{B}}{\partial t}=\mathrm{\nabla }\times (\mathbf{u}\times \mathbf{B}),\hfill \end{array}$$(4)

$$\begin{array}{cc}& \mathrm{\nabla }·\mathbf{B}=0.\hfill \end{array}$$(5)

Drury & Falle (1986) modeled the evolution of the CR pressure through the CR transport equation with diffusion coefficient D, and did not consider the effect of the magnetic field. By simultaneously perturbing ρ, u, P, and P_CR, using a rigorous perturbative approach, they found that these perturbations grow when the length scale of the CR pressure gradient, P_CR/∇P_CR, is smaller than D/c_s, where $c_{\mathrm{s}}=\sqrt{\gamma P_0/{\rho }_0}$ is the sound speed. If this condition is not met, the perturbations are damped due to friction with the CRs (Ptuskin 1981). In shock precursors, damping effects can be safely neglected because P_CR/∇P_CR ∼ D/v_sh ≪ D/c_s, where v_sh ≫ c_s is the shock speed. As such, for the case of a shock precursor, it is not necessary to go through the elegant but cumbersome calculations of DF, and the underlying physics describing the instability can be captured by perturbing only the plasma quantities ρ, u, and P, while assuming that ∇P_CR is constant. Although this neglects the back-reaction of the instability onto CRs, it also reduces the complexity of the problem to a purely MHD one. We later comment on the implications of this assumption and on the reaction of the system to the enhanced scattering that may follow the excitation of the instability.

A clear advantage of the MHD formulation is that the dispersion relation and the growth rate of the instability can be derived straightforwardly. In the reference frame of the plasma, we can assume that the background plasma is uniform¹ for Eqs. (1)–(5), with ρ = ρ₀, u = 0, P = P₀, and B = B₀ while we introduce small perturbations–δρ, δu, δP, and δB–into the system. This approximation is valid for perturbations on scales much smaller than that of the CR pressure gradient. We assume that the perturbations are proportional to exp[i(k ⋅ x − ωt)]. The dispersion relation is

$$\begin{array}{cc}\hfill ({\omega }^2-v_{\mathrm{A}}^2{k}_{\Vert }^2)& [{\omega }^4-{\omega }^2(k^2{c}_{\mathrm{ms}}^2+\frac{i\mathbf{k}·\mathrm{\nabla }{P}_{\mathrm{CR}}}{{\rho }_0})\hfill \\ \hfill & +k^2{k}_{\Vert }^2{v}_{\mathrm{A}}^2{c}_{\mathrm{s}}^2+i{v}_{\mathrm{A}}^2{k}^2{k}_{\Vert }\frac{{\partial }_{\Vert }{P}_{\mathrm{CR}}}{{\rho }_0}]=0,\hfill \end{array}$$(6)

where $v_{\mathrm{A}}=B_0/\sqrt{4\pi {\rho }_0}$ is the Alfvén speed, $c_{\mathrm{ms}}=\sqrt{c_{\mathrm{s}}^2+v_{\mathrm{A}}^2}$ is the magnetosonic speed, and the ∥ subscript refers to components parallel to B₀.

In addition to the usual Alfvén waves (ω² = v_A²k_∥²), the dispersion relation shows that there are unstable modes obtained from roots of the expression inside the square brackets. When k_∥ = 0, their dispersion relation is

$$\begin{array}{c}\hfill {\omega }^2=c_{\mathrm{ms}}^2{k}^2+\frac{i\mathbf{k}·\mathrm{\nabla }{P}_{\mathrm{CR}}}{{\rho }_0},\end{array}$$(7)

while for |k|=k_∥ it is the same as in the absence of a magnetic field, namely

$$\begin{array}{c}\hfill {\omega }^2=c_{\mathrm{s}}^2{k}^2+\frac{i\mathbf{k}·\mathrm{\nabla }{P}_{\mathrm{CR}}}{{\rho }_0}·\end{array}$$(8)

The growth rate of the instability, Γ, is characterized by two distinct regimes, which depend on the ratio between the real and imaginary parts of ω². For example, Eq. (8) with k ∥ ∇P_CR yields

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\{\begin{array}{cc}\sqrt{\frac{k|\mathrm{\nabla }{P}_{\mathrm{CR}}|}{2{\rho }_0}},\hfill & k\ll \frac{|\mathrm{\nabla }{P}_{\mathrm{CR}}|}{{\rho }_0{c}_{\mathrm{s}}^2}\hfill \\ \frac{|\mathrm{\nabla }{P}_{\mathrm{CR}}|}{2{\rho }_0{c}_{\mathrm{s}}},\hfill & k\gg \frac{|\mathrm{\nabla }{P}_{\mathrm{CR}}|}{{\rho }_0{c}_{\mathrm{s}}^2}\hfill \end{array}.\end{array}$$(9)

The real part of ω, which determines the phase velocity of the perturbations, is equal to Γ for long wavelengths (k ≪ |∇P_CR|/ρ₀c_s²). Instead, the real part of ω is equal to c_sk for short wavelengths (k ≫ |∇P_CR|/ρ₀c_s²). In the latter regime, the perturbations are modified sound waves, because Γ ≪ c_sk. In the other regime, the condition that the perturbations must occur on scales smaller than the size of the precursor, L ≈ P_CR/∇P_CR, implies that the sonic Mach number must be $M_s=u_0/c_{\mathrm{s}}\gg {\xi }_{\mathrm{CR}}^{-1/2}$, where ξ_CR = P_CR/ρ₀u₀². This condition is certainly satisfied for shocks of astrophysical relevance.

To better understand the nature of the instability, we considered only modes with k ∥ ∇P_CR and computed the eigenmodes of the system. If the density perturbation is given by δρ = δρ₀ e^{i(k ⋅ x − ωt)}, one can show that pressure perturbations are δP = δρ₀c_s² e^{i(k ⋅ x − ωt)}, just as in regular sound waves, while velocity perturbations become

$$\begin{array}{c}\hfill \delta \mathbf{u}=\pm \frac{\mathrm{\Omega }}{k}\frac{\delta {\rho }_0}{{\rho }_0}{e}^{i(\mathbf{k}·\mathbf{x}-\omega t+\varphi /2)}\widehat{\mathbf{k}},\end{array}$$(10)

where $\mathrm{\Omega }=|\omega |=\sqrt[4]{c_s^4{k}^4+{(\mathbf{k}·\mathrm{\nabla }{P}_{\mathrm{CR}}/{\rho }_0)}^2}$ and ϕ = arg(ω²) = arctan(k ⋅ ∇P_CR/ρ₀c_s²k²) is the principal value (−π/2 < ϕ ≤ π/2). The key point behind the instability is that velocity and pressure perturbations are out of phase, as was also noted in previous work (e.g., Begelman & Zweibel 1994). Consider only the positive sign of Eq. (10). Assuming k ⋅ ∇P_CR = k|∇P_CR| (i.e., ω_R, ω_I > 0), we have 0 < ϕ < π/2, and therefore the pressure lags the velocity in phase, which leads to instability. This is similar to pushing a swing soon after it starts moving away, increasing its amplitude of oscillation. Meanwhile, we find −π/2 < ϕ < 0 if we take waves in the opposite direction (k ⋅ ∇P_CR = −k|∇P_CR|, or equivalently ω_R, ω_I < 0), making pressure precede velocity in phase. The result is akin to pushing a swing too early, damping its motion. The equivalent analysis taking the negative sign for δu in Eq. (10) reverses the effects of leading/lagging pressure and is left for the reader to carry out.

Figure 1 shows the real and complex parts of ω as functions of k, as predicted by Eqs. (7) and (8), along with their asymptotic behaviors. We assumed that k ∥ ∇P_CR and the CR pressure gradient is constant, |∇P_CR|=ξ_CRρ₀u₀²/L, where L is the length of the precursor and ξ_CR is the fraction of the shock ram pressure converted into CRs. The number of e-folds available for the instability to grow on the precursor scale can be estimated as Γt_adv, where t_adv ≈ L/u₀ is the advection timescale. In the absence of magnetic fields, Eq. (9) gives

$$\begin{array}{c}\hfill \mathrm{\Gamma }{t}_{\mathrm{adv}}=\{\begin{array}{cc}\sqrt{\frac{kL{\xi }_{\mathrm{CR}}}{2}},\hfill & kL\ll {\xi }_{\mathrm{CR}}{M}_s^2\hfill \\ \frac{{\xi }_{\mathrm{CR}}{M}_s}{2},\hfill & kL\gg {\xi }_{\mathrm{CR}}{M}_s^2\hfill \end{array}.\end{array}$$(11)

Fig. 1. Real and imaginary parts of ω. The CR pressure gradient is constant, |∇PCR|=ξCRρ0u02/L, and the perturbation wave vector, k, is aligned with ∇PCR. We assume ξCR = 0.1, Ms = 100, and $M_{\mathrm{ms}}=100/\sqrt{5}$. The dashed and dotted lines show the asymptotic behavior of ω.

The instability develops when Γt_adv> a few. If the magnetic field is perpendicular to ∇P_CR (and consequently to k), M_s should be replaced by the magnetosonic Mach number, M_ms = u₀/c_ms.

3. Simulations

In order to investigate the linear evolution of AI discussed in Sect. 2 and, more importantly, in order to study its nonlinear evolution, we performed MHD simulations using the PLUTO code (Mignone et al. 2007). Although most of our simulations are in 2D, we also compared the results of selected cases of physical interest with 3D simulations. Inspired by DD, our simulation box consists of a uniform 2D grid of dimensions N_x × N_y embedded in a Cartesian coordinate system with dimensions [0, L]×[0, L/8]. This region represents the precursor of a shockwave in the shock rest frame, with the shock itself lying just outside the right boundary at x = L. The x = 0 boundary represents upstream infinity (i.e., the ISM), where plasma of density ρ₀, pressure P₀, magnetic field B₀, and adiabatic index γ = 5/3 enters the simulation box at the shock speed ${\mathbf{u}}_0=u_0\widehat{\mathbf{x}}$.

The CR pressure gradient is introduced as a body force using the “VECTOR” prescription in PLUTO. For simplicity, we assumed that the CR pressure increases linearly from zero at x = 0 to a fraction, ξ_CR, of the shock ram pressure at x = L:

$$\begin{array}{c}\hfill P_{\mathrm{CR}}(x)={\xi }_{\mathrm{CR}}{\rho }_0{u}_0^2\frac{x}{L},\mathrm{\nabla }{P}_{\mathrm{CR}}(x)={\xi }_{\mathrm{CR}}\frac{{\rho }_0{u}_0^2}{L}\widehat{\mathbf{x}}.\end{array}$$(12)

The presence of ∇P_CR induces nonuniform profiles along x in the density, velocity, pressure, and magnetic field of the background plasma. The steady-state profiles can be found by solving Eqs. (1)–(5) in their conservative form:

$$\begin{array}{c}\hfill U(x)\equiv \frac{u(x)}{u_0}=\frac{{\rho }_0}{\rho (x)}={\left[\frac{P_0}{P(x)}\right]}^{1/\gamma }=\frac{B_{y,0}}{B_y(x)}\end{array}$$(13)

and B_x(x) = B_x, 0, where U(x) satisfies

$$\begin{array}{c}\hfill U+\frac{U^{-\gamma }}{\gamma M_s^2}+\frac{U^{-2}}{2M_{\mathrm{A},y}^2}+{\xi }_{\mathrm{CR}}\frac{x}{L}=1+\frac{1}{\gamma M_s^2}+\frac{1}{2M_{\mathrm{A},y}^2},\end{array}$$(14)

where M_A, y = u₀/v_A, y and $v_{\mathrm{A},y}=B_{y,0}/\sqrt{4\pi {\rho }_0}$. Note that this is different from the Alfvénic Mach number M_A = u₀/v_A.

Equation (14) can be solved numerically to find U(x) exactly. In the limit of high Mach numbers, one can approximate U(x)≈1 − ξ_CRx/L. The resulting density, velocity, pressure, and magnetic-field profiles are used as the initial conditions for the simulations. The time it takes for a fluid element entering from the left at x = 0 to reach a position x inside the box in the presence of a shock precursor is

$$\begin{array}{c}\hfill t(x)={\displaystyle {\int }_0^x}\frac{dx\prime }{u(x\prime )}\approx \frac{L}{u_0}\frac{1}{{\xi }_{\mathrm{CR}}}ln\left(\frac{1}{1-{\xi }_{\mathrm{CR}}x/L}\right)\gtrsim \frac{x}{u_0}·\end{array}$$(15)

PLUTO implements boundary conditions (BCs) in each timestep by updating the value of physical quantities in inactive “ghost” cells surrounding the active N_x × N_y simulation grid. We used periodic BCs at the y = 0 and y = L/8 boundaries and an outflow (i.e., zero-gradient) BC at the x = L boundary, so that the plasma cannot flow back from the shock into the precursor. We extrapolated the variables into the ghost cells at x = 0 to avoid discontinuities of the derivatives. To avoid unphysical pressure waves that propagate back into the simulation box, we set the pressure at the right boundary to a very low value.

In Sect. 3.1, we used the “HLLC” Riemann solver with “MP5” reconstruction and third-order Runge–Kutta time-stepping (with Courant number CFL = 0.4 for 2D simulations and CFL = 0.3 for 3D simulations), while in Sect. 3.2 the solver and reconstruction were respectively downgraded to “HLL” and “PARABOLIC” for numerical stability purposes. The ∇ ⋅ B = 0 condition was controlled using the “CONSTRAINED_TRANSPORT” method under the “UCT_HLL” scheme. “ENTROPY_SWITCH” was always kept off.

3.1. The linear regime

We considered the linear regime where the perturbations produced by AI are much smaller than the unperturbed quantities in the entire simulation box. At x = 0, the BC was set to

$$\begin{array}{c}\hfill \rho (x=0,y,t)={\rho }_0+\delta {\rho }_{0}sin(k_x{u}_{0}t+k_{y}y+\varphi ),\end{array}$$(16)

where 0 ≤ ϕ < 2π is a random phase, $\mathbf{u}=u_0\widehat{\mathbf{x}}$, P = P₀, and $\mathbf{B}=B_{x,0}\widehat{\mathbf{x}}$ or $\mathbf{B}=B_{y,0}\widehat{\mathbf{y}}$. Eq. (16) corresponds to a continuous inflow of density fluctuations with wavelengths λ_x = 2π/k_x and λ_y = 2π/k_y along the x and y directions. After one advection time, all transients due to initial conditions and boundary discontinuities at t = 0 disappear and the system reaches a quasi-stationary state. In Sects. 3.1.1 and 3.1.2 we study the behavior of k_x modes (where k_y = 0) and k_y modes (where k_x = 0) and test the growth rates of AI derived in Sect. 2. In Sect. 3.1.3 we investigate the more realistic scenario with many k_x and k_y modes.

3.1.1. Single k_x mode

The simulations in this section were performed on a default grid resolution of N_x = 4000 and N_y = 4. The chosen value of N_x is large enough to avoid significant numerical damping. Since ∇P_CR is oriented along the x-axis, density perturbations with k_y = 0 are expected to induce δu_x, δP, and δB_y perturbations that grow due to AI. In order to study the evolution of perturbations, we subtracted the background profiles from the total quantities provided as output by PLUTO, for example

$$\begin{array}{cc}\hfill \delta \rho (x,y,t)& =\rho (x,y,t)-\frac{{\rho }_0}{U(x)},\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill \delta u_x(x,y,t)& =u_x(x,y,t)-u_{0}U(x).\hfill \end{array}$$(18)

To assess the validity of the analytic results derived in Sect. 2, we performed a scan over many different values of the parameters ξ_CR, M_s, δρ₀/ρ₀, v_A/c_s, and k_xL = 2πλ_x⁻¹L. In Fig. 2, we show some selected simulation results at steady state. Each subplot depicts constant-y slices of density (black), pressure (red), velocity (blue), and magnetic-field (green) perturbations, along with the corresponding expected growths ∝exp[Γt(x)] (dashed lines), which are all consistent with the simulation results². In the top left panel, we adopted ξ_CR = 0.1, M_s = 100, δρ₀/ρ₀ = 10⁻⁴, k_xL = 80π, and no background magnetic field, while in the bottom left panel we add a magnetic field perpendicular to the shock normal. As expected, a perpendicular magnetic field suppresses the growth rate of the instability. We also confirmed that a parallel field ${\mathbf{B}}_0\Vert \widehat{\mathbf{x}}$ reproduces the zero-field result. In the top right panel, we probe the small-wavelength limit (k ≳ |∇P_CR|/ρ₀c_s² = ξ_CRM_s²/L), which requires a factor-ten increase in grid resolution (N_x = 40 000) to avoid numerical dissipation. Finally, in the bottom right panel, we consider larger values of δρ₀/ρ₀ and ξ_CR. In this case, the amplitude of the perturbations approaches the nonlinear regime close to the shock. As we discuss in Sect. 3.2, in this regime the perturbations steepen and eventually lead to the formation of shocklets.

Fig. 2. Quasi-stationary x-profiles of the dimensionless density, pressure, velocity, and magnetic field perturbations from kx modes. We show the simulation results as solid lines and the expected growths for AI derived in Sect. 2 as dashed lines. In the bottom left panel, the initial magnetic field is pointed along y, and the lighter gray line corresponds to the expected growth in the absence of magnetic fields.

3.1.2. Single k_y mode

The dispersion relations given by Eqs. (7) and (8) show that only perturbation wave vectors parallel to ∇P_CR lead to instability. However, perturbations with k_x = 0 are also affected by the CR pressure gradient because the acceleration −∇P_CR/ρ of a fluid element depends on the local density (overdensities decelerate slightly less than underdensities). For density perturbations along y, one has

$$\begin{array}{c}\hfill -\frac{\mathrm{\nabla }{P}_{\mathrm{CR}}}{{\rho }_0+\delta {\rho }_{0}sin(k_{y}y)}\approx -\frac{\mathrm{\nabla }{P}_{\mathrm{CR}}}{{\rho }_0}+\frac{\mathrm{\nabla }{P}_{\mathrm{CR}}}{{\rho }_0}\frac{\delta {\rho }_0}{{\rho }_0}sin(k_{y}y).\end{array}$$(19)

The first term corresponds to an overall deceleration of the fluid, while the second term leads to a relative acceleration of overdensities with respect to underdensities. The evolution of δu_x is given by (see Appendix A)

$$\begin{array}{c}\hfill \frac{\delta u_x(x,y)}{u_0}\approx \frac{\mathrm{\nabla }{P}_{\mathrm{CR}}}{{\rho }_0{u}_0^2}\frac{\delta {\rho }_0}{{\rho }_0}sin(k_{y}y)\frac{sin(k{\prime }_{x}x)}{k{\prime }_x},\end{array}$$(20)

where

$$\begin{array}{c}\hfill k{\prime }_x=k_y\sqrt{\frac{B_{y,0}^2}{4\pi {\rho }_0{u}_0^2{U}^3}}=\frac{k_y}{M_{A,y}{U}^{3/2}}·\end{array}$$(21)

Equation (20) can be interpreted as follows. Underdensities are more decelerated by the CR pressure gradient with respect to overdensities. Magnetic-field lines bend because they are frozen in the fluid, and magnetic tension acts as a restoring force. Then, δu_x oscillates while the fluid moves toward the shock. When there is no magnetic field, δu_x increases linearly.

This effect can be seen in Fig. 3, which shows the results of simulations performed on a grid with N_x = 2000 and N_y = 250. The top panels show the steady-state result of a simulation with the same parameters as in the top left panel of Fig. 2, but with a k_y = 40π mode instead of a k_x = 80π mode. On the left, we show δρ/ρ₀. Sinusoidal oscillations along the y-axis do not grow exponentially, as expected. On the right, we show δu_x. We considered different y slices: y = y_peak (red), which is a peak of δρ/ρ₀; y = y_trough (blue), which is a trough; and y = y_node (black), which is a node. The evolution also closely matches the prediction of Eq. (20) when B_y, 0 = 0; this is shown as dotted lines.

Fig. 3. Stationary profiles of dimensionless density (left) and velocity (right) perturbations from ky modes. On the right, red, blue, and black lines show the simulation results for δux/cs at a fixed y coordinate corresponding respectively to the peak, trough, and node of δρ/ρ0 oscillations. Dotted and dashed lines correspond to the analytical predictions of Eq. (20). Top: No magnetic fields. Bottom: Magnetic field perpendicular to the shock normal.

The same quantities are shown in the bottom panels of Fig. 3, but for a different choice of parameters–notably with the addition of a perpendicular magnetic field such that v_A = 3c_s. The density profile looks qualitatively similar to the one in the top panel, indicating that no significant forces arise along the y direction. However, the velocity perturbations oscillate around zero, rather than growing linearly, because of magnetic tension. We verified the appearance of δB_x in this simulation, thereby confirming our interpretation.

3.1.3. Many modes

We now address the more realistic case of ISM inhomogeneities across multiple scales and along all directions. The injected density profile is

$$\begin{array}{cc}\hfill \rho (x=0,y,t)& ={\rho }_0[1+\frac{\delta {\rho }_0}{{\rho }_0}(u_{0}t,y)],\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill \frac{\delta {\rho }_0}{{\rho }_0}(x,y)& ={\displaystyle \sum _{\begin{array}{c}a,b=0\end{array}}^{N_x,N_y}}{A}_{\mathit{ab}}cos[\frac{2\pi }{L}(ax+by)+{\varphi }_{\mathit{ab}}],\hfill \end{array}$$(23)

where A_ab and ϕ_ab are random amplitudes and phases of each perturbation. The latter are drawn uniformly from 0 ≤ ϕ_ab < 2π, while the former were chosen so as to produce a flat omnidirectional power spectrum of ISM density fluctuations, E_δρ0/ρ0(k)∝k⁻¹, with a root-mean-square (RMS) amplitude of δρ_0, rms/ρ₀ = 10⁻⁴. This procedure is described in Appendix B.1. Perturbations with wavelengths of > L/8 are set to have zero amplitude, while those with wavelengths down to grid-resolution scales are still injected into the system.

All power spectra in this work are calculated in the last eighth of the box (i.e., 7L/8 < x < L), since the amplitude of fluctuations is largest and nonlinear structures are most prominent in the part of the precursor that is closest to the shock. We also time-averaged the power spectra over several snapshots spanning one advection time, making them stationary.

Figure 4 shows the steady-state δρ/ρ₀, δu_x, y, and δB_x, y profiles (left) and the corresponding δρ/ρ₀ power spectra (right), from a simulation performed in a 2048 × 256 grid using ξ_CR = 0.1, M_s = 100, and B₀ directed along $\widehat{\mathbf{y}}$ such that v_A = c_s. The perturbations remain small before reaching the shock, while still growing exponentially at the nominal AI rate. Tracking the growth of each individual mode along the box becomes impractical, as it requires the use of sophisticated techniques such as wavelet transforms (Farge & Schneider 2015). Here, we adopted the simple Fourier transform (FT), whose formalism is outlined in Appendix B, sacrificing spatial resolution in order to obtain the power of each k mode averaged over the whole box.

Fig. 4. Linear-regime simulation results for initial density fluctuations following Eq. (23) with a flat power spectrum in a 2048 × 256 grid. Left: Perturbation profiles for (from top to bottom) density, velocity x and y components, and magnetic-field x and y components. Note that δux/cs has been multiplied by 1/3 and δuy/cs by two, so all quantities can be visualized using the same color-bar scale. Right: Dimensionless 2D (top) and omnidirectional (bottom) density fluctuations’ power spectra, k2Pδρ/ρ0(kx, ky) and kEδρ/ρ0(k), respectively.

The top right panel of Fig. 4 shows the 2D density power spectrum P_δρ/ρ0(k_x, k_y), multiplied by k² to make it dimensionless. For k ≲ 10³, the observed power spectrum resembles the expected AI growth rate. In other words, the amplitude of each (a, b) mode in Eq. (23) grows along the box as A_ab(t) = A_abexp[Γ(k_a, k_b) t], where Γ(k_a, k_b) is the AI growth rate obtained from Eq. (6). Numerical damping significantly affects the simulation results for k ≳ 10³, evidenced by a clear break at k ∼ 10³ in the time-averaged omnidirectional power spectrum E_δρ/ρ0(k) shown in the bottom right panel. This is an unavoidable and critical feature of the simulations, which was also pointed out by DD and VLS.

3.2. Nonlinear regime

We considered the nonlinear regime where large perturbations are present near the shock surface (δρ, δP, δu, δB ≳ ρ₀, P₀, u₀, B₀). Large δB perturbations are expected to efficiently scatter particles, thus increasing their maximum energy. We show that AI can lead to large amplitude fluctuations even when the perturbations are small at upstream infinity.

We employed the simulation framework described in Sect. 3.1.3, injecting a superposition of perturbations down to the grid resolution scale. While the initial growth of preexisting small perturbations at upstream infinity is due to the onset of AI in the CR precursor pressure gradient, their nonlinear evolution is connected with the onset of turbulence. A similar investigation was carried out by DD using 2D MHD simulations and by VLS using 3D simulations. Both DD and VLS adopted parameters aimed at maximizing MFA, but that were not viable on physical grounds; for instance, the efficiency of particle acceleration was ξ_CR = 0.6, which is now considered disproportionately large with respect to the results of hybrid simulations, as well as phenomenological considerations, suggesting ξ_CR ≃ 0.1. In addition, previous studies considered Mach numbers of M_s ∼ 100, while shocks at the beginning of the Sedov phase of a SNR are expected to have M_s ∼ 500 and be even larger at earlier times.

The adoption of more realistic values of these parameters has profound consequences for MFA and for the onset of turbulence. In Fig. 5, we show 2D profiles of ρ/ρ₀, u_x, y/c_s, B_x, y/B₀, and (B/B₀)² for our benchmark choice of parameters (ξ_CR = 0.1, M_s = 500, δρ_0, rms/ρ₀ = 0.1, and v_A = c_s; left column), as well as for the parameters adopted in previous literature (ξ_CR = 0.6, M_s = 100, δρ_0,rms/ρ₀ = 0.1, and v_A = c_s; right column). The main qualitative differences are as follows: 1) the left profiles contain fluctuating structures predominantly on smaller scales with respect to the right profiles; 2) the amplitude of (δB/B₀)² fluctuations is smaller on the left, indicating a lower level of MFA; and 3) the combination of a longer advection time and a larger ξ_CR allows for early eddy formation in the right panels. Meanwhile, in the left panels, eddies form only very close to the shock, and most of the precursor is dominated by small-scale shocklets. Kang et al. (1992) suggested that these shocklets can play a role in energizing particles, even before they cross the shock surface. In both cases, the morphology of δB_x is elongated along x, while δB_y is elongated along y since the latter grows mainly due to plasma compression at the location of the shocklets.

Fig. 5. Two-dimensional snapshots of quasi-stationary density, velocity, magnetic field, and magnetic energy density profiles, for different parameter choices, as indicated above each column, in a 4096 × 512 simulation grid. All quantities have been normalized so as to become dimensionless. For better contrast, color bars use linear scaling for ρ/ρ0 and ux/cs, logarithmic scaling for (B/B0)2, and symmetric logarithmic scaling (linear near zero and logarithmic at large amplitudes) for other quantities.

To quantitatively analyze our results, we calculated both kinetic and magnetic omnidirectional turbulent-energy-density spectra, E_turb, kin(k) and E_turb, mag(k), respectively, as defined in Appendix B. We first obtain the density-weighted velocity, $\mathbf{w}(x,y)=\sqrt{\rho (x,y)}\mathbf{u}(x,y)$ (as is customary in studies of compressible MHD turbulence; e.g., Kida & Orszag 1990; Grete et al. 2017), and the magnetic field, B(x, y), from the simulation output, and then subtracted their mean background profiles³ $\overline{\mathbf{w}}(x)$ and $\overline{\mathbf{B}}(x)$ to obtain the turbulent components:

$$\begin{array}{cc}\hfill \delta \mathbf{w}(x,y)& =\mathbf{w}(x,y)-\overline{\mathbf{w}}(x),\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill \delta \mathbf{B}(x,y)& =\mathbf{B}(x,y)-\overline{\mathbf{B}}(x).\hfill \end{array}$$(25)

The turbulent energy spectra are defined as the RMS energy density in turbulent fluctuations on scales between k and k + dk:

$$\begin{array}{cc}\hfill E_{\mathrm{turb},\mathrm{kin}}(k)dk& =\frac{\delta w_{\mathrm{rms}}^2(k)}{2},\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill E_{\mathrm{turb},\mathrm{mag}}(k)dk& =\frac{\delta B_{\mathrm{rms}}^2(k)}{8\pi }·\hfill \end{array}$$(27)

We also calculated the strength of MFA as a function of the distance from the shock surface:

$$\begin{array}{c}\hfill \frac{\delta B_{\mathrm{rms}}^2(x)}{B_0^2},\delta B_{\mathrm{rms}}^2(x)=\frac{1}{L/8}{\displaystyle {\int }_0^{L/8}}{|\delta \mathbf{B}(x,y)|}^{2}dy.\end{array}$$(28)

In Fig. 6, we show the turbulent power spectra and MFA for different choices of parameters (M_s, ξ_CR, v_A/c_s, δρ_0,rms/ρ₀, and magnetic-field orientation). In the top panels, we fix δρ_0,rms/ρ₀ = 0.1, v_A = c_s, and ${\mathbf{B}}_0\Vert \widehat{\mathbf{y}}$ while varying the remaining parameters. We compare the results for the benchmark values of the parameters (blue solid) with the DD/VLS parameters (dashed yellow) and with the more conservative case M_s = 100, ξ_CR = 0.1 (dotted green). It is worth recalling that for a typical SNR shock, the Mach number in the phase that is expected to be most important for particle acceleration is typically ∼100 − 1000 and the efficiency is ξ_CR ∼ 0.1. From the top-right panel we immediately see that the scenario with M_s = 100 and ξ_CR ∼ 0.1 does not yield any appreciable MFA; the background field still carries more energy than the fluctuations in this case. Note that in our benchmark case the turbulent magnetic-energy spectrum is very hard, with most of the power concentrated at small scales, which could translate into enhanced scattering of particles at the energies of interest (see discussion in Sect. 4). Meanwhile, strongly modified shocks with ξ_CR = 0.6 can easily produce δB_rms/B₀ ≳ 10, but their turbulent magnetic-energy spectra are somewhat steep, with most power concentrated on large scales. This is the only scenario in which we are anywhere close to reaching energy equipartition:

$$\begin{array}{c}\hfill E_{\mathrm{turb},\mathrm{kin}}(k)\approx E_{\mathrm{turb},\mathrm{mag}}(k),\mathrm{for}k>k_{\mathrm{eq}},\end{array}$$(29)

Fig. 6. Turbulent power spectra in units of ρ0cs2 (left) and MFA (right), calculated via Eq. (28), for different parameter values and a 4096 × 512 grid. Thick/thin lines in the left panels are kinetic/magnetic power spectra. The right panels’ legends indicate the colors/line styles corresponding to each choice of parameters. Top: We fix δρ0,rms/ρ0 = 0.1, vA = cs, and ${\mathbf{B}}_0\Vert \widehat{\mathbf{y}}$, and we vary Ms and ξCR. Middle: We fix Ms = 500, ξCR = 0.1, and ${\mathbf{B}}_0\Vert \widehat{\mathbf{y}}$, and we vary δρ0,rms/ρ0 and vA/cs. Bottom: We fix Ms = 500, ξCR = 0.1, δρ0,rms/ρ0 = 0.1, and vA = cs, and we vary the background field orientation.

where $k_{\mathrm{eq}}^{-1}$ is the equipartition length scale. For smaller ξ_CR, the Mach number must be large to obtain any meaningful MFA. However, increasing the Mach number also increases E_{turb,  kin} and k_eq. Ultimately, it is hard to estimate k_eq from simulations because this scale is never numerically resolved and E_{turb,  mag}(k) does not follow the ∼k^−1/2 scaling predicted by Beresnyak et al. (2009).

In the middle panels of Fig. 6 we instead fix M_s = 500 and ξ_CR = 0.1 and vary the remaining parameters while keeping ${\mathbf{B}}_0\Vert \widehat{\mathbf{y}}$. For the dashed yellow curves, we set v_A = 0.1c_s, while for the dotted green curves, we doubled the RMS amplitude of the initial fluctuations, δρ_0,rms/ρ₀ = 0.2. Inspecting the low-k region of the power spectra, where numerical damping is negligible, we find that increasing v_A by a factor of ten increases E_{turb,  mag} by a factor of ≈100. Similarly, doubling δρ_0,rms/ρ₀ makes E_{turb,  mag} larger by a factor of ≈4. This is consistent with the natural expectation that $E_{\mathrm{turb},\mathrm{mag}}\propto B_0^2{(\delta {\rho }_{0,\mathrm{rms}}/{\rho }_0)}^2$.

In the bottom panels we keep M_s, ξ, δρ_0,rms/ρ₀, and v_A/c_s fixed to their benchmark values to isolate the effect of field orientation. The results are somewhat consistent with the interpretation put forward by Downes & Drury (2014), the authors of which postulated that amplification occurs primarily to the component δB_y. If the shock is parallel (${\mathbf{B}}_0\Vert \widehat{\mathbf{x}}$, dotted green lines), magnetic fluctuations along y arise only after the system turns nonlinear. Instead, if ${\mathbf{B}}_0\Vert \widehat{\mathbf{y}}$ (solid blue lines), δB_y is immediately generated in AI eigenmodes. The latter is the most optimistic scenario for MFA. However, more generally, SNR shocks expand onto oblique fields and MFA saturates between the two extreme cases discussed above. As an example, we show the case in which B₀ is inclined by 30° with respect to the shock normal (dashed yellow lines). At low-k, its magnetic power spectrum is ∼4 times lower than the perpendicular shock case. This indicates that $\delta B_y\propto {\mathbf{B}}_0·\widehat{\mathbf{y}}=B_{0}sin30°=B_0/2$ indeed dominates the amplified field.

To meaningfully interpret these results, in Fig. 7 we show a numerical convergence test, adopting the benchmark choice of parameters. In the left panel, we show the turbulent kinetic (thick lines) and magnetic (thin lines) power spectra, for different grid resolutions, normalized to ρ₀c_s². The magnetic power peaks at larger k for higher resolution, but it is always far from reaching equipartition with the kinetic power. Note that doubling the grid resolution doubles the wavenumber, k_d, above which damping is important (k_dL ≈ 1.3 × 10³ for the 1024 × 128 resolution) and changes the shape of E_{turb,  mag}(k) below k_d, making it harder. This is not only due to the reduced numerical damping for k ≲ k_d, it also happens because less power can be generated and transferred from higher to smaller k via inverse cascading. In the right panel, we show (δB_rms/B₀)² as a function of x. As we can see, MFA increases with resolution, without converging. Our δB/B₀ of order a few can only be interpreted as a conservative lower limit, because most of the magnetic energy is expected to be located at damping-dominated scales of kL ≳ 2000. As we show in Sect. 4, these scales play a crucial role in particle scattering.

Fig. 7. Left: Turbulent kinetic (thick lines) and magnetic (thin lines) power spectra in units of ρ0cs2 for different grid resolutions (see right panel legend). Right: MFA as a function of x for different grid resolutions.

In this work, we focused on 2D simulations in order to access smaller scales, which are the most important because they carry most of the turbulent magnetic energy. Both DD and VLS have shown that simulations in 2D yield less power on small scales than in 3D. This is likely because turbulent kinetic energy predominantly undergoes an inverse cascade in 2D and a forward cascade in 3D. However, for this difference to be appreciable, turbulence should develop and have enough time to transfer kinetic energy between scales. Energy can safely cascade from k₁ → k₂ when the eddy turnover time t_eddy(k)∼1/k δu(k) at all intermediate wavenumbers k₁ < k < k₂ is much shorter than the advection time in the precursor:

$$\begin{array}{c}\hfill k\delta u_{\mathrm{rms}}(k)\gg t_{\mathrm{adv}}^{-1}=u_0/L,\mathrm{for}{k}_1<k<k_2.\end{array}$$(30)

The RMS amplitude of velocity fluctuations at a given k, δu_rms(k), is defined via E_δu(k) dk = ρ₀δu_rms²(k)/2, where the velocity power spectrum, E_δu, is defined in Appendix B. Note that E_δu is different from E_{turb,  kin}, which is defined in terms of $\mathbf{w}=\sqrt{\rho }\mathbf{u}$.

In Fig. 8, we compare power spectra (left panel), MFA (middle panel), and characteristic timescales (right panel) for our benchmark set of parameters and 1024 × 128 (×128) resolution, in 2D (blue solid lines) and 3D (yellow dashed lines) simulations, respectively. Indeed, we find slightly more kinetic power at higher k in the 3D scenario, as well as comparable MFA. However, numerical damping once again prevents us from resolving the scales in which this difference could become appreciable. By comparing the eddy turnover times in the right panel with the advection time (dotted green line), we see that condition (30) is only marginally satisfied, and significant energy cascading is not expected to take place. In contrast, using the overly optimistic ξ_CR = 0.6 and M_s = 100 adopted in DD and VLS, one expects stronger cascading simply because Γ t_adv is larger in the kL ≪ ξ_CRM_ms² limit, thus reaching the nonlinear/turbulent regime earlier. In summary, while 2D simulations allowed us to access perturbations with large wavenumbers, 3D simulations may provide a better description of the nonlinear phase if small scales are resolved. On the other hand, it is clear that accessing such scales in 3D simulations quickly becomes unbearably expensive from the computational point of view.

Fig. 8. Left: Turbulent kinetic (thick lines) and magnetic (thin lines) power spectra in units of ρ0cs2 for 2D (1024 × 128 grid; blue solid) and 3D (1024 × 128 × 128 grid; dashed yellow) simulations. Center: MFA as a function of x for 2D and 3D simulations. In 3D, we performed a yz averaging rather than just a y averaging. Right: Inverse eddy turnover timescales for the 2D (solid blue) and 3D (dashed yellow) runs, compared to the inverse advection timescale $t_{\mathrm{adv}}^{-1}=u_0/L$ (dotted green), in units of cs/L.

4. Competition and cooperation between acoustic and nonresonant instability

Particle acceleration at SNR shocks must be accompanied by substantial MFA, not only to reach a higher value of the maximum energy, but also to explain the thin filaments observed in X-rays (Vink 2012). It is, however, important to keep in mind that while the former point requires MFA upstream of the shock, namely in the CR precursor, to reduce the acceleration time, the latter only requires magnetic-field enhancement behind the shock.

All processes that can in principle lead to MFA are associated with a gradient in the CR distribution upstream: in the case of the nonresonant instability (Bell 2004), it is the current of particles escaping the upstream region that leads to the growth of the instability. In general, this phenomenon occurs on scales larger than the precursor size, which is L ∼ D(E_max)/v_sh, where D is the diffusion coefficient and E_max is the maximum energy of the accelerated particles (for convenience, we use v_sh rather than u₀ in this section). Such current is directly related to the gradient of CR number density, and the process may lead to δB/B ≫ 1. For the resonant streaming instability (Lagage & Cesarsky 1983), it is the gradient of the CR density upstream to produce Alfvén waves at wavenumbers resonant with the CR momentum (r_L(E)≃k⁻¹). This latter process generally saturates at δB/B ≲ 1, especially if some type of damping is effective. Finally, AI discussed in this article grows due to the dynamical action of the CR pressure gradient upstream of the shock on density fluctuations present in the ISM and can lead to δB/B ≫ 1, depending on whether the instability has enough time to grow in one advection time. Below we comment on the location where these instabilities can be excited, on the scales involved and which one is expected to dominate as a function of the parameters.

If, for simplicity, we assume that the spectrum of accelerated particles is a power law in momentum space, f(p)∝p⁻⁴, which is appropriate for strong shocks, then the condition for the growth of the nonresonant instability can be written as

$$\begin{array}{c}\hfill \frac{3{\xi }_{\mathrm{CR}}{\rho }_0{v}_{\mathrm{sh}}^2}{\mathrm{\Lambda }}\frac{v_{\mathrm{sh}}}{c}\ge \frac{B_0^2}{4\pi },\end{array}$$(31)

where Λ ≈ ln(E_max/m_pc²)∼10. This condition is equivalent to requiring that the energy density carried by the current is larger than the energy density in the preexisting magnetic field (Bell 2004; Amato & Blasi 2009). If one introduces the Alfvénic Mach number in the preexisting field, then this condition can be rewritten as

$$\begin{array}{c}\hfill M_A\ge {\left(\frac{c\mathrm{\Lambda }}{3{\xi }_{\mathrm{CR}}{v}_A}\right)}^{1/3}\approx 100{\left(\frac{{\xi }_{\mathrm{CR}}}{0.1}\right)}^{-1/3}{\left(\frac{v_A}{10\mathrm{km}/\mathrm{s}}\right)}^{-1/3}.\end{array}$$(32)

Typically, when the shock velocity drops below ∼1000 km/s, the nonresonant instability stops being excited and one has only the resonant streaming instability to rely upon.

At saturation, the magnetic field, B_sat, can be estimated from Eq. (31) by replacing B₀ with B_sat and interpreting it as an equality. Following Schure & Bell (2013), Cristofari et al. (2020), the maximum energy is obtained by requiring that the instability grows for n e-folds, where typically it is assumed that n ≈ 5, namely γ_maxt ≈ 5, where γ_max is the maximum growth rate of the instability. When the nonresonant instability is excited, it is therefore driven by the current of escaping particles over a region with size $L_{\mathit{MFA}}\sim 5c{\gamma }_{\max }^{-1}$, where γ_max = k_maxv_A and

$$\begin{array}{c}\hfill k_{\max }=\frac{4\pi }{c{B}_0}{J}_{\mathrm{CR}}(>E)=\frac{4\pi }{c{B}_0}\frac{3{\xi }_{\mathrm{CR}}{\rho }_0{v}_{\mathrm{sh}}^2}{\mathrm{\Lambda }E}e{v}_{\mathrm{sh}}\end{array}$$(33)

is the scale where the growth is the fastest (Bell 2004). Here, J_CR(> E) is the current carried by particles with energy > E. It can be easily shown that the condition in Eq. (31) is also equivalent to requiring k_max ≥ r_L(E)⁻¹, where r_L(E) is the Larmor radius of the particles with energy E. The region affected by MFA, L_MFA, is typically much larger than the size of the precursor, and the amplification of the magnetic field occurs on such a larger region. It is also important, for the considerations below, to keep in mind that the instability also results in overdensities of δρ/ρ₀ ≳ 1 in the same region where the field is amplified.

At saturation, the nonresonant instability produces magnetic structures on the scale of the Larmor radius of the escaping particles in the amplified magnetic field. Using the considerations above, one finds that at the highest energy, E_max, the Larmor radius in the amplified field is

$$\begin{array}{c}\hfill r_L(E_{\max })=\frac{1}{5}{\left(\frac{3{\xi }_{\mathrm{CR}}}{\mathrm{\Lambda }}\right)}^{1/2}{\left(\frac{v_{\mathrm{sh}}}{c}\right)}^{1/2}{v}_{\mathrm{sh}}t.\end{array}$$(34)

At this point, the size of the precursor would be of the order of L ≃ D(E_max)/v_sh, corresponding to a scale of $L{r}_{\mathrm{L}}^{-1}\simeq c/3v_{\mathrm{sh}}$ in the Bohm limit. For v_sh = 10³ − 10⁴ km/s this corresponds to $L{r}_{\mathrm{L}}^{-1}\simeq 100-10$. Notice that these latter considerations apply to the nonresonant instability at saturation, while the wavenumber where the instability initially (in the linear regime) grows the fastest, k_max, corresponds to much smaller scales.

It is also useful to estimate L/R_sh, where R_sh is the radius of the shock, which is obtained by approximating v_sht ≈ R_sh:

$$\begin{array}{c}\hfill \frac{L}{R_{\mathrm{sh}}}\approx \frac{1}{15}{\left(\frac{3{\xi }_{\mathrm{CR}}}{\mathrm{\Lambda }}\right)}^{1/2}{\left(\frac{v_{\mathrm{sh}}}{c}\right)}^{-1/2}.\end{array}$$(35)

For v_sh = 10³ − 10⁴ km/s this corresponds to $L{R}_{\mathrm{sh}}^{-1}\simeq 0.2-0.06$, which is consistent with the typical approximation that D(E_max)/v_sh ∼ 0.1R_sh.

The simple estimates presented here are formally valid for a parallel shock, although Bell (2005) showed that similar results are obtained for oblique shocks. Hence, we adopted the estimate of L provided above for oblique shocks as well. For shocks where the inclination to the shock normal is ≳45°, injection has been shown to be suppressed (Caprioli & Spitkovsky 2014), and they are of lesser interest here.

Let us now estimate these same quantities for AI. The fastest growing modes grow at a rate of ∼ξ_CRM_msv_sh/2L, which is larger than the nonresonant growth rate at E_max when M_ms ≳ (10/ξ_CR)(L/R_sh)∼10, where the last inequality holds if we adopt ξ_CR and L/R_sh ∼ 0.1. Based on these considerations, AI seems to grow faster for the high Mach number shocks that are typical of young SNRs.

As discussed above, the saturation of the instability is not easily accessible to current simulations; it requires the small-scale magnetic modes to reach some sort of equipartition with kinetic perturbations and initiate a nonlinear cascade toward larger scales. If trusting the qualitative expectations of DD, one could suggest that at saturation

$$\begin{array}{c}\hfill \frac{\delta B^2}{8\pi }\approx {\xi }_{\mathrm{CR}}^2{\rho }_0{v}_{\mathrm{sh}}^2{\left(\frac{\delta \rho }{{\rho }_0}\right)}^2,\end{array}$$(36)

and assuming δρ/ρ₀ ∼ 1 (which is the result of the nonlinear growth of the instability) we obtain the following ratio between the saturation values of the acoustic/turbulent to nonresonant instability:

$$\begin{array}{c}\hfill \frac{{(\delta B^2)}_{\mathrm{ai}}}{{(\delta B^2)}_{\mathrm{nr}}}\approx \frac{2{\xi }_{\mathrm{CR}}\mathrm{\Lambda }}{3}\frac{c}{v_{\mathrm{sh}}}·\end{array}$$(37)

For typical values of the parameters, this ratio exceeds unity, which suggests that AI may prevail. As stressed above, numerical simulations of AI are not yet suitable to really test this regime and check the viability of Eq. (36).

We want to stress that the two CR-induced instabilities studied here may operate at the same time and in fact cooperate rather than compete with each other: the fact that the Bell instability operates due to the current of particles escaping the acceleration region implies that this instability produces density perturbations of δρ/ρ₀ ≳ 1 upstream of a shock, over a region much larger than the size of the precursor. When these regions enter the shock precursor, they can undergo additional MFA due to the onset of AI and its nonlinear counterpart. However, this conclusion requires some caution; AI is expected to amplify the perpendicular component of the magnetic field, while, as discussed above, MFA is reduced in the case of parallel shocks. On the other hand, it is well known that the nonresonant instability creates large perpendicular magnetic fields away from the shock. Hence, by the time that these perturbations enter the precursor they can further be amplified by AI and the considerations above would apply.

While this synergic relation between the two processes can only be proven with dedicated simulations in which CRs are active components (hybrid or MHD+PIC simulations), it is worth stressing that hybrid simulations have already shown evidence for episodes in which particle trapping in large magnetic structures occurs upstream before these structures are eventually advected downstream (Zeković et al. 2025). It is possible that such structures may be the result of overdense regions being unstable in the CR pressure gradient.

5. Conclusions

Cosmic-ray-induced MFA is an integral part of the process of particle acceleration at shocks, especially in the case of SNR shocks. In this work we investigated the possibility that the amplification may be due to the reaction of small density perturbations to the CR pressure gradient in the upstream region of the shock, as originally suggested by DD. We expanded on previous investigations carried out using MHD simulations with an assigned pressure gradient (DD, VLS) by improving resolution, deepening the discussion around numerical damping and the role of dimensionality, and adopting parameter values that are closer to the ones expected for a realistic SNR shock.

Our conclusions can be summarized as follows. 1) AI can lead to the nonlinear growth of small perturbations even for initial perturbations of δρ/ρ₀ ∼ 10⁻¹. This phenomenon leads to considerable MFA, especially when the original field is perpendicular to the shock normal. In contrast, previous works started with large initial perturbations, so that the role of AI was somehow reduced and nonlinear effects started immediately. We provided expressions for the growth rate and the scales involved. 2) The detailed study of the dependence of our results on resolution showed that by increasing the resolution, more power appears at large wavenumbers, as expected based on AI, but even the highest resolution runs were incapable of properly resolving the small scales where equipartition of kinetic and magnetic power is expected, for the realistic values of ξ_CR and Mach number, M_s, adopted here. This consideration applies even more to previous investigations. As a consequence, we were only able to impose a lower limit to the magnitude of MFA. 3) When AI leads to δρ/ρ₀ ≫ 1, we occasionally observe the formation of shocklets that temporarily stall and are later advected with the fluid across the shock. These structures look similar to the SLAMS known to the space plasma community (Schwartz et al. 1992) and may play an important role in particle acceleration. 4) For values of the Mach number appropriate for SNRs, we showed that AI may grow faster than the nonresonant instability and lead to a larger MFA. We comment on the possibility that Bell instability may cooperate by first creating regions with δρ/ρ₀ ≳ 1 far upstream, and that these regions may induce further amplification due to the CR pressure gradient. Future investigation, possibly using an MHD+PIC approach, is planned to better study the nonlinear reaction of accelerated particles to this instability.

Acknowledgments

The authors are grateful to an anonymous referee for their comments that helped us improve the manuscript. They also acknowledge precious conversations with E. Amato, V. Berta, N. Bucciantini, D. Caprioli, B. Olmi. They are also thankful to L. Drury and T.P. Downes for clarifications on their earlier work. The work of PB and AC was partially funded by the European Union – Next Generation EU under PRIN-MUR 2022TJW4EJ “Unveiling the footprints of the cosmic ray journey through the Galaxy and beyond”. The work of PB was also partially funded under MUR National Innovation Ecosystem grant ECS00000041 – VITALITY/ASTRA – CUP D13C21000430001. The work of ES was supported by a Rita Levi Montalcini fellowship.

References

Appendix A: Velocity perturbations for k_y modes

In this Appendix, we derive Eq. (20). We consider small density perturbations along y,

$$\begin{array}{c}\hfill \rho (x,y)=\frac{{\rho }_0+\delta \rho (x)sin(k_{y}y)}{U(x)}.\end{array}$$(A.1)

We assume that the density perturbations are nearly independent of x, namely δρ ≃ δρ₀ (see Fig. 3). The perturbed magnetic field is $\mathbf{B}=B_{y,0}{U}^{-1}\widehat{\mathbf{y}}+\delta \mathbf{B}$, and the perturbed velocity is $\mathbf{u}=u_{0}U(x)\widehat{\mathbf{x}}+\delta \mathbf{u}$.

Since ξ_CR ≪ 1, the derivatives of U can be neglected with respect to the derivatives of the perturbations. Using ∇ ⋅ δB = 0, the steady-state induction equation, ∇ × (u × B) = 0, can be approximated as

$$\begin{array}{cc}& \frac{\partial \delta B_y}{\partial x}=-\frac{B_{y,0}}{u_0{U}^2}\frac{\partial \delta u_x}{\partial x},\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& \frac{\partial \delta B_x}{\partial x}=\frac{B_{y,0}}{u_0{U}^2}\frac{\partial \delta u_x}{\partial y}.\hfill \end{array}$$(A.3)

Assuming v_A, y² = B_y, 0²/4πρ₀ ≪ u₀², and using Eq. (A.2), the momentum equation can be approximated as

$$\begin{array}{c}\hfill u_{0}U\frac{\partial \delta u_x}{\partial x}=U\frac{\mathrm{\nabla }{P}_{\mathrm{CR}}}{{\rho }_0}\frac{\delta {\rho }_0}{{\rho }_0}sin(k_{y}y)+\frac{B_{y,0}}{4\pi {\rho }_0}\frac{\partial \delta B_x}{\partial y},\end{array}$$(A.4)

Making the ansatz δu_x = a(x)sin(k_yy), Eq. (A.3) gives

$$\begin{array}{c}\hfill \frac{\partial \delta B_x}{\partial x}=a(x)\frac{k_y{B}_{y,0}}{u_0{U}^2}cos(k_{y}y).\end{array}$$(A.5)

Taking the partial derivative of (A.4) with respect to x, we get

$$\begin{array}{c}\hfill a^{″}(x)=-\frac{k_y^2{B}_{y,0}^2}{4\pi {\rho }_0{u}_0^2{U}^3}a(x).\end{array}$$(A.6)

Solving Eq. (A.6) with BCs a(0) = 0 and a′(0) = (1/u₀)(∇P_CR/ρ₀)(δρ₀/ρ₀) yields

$$\begin{array}{c}\hfill \delta u_x(x,y)=\frac{1}{u_0}\frac{\mathrm{\nabla }{P}_{\mathrm{CR}}}{{\rho }_0}\frac{\delta {\rho }_0}{{\rho }_0}\frac{sin(k_x^{\prime }x)}{k_x^{\prime }}sin(k_{y}y),\end{array}$$(A.7)

where k_x′=k_y(v_A/u₀)U^−3/2.

Appendix B: Fourier transforms and power spectra

This Appendix outlines the formalism behind the discrete FT (DFT) and the power spectrum definitions used in Sect. 3. Below, f(x) denotes a generic field, representing density or a single velocity/magnetic field component. When specifically dealing with vector field components, we use Greek letter subscripts, f = (f_α).

Our 2D simulation box has area A = L²/8 embedded within a grid of N_tot = N_xN_y points with uniform spacings along each direction, Δx = L/N_x and Δy = L/8N_y, and spatial coordinates

$$\begin{array}{c}\hfill {\mathbf{x}}_{\mathit{ij}}=(x_i,y_j)=(i\mathrm{\Delta }x,j\mathrm{\Delta }y),\{\begin{array}{c}i=0,...,N_x-1\hfill \\ j=0,...,N_y-1\hfill \end{array}.\end{array}$$(B.1)

This grid allows for specific wavenumbers along each direction:

$$\begin{array}{c}\hfill {\mathbf{k}}_{\mathit{ab}}=(k_a,k_b)=\frac{2\pi }{L}(a,8b),\{\begin{array}{c}a=-\frac{N_x}{2}+1,...,\frac{N_x}{2}\hfill \\ b=-\frac{N_y}{2}+1,...,\frac{N_y}{2}\hfill \end{array}.\end{array}$$(B.2)

Grid cell areas are Δ²x = ΔxΔy and Δ²k = (2π)²/A in physical and reciprocal spaces, respectively. Note that, along each direction, the maximum allowed k corresponds to the Nyquist limit, k_x, max = 2π/2Δx and k_y, max = 2π/2Δy, while the minimum allowed k (besides the zero mode) corresponds to a wavelength the size of each box dimension, k_x, min = 2π/L and k_y, min = 2π/(L/8) = 16π/L.

We adopt the convention relating a two-dimensional field f(x) to its FT $\stackrel{\sim }{f}(\mathbf{k})$ via

$$\begin{array}{cc}& \stackrel{\sim }{f}(\mathbf{k})=\frac{1}{{(2\pi )}^2}{\displaystyle \int f}(\mathbf{x})e^{-i\mathbf{k}·\mathbf{x}}{d}^2\mathbf{x},\hfill \end{array}$$(B.3)

$$\begin{array}{cc}& f(\mathbf{x})={\displaystyle \int \stackrel{\sim }{f}}(\mathbf{k})e^{i\mathbf{k}·\mathbf{x}}{d}^2\mathbf{k}.\hfill \end{array}$$(B.4)

Upon discretization, this yields the DFT

$$\begin{array}{c}\hfill {\stackrel{\sim }{f}}_{\mathit{ab}}=\frac{1}{N_{\mathrm{tot}}}{\displaystyle \sum _{\mathit{ij}}}{f}_{\mathit{ij}}{e}^{-i{\mathbf{k}}_{\mathit{ab}}·{\mathbf{x}}_{\mathit{ij}}},f_{\mathit{ij}}={\displaystyle \sum _{\mathit{ab}}}{\stackrel{\sim }{f}}_{\mathit{ab}}{e}^{i{\mathbf{k}}_{\mathit{ab}}·{\mathbf{x}}_{\mathit{ij}}},\end{array}$$(B.5)

where one defines f_ij ≡ f(x_ij) and ${\stackrel{\sim }{f}}_{\mathit{ab}}\equiv {\mathrm{\Delta }}^{2}k\stackrel{\sim }{f}({\mathbf{k}}_{\mathit{ab}})$, obeying Parseval’s identity

$$\begin{array}{c}\hfill {\displaystyle \sum _{\mathit{ij}}}|f_{\mathit{ij}}{|}^2=N{\displaystyle \sum _{\mathit{ab}}}{|{\stackrel{\sim }{f}}_{\mathit{ab}}|}^2.\end{array}$$(B.6)

These DFTs can be efficiently calculated using NumPy’s ‘fft’ package (Harris et al. 2020) under the ‘forward’ normalization convention.

Under statistical homogeneity, the power spectrum tensor P_αβ(k_ab) is defined for vector fields via the ensemble average

$$\begin{array}{c}\hfill {\mathrm{\Delta }}^{2}k⟨{\stackrel{\sim }{f}}_{\alpha }({\mathbf{k}}_{\mathit{ab}}){\stackrel{\sim }{f}}_{\beta }^{\ast }({\mathbf{k}}_{a^{\prime }{b}^{\prime }})⟩={\delta }_{a{a}^{\prime }}{\delta }_{b{b}^{\prime }}{P}_{\alpha \beta }({\mathbf{k}}_{\mathit{ab}}),\end{array}$$(B.7)

or, equivalently, as the FT of the correlation tensor

$$\begin{array}{c}\hfill {\xi }_{\alpha \beta }({\mathbf{r}}_{\mathit{ij}})=⟨f_{\alpha }(\mathbf{x})f_{\beta }^{\ast }(\mathbf{x}-{\mathbf{r}}_{\mathit{ij}})⟩.\end{array}$$(B.8)

Scalar fields yield a scalar power spectrum, P(k_ab), with an equivalent definition dropping the component indices. Rigorously performing ensemble averages requires many simulations, which is too computationally expensive. The best one can do from a single simulation is compute the estimator ${\widehat{P}}_{\alpha \beta }({\mathbf{k}}_{\mathit{ab}})={\mathrm{\Delta }}^{2}k{\stackrel{\sim }{f}}_{\alpha }({\mathbf{k}}_{\mathit{ab}}){\stackrel{\sim }{f}}_{\beta }^{\ast }({\mathbf{k}}_{\mathit{ab}})={\stackrel{\sim }{f}}_{\alpha ,ab}{\stackrel{\sim }{f}}_{\beta ,ab}^{\ast }/{\mathrm{\Delta }}^{2}k$, which is the FT of ${\widehat{\xi }}_{\alpha \beta }({\mathbf{r}}_{\mathit{ij}})=(1/N){\displaystyle \sum _{\mathit{kl}}}{f}_{\alpha }({\mathbf{x}}_{\mathit{kl}})f_{\beta }^{\ast }({\mathbf{x}}_{\mathit{kl}}-{\mathbf{r}}_{\mathit{ij}})$. Note that its trace at zero separation is

$$\begin{array}{c}\hfill {\widehat{\xi }}_{\alpha \alpha }(\mathbf{0})=\frac{1}{N}{\displaystyle \sum _{\mathit{ij}}}{|\mathbf{f}({\mathbf{x}}_{\mathit{ij}})|}^2={\displaystyle \sum _{\mathit{ab}}}{\mathrm{\Delta }}^{2}k{\widehat{P}}_{\alpha \alpha }({\mathbf{k}}_{\mathit{ab}}).\end{array}$$(B.9)

We shall drop the hats from now on, always working with the estimators. For a statistically homogeneous and isotropic two-dimensional field satisfying ∇ ⋅ f = 0, one can write (Batchelor 1953)

$$\begin{array}{c}\hfill P_{\alpha \beta }(\mathbf{k})=P(k)({\delta }_{\alpha \beta }-\frac{k_{\alpha }{k}_{\beta }}{k^2}),\end{array}$$(B.10)

where $P(k)=P_{\alpha \alpha }(\mathbf{k})={\mathrm{\Delta }}^{2}k{|\stackrel{\sim }{\mathbf{f}}(\mathbf{k})|}^2$ is the tensor’s trace, which depends only on k = |k|. For a scalar field, isotropy also implies that P(k) = P(k).

From P(k_ab) computed for either a scalar or a vector field f, one can define an omnidirectional power spectrum E(k) (also commonly called an “energy spectrum”) by binning k and adding all the power within each bin k_ab = |k_ab|∈[k_n, k_n + 1 = k_n + Δk_n):

$$\begin{array}{c}\hfill E({\overline{k}}_n)\mathrm{\Delta }{k}_n=C{\displaystyle \sum _{\begin{array}{c}\mathit{ab}\\ k_{\mathit{ab}}\in [k_n,k_{n+1})\end{array}}}{\mathrm{\Delta }}^{2}kP(k_{\mathit{ab}}),\end{array}$$(B.11)

where ${\overline{k}}_n=\sqrt{k_n{k}_{n+1}}$ is a representative k in bin n (i.e. the geometric mean of its edges) and C is a constant prefactor depending on the specific field. This prefactor is introduced such that ${\displaystyle \sum _n}E({\overline{k}}_n)\mathrm{\Delta }{k}_n$ is an energy density, or some other quantity of interest, with $E({\overline{k}}_n)\mathrm{\Delta }{k}_n$ being the contribution coming from bandwidth [k_n, k_n + Δk_n]. In this work, we consider 4 types of omnidirectional power spectra:

1. For density fluctuations, the field is f = δρ/ρ₀, and we choose C = 1 such that

   $$\begin{array}{cc}\hfill {\displaystyle \sum _n}& E_{\delta \rho /{\rho }_0}({\overline{k}}_n)\mathrm{\Delta }{k}_n={\displaystyle \sum _{\mathit{ab}}}{\mathrm{\Delta }}^{2}kP(k_{\mathit{ab}})=\xi (\mathbf{0})\hfill \\ \hfill & =\frac{1}{A}{\displaystyle \sum _{\mathit{ij}}}\mathrm{\Delta }x\mathrm{\Delta }y{\left[\frac{\delta \rho }{{\rho }_0}({\mathbf{x}}_{\mathit{ij}})\right]}^2={\left(\frac{\delta {\rho }_{\mathit{rms}}}{{\rho }_0}\right)}^2,\hfill \end{array}$$(B.12)

   where δρ_rms/ρ₀ is the RMS δρ/ρ₀ over all scales in the box.

2. For magnetic fields fluctuations, f = δB, we choose C = 1/8π whereby the turbulent magnetic energy density power spectrum E_δB(k)≡E_turb, mag(k) satisfies

   $$\begin{array}{c}\hfill {\displaystyle \sum _n}{E}_{\mathrm{turb},\mathrm{mag}}({\overline{k}}_n)\mathrm{\Delta }{k}_n=\frac{1}{A}{\displaystyle \sum _{\mathit{ij}}}\mathrm{\Delta }x\mathrm{\Delta }y\frac{|\delta \mathbf{B}({\mathbf{x}}_{\mathit{ij}}){|}^2}{8\pi }=\frac{\delta B_{\mathit{rms}}^2}{8\pi }\end{array}$$(B.13)

   is the average magnetic energy density in perturbations at all scales, while $E_{\mathrm{turb},\mathrm{mag}}({\overline{k}}_n)\mathrm{\Delta }{k}_n=\delta B_{\mathit{rms}}^2({\overline{k}}_n)/8\pi$ is that coming from perturbations with wavenumbers k ∈ [k_n, k_n + Δk_n]. In other words,

   $$\begin{array}{c}\hfill \delta B_{\mathit{rms}}({\overline{k}}_n)=\sqrt{8\pi E_{\mathrm{turb},\mathrm{mag}}({\overline{k}}_n)\mathrm{\Delta }{k}_n}\end{array}$$(B.14)

   is the RMS value of magnetic field fluctuations with characteristic scale $\ell =2\pi /{\overline{k}}_n$.

3. For velocity fields, e.g.δu, we choose C = ρ₀/2 such that

   $$\begin{array}{c}\hfill {\displaystyle \sum _n}{E}_{\delta u}({\overline{k}}_n)\mathrm{\Delta }{k}_n=\frac{1}{A}{\displaystyle \sum _{\mathit{ij}}}\mathrm{\Delta }x\mathrm{\Delta }y\frac{{\rho }_0{|\delta \mathbf{u}({\mathbf{x}}_{\mathit{ij}})|}^2}{2}=\frac{{\rho }_0\delta u_{\mathit{rms}}^2}{2}.\end{array}$$(B.15)

   Analogously to the magnetic field case, the RMS amplitude of velocity perturbations at a characteristic scale $\ell =2\pi /{\overline{k}}_n$ can be found through

   $$\begin{array}{c}\hfill \delta u_{\mathit{rms}}({\overline{k}}_n)=\sqrt{\frac{2E_{\delta u}({\overline{k}}_n)\mathrm{\Delta }{k}_n}{{\rho }_0}}.\end{array}$$(B.16)

4. For studies of compressible turbulence, it is common to consider the density-weighted velocity $\mathbf{w}=\sqrt{\rho }\mathbf{u}$ to define kinetic energy spectra (Kida & Orszag 1990; Grete et al. 2017). Let δw be its fluctuations around a mean background profile. We can then choose C = 1/2 such that its corresponding power spectrum E_δw(k)≡E_turb, kin(k) is the turbulent kinetic energy spectrum:

   $$\begin{array}{c}\hfill {\displaystyle \sum _n}{E}_{\mathrm{turb},\mathrm{kin}}({\overline{k}}_n)\mathrm{\Delta }{k}_n=\frac{1}{A}{\displaystyle \sum _{\mathit{ij}}}\mathrm{\Delta }x\mathrm{\Delta }y\frac{|\delta w({\mathbf{x}}_{\mathit{ij}}){|}^2}{2}.\end{array}$$(B.17)

In practice, we calculate each E(k) by defining 50 logarithmically spaced k-bins between k_min = max(k_x, min, k_y, min) and k_max = min(k_x, max, k_y, max), and add their power following Eq. (B.11). This leads to

$$\begin{array}{c}\hfill E(k)\mathrm{\Delta }k\approx kE(k)\mathrm{\Delta }lnk=kE(k)\frac{ln(k_{\mathit{max}}/k_{\mathit{min}})}{50},\end{array}$$(B.18)

with ln(k_max/k_min)/50 ≈ 0.1 in our simulations. This means that δB_rms(k) and δu_rms(k) can be estimated, within factors of order unity, from the logarithmic energy spectrum kE(k) via:

$$\begin{array}{cc}\hfill \delta B_{\mathit{rms}}(k)& \sim \sqrt{k{E}_{\delta B}(k)},\hfill \end{array}$$(B.19)

$$\begin{array}{cc}\hfill \delta u_{\mathit{rms}}(k)& \sim \sqrt{k{E}_{\delta u}(k)}.\hfill \end{array}$$(B.20)

Generating Initial Density Perturbations

In Sects. 3.1.3 and 3.2, we inject a density perturbations field δρ₀(x)/ρ₀ which is superposition of modes with different wave vectors, according to Eq. (23), following a predetermined omnidirectional power spectrum (with C = 1) between max(k_x, min, k_y, min) < k < min(k_x, max, k_y, max),

$$\begin{array}{c}\hfill E_{\delta \rho /{\rho }_0}(k)=E^{\ast }{\left(\frac{k}{k^{\ast }}\right)}^{\alpha }\end{array}$$(B.21)

and zero outside this range, where E^* is a normalization setting the power in fluctuations at a reference scale 2π/k^*. We fix k^* = min(k_x, max, k_y, max). Using the formalism presented above, we now describe our method for generating such a field⁴. This is done mode by mode in Fourier space, rather than directly in physical space, since it makes the connection with the desired power spectrum more direct.

We consider density perturbations to be Gaussian in nature, meaning that $\delta {\stackrel{\sim }{\rho }}_{\mathit{ab}}/{\rho }_0={\mathrm{\Delta }}^{2}k\delta \stackrel{\sim }{\rho }({\mathbf{k}}_{\mathit{ab}})/{\rho }_0=X_{\mathit{ab}}+i{Y}_{\mathit{ab}}$, with X_ab and Y_ab randomly sampled from a Gaussian distribution with zero mean and variance σ_ab². Writing this in polar form,

$$\begin{array}{c}\hfill \frac{\delta {\stackrel{\sim }{\rho }}_{\mathit{ab}}}{{\rho }_0}=A_{\mathit{ab}}{e}^{i{\varphi }_{\mathit{ab}}},\end{array}$$(B.22)

it can be easily shown that A_ab follows a Rayleigh distribution with scale parameter σ_ab, implying $A_{ab,rms}=\sqrt{2}{\sigma }_{\mathit{ab}}$, while ϕ_ab is uniformly distributed between [0, 2π). Recall that producing a real δρ(x_ij)/ρ₀ requires $\delta \stackrel{\sim }{\rho }(-{\mathbf{k}}_{\mathit{ab}})=\delta {\stackrel{\sim }{\rho }}^{\ast }({\mathbf{k}}_{\mathit{ab}})$. This allows us to generate only half of the k-modes; once we have $\delta {\stackrel{\sim }{\rho }}_{\mathit{ab}}/{\rho }_0$ for all a, b ≥ 0, the modes with a, b < 0 are all set by the reality condition⁵: A_−a, −b = A_ab and ϕ_−a, −b = −ϕ_ab.

The power spectrum resulting from Eq. (B.22) is P(k_ab) = A_ab²/Δ²k, while its inverse FT yields

$$\begin{array}{c}\hfill \frac{\delta \rho }{{\rho }_0}({\mathbf{x}}_{\mathit{ij}})={\displaystyle \sum _{\mathit{ab}}}{A}_{\mathit{ab}}cos({\mathbf{k}}_{\mathit{ab}}·{\mathbf{x}}_{\mathit{ij}}+{\varphi }_{\mathit{ab}}),\end{array}$$(B.23)

in agreement with Eq. (23). The RMS density perturbation over all scales can be be shown to match our expectation:

$$\begin{array}{c}\hfill {\left(\frac{\delta {\rho }_{\mathrm{rms}}}{{\rho }_0}\right)}^2=\frac{1}{N}{\displaystyle \sum _{\mathit{ij}}}{\left[\frac{\delta \rho }{{\rho }_0}({\mathbf{x}}_{\mathit{ij}})\right]}^2={\displaystyle \sum _{\mathit{ab}}}{A}_{\mathit{ab}}^2=E^{\ast }{\displaystyle \sum _n}{\left(\frac{{\overline{k}}_n}{k^{\ast }}\right)}^{\alpha }\mathrm{\Delta }{k}_n.\end{array}$$(B.24)

Crucially, Eq. (B.24) allows us to find E^* for any desired value of δρ_rms/ρ₀. Meanwhile, the RMS value of perturbations on a scale k_ab is simply A_ab, rms. The power inside each k-bin is

$$\begin{array}{c}\hfill E_{\delta \rho /{\rho }_0}({\overline{k}}_n)\mathrm{\Delta }{k}_n={\displaystyle \sum _{\begin{array}{c}\mathit{ab}\\ k_{\mathit{ab}}\in [k_n,k_{n+1})\end{array}}}{A}_{\mathit{ab}}^2=N_n\overline{A_{\mathit{ab}}^2},\end{array}$$(B.25)

where $N_n\simeq 2\pi {\overline{k}}_n\mathrm{\Delta }{k}_n/{\mathrm{\Delta }}^{2}k$ is the number of grid points in the bin and $\overline{A_{\mathit{ab}}^2}$ is the bin-averaged power per mode, which we set to A_ab, rms. In the Δk_n → 0 limit, this approaches

$$\begin{array}{c}\hfill A_{ab,\mathrm{rms}}=\sqrt{2}{\sigma }_{\mathit{ab}}=\sqrt{\frac{E_{\delta \rho /{\rho }_0}(k_{\mathit{ab}})}{2\pi k_{\mathit{ab}}}{\mathrm{\Delta }}^{2}k}.\end{array}$$(B.26)

This relation determines σ_ab used to randomly sample each A_ab mode, guaranteeing that the resulting spectrum will be properly normalized and follow the desired scaling with k. Once we have sampled every $\delta \stackrel{\sim }{\rho }({\mathbf{k}}_{\mathit{ab}})/{\rho }_0$, we simply take its inverse FT to get the final δρ(x_ij)/ρ₀ entering the precursor.

All Figures

	Fig. 1. Real and imaginary parts of ω. The CR pressure gradient is constant, |∇PCR|=ξCRρ0u02/L, and the perturbation wave vector, k, is aligned with ∇PCR. We assume ξCR = 0.1, Ms = 100, and $M_{\mathrm{ms}}=100/\sqrt{5}$. The dashed and dotted lines show the asymptotic behavior of ω.
In the text