Spatio-temporal evolution of the nonresonant instability in shock precursors of young supernova remnants

TL;DR

This study uses advanced kinetic simulations with open boundaries to explore the spatial and temporal evolution of nonresonant cosmic-ray-driven instabilities in supernova remnant shock precursors, revealing shock formation, turbulence, and particle heating.

Contribution

First fully kinetic simulations with open boundaries to investigate both spatial and temporal development of the instability in shock precursors.

Findings

Magnetic-field amplification similar to previous studies.

Formation of shock-like compression structures with ion reflection.

Turbulent electric fields cause nonadiabatic plasma heating and modify cosmic-ray energy distribution.

Abstract

A nonresonant cosmic-ray-current-driven instability may operate in the shock precursors of young supernova remnants and be responsible for magnetic-field amplification, plasma heating and turbulence. Earlier simulations demonstrated magnetic-field amplification, and in kinetic studies a reduction of the relative drift between cosmic rays and thermal plasma was observed as backreaction. However, all published simulations used periodic boundary conditions, which do not account for mass conservation in decelerating flows and only allow the temporal development to be studied. Here we report results of fully kinetic Particle-In-Cell simulations with open boundaries that permit inflow of plasma on one side of the simulation box and outflow at the other end, hence allowing an investigation of both the temporal and the spatial development of the instability. Magnetic-field amplification proceeds as in studies with periodic boundaries and, observed here for the first time, the reduction of relative drifts causes the formation of a shock-like compression structure at which a fraction of the plasma ions are reflected. Turbulent electric field generated by the nonresonant instability inelastically scatters cosmic rays, modifying and anisotropizing their energy distribution. Spatial CR scattering is compatible with Bohm diffusion. Electromagnetic turbulence leads to significant nonadiabatic heating of the background plasma maintaining bulk equipartition between ions and electrons. The highest temperatures are reached at sites of large-amplitude electrostatic fields. Ion spectra show supra-thermal tails resulting from stochastic scattering in the turbulent electric field. Together, these modifications in the plasma flow will affect the properties of the shock and particle acceleration there.

Figures

Equations

$$\mathrm{Im}\omega ={\gamma }_{\max }=v_{\mathrm{A}}{k}_{\parallel \max }$$

$${\lambda }_{\max }=\frac{4\pi B_{\parallel 0}}{{\mu }_0{j}_{\parallel \mathrm{CR}}},$$

$$K(t)=\frac{⟨(x-x_0{)}^2⟩+⟨(y-y_0{)}^2⟩}{4t},$$

$$\frac{K(t)}{{\lambda }_{\max }^2{\gamma }_{\max }}=20.8{\gamma }_{\max }t+\frac{10^5}{3{\gamma }_{\max }t}{\sin }^2\left(\frac{{\gamma }_{\max }t}{40}\right),$$

$$D_{\mathrm{B}}=\frac{r_{\mathrm{CRg},\mathrm{f}}^2}{\tau },$$

$$D_{\mathrm{B}}\approx 740{\lambda }_{\max }^2{\gamma }_{\max },$$

$$\delta {\Pi }_{\mathrm{CR}}=-{\rho }_{\mathrm{up}}{V}_{\mathrm{up}}\delta V\Rightarrow \delta V=\frac{V_{\mathrm{up}}}{6}\frac{U_{\mathrm{CR}}}{U_{\mathrm{up}}}\simeq 2v_{\mathrm{A}}{N}_{\mathrm{e}-\mathrm{fold}},$$

$$L\simeq \frac{D}{v_{\mathrm{up}}}=\frac{c}{3v_{\mathrm{up}}}{\lambda }_{\mathrm{mfp}},$$

Full text

††pubyear: 2017††pagerange: Spatio-temporal evolution of the nonresonant instability in shock precursors of young supernova remnants–References

1 Introduction

Galactic cosmic rays (CRs) with energies up to a few peta-electronvolt (PeV) are widely believed to originate from young shell-type supernova remnants (SNRs), at whose forward shock waves particles can be energized in the diffusive shock-acceleration (DSA) process (Reynolds, 2008). The confinement of high-energy particles near the shock and the inferred rate of shock crossings require magnetic turbulence amplified to levels much larger than those typically found in the interstellar medium (ISM). Observations of X-ray synchrotron emission from young SNRs, that is partially produced in narrow filaments (e.g., Vink & Laming, 2003; Bamba et al., 2005) and can be variable on short time scales (Uchiyama et al., 2007), provide evidence for magnetic-field amplification by factors $10-100$ compared to the ISM level, even if turbulence damping accounts for the filaments (Pohl et al., 2005; Rettig & Pohl, 2012; Tran et al., 2015) and dynamic turbulence generates the variability (Bykov et al., 2008). Detection of TeV-energy $\gamma$-ray emission from the same sources (e.g., Acciari et al., 2011, 2010; Aharonian et al., 2007, 2009) proves electron and/or proton acceleration to very high energies and suggests a direct relation between energetic particle production and non-linear magnetic-field amplification. The latter is also needed in systems hosting relativistic shocks, such as jets of radio galaxies (e.g., Araudo et al., 2015) or gamma-ray bursts (GRB, see, e.g., Milosavljević & Nakar, 2006).

It has been long recognized on theoretical grounds that upstream of the shock magnetic fields can be amplified by the CRs themselves. Diffusing streaming of energetic particles in a quasi-parallel shock precursor can excite plasma instabilities generating magnetohydrodynamic (MHD) waves that are either resonant (Lerche, 1967; Kulsrud & Pearce, 1969; Wentzel, 1974; Skilling, 1975a, b, c; Achterberg, 1983) or nonresonant with CRs. It was shown by Bell (2004) that in young SNRs with efficient particle production the nonresonant CR streaming instability is the fastest one and generates short-scale magnetic turbulence with wavelengths much smaller than the gyroradii of the CRs driving it. The instability appears due to the current density, $\mn@boldsymbol{j}_{\rm CR}$, carried by CR protons drifting with respect to the magnetized upstream plasma. The most rapid growth occurs for a wavevector, $\mn@boldsymbol{k}$, aligned with the mean homogeneous magnetic field $\mn@boldsymbol{B}_{0}$ (Bell, 2005). In this case, and for $\mn@boldsymbol{j}_{\rm CR}$ parallel to $\mn@boldsymbol{B}_{0}$, the magnetic fluctuations take a form of a right-hand circularly polarized wave with the complex frequency $\omega$ and maximum growth rate

[TABLE]

at the wavelength

[TABLE]

where the subscript ‘max’ stands for the most unstable mode, $\mu_{0}$ is the permeability of free space and $v_{\rm A}$ is the Alfvén velocity. The CR current is balanced by the return current, $\mn@boldsymbol{j}_{\rm ret}=-\mn@boldsymbol{j}_{\rm CR}$, carried by the background plasma. The resulting Lorentz force, $\mn@boldsymbol{F}=\mn@boldsymbol{j}_{\rm ret}\times\mn@boldsymbol{B}$, acts on the plasma increasing the amplitude of the perturbed magnetic field (Zirakashvili et al., 2008). As the direction of the force is always perpendicular to the CR current, a purely growing (${\rm Re}\,\omega\approx 0$) mode results and transverse magnetic-field components are amplified. A plasma magnetization condition, $\omega\ll\Omega_{i}$, must apply, where $\Omega_{i}$ is the ion gyrofrequency. The original MHD treatment of the instability was confirmed by kinetic studies (e.g. Reville et al., 2006; Amato & Blasi, 2009).

The nonlinear evolution of the nonresonant instability has been numerically studied with MHD (Bell, 2004; Zirakashvili et al., 2008; Schure & Bell, 2013; Beresnyak & Li, 2014), hybrid kinetic (Lucek & Bell, 2000; Gargaté et al., 2010) and fully kinetic particle-in-cell (PIC) simulations (Niemiec et al., 2008; Riquelme & Spitkovsky, 2009; Stroman et al., 2009). These studies naturally model only a small portion of the far-upstream shock precursor and all of them employed periodic boundary conditions, implying the spatial impact of the instability is ignored. The spatial structure of the perpendicular case of the instability, in which $\mn@boldsymbol{j}_{\rm CR}\perp\mn@boldsymbol{B}_{0}$, was studied by Riquelme & Spitkovsky (2010). The primary effect of $\mn@boldsymbol{j}_{\rm CR}\perp\mn@boldsymbol{B}_{0}$ is a homogeneous linear acceleration of the upstream plasma with $\delta v=2\,V_{\mathrm{A}}$ for each growth time of the instability. One may also envision localized regions with $\mn@boldsymbol{j}_{\rm CR}\perp\mn@boldsymbol{B}_{0}$ which would then be differentially accelerated, inducing vorticity and further magnetic-field amplification by turbulent dynamo action. This situation may arise for lower-energy CRs which are confined closer to the shock and can produce a transverse CR current by interacting with the turbulent field structures generated already far upstream by high-energy CRs via the nonresonant instability with $\mn@boldsymbol{j}_{\rm CR}\parallel\mn@boldsymbol{B}_{0}$, that is presumed to have provided large chunks of plasma with magnetic field perpendicular to the shock normal.

