Monte Carlo modelling of particle acceleration in collisionless shocks with an effective mean electric field

TL;DR

This paper introduces a nonlinear Monte Carlo model for particle acceleration in collisionless shocks, incorporating an effective mean electric field and backreaction effects, to better understand cosmic ray acceleration in supernova remnants.

Contribution

It presents a novel Monte Carlo simulation that explicitly includes an effective mean electric field and backreaction, advancing the modeling of shock acceleration processes.

Findings

Demonstrates the impact of mean electric field on particle acceleration efficiency

Shows nonlinear feedback effects modify shock structure and acceleration dynamics

Provides insights into magnetic field amplification mechanisms in supernova remnants

Abstract

Relativistic particle acceleration in collisionless shocks of supernova remnants is accompanied by magnetic field amplification from cosmic ray (CR) driven plasma instabilities. Bell's fast CR-current instability is predicted to produce turbulence with a non-zero mean electric field in the shock precursor. We present a Monte Carlo model of Fermi shock acceleration explicitly taking into account an effective mean upstream electric field. Our model is nonlinear and includes the backreaction effects of efficient Fermi acceleration on the shock structure.

Tables

Table 1. Table 1: Results of the self-consistence calculations.

\br$u_{\mathrm{sh}}$ [km s-1]	$n_0$ [cm-3]	$\text{0}{R}_{\mathrm{tot}}$	\m$B_{\mathrm{eff},d}$[mG]	\m$\phi$ [TeV]
\mr\0\0\0 5000	\0\0 0.3	9.4 (9.2)	\0 0.63 (0.61)	\0\0\0 2.4
\0\0\0 5000	\0\0 1.0	8.9 (8.7)	\0\0 1.2 (1.1)	\0\0\0 3.9
\0\0 10000	\0\0 0.3	7.8 (7.7)	\0\0 1.4 (1.3)	\0\0 10.7
\0\0 10000	\0\0 1.0	7.5 (7.4)	\0\0 2.6 (2.5)	\0\0 15.8
\0\0 20000	\0\0 0.3	6.7 (6.7)	\0\0 3.1 (2.8)	\0\0 37.6
\0\0 20000	\0\0 1.0	6.6 (6.6)	\0\0 5.5 (5.1)	\0\0 51.6
\br

Equations

$$\lambda (x,p)=\frac{1}{{\lambda }_{B,\mathrm{st}}(x,p{)}^{-1}+{\lambda }_{\mathrm{ss}}(x,p{)}^{-1}},$$

$$v_{\mathrm{scat}}(x)=\frac{-{\int }_0^{\infty }\Gamma (x,k)W(x,k)dk}{d{P}_{\mathrm{cr}}(x)/dx}.$$

$$E_0(x)=\frac{-{\int }_0^{\infty }\Gamma (x,k)W(x,k)dk}{j_{\mathrm{cr}}(x)}.$$

$$E_{\mathrm{est}}=-v_{\mathrm{scat}}{B}_0/c.$$

Full text

To be published in Journal of Physics: Conference Series (in print, 2019).

1 Introduction

The converging plasmas on either side of a collisionless shock allow particle acceleration as particles scatter nearly elastically between the upstream and downstream regions. The scattering centers are magnetic field fluctuations moving with the background plasma and multiple crossing of the collisionless shock front can result in significant energy gains for individual particles.

This first-order Fermi mechanism [also called diffusive shock acceleration (DSA)] is expected to be efficient at the strong shocks in young supernova remnants (SNRs) (e.g. [1]). There is strong evidence to show that anisotropic CR distributions in the shock precursor leads to the development of plasma instabilities and magnetic field amplification (MFA) (e.g. [2]). Furthermore, if the acceleration is efficient, the gradient of the CR pressure will cause a modification of the bulk flow of the background plasma in the upstream region. The anisotropy of the CR distribution determines the growth of magnetic fluctuations in the precursor and the energy dependent scattering rate of fast particles. Since both of these depend critically on the bulk flow, DSA, with self-consistent MFA, is an essentially nonlinear problem which can only be fully addressed with computer simulations (see e.g. [3, 4]).

Several nonlinear models of efficient DSA in SNRs were constructed to account for the backreaction of accelerated particles on the dynamics of the background plasma which carries the scattering centers (see, for example, [5, 6, 7]). Generally, the kinetic description of fast particle transport and acceleration is based on macroscopic recipes of particle scattering by magnetic turbulence. Microscopic plasma instabilities amplify magnetic fluctuations which are described after an appropriate averaging procedure on the ensemble of scattering centers. In the case of the resonant CR streaming instability [8, 9], the scattering centers can move relative to the background plasma with some speed.