The studies using periodic boxes and field-aligned CR currents demonstrated significant magnetic-field amplification by an initially parallel wave with preferred wavelength $\lambda_{\rm max}$ as predicted by the linear theory. In the nonlinear phase, density fluctuations are imposed in the quasi-thermal plasma, and the plasma acquires strong transverse velocity fluctuations. The net $\mn@boldsymbol{j}\times\mn@boldsymbol{B}$ force pushes on the plasma that eventually forms low-density cavities surrounded by filaments of compressed material. As the cavities expand and merge with each other, the characteristic scale of magnetic-field fluctuations increases, but the turbulence growth eventually slows down or saturates. The exact saturation mechanism differs with the choice of the simulation technique, and the saturation levels vary considerably.

MHD simulations use a fluid model for the background plasma and represent CRs as a constant external current, which means that backreaction on CRs is suppressed. The saturation is due to magnetic-field tension opposing the $\mn@boldsymbol{j}\times\mn@boldsymbol{B}$ force and scales with wavenumber to yield a saturation level $B_{\rm sat}\sim\mu_{0}j_{\parallel\rm CR}/k$. Since the growth rate increases with both the wavenumber and the CR current, $\gamma\sim(kj_{\parallel\rm CR})^{1/2}$, the faster growing short-scale modes saturate earlier and at lower amplitudes than the slowly-growing longer-wavelength modes. Consequently, as the turbulence evolves its characteristic wavelength also increases until it becomes commensurate with the CR Larmor radius, $kr_{\mathrm{CRg}}\approx 1$, if the CR current is assumed constant. However, in reality the turbulence growth slows down once the magnetic-field amplitude increases to $\delta B/B_{0}\ga 1$ and density fluctuations appear, grow in size, acquire transverse motions and start to collide with each other. The rms field amplitude will still slowly grow, and a full saturation is never observed in MHD simulations (Zirakashvili et al., 2008), unless imposed by the boundary conditions and the finite size of the simulation box (Bell, 2004). The final unsaturated field amplitudes obtained with the MHD approach reach $\delta B/B_{0}\simeq 30-50$.

Fully kinetic PIC simulations allow treating both the electron and ion dynamics self-consistently and can thus account for the feedback on CRs of the amplified magnetic field. The backreaction causes a net transfer of momentum to the plasma. CRs slow down in bulk, and the background plasma is accelerated in the direction of the CR drift, reducing the relative speed and removing the CR current seen by the plasma. The resulting magnetic-field saturation level, $\delta B/B_{0}\simeq 10-20$, is consistently lower than that obtained in MHD simulations (Riquelme & Spitkovsky, 2009; Stroman et al., 2009).

Similar saturation mechanisms and levels were observed in PIC simulations of relativistic ion beams studied in the limit of magnetized ambient plasma (Niemiec et al., 2010), in which the turbulence arises from growth of the nonresonant modes (Reville et al., 2006). Relativistic ion beams in this regime can be responsible for magnetic-field amplification upstream of external GRB shocks (Milosavljević & Nakar, 2006) and at relativistic shocks associated with jets of active galactic nuclei. However, in such systems turbulence driven by the CR current can be accompanied with electromagnetic filamentation and electrostatic Buneman modes (Niemiec et al., 2010). The saturation physics is also consistent with quasi-linear predictions for nonresonant modes derived for nonrelativistic beams (Winske & Leroy, 1984) and for CR streaming distributions with power-law momentum dependence (Luo & Melrose, 2009). Hybrid kinetic (fluid electrons and PIC model for ions) simulations by Gargaté et al. (2010) yielded comparable amplification levels of $\delta B/B_{0}\simeq 10$ and noted the role of effective de-magnetization of the plasma through its rapid transverse motions. The importance of kinetic plasma treatment was further emphasized in PIC simulations that assumed a constant external CR current as in MHD studies (Ohira et al., 2009; Riquelme & Spitkovsky, 2009), in which case the saturation is still caused by the reduction of the relative speed between CRs and plasma, but now the plasma is accelerated to match the constant CR bulk drift (see also Niemiec et al., 2010).

All previously published simulations used periodic boundary conditions, which suppress the spatial impact of the bulk acceleration of CRs and quasi-thermal plasma. In reality, continuity mandates compression or expansion to go along with bulk acceleration (Ohira et al., 2009), and CRs will continuously stream into fresh, unaffected plasma in which the instability must build up. The aim of this work is to study the nonlinear development and saturation processes of the nonresonant CR current-driven instability with open boundaries that permit inflow of plasma on one side of the simulation box and outflow at the other end. In this way the mass is conserved in the system and compression is naturally accounted for. This simulation represents the next step toward realistic modeling of nonresonant magnetic-field amplification in SNR shock precursors. Our study reproduces some of the earlier results and also finds new effects specific to our new setup, that allows us to investigate both the temporal and the spatial development of the nonresonant instability. Particular attention is given to the microphysics of the saturation process, CR scattering and plasma heating.

Although our approach offers a global model of the far-upstream region of the shock precursor, we do not directly include the shock that accelerates the CRs. Therefore, we cannot self-consistently determine the shock structure as modified by the CRs. Nonlinear steady-state shock models (e.g., Zirakashvili et al., 2008; Schure & Bell, 2013; Bykov et al., 2014) seek balance between magnetic-field amplification, that allows for CR confinement and enables rapid acceleration by the DSA process and CR escape that provides the current driving the field generation (Bell et al., 2013). Our study can define conditions at which stationary solutions are possible and aid in the validation of the assumptions on particle scattering and heating which are important ingredients of these models.

The numerical PIC model is described in Section 2 and results of our modeling are presented in Section 3. Section 3.1 describes the main features of the nonresonant instability. The microphysics of the saturation process is discussed in Section 3.2, that treats CR scattering and plasma heating. We present a discussion in Section 4 and summarize our results in Section 5.

2 Simulation setup

The nonresonant instability arises in the precursor of a strong nonrelativistic parallel shock when shock-accelerated CRs drift against the cold upstream electron-ion plasma. Our earlier simulations of this instability with periodic boundary conditions were conducted in the plasma-ion rest frame (Niemiec et al., 2008; Stroman et al., 2009). In this reference frame, an isotropic population of relativistic CR ions with number density $N_{\mathrm{CR}}\ll N_{e}$ moves with the shock velocity, $v_{\mathrm{sh}}$, along a homogeneous magnetic field. The current density carried by CR ions, $j_{\mathrm{CR}}=eN_{\mathrm{CR}}v_{\mathrm{sh}}$, where $e$ is the electric charge, is compensated by the return current density, $j_{\mathrm{ret}}=-eN_{\mathrm{e}}v_{\mathrm{d}}=-j_{\mathrm{CR}}$, provided by electrons of number density $N_{\mathrm{e}}$ that slowly drift relative to background ions with $v_{\mathrm{d}}=v_{\mathrm{sh}}N_{\mathrm{CR}}/N_{\mathrm{e}}$. We add an excess electron population with number density $\delta N_{\mathrm{e}}=N_{\mathrm{e}}-N_{\mathrm{i}}=N_{\mathrm{CR}}$ to maintain charge neutrality.

Maintaining stable conditions with open boundaries requires that the calculations in this work be performed in the cosmic-ray rest frame. As illustrated in Fig. 1, in the CR rest frame the electron-ion plasma forms a beam with velocity $v_{sh}$. Charge compensation is established with excess electrons that are initiated at rest to maintain current balance. With this method stable injection of the electron-ion beam into the numerical box is simpler to achieve than with, e.g., adding the excess electrons to the electron-ion beam.