On the other hand, the short-scale modes produced by the non-resonant Bell instability, which grow very fast in the shock precursor, have zero or small phase speed in comparison with the Alfvén velocity [10]. However, it was shown (see e.g. [10, 11]) that there is a mean, large-scale electric field directed opposite to the CR current producing Bell’s instability. This electric field is directed opposite to the CR current and modifies the energy flux of Fermi accelerated particles in the shock precursor. The effect of the mean electric field on the spectrum of shock accelerated particles was discussed in [12] with a test-particle approach using the advection-diffusive equation with a mean electric field [13]. Within this test-particle model, it was shown that including the mean electric field resulted in a decrease in the CR proton acceleration efficiency. Here, we study the effect of the mean electric field on DSA with a nonlinear Monte Carlo model which was developed in [7]. We show below that nonlinear effects are essential for modeling CR ions with a self-consistent mean electric field. We don’t consider CR electrons here since their precursor pressure is normally much less than that of ions and doesn’t affect the flow dynamics.

2 Monte Carlo model of DSA with electric fields

We construct a steady-state model of a plane-parallel, nonrelativistic collisionless shock where the nonlinear shock structure is determined iteratively. The model convergence is determined by satisfying the conservation of mass, momentum and energy fluxes to some limit of accuracy. Particles diffuse by pitch-angle scattering off magnetic field fluctuations (see [7, 14]). The scattering is assumed to be isotropic and elastic in the local scattering center frame.

For CRs, the speed of the scattering centers, in the shock rest frame, is equal to the sum of the speed of the background plasma, $u(x)$, and the speed of the scattering centers relative to the background plasma $v_{\mathrm{scat}}(x)$. For the thermal background plasma, the scattering center speed is equal to the bulk plasma speed, $u(x)$ (see [15] for a full description of the pitch-angle diffusion process). An important property of the Monte Carlo technique is that particle injection is included self-consistently once the scattering assumptions are made. A particle is “injected” into the Fermi process if it is able to scatter back upstream having at least once crossed into the downstream region. Superthermal particles emerge smoothly from the background population and no artificial distinction between thermal particles and CRs is made. There is no fluid approximation in the Monte Carlo method; the entire plasma, thermal and CR, is described as particles moving with a single set of scattering assumptions.

We begin our simulation with a mean electric field, $\mathbf{E}_{0}(x)$, at each spatial grid point. As particles scatter during time $dt$ (short compared to a collision time), they experience a momentum change from this field $e\mathbf{E}_{0}dt$, where $e$ is the elementary charge (we consider the plasma consisting of protons only). This additional momentum change is included in our iterative scheme where we adjust the bulk flow speed $u(x)$, the background magnetic field turbulence, the overall compression ratio $R_{\mathrm{tot}}$, and the spatial dependence of the mean field $E_{0}(x)$, until mass, momentum, and energy fluxes are conserved throughout the shock. The momentum and energy fluxes are calculated from the particle distribution functions and the magnetic field amplification is calculated self-consistently from the resonant and Bell’s instabilities as in [7]. For simplicity, we do not include turbulent cascade here.

Ours is a plane-parallel model where all quantities depend only on the $x$-coordinate. The viscous shock is at $x=0$ and negative values of $x$ correspond to upstream values, i.e., the shock is directed along the $-x$-axis. The nonperturbed far upstream magnetic field, $B_{0}$, and the average electric field are both directed parallel to the shock normal. Particles leave the system by either convecting far downstream or by escaping at an upstream free escape boundary (FEB), $x_{\mathrm{FEB}}$. For our results here, $B_{0}=3\,\mu$G, $x_{\mathrm{FEB}}=-10^{8}\,r_{g0}$, $r_{g0}\equiv m_{p}u_{0}c/(eB_{0})$, $u_{0}=5000$ km s*-1*, $m_{p}$ is the proton mass, and $c$ is the speed of light.

The turbulent magnetic field is $\displaystyle B_{w}\left(x\right)=\sqrt{4\pi\int_{0}^{\infty}W\left(x,k\right)dk}$, where $W\left(x,k\right)$ is the spectral energy density of the turbulence for wavenumber $k$, $B_{w}(x_{\mathrm{FEB}})=B_{0}$, and $\displaystyle W(x_{\mathrm{FEB}},k)\sim k^{-5/3}$. The energy-containing scale is equal to 10 pc. The effective magnetic field is given by $\displaystyle B_{\mathrm{eff}}(x)=\sqrt{B_{w}^{2}\left(x\right)+B_{0}^{2}}$ (see [16] for a fuller discussion of the turbulence calculation).

The position- and momentum-dependent mean free path of particles in the local scattering center frame is

[TABLE]