In the linear phase the nonresonant mode has its peak growth rate at a wavelength smaller than the CR gyroradius, $\lambda_{\mathrm{max}}\ll r_{\mathrm{CRg}}$. The simulation results presented here are based on a PIC experiment utilizing an extremely large two-dimensional computational grid of size $L_{x}\times L_{y}=130,000\times 12,000$ cells and a total of $5.6\cdot 10^{10}$ macro-particles. In terms of $\lambda_{\mathrm{max}}$, the size of the numerical box is $L_{x}\times L_{y}\simeq(293\times 27)\,\lambda_{\mathrm{max}}^{2}$. All system constituents – CR ions, electrons and ions of the ambient plasma beam as well as the electrons of the excess population – initially fill up the whole computational box. Periodic boundary conditions for particles and electromagnetic fields are applied in the transverse ($y$) direction. The right and left box boundaries are open for the ambient beam particles and fields. Fresh particles of the electron-ion beam are injected at the left side of the box, where absorbing boundary conditions are applied for fields. Outflowing beam particles are removed at the right side of the box, 50 cells off the grid boundary. Here the mean fluxes of electrons and ions are nearly identical, and small-scale electromagnetic impulses generated by local charge imbalances propagate freely towards the open right boundary for fields, at which they are effectively absorbed. CR ions are kept inside the box and reflected, if they encounter the open boundaries for the ambient plasma. The electron-ion plasma beam pushes the excess electrons to the right, leaving the much heavier CRs initially unaffected. The simulated system quickly adjusts itself to the new situation in a numerically stable way, producing a suitable return current and compensating for the small charge imbalance, exactly as one would envision a real plasma to respond. The excess electrons concentrate towards the right boundary for particles. They are not removed from the system in order to maintain charge conservation. On account of the high ambient plasma-to-CR density ratio, the concentrations in excess electrons are always a small perturbation to the total electron distribution, even in the nonlinear stage of the instability, during which density and field structures are produced that focus the plasma into compact islands. Therefore, the excess electrons do not play a significant role in the dynamics of the system. Numerical artifacts at the right boundary cause field fluctuations with much smaller amplitude than those of the turbulence in the region of interest. These fluctuations do not radiate into the box in any significant form. Therefore, the plasma outflow boundary does not have a noticeable influence on the evolution of Bell’s modes in our numerical experiment.

The numerical parameters are similar to those used in our earlier simulation with periodic boundaries (Stroman et al., 2009), to permit a direct comparison. As in that study, we also choose $\mn@boldsymbol{j}_{\rm CR}$ to be aligned with $\mn@boldsymbol{B}_{0}$, in which case the turbulence takes on the simplest form (Bell, 2005). The chosen parameters were shown in Stroman et al. (2009) to well resolve all physical characteristics of the nonresonant instability. The key issue in properly reproducing the instability characteristics is to assure the so-called magnetization condition, $\gamma_{\rm max}/\Omega_{i}\ll 1$. This condition does not explicitly depend on specific values of $\gamma_{\rm CR}$ or $N_{\rm CR}/N_{i}$, but these quantities must be chosen so that the instability condition is fulfilled. It is thus justified to say that our simulation probes the far-upstream structure of the precursor with high-energy cosmic rays. The rationale for selecting specific parameters was extensively given in Niemiec et al. (2008) and in Stroman et al. (2009), and among the drivers is the low growth rate for high-energy CRs of realistic density, that does not permit PIC simulations even in 2D. This rationale is supported with numerical tests, in particular for the range of $N_{\rm CR}/N_{i}$ in Niemiec et al. (2008) (also in Niemiec et al., 2010) and $\gamma_{\rm CR}$ in Stroman et al. (2009). The electron skin depth $\lambda_{\rm se}=c/\omega_{\rm pe}=4\Delta$, where $\Delta$ is the grid cell size, $c$ is the speed of light, $\omega_{\rm pe}=\sqrt{e^{2}N_{\mathrm{e}}/\epsilon_{0}m_{\mathrm{e}}}$ the electron plasma frequency and $m_{\mathrm{e}}$ is the electron mass. CR ions are represented by a population of isotropic and monoenergetic particles with Lorentz factor $\Gamma_{\rm CR}=50$. The number density ratio between plasma ions and CR ions is $N_{\mathrm{i}}/N_{\mathrm{CR}}=50$, and the assumed plasma magnetization, given as the ratio of the peak growth rate to the nonrelativistic ion gyrofrequency, is $\gamma_{\rm max}/\Omega_{\mathrm{i}}=0.4$. The homogeneous magnetic field, $B_{0}=B_{\parallel 0}$, is chosen for a CR gyroradius $r_{\rm CRg}\simeq 140,000\Delta$ and an Alfvén speed $v_{\rm A}=0.01c$, with $v_{\rm A}=[B_{\parallel 0}^{2}/\mu_{0}(N_{\mathrm{e}}m_{\mathrm{e}}+N_{\mathrm{i}}m_{\mathrm{i}}]^{1/2}$. The ambient plasma beam moves in the $x$-direction, parallel to the uniform magnetic field, with velocity $v_{\mathrm{sh}}=0.4c$. The initial thermal speed of electrons, $v_{\mathrm{e,th}}=0.01c$, is the same for both plasma and excess electrons. The ions are in thermal equilibrium with the electrons. Computational constraints force us to use a reduced ion-to-electron mass ratio, $m_{\mathrm{i}}=50\,m_{\mathrm{e}}$, for which the spatial and temporal scales are still clearly separated. The wavelength of the most unstable mode is $\lambda_{\rm max}\simeq 444.3\Delta$ and the inverse growth rate is $\gamma_{\rm max}^{-1}\simeq{\textrm{14,142}}\,\delta t$, where the timestep $\delta t=0.125/\omega_{\rm pe}$.

We follow the system evolution through the strongly nonlinear phase up to $t\,\gamma_{\rm max}=26.9$. Beyond this time the largest scales of the turbulence structure start to become comparable to the transverse size of our numerical box, $L_{y}$, and cannot thus be properly resolved. The simulation uses $N_{\rm ppc}=9$ particles per cell per species, and statistical weights, $w_{g}=0.02$, are applied to CRs and excess electrons to establish the desired density ratio. Pairs of ions and electrons are initiated at the same locations, so that the initial charge and current vanish. To stabilize the plasma against nonphysical effects on very long time scales and with a moderate number of particles per cell we apply the numerical PIC model used recently in Niemiec et al. (2012). The model involves second-order particle shapes (TSC – Triangular-Shaped Cloud) and uses a second-order FDTD field-solver with a weak Friedman filter (Greenwood et al., 2004) to suppress small-scale noise. The simulations were performed with a modified version of the TRISTAN code (Buneman, 1993) parallelized using MPI (Niemiec et al., 2008). We follow all three components of the particle velocities and electromagnetic fields. Note that this 2D3V model was demonstrated to well reproduce the results of fully three-dimensional kinetic simulations, including both the magnetic field amplification level and the turbulence structure (Niemiec et al., 2008; Riquelme & Spitkovsky, 2009). For the unstable wavevector aligned with $\mn@boldsymbol{B}_{0}$ and contained in the simulation plane, spatial distributions of particle densities and electromagnetic fields are projections onto a 2D plane of the turbulent structures observed in three-dimensions with kinetic and also MHD simulations. The neglect of gradients in z-direction inherent to 2D3V simulations may somewhat impact the development of 3D turbulence structures in the deep non-linear phase of the instability, e.g. concerning the size and lifetime of plasma voids, but that is more than balanced by the large benefit of our 2D3V numerical experiment covering a much wider range of dynamical scales of nonlinearly evolving turbulence than can be possibly afforded with a 3D simulation under current computing capabilities.

The numerical model used in this study has been extensively tested. Test simulations with periodic boundary conditions in all directions verified that the physical features of the nonresonant instability (the growth rate, the wavelength of the most unstable mode, the nonlinear response) in the cosmic-ray rest frame match exactly those obtained in the plasma rest frame with the same parameters (Stroman et al., 2009). The results do not significantly change if the system resolution (the electron skindepth) or the number of particles per cell, $N_{\rm ppc}$, are increased. The numerical stability of plasma injection and escape at the open boundary has been verified in test simulations with an electron-ion beam propagating in an otherwise empty box. Stability at the left boundary is ensured by the charge-neutral injection of particles. Radiation-absorbing boundary conditions are used as originally implemented in the TRISTAN code (Buneman, 1993).

3 Results

The open boundaries used in this numerical study allow us to investigate both the temporal and the spatial development of Bell’s nonresonant instability, and they permit maintaining mass conservation in the system. In Section 3.1 we demonstrate magnetic-field amplification through the nonresonant instability, its saturation and the main features of its backreaction on CRs that re-confirm our earlier results obtained in PIC simulations with periodic boundaries. We also describe new effects specific to our current, more realistic simulation setup. Section 3.2 presents the microphysical picture of the nonlinear particle-wave interactions during the saturation phase that result in kinetic-energy transfer, CR scattering and plasma heating.

3.1 Nonresonant instability: main features

In this section we describe the state of the system at time $t\gamma_{\rm max}=13$, when the magnetic-field growth reaches saturation, and at time $t\gamma_{\rm max}=26.9$, when the system is in the strongly nonlinear phase. We present a global spatial structure of the far-upstream precursor in Section 3.1.1. Spatio-temporal development of magnetic field turbulence is described in Section 3.1.2, the effects of the nonlinear plasma response in Section 3.1.3 and the electric field turbulence in Section 3.1.4. Finally, in Section 3.1.5 we demonstrate the effects of a compression wave – a new feature revealed through our modeling with open boundaries – on the field turbulence and particles in the precursor.

3.1.1 Spatial structure of far-upstream precursor

Figs. 2 and 3 show the distributions of the transverse magnetic-field component, $B_{z}$, and the number density of ambient ions at the two points in time, respectively.

To be noted from Figs. 2 and 3 is the spatial imprint of the time development of the instability that arises from the motion of the wave-carrying plasma from left to right. Growth of magnetic turbulence becomes evident at $x/\lambda_{\rm max}\approx 22$ and reaches its saturation at $x/\lambda_{\rm max}\approx 75$. In studies with periodic boundaries only the temporal evolution was observable. As in the present setup the entire simulation box is filled with CRs and with streaming plasma at the beginning, the instability development is mainly temporal in the right part of the box in which the beam particles have travelled the same distance on the CR background from the time of their initial injection. Wave growth begins only close to the left boundary, where fresh plasma is injected, and as the instability needs a few growth times to reach a significant amplitude, the wave-carrying plasma will have streamed a certain distance. Later phases of turbulence development need more time, and hence the plasma will have propagated a larger distance. In the end, the temporal stages of the instability – exponential growth, nonlinear response, saturation and cascading – are mapped to certain locations in $x$ that slowly move to the right on account of the displacement of the CRs, that are pushed on by the streaming plasma in the nonlinear phase. Travel time is thus unimportant for $x/\lambda_{\rm max}\gtrsim 75$ at $t\gamma_{\rm max}=13$ (Fig. 2) and for $x/\lambda_{\rm max}\gtrsim 100$ at the end of the simulation (Fig. 3), where we primarily observe the temporal development.

One feature of nonlinear feedback in periodic simulations using plasma-ion rest frame was bulk deceleration of cosmic rays and acceleration of the thermal plasma until the remaining streaming was insufficient to further drive turbulence. As we discuss in Section 3.1.2, the same effects are observed in our current study performed in the CR rest frame. One consequence of this backreaction in our realistic modeling with open boundaries is the appearance of a strong density enhancement, clearly visible in Fig. 3a at $x/\lambda_{\rm max}\approx 100-130$. This compression is caused by the collision of the incoming electron-ion plasma with the decelerated plasma in the turbulence zone. The compression factor still increases at the end of the simulation and approaches the value of 4, as in a strong nonrelativistic shock. Detailed features of this shock-like structure will be described in the sections that follow.

3.1.2 Spatio-temporal evolution of magnetic turbulence

The spatial and the temporal development of the instability are connected through the streaming of the plasma, and proper interpretation requires knowledge of its bulk velocity. We have averaged over coordinate $y$ various plasma variables to produce their spatial profiles in coordinate $x$. Fig. 4 displays such profiles of the ion number density, the magnetic and electric-field amplitudes, the $E_{x}$ component of the electric field and particle velocities.

To be noted from Fig. 4b is the magnetic-field amplification to about $20\,B_{\parallel 0}$ at $t\gamma_{\rm max}=26.9$, slightly exceeding that seen in simulations with periodic boundaries (Stroman et al., 2009) on account of compression. The bulk acceleration/deceleration observed in earlier studies is also evident (Fig. 4d), and at the location of the peak amplitude of the mean magnetic field the bulk velocities of plasma and CRs have converged to $\delta V\lesssim 0.05\,$c. Associated plasma compression by up to a factor 4, identified already in Section 3.1.1, is seen in Fig. 4a for $x/\lambda_{\rm max}\approx 90-130$. Essentially all of it arises late in the evolution, starting from $t\gamma_{\rm max}\approx 13.0$, close to the instability saturation time, and continuing until the end of the simulation at $t\gamma_{\rm max}=26.9$. The onset of compression moves with the plasma flow that pushes on the turbulent zone, and at times $t\gamma_{\rm max}=13$ and $26.9$ it is contained in the white rectangular boxes in Figs. 2 and 3, respectively.

Knowledge of the plasma bulk speed permits decoupling of the temporal and the spatial evolution of the instability, at least in the left part of the simulation box. For that purpose we have defined 4 locations labeled $x_{1}$, $x_{2}$, $x_{3}$ and $x_{4}$ and marked them in Fig. 2. We also marked in white a region around $x_{4}$. These 4 locations correspond to the first emergence of waves (1), their exponential growth (2), entering the regime $\delta B\gtrsim B_{\parallel 0}$ (3) and saturation (4). As wave-particle interactions can begin only at the left boundary of the simulation box and a certain distance implies a specific streaming time, the four locations correspond to points $t_{1}$, $t_{2}$, $t_{3}$ and $t_{4}$ in the temporal evolution of the instability.

Recalling that the streaming time for $135\lesssim x/\lambda_{\rm max}\lesssim 180$ is not smaller than the simulation time, implying we primarily observe the temporal evolution, we can directly compare our new results with the earlier investigations. In Fig. 5 we plot with solid lines the time evolution of the magnetic-field amplitudes spatially averaged in the region $135\lesssim x/\lambda_{\rm max}\lesssim 180$, that does not contain the compression wave. We add the time markers $t_{1}$ – $t_{4}$ and can now compare the state of the system with that observed at the corresponding locations $x_{1}$ – $x_{4}$ in Fig. 2.

The waves that start to emerge at $x/\lambda_{\rm max}=x_{1}\gtrsim 22.5$ (Fig. 2) move with the beam (i.e., they are stationary in the ambient plasma rest frame) and have wave vectors parallel to the plasma drift with wavelength $\sim\lambda_{\rm max}$. At this location the electron-ion beam has propagated through the CR plasma for $t\gamma_{\rm max}=t_{1}\approx 3.5$, and we see in Fig. 5 that at this time the first sign of growth in $B_{\perp}$ appears. At $x/\lambda_{\rm max}=x_{2}\approx 32$ the waves are clearly visible in Fig. 2 and in the temporal development depicted in Fig. 5 the exponential growth with rate $\gtrsim 0.9\,\gamma_{\rm max}$ begins. At $t\gamma_{\rm max}=t_{3}\approx 8$ the turbulent field becomes comparable to the strength of the homogeneous magnetic field, $\langle\delta B/B_{\parallel 0}\rangle\ga 1$. We observe at the corresponding location $x_{3}$ that the magnetic fluctuations start to deviate from a purely parallel mode, and the dominant length scale begins to exceed $\lambda_{\rm max}$. The subsequent evolution leads to isotropic turbulence beyond $x/\lambda_{\rm max}=x_{4}\approx 75$, corresponding to the saturation at $t\gamma_{\rm max}\gtrsim t_{4}=13$ at the level $\langle\delta B/B_{\parallel 0}\rangle\approx 15$. At later times the mean field amplitude stays at roughly the same level but the turbulence scale quickly evolves to very long wavelengths, as is seen in Fig. 3 for $x/\lambda_{\rm max}\ga 110$.

The magnetic-field is amplified to essentially the same amplitude as in our earlier simulations using periodic boundaries (Stroman et al., 2009). Fig. 5 also displays with dashed lines the magnetic-field amplitudes in the compression region where, as noted before, larger amplification is reached. One can note that magnetic-field amplification in the compression wave grows steadily, and by the end of our numerical experiment the field has an amplitude $\sim 30\%$ higher than that provided through the nonresonant instability (compare Fig. 4b). It does not reach saturation, which suggest that the compression structure has not yet attained a steady state, as shown by the blue dashed curve in Fig. 5. The total magnetic-field amplification in a portion of the precursor thus results from compression, and our simulation is uniquely suited to recover it.

3.1.3 Nonlinear plasma response

As the magnetic field grows to amplitudes $\langle\delta B/B_{\parallel 0}\rangle\ga 1$, plasma-density perturbations appear (Fig. 2a for $x\ga x_{3}$), that are accompanied by velocity fluctuations. In this strongly nonlinear phase, plasma parcels begin to move, collide with each other and collapse into turbulently moving filaments, causing distortions of the initially parallel magnetic modes. The plasma also slows down in bulk, in particular during the late phase. Whereas at time $t\gamma_{\rm max}=13$ the electron-ion beam decelerated only from the initial drift speed of $0.4c$ to $\sim 0.35c$ in the region of strong turbulence ($x/\lambda_{\rm max}\ga 60$), by $t\gamma_{\rm max}=26.9$ the beam drifts slowly with $v\sim 0.15c$ (red lines in Fig. 4d). Bulk CR motion with mean velocity of $\sim 0.1c$ is likewise observed, and the bulk velocities of the three plasma species converge (the slight reduction in CR bulk speed toward $x\approx 180\,\lambda_{\rm max}$ results from the reflecting boundary condition for CRs). Further growth of the turbulence is thus limited in this region and the instability saturates.

If we envision Bell’s instability to be an agent of magnetic-field amplification in the precursor region of SNR shocks, we must also consider it as a means of shock modification by deceleration and compression of the upstream plasma (for a review see, e.g., Blandford & Eichler, 1987). We initiate our simulations with a constant CR density, and so the backreaction is not driven by a gradient in CR pressure, but generates it. In any case, the feedback on the upstream plasma will also modify the thermal sub-shock, and a new balance between particle acceleration and its feedback will be established. We do not explicitly model the shock that provides the CRs, and so we cannot determine that statistical equilibrium between shock acceleration of CRs, the build-up of magnetic turbulence in the precursor, and the backreaction of CRs and turbulence on the inflowing plasma.

3.1.4 Electric field turbulence

Associated with the magnetic turbulence in the nonlinear stage is a turbulent electric field whose amplitude is an order of magnitude smaller than that of the magnetic field (see Fig. 4b). This electric field initially arises as $\mn@boldsymbol{v}\times\mn@boldsymbol{B}$ motional term of the bulk flow with the mainly transverse magnetic fluctuations, and hence its dominant components are $E_{y}$ and $E_{z}$. The electric fluctuations appear when the magnetic modes have grown to nonlinear amplitude (e.g., at $x/\lambda_{\rm max}\approx 45$ at $t\gamma_{\rm max}=13$ in Fig. 4b).

Once the plasma attains turbulent motions, an $E_{x}$ component appears as well (for $x/\lambda_{\rm max}\approx 55$ at $t\gamma_{\rm max}=13$, Fig. 4b). In contrast to the perpendicular electric-field components that oscillate around zero, the average longitudinal component in the turbulent zone is essentially always positive, as demonstrated in Fig. 4c. Most of the $E_{x}$ field results from $v_{y}B_{z}$ and $v_{z}B_{y}$ terms. The presence of this mean parallel electric field is a response of the system to the CR electric current (e.g., Bell, 2005; Zirakashvili et al., 2008). Its role in the microscopic particle motions will be discussed in Section 3.2.

The spatial structure of the plasma velocity profile can be traced in Fig. 6 that shows the longitudinal phase-space distribution at $t\gamma_{\rm max}=26.9$ (see also discussion on Fig. 9 in Sec. 3.2). Evident in Fig. 6 for $x/\lambda_{\rm max}\ga 90$ is the significant plasma heating in the region of the strong density and velocity fluctuations. The irregular phase-space structure for electrons visible in Fig. 6a reflects localized heating, the details of which will be discussed in the next section. Here we only note that although the dominant part of the electric fields arises from magnetic-field transport, there also exists an electrostatic field component. It can be calculated by either performing a Lorentz transformation of the fields into the local plasma rest frame or by directly solving the Poisson equation. Both methods reveal electric-field fluctuations that can inelastically scatter particles and modify their distribution.

3.1.5 Phase-space structure of the compression wave

As noted above, our realistic precursor modeling reveals the compression wave that can be observed in Figs. 3 and 4 for $x/\lambda_{\rm max}\approx 90-130$. The phase-space structure in the vicinity of this feature is shown in Fig. 6. Beside the density increase by a factor of four (Fig. 4a), the plasma is strongly heated, as observed in nonrelativistic shock waves. The structure also exhibits a distinctive feature of a supercritical shock – in the upstream region a population of the reflected plasma ions is observed that it is clearly visible in Fig. 6b as a diffuse particle beam with negative momenta for $x/\lambda_{\rm max}\la 90$, extending almost to the plasma injection boundary at the left side of the simulation box. Although the phase-space density of the reflected ions is very low, they cause substantial heating of the incoming plasma electrons. They possibly also modify the magnetic-field structure through, e.g., ion-beam streaming instabilities (e.g., Niemiec et al., 2012). Any distortions of the magnetic turbulence resulting from such instabilities would appear as magnetic filaments and are thus difficult to identify in the fluctuations that are already organized in the form of filaments upstream of the shock-like structure. Additional effects of the compression wave will be discussed in particle spectra presented in the section that follows.

3.2 Microphysics of the saturation process

The energy and momentum transfer between CRs and the plasma that leads to the saturation of the nonresonant instability is mediated by highly turbulent large-amplitude electromagnetic fields. These fields scatter the particles and significantly modify their distributions. In this section we describe these processes at the micro-physical level. CR scattering in the Bell’s turbulence is presented in Section 3.2.1 and discussed in the context of their spatial diffusion in Section 3.2.2. The background electron-ion plasma heating and scattering is described in Sections 3.2.3, 3.2.4 and 3.2.5.

3.2.1 Cosmic-ray scattering

Figs. 7 and 8 show CR kinetic-energy spectra and phase-space distributions at the end of our simulation at $t\gamma_{\rm max}=26.9$. Results for the four physically distinct regions marked in Fig. 6 are presented with different colours. The degree of CR scattering and anisotropy is apparent. It is also evident that most of the scattering is inelastic and results in either acceleration or deceleration of CR particles as they probe the turbulent precursor regions.

Deeper insight into CR scattering is provided by Fig. 9, that shows the kinetic-energy evolution and the $x$-component of particle trajectories for typical CRs. Scattering in energy commences at time $t\gamma_{\rm max}\approx t_{3}=8$ (see Fig. 5), when the amplitude of magnetic turbulence becomes comparable to the strength of the homogeneous field, $\langle\delta B/B_{\parallel 0}\rangle\ga 1$. The mean energy increases from the initial $\Gamma-1=49$ to approximately $50$ by $t\gamma_{\rm max}\approx 18$ and shows little change at later times. Much of the mean energy gain arises from the bulk energy transfer from the ambient plasma beam that is demonstrated in Fig. 4d. Correspondingly, the most efficient scattering takes place before $t\gamma_{\rm max}\approx 18$. Inspection of individual particle trajectories shows that CRs initially move on almost rectilinear paths on account of their large Larmor radii and then become deflected in pitch angle as they pass through regions of strong magnetic field of varying polarity. Part of the scattering may be resonant with the waves generated via the instability (Stroman et al., 2009; Luo & Melrose, 2009), but most of it results from inelastic scattering off electromagnetic turbulence. Consequently the time histories of particle energy in Fig. 9a show corrugated structures. Scattering in energy occurs in the electric fields associated with the motion of turbulent magnetic field.

An interesting feature can be observed in Fig. 9: CRs that travel along with the plasma flow in the positive $x$-direction are accelerated (e.g., particles No. 1 and 2, shown with black and magenta lines), whereas those moving against the flow become decelerated (e.g., CR No. 6, shown with the yellow line). Particles that change the $x$-component of their motions upon scattering seem to gain or lose their energy accordingly (e.g., CRs No. 3, 4, and 5, shown with blue, green and red lines). To be scattered, a particle must reside in the region of strongly nonlinear turbulence, and hence fluctuations in energy are seen predominantly at times $t\gamma_{\mathrm{max}}\gtrsim 10$ and particles at locations $x\gtrsim 80\lambda_{\mathrm{max}}$. These characteristics of CR scattering result from the structure of turbulent electric fields described in Section 3.1.4.

The $E_{x}$ field component in the turbulence region is on average nonzero and directed in the positive $x$-direction. It is induced in the system by the CR electric current and allows bulk-energy transfer between the CRs, the ambient plasma and the turbulence. In the initial CR rest frame of our simulation the mean electric field will extract energy from the ambient plasma beam at a rate $-E_{x}j_{\rm ret}$, causing a bulk acceleration of CRs and a bulk slow-down of the ambient beam flow velocity, as observed in Figs. 4d and 9a. At the level of individual particle trajectories the same field causes a net acceleration for CR ions traveling along the flow and a net deceleration for CR particles moving in the negative $x$-direction. Momentum change by the dominant $E_{y}$ and $E_{z}$ fields averages to zero for CRs. The thermal speed of ions is small, as is their Larmor radius, and the ions are more susceptible to the effects of the turbulent fields than are CR particles. We will describe in detail the microphysics of ion scattering and their heating in Sections 3.2.3 – 3.2.5.

The yellow line in Fig. 7 indicates a high-energy bump resulting from efficient CR scattering, which requires a large region of strong turbulence. Evidently our simulation setup provides a sufficient volume filled with turbulence, but in reality, this region may be limited in size by the presence of a shock and also by lower-energy CRs that reside closer to the shock and modify the turbulence structure there. Hence, the CR spectra formed in this region are presented here to illustrate the main features of the physical processes in the precursor and not to provide quantitative predictions. The latter would require a full shock simulation, including the electron dynamics. A similar remark applies to the spectrum and phase-space distribution of CRs calculated for the very far upstream region ($x/\lambda_{\rm max}=10-30$, blue line/points in Fig. 7/Fig. 8a). In this zone CRs scattered off their initial distribution might be overabundant on account of the reflecting boundary that prohibits leakage of these particles away from the strongly turbulent region in which they originated. A similar caveat applies to the region $x/\lambda_{\rm max}=160-180$ where reflected CRs from the rear boundary lower the bulk speed in the simulation. The distributions calculated in the zone of small-amplitude turbulence ($x/\lambda_{\rm max}=70-90$, green line/points) and at the compression wave ($x/\lambda_{\rm max}=90-110$, red line/points) are least influenced by boundary conditions. Some anisotropy in these distributions comes from particles reflected off the compressed field of the shock-like structure.

3.2.2 Spatial diffusion of cosmic rays

The motion of CRs in the region of strong magnetic turbulence resembles diffusion in space. We can calculate a running diffusion coefficient as

[TABLE]

where $x_{0}$ and $y_{0}$ are the initial particle locations. The running diffusion coefficient, $K(t)$, and its time derivative are shown in Fig. 10 with solid red lines. The coefficient is calculated from time $t\gamma_{\rm max}=t_{4}=13$ until the end of the simulation and for all CRs that continuously reside in the strongly turbulent zone at $x/\lambda_{\rm max}\ga 110$. Particles diffusing in from the small-amplitude region, as well as those that leak into the latter zone, are ignored.

One notes that CR motion is not truly diffusive, as the running diffusion coefficient does not reach constancy, $K(t)={\rm const}$. However, it clearly deviates from helical motion. If particles would only gyrate in the initial homogeneous magnetic field, then the running diffusion coefficient would behave as

[TABLE]

which leads to the dashed green lines in Fig. 10. $K(t)$ clearly deviates from that but does not reach a constant level during the simulation. As Bell’s instability operates at $k\,r_{\mathrm{CRg}}\gg 1$, the magnetic-field fluctuations can only impose scattering, and the time derivative of $K(t)$ in Fig. 10b indicates a linear decrease toward $dK(t)/dt=0$, whereas gyration would lead to a parabolic decline. One can estimate the time, $t_{\rm diff}$, needed for CRs to transition into the diffusion regime by extrapolating the derivative down to $dK/dt=0$, i.e., $K(t_{\rm diff})={\rm const}$. The required time is $t_{\rm diff}\gamma_{\rm max}\approx 36.5\pm 1$, and the spatial diffusion coefficient would reach $D=K(t_{\rm diff})=(770\pm 20)~{}\lambda^{2}_{\rm max}\gamma_{\rm max}$. The time $t_{\rm diff}$ considerably extends beyond the duration of our numerical experiment.

Let us compare our approximate value of the diffusion coefficient, $D$, with that allowed in the Bohm limit. For that purpose we shall write the Bohm diffusion coefficient as

[TABLE]

where $r_{\rm CRg,f}$ is the final CR gyro-radius in the amplified magnetic field and $\tau$ is the mean time between scattering events. In the root-mean square field, the final CR gyro-radius is $r_{\rm CRg,f}\approx 18\lambda_{\rm max}$. The typical coherence length of the magnetic turbulence is about $7\lambda_{\rm max}$, and the time needed for CRs to propagate one coherence length may be taken as proxy of the scattering time, $\tau\gamma_{\rm max}\approx 0.44$. Our estimate of the Bohm diffusion coefficient hence is

[TABLE]

which is numerically similar to our estimated value of $D$. One can thus conclude that the spatial CR scattering observed in the strongly turbulent region generated due to the nonresonant instability is compatible with Bohm diffusion in the amplified magnetic field. Note, that similar conclusion was achieved for pitch-angle scattering of a mildly relativistic CR beam interacting with Bell’s turbulence (Niemiec et al., 2010).

3.2.3 Spectra of ambient plasma

We extract particle data for the four regions along the flow, that are marked in Fig. 6, and compute kinetic-energy spectra for electrons and ions at the end of the simulation run at time $t\gamma_{\rm max}=26.9$, that are shown in Figs. 11 and 12, respectively. The spectra are calculated in the local plasma rest frame, i.e., accounting for turbulent plasma motions.

The plasma particles are successively heated as they stream through the background of CRs, exciting the nonresonant instability that nonlinearly grows in the region of strong turbulence. Effective plasma heating is evident in the band $x/\lambda_{\rm max}=90-110$ (red lines, No. 3), and the electron spectra have largely relaxed to a Maxwellian in the area $x/\lambda_{\rm max}=160-180$ (yellow lines, No. 4). The ion spectra still clearly deviate from Maxwellians and show supra-thermal tails. The two species remain close to equipartition in bulk, but a substantial fraction of ions can achieve more than an order of magnitude higher energies than do electrons. This suggests that although the heating mechanism for both plasma species appears similar, either thermal relaxation of ions is slow or there exists another process that enables further energization of ions beyond the thermal pool.

The spectra calculated in the region of the compression wave ($x/\lambda_{\rm max}=90-110$, red lines, No. 3), through which the system transitions from small-amplitude fluctuations to highly nonlinear turbulence, show a complex structure, as expected. As discussed already in Section 3.1.5, plasma compression causes reflection of ions back to the upstream region of the compression zone (compare Fig. 6). The reflected population can be observed as a high-energy bump in the ion spectrum for the region $x/\lambda_{\rm max}=70-90$ (green line No. 2 in Fig. 12). These particles interact with the incoming plasma and cause a tail in the otherwise thermal ion distribution that covers the energy range $5\cdot 10^{-3}\ga\Gamma_{i}-1\la 10^{-1}$. The reflected particles reach far upstream and are still visible in the spectra for the area $x/\lambda_{\rm max}=10-30$ (blue line No. 1 in Fig. 12). The counterstreaming ions also influence the electrons that react with the formation of spectral tails (blue and green lines in Fig. 11).

3.2.4 Micro-physics of plasma heating

To investigate the micro-physics of plasma heating we analyze the turbulence structure in a small portion of the computational box located in the region of strongly nonlinear field fluctuations. Fig. 13 shows maps of the $B_{z}$ magnetic-field amplitude and the ion number density, overlaid on the in-plane magnetic and electric field lines, respectively (Fig. 13a-b). We also display maps of the mean kinetic energy of ions and electrons (Fig. 13c-d) that serve as proxy for their temperatures. The maps are calculated at time $t\gamma_{\rm max}=15$ and are representative of the strongly nonlinear stage of turbulence development. To aid the cross-correlation of features in the maps we add crosses as marks.

To be noted from the figure is the organization of the turbulent plasma into filaments and islands of enhanced density that surround plasma voids of very low density (Fig. 13b). Inspection of the plasma velocity fluctuations reveals that these structures are vortices that are filled with very strong perpendicular magnetic field ($B_{z}$) and wrapped in quasi-circular in-plane magnetic field ($B_{x}$ and $B_{y}$). These features are known from earlier studies of the nonresonant instability and have been explained as consequence of the $\mn@boldsymbol{j}_{\rm ret}\times\mn@boldsymbol{B}$ force, that accelerates the plasma away from the cavities, due to a local return current, $\mn@boldsymbol{j}_{\rm ret}$, that is unbalanced by the homogeneously distributed CRs (see, e.g., Niemiec et al., 2008). The ion temperature (Fig. 13c) is particularly high in regions of very low plasma density, an example of which are the cavity structures at $(x/\lambda_{\rm max},y/\lambda_{\rm max})\approx(125.5,16-20),(129,23)$ and $(130.5,20.5)$, as well as the filaments at $(x/\lambda_{\rm max},y/\lambda_{\rm max})\approx(127,19-21)$ and $(138.5,19-21)$. Likewise, the ions are cold in filaments and islands of high plasma density, for example at $(x/\lambda_{\rm max},y/\lambda_{\rm max})\approx(134.5,20-24)$. The electron temperature largely follows that of the ions, although small differences in the locations of the most efficient heating can be noticed (Fig. 13d). The anticorrelation between plasma density and temperature indicates that plasma heating does not result from adiabatic compression in turbulently moving plasma parcels. In fact, the extent of heating depends on the strength of the electric fields that are associated with the magnetic turbulence.

As discussed in Section 3.1, the vorticity present in the region of strongly nonlinear turbulence induces motional electric field in all directions. In particular, the $E_{x}$ field component is contributed by the terms $v_{y}B_{z}$ and $v_{z}B_{y}$. One can note, that the electric field is weak or disappears completely in most of the density filaments. This is because the magnetic-field amplitude is in most locations dominated by the $B_{z}$ component, whose change of sign between magnetic vortices leads to the cancellation of the motional electric fields. The plasma heating can therefore occur in neighboring regions with weaker magnetic field, in which the electric field is present and can accelerate low-density particles. Correspondingly, efficient heating arises in plasma cavities where large-amplitude electric field is present. Such regions can be clearly identified in Fig. 13 around $(x/\lambda_{\rm max},y/\lambda_{\rm max})\approx(125.5,19.5)$ or $(130.5,16.5)$.

The amount of heat acquired by electrons or ions at specific locations depends on the configuration of the electromagnetic fields and the local plasma flow. Motional electric field by definition vanishes in the local flow frame, and so particles that by their kinetic energy are indistinguishable from the bulk will follow the plasma flow without significant energization. Additional electric field is required to heat such particles. In Fig. 14 we display the charge density of the plasma and the corresponding electrostatic field in the same region as chosen for Fig. 13. The Poisson equation, ${\nabla\cdot}\mn@boldsymbol{E}=\frac{e}{\epsilon_{0}}(N_{\mathrm{i}}-N_{\mathrm{e}})$, provides virtually the same field as that calculated in the local plasma rest frame, but details are better visible because the local plasma frame does not need to be resolved at small scales. To be noted from the figure is that many but not all locations with high plasma temperature coincide with regions of large-amplitude electrostatic field. Even in the plasma cavities one finds electrostatic field, and the potential energy is transferred to few particles, leading to a high kinetic energy. The charge separation that provides electrostatic field must be driven, and plasma heating damps the fields. A direct comparison of field strength and plasma temperature is impeded by the electrostatic field being a snapshot and the temperature representing the heating history of a plasma parcel. As electrons are far more mobile along the magnetic field than are ions, heat is redistributed efficiently. Hence the spatial variation in electron temperature appears less dramatic than that of ions.

3.2.5 Plasma ion and electron scattering

The electrostatic heating applies to low-energy (thermal) electrons and ions that are tight to the bulk flow. A few particles, that manage to gain higher energies, are pulled out of the thermal pool and further accelerated. They sample the turbulent convective fields, undergoing stochastic scattering. The trajectory of such an ion is shown in Fig. 13, the details of its energization history are presented in Fig. 15. A particle traveling in the turbulent field is deflected by magnetic irregularities and scattered upon passage through regions of opposite polarity. The motional electric field associated with the magnetic turbulence then accelerates or decelerates a particle depending on its local direction of motion, which is demonstrated in Fig. 15. At time $t\gamma_{\rm max}=15$, for example, the ion is instantaneously accelerated on account of work in the $E_{x}$ and $E_{z}$ electric-field components that provide a net positive product $q(\mn@boldsymbol{v}\cdot\mn@boldsymbol{E})$ (see Fig. 15b-c). Shortly before that, at $t\gamma_{\rm max}=14.8$, the particle was slightly decelerated when it moved against the then dominant $E_{y}$ electric field component. The excellent correlation between the total rate of energization and that of work in the electric field suggests that all gains and losses in ion energy measured in the simulation arise from work in the turbulent electric field. To be noted from Fig. 15c is that while the work in the $E_{y}$ and $E_{z}$ fields averages to zero and would result in scattering in momentum space, the work in $E_{x}$ is on average positive, implying quasi-continuous gain in energy. This behaviour resembles that observed for CR ions. As noted in Section 3.2.1, CR ions tend to gain energy when moving in positive $x$-direction and to lose energy otherwise.

Scattering of supra-thermal electrons proceeds in a similar way as for the ions. However, due to their small mass compared to the ions, the electrons tend to experience locally strong acceleration or deceleration, and abrupt energy gains and losses can occur as a particle moves through different turbulence regions. Consequently, the work of all electric-field components can contribute. Beside scattering in the convective electric field, some electrons can be accelerated directly in the $E_{z}$ field at the magnetic-reconnection sites that may be formed between plasma parcels. Since these sites are very few, acceleration in the reconnection structures is much less efficient than that resulting from motional electric fields.

4 Discussion

Our new realistic setup with open boundaries permits us to study the spatial structure of a far-upstream SNR precursor region in which the nonresonant instability generates strong magnetic turbulence. One of the new features revealed in our study is the compression structure, reminiscent of a strong nonrelativistic supercritical shock with almost four-fold density increase. It arises from the difference in the flow speed between the incoming quasi-laminar and weakly-perturbed plasma and the plasma in the region of strong turbulence that has already slowed down due to backreaction effects. In comparison to the energy density and inertia of the background plasma, the energy density and bulk momentum of CRs in real shock precursors of young SNRs are likely a bit smaller than they are in this simulation. The upstream plasma will always suffer some deceleration, although that might be weaker in real precursors than is observed in this study. Also, the CRs will always be pushed back, and that may be the dominant effect for realistic energy-density ratios. Hence, the level of compression will clearly depend on the acceleration efficiency of the shock. In any case, Bell’s nonresonant mode is an efficient means to transfer bulk momentum from the CRs to the upstream plasma. Ignoring escape and assuming isotropy at the shock, the CR pressure in the precursor is $\Pi_{\mathrm{CR}}\approx U_{\mathrm{CR}}/3$, where $U_{\mathrm{CR}}$ is the cosmic-ray energy density. In the fluid picture the gradient of the CR pressure is compensated by the gradient of the plasma momentum flux, yielding for the total change in plasma flow velocity

[TABLE]

where $U_{\mathrm{up}}=\rho_{\mathrm{up}}\,V_{\mathrm{up}}^{2}/2$ is the far-upstream bulk-flow energy density of the unshocked plasma. For the last equality we used Eqs. 10-12 of Niemiec et al. (2008) to relate the energy-density ratio to the Alfvén speed upstream of the shock and the number of e-foldings that is available for linear growth of Bell’s instability before the plasma is captured by the shock. For standard parameters ($B_{0}\approx 10\ \mu$G and $n_{\mathrm{up}}\simeq 0.5\ \mathrm{cm}^{-3}$) we find $v_{\mathrm{A}}\simeq 30\ \mathrm{km/s}$, and reaching the saturation level of Bell’s mode should require at least the 15 e-foldings that it needs in our simulations. The expected velocity decrement of the upstream flow is then \delta V\approx 900\km/s, much larger than the sound speed, and so it is quite conceivable that also under realistic conditions a compression front develops in the shock precursor, provided Bell’s instability accounts for magnetic-field amplification to saturation level. In any case, the modification to the upstream flow would be large.

Electromagnetic turbulence is accompanied by significant heating of background electrons and ions. We show that this heating does not result from adiabatic compression in density filaments formed by turbulently moving plasma parcels, as was deduced from MHD studies. Instead, the highest plasma temperatures are reached in regions of very low plasma density that harbor large-amplitude electrostatic fields. We expect that plasma heating in the precursor has impact on the properties of the shock at which CRs are accelerated, and the system will likely approach a new balance between shock acceleration, driving of the nonresonant mode and its feedback in the precursor.

The scattering of CRs in space containing the turbulence amplified through the nonresonant instability is shown here to be compatible with Bohm diffusion, when the turbulence coherence scale is taken into account. As the precursor size scales with the mean free path for scattering,

[TABLE]

the system is about $25\,\lambda_{\mathrm{mfp}}$ in size, and CRs should reach the diffusive regime. Our estimate of the spatial diffusion rate being close to the Bohm limit is compatible with recent studies of nonresonant turbulence. Reville et al. (2008) investigated particle transport by following CR trajectories in a static snapshot of the magnetic field amplified through the nonresonant instability to nonlinear amplitude of $\delta B/B_{\parallel 0}\sim 5$, i.e., before the turbulence growth has started to slow down and saturate. They found the diffusion rate below the Bohm limit in the pre-existing uniform field for particles with an energy considerably lower than that of CRs driving the turbulence, while $k_{\rm max}r_{\rm CRg}\gtrsim 1$ is still satisfied. Reville & Bell (2013) ran coupled MHD-kinetic simulations that allowed for the CR filamentation instability and included the full CR dynamics, but regulated the driving to maintain an approximately constant CR current. They demonstrated sub-Bohm diffusion in the pre-existing field for CR driving the field growth, but estimated the diffusion coefficient to be comparable with the Bohm limit in the root-mean-squared magnetic-field amplitude. Similar results were reported for full shock hybrid-kinetic simulations of strong shocks by Caprioli & Spitkovsky (2014) for particles with energies close to the maximum energy, $E_{\rm max}$, propagating in strongly nonlinear turbulence generated via nonresonant instability by escaping CRs with $E\gtrsim E_{\rm max}$.

As a final note, the nonresonant instability, stimulating the magnetic field growth on scales small compared with the gyroradii of CRs driving it, should be accompanied by CR filamentation (Bell, 2005; Reville & Bell, 2012). The growth rate of the instability is proportional to the root mean square of the transverse component of the magnetic field, $\langle\delta B_{\perp}\rangle$, and so filamentation can efficiently grow only after the small-scale field has been significantly amplified. The scale above which filamentation dominates over the nonresonant instability can be estimated as $kr_{\rm CRg}=\langle\delta B_{\perp}/B_{\parallel 0}\rangle^{2}v_{\rm sh}/c$ (Reville & Bell, 2012). For our simulation the wavelength, $\lambda\simeq 10,000\Delta$, is comparable with the transverse size of the simulation box, and so CR filamentation cannot be captured. However, at late times we observe nonuniform structures in the CR number density, suggesting that the filamentation instability partially operates, although extended CR filamentary structures and the accompanying magnetic fields are not observed. Likewise, secondary MHD instabilities operating on scales larger than the CR Larmor radius cannot be described by our simulation (e.g., Bykov et al., 2011; Rogachevskii et al., 2012).

5 Summary

We present results of a 2.5-dimensional PIC simulation of a nonresonant cosmic-ray-current-driven instability that is believed to operate at parallel shocks of young supernova remnants and provide magnetic-field amplification there. Earlier numerical investigations of this instability used periodic boundary conditions that allow to follow only the temporal development of the system. They also did not account for mass conservation, and hence did not properly describe the nonlinear backreaction of the turbulence that involves bulk deceleration of CRs and background plasma. Our current study should be more realistic, as it for the first time employs open boundaries that permit plasma inflow on one side of the simulation box and outflow at the other end. This alleviates the continuity issue and also allows studying both the spatial and temporal evolution of the system. Our main results can be summarized as follows:

• •

We demonstrate magnetic-field amplification consistent with that found in studies with periodic simulation boxes. The evolutionary stages of the nonresonant instability, i.e., the first emergence of waves, their initial exponential growth, their attaining nonlinear wave amplitudes, their saturation and wave cascading, that have been observed in the earlier studies only as temporal evolution, can now be mapped to specific locations in the shock precursor, and the relation between the spatial and the temporal development becomes visible.

• •

We confirm the nonlinear backreaction of the amplified magnetic turbulence as modification of the bulk motion of CRs and plasma. Saturation of the instability is provided by an effective convergence of the bulk velocities of CRs and ambient plasma.

• •

The magnetic field is amplified through the nonresonant instability to $\langle\delta B/B_{\parallel 0}\rangle\approx 15$, i.e., essentially the same level as observed in our earlier simulations using periodic boxes (Stroman et al., 2009). Stronger amplification can be achieved through congestion in the region where incoming ISM plasma is decelerated to form a compression region.

• •

The formation of this compression region is revealed for the first time. The plasma density abruptly increases there by a factor of $4$, and a population of ions reflected from this shock-like structure appears. Streaming of the reflected ions against the ambient plasma flow may cause further magnetic-field generation via filamentation instabilities.

• •

In the nonlinear phase the magnetic turbulence is accompanied by strong density and velocity fluctuations. Turbulent electric field also appears at amplitudes an order of magnitude smaller than those of the magnetic field.

• •

Turbulent electric field mainly arises from magnetic-field transport, but there is also an electrostatic field component resulting from charge separation in the background plasma. On account of significant transverse plasma motions, an electric-field component in plasma drift direction with nonzero mean amplitude is present in the turbulence region. Bulk energy is transferred between CRs, the ambient plasma and turbulence through work in this field.

• •

The development of nonlinear turbulence is accompanied by significant heating of the background plasma. Ambient electrons and ions remain close to equipartition in bulk. While electrons are largely thermalized, the ions strongly deviate from a Maxwellian distribution and show supra-thermal spectral tails. The plasma heating is not adiabatic, is most efficient in low-density regions, and results from large-amplitude electrostatic fields.

• •

Turbulent electric field inelastically scatters CRs, introducing significant anisotropy and modifying their energy distribution. Stochastic scattering off the turbulent electric field is also responsible for the generation of supra-thermal tails in the spectra of plasma ions.

• •

Spatial of scattering of CRs in the turbulence generated via the nonresonant instability is compatible with Bohm diffusion, when corrected for the coherence scale of the turbulence.

In conclusion, our results support the notion that the nonresonant instability can produce significant magnetic-field amplification far upstream in the precursors of young SNR shocks. This process, when combined with additional field generation through instabilities that operate on larger scales, in close proximity to the shock front and/or at the shock itself, may be able to account for the level of magnetic-field amplification inferred from SNR observations. Similar conclusions may apply to shocks in relativistic jets and GRBs. The amplitude of the amplified field and CR scattering by the magnetic turbulence can both provide suitable conditions for particle acceleration to PeV energies at shocks in young SNRs. The backreaction of the magnetic-field amplification through Bell’s non-resonant mode involves plasma deceleration, compression, and heating, all of which will affect the sub-shock of the thermal plasma, at which CRs are accelerated. Continuous plasma deceleration is a known effect of CRs in shock precursors (Blandford & Eichler, 1987), but the shock-like compression feature that we observe in our simulation is a new effect whose impact on the shock structure and particle acceleration in general has yet to be explored.

Acknowledgements

This work has been supported by Narodowe Centrum Nauki through research projects DEC-2011/01/B/ST9/03183 (J.N.) and DEC-2013/10/E/ST9/00662 (O.K., J.N. and A.B.). M.P. acknowledges support through grant PO 1508/1-2 of the Deutsche Forschungsgemeinschaft. The numerical experiment was possible through a 7 Mcore-hour allocation on the 2.399 PFlop Prometheus system at ACC Cyfronet AGH. Part of the numerical work was conducted on resources provided by The North-German Supercomputing Alliance (HLRN) under project bbp00003 and also at the Pleiades facility at the NASA Advanced Supercomputing (NAS).