where $\displaystyle\lambda_{B,\mathrm{st}}(x,p)=\frac{pc}{eB_{\mathrm{ls,st}}(x,k_{\mathrm{res}})}$, $\displaystyle\lambda_{\mathrm{ss}}(x,p)=\left(\frac{pc}{\pi e}\right)^{2}\left[{\int_{k_{\mathrm{res}}}^{\infty}\frac{W\left(x,k\right)}{k}dk}\right]^{-1}$, $\displaystyle B_{\mathrm{ls}}(x,k)=\sqrt{4\pi\int_{0}^{k}W(x,k)dk+B_{0}^{2}}$, $\displaystyle B_{\mathrm{ls,st}}(x,k)=\sqrt{4\pi\int_{0}^{k}W(x,k)dk}$, $p$ is the particle momentum in the scattering center frame, and $\displaystyle\frac{k_{\mathrm{res}}pc}{eB_{\mathrm{ls}}(x,k_{\mathrm{res}})}=1$.

The rate of increase of the turbulent energy flux due to the CR driven plasma instabilities is equal to $\displaystyle\int_{0}^{\infty}\Gamma\left(x,k\right)W\left(x,k\right)dk$ , where $\displaystyle\Gamma\left(x,k\right)$ is the growth-rate of the plasma instability (see [7]). The change of the CR energy flux $F_{\mathrm{cr}}(x)$ due to CR interactions with the microscopic background plasma fields can be expressed as: $\displaystyle\frac{dF_{\mathrm{cr}}(x)}{dx}=<\mathbf{J_{cr}}\mathbf{E}>$, where the averaging of the product of microscopic CR current $\mathbf{J_{cr}}$, and the electric field $\mathbf{E}$, is over scales larger than the CR scattering length. In the macroscopic approach used in [7, 12], as well as here, $\displaystyle<\mathbf{J_{cr}}\mathbf{E}>$ can be expressed in two alternative ways.

First, it can be presented as $\displaystyle<\mathbf{J_{cr}}\mathbf{E}>=\left[u\left(x\right)+v_{\mathrm{scat}}(x)\right]\frac{dP_{\mathrm{cr}}(x)}{dx}$, where $P_{\mathrm{cr}}(x)$ is the CR pressure and the scattering center speed $v_{\mathrm{scat}}$ can be derived self-consistently as described in [7]. Alternatively, the average can be represented as in [12]: $\displaystyle<\mathbf{J_{cr}}\mathbf{E}>=u\left(x\right)\frac{dP_{\mathrm{cr}}(x)}{dx}+j_{\mathrm{cr}}(x)E_{0}\left(x\right)$, where $j_{\mathrm{cr}}(x)$ is the mean CR electric current in the local scattering center frame. The first term on the right hand sides of these equations describes the adiabatic change of the CR energy flux.

The work done by CRs to increase the turbulent energy flux due to the CR-driven plasma instabilities is represented in the macroscopic iterative model in two alternative ways. One uses $v_{\mathrm{scat}}$ as a macroscopic parameter for the closure of the energy conservation equation, i.e.,

[TABLE]

The other method uses the mean effective electric field, i.e.,

[TABLE]

Next we compare these two approaches in the nonlinear Monte Carlo DSA model.

3 Results

We first perform test-particle calculations to compare the effects of an average electric field (Eq. 3) against one where the scattering center speed is determined with Eq. (2). In test-particle models, in general, there is no modification of the upstream flow, i.e., $u(x)=u_{\mathrm{sh}}$, where $u_{\mathrm{sh}}$ is the shock speed measured in the frame of the far upstream plasma. We further assume constant values of the average electric field, $E_{0}$, the upstream scattering center speed, $v_{\mathrm{scat}}$, and we ignore MFA. We estimate values for $E_{0}$ and $v_{\mathrm{scat}}$ that influence the CR spectrum in a similar fashion. For this purpose, we equate changes of CR energy fluxes, i.e., $j_{\mathrm{cr}}E_{0}=v_{\mathrm{scat}}dP_{\mathrm{cr}}(x)/dx,$ for models with $E_{0}\neq 0$ or $v_{\mathrm{scat}}\neq 0$. As definition of the CR current and pressure gradient, we use values from the model with $v_{\mathrm{scat}}\neq 0$ and $E_{0}=0$.

This model is equivalent to standard unmodified DSA [17] with the scattering center speed in the upstream $u^{u}_{sc}=u_{\mathrm{sh}}+v_{\mathrm{scat}}$ and in the downstream $u_{\mathrm{sc}}^{d}=u_{\mathrm{sh}}/R_{\mathrm{tot}}$. In this case, the upstream CR distribution function is $\displaystyle f(x,p)=f_{0}(p)\exp{\left[u_{\mathrm{sc}}^{u}x/D(p)\right]}$, where $D(p)=\lambda(p)v/3$, $v$ is the particle speed, and the upstream CR current is given by $\displaystyle j_{\mathrm{cr}}(x)=-eu_{\mathrm{sc}}^{u}\int{f(x,p)d^{3}p}$. By definition, $\displaystyle P_{\mathrm{cr}}(x)=\frac{1}{3}\int{vpf(x,p)d^{3}p}$, and using the expression for $f(x,p)$, we obtain $\displaystyle\frac{dP_{\mathrm{cr}}(x)}{dx}=u_{\mathrm{sc}}^{u}\int{\frac{p}{\lambda(p)}f(x,p)d^{3}p}$.

In our test-particle calculations we assume $\lambda(p)\approx r_{g}$, where $r_{g}=pc/(eB_{0})$ is the gyro-radius of a particle with momentum $p$ in a magnetic field $B_{0}$. We also assume that the CR current doesn’t change due to the action of the electric field. With these approximations, the estimated average electric field is

[TABLE]

The particle spectra obtained with our test-particle calculations are shown in the left panel of Fig. 1. In the model with $v_{\mathrm{scat}}\neq 0$ and $E_{0}=0$ the distribution function at superthermal energies is a power-law $f_{0}(p)\sim p^{-\alpha}$, where $\alpha=3R_{\mathrm{eff}}/(R_{\mathrm{eff}}-1)$ and $R_{\mathrm{eff}}$ is the shock compression ratio, i.e., $R_{\mathrm{eff}}=u_{\mathrm{sc}}^{u}/u_{\mathrm{sc}}^{d}$. The subscript on $R_{\mathrm{eff}}$ indicates “effective”.

For our self-consistent, nonlinear models, we use two values for the background ambient number density, $n_{0}$, and three values for the shock speed, $u_{\mathrm{sh}}$ (see Table 1). All other parameters are the same except for the two closure models indicated by eqs. (2) and (3). The parameters are typical for forward shocks in SNRs interacting with the interstellar medium. As indicated by the red and blue curves in the right panel of Fig. 1, these two closure models give nearly identical results for our nonlinear shocks. The nonlinear profiles for $E_{0}(x)$ and $v_{\mathrm{scat}}(x)$ are presented in Fig. 2. Other quantities are given in Table 1 [values in parentheses indicate $v_{\mathrm{scat}}(x)\neq 0$], where $R_{\mathrm{tot}}$ is the full compression, $B_{\mathrm{eff,d}}$ is the effective downstream field, and $\varphi$ is the energy a proton acquires by traversing the upstream region through the average electric field, i.e., $\varphi=e\int_{x_{\mathrm{FEB}}}^{0}E_{0}(x)dx$.

4 Summary

The self-generation of magnetic turbulence is an essential part of nonlinear Fermi acceleration. However, the production and amplification of turbulence is a complicated process, with an extremely wide dynamic range if SNRs are to be modeled. While particle-in-cell simulations give important information on small scales, none are yet large enough to give a consistent picture of Fermi acceleration and MFA for shocks typical of those seen in young SNRs. With important approximations, the Monte Carlo techniques we discuss can model critical aspects of nonlinear Fermi acceleration in shocks with large dynamic ranges. We have generalized previous results to include descriptions of the speed of scattering centers through the self-generated turbulence and the dynamic effects of non-zero local electric fields generated by Bell’s instability.

While further study is required to determine which, if either, of these effects dominate the Fermi acceleration of ions, we have shown that both processes, expressed in Eqs. (2) and (3), influence the accelerated proton spectrum in very similar ways (right panel in Fig. 1). It is also important to recognize that radiation from SNRs can’t be modeled properly unless electrons are included consistently (see, for example, [18]). As noted by [12], the fact that electrons have the opposite sign from protons may be important for the processes discussed here. However, since electrons carry much less momentum than protons, they will have little influence on the shock dynamics and we don’t include them in our simulations.

In general, an average electric field $E_{0}(x)>0$ removes energy from protons and a scattering center speed $v_{\mathrm{scat}}(x)<0$, i.e., one where the scattering centers move away from the viscous subshock in the upstream plasma frame, results in a weaker effective compression ratio. While both of these effects tend to soften the proton spectrum, the nonlinear shock modification, along with the self-consistent injection included in the Monte Carlo model, still produce a concave spectral shape with a slope harder than $p^{-4}$ before the cutoff at the highest energies (red and blue curves in the right panel of Fig. 2). We note that $E_{0}(x)>0$ will tend to harden an electron spectrum while $v_{\mathrm{scat}}(x)<0$ will tend to soften spectra regardless of sign.

\ack

The authors acknowledge support from RAS Presidium program No. 12. The results of this work were obtained using computational resources of the Peter the Great Saint-Petersburg Polytechnic University Supercomputing Center (http://www.spbstu.ru).