Bell Instability–mediated Diffusive Shock Acceleration at Supernova Blast Wave Shock Propagating in the Interstellar Medium

Supernova blast wave shock is a very important site of cosmic-ray acceleration. However, the detailed physical process of acceleration, in particular, nonlinear interplay between cosmic-ray streaming and magnetic field amplification, has not been studied under a realistic environment. In this paper, using a unique and novel numerical method, we study cosmic-ray acceleration at supernova blast wave shock propagating in the interstellar medium with well-resolved magnetic field amplification by nonresonant hybrid instability (or Bell instability). We find that the magnetic field is mildly amplified under typical interstellar medium conditions that leads to maximum cosmic-ray energy ≃30 TeV for supernova remnants with age ≃1000 yr consistent with gamma-ray observations. The strength of the amplified magnetic field does not reach the so-called saturation level because the cosmic-ray electric current toward the shock upstream has a finite spatial extent, by which Bell instability cannot experience many e-folding times.


Introduction
Supernova blast wave shock is known to be an active site of particle acceleration (Koyama et al. 1995;Aharonian et al. 2008).It is widely accepted that the diffusive shock acceleration (DSA) mechanism accounts for the acceleration of charged particles with power-law momentum distribution (Drury 1983;Blandford & Eichler 1987).Since the rate of galactic supernovae can provide sufficient energy density for cosmic rays (CRs) below the so-called knee energy (;3 PeV), it has been believed that historical young supernova remnants (SNRs) with age 1000 yr can be PeVatrons.In order for the young SNRs to be PeVatrons, magnetic field must be amplified at least by 2 orders of magnitude in upstream of the blast wave shock.Nonresonant hybrid instability (NRHI or Bell instability) is known to be an effective mechanism of the magnetic field amplification (Bell 2004;Marcowith et al. 2016).Near the SNR forward shock, a stream of accelerated CR nuclei (mostly protons) makes an electric current toward upstream that induces a return current of thermal electrons to keep charge neutrality.The induced return current exerts the Lorentz force on background plasma that amplifies the magnetic field by exciting circularly polarized Alfvén-like waves.
According to particle-in-cell (PIC) simulations by Caprioli & Spitkovsky (2014a), the upstream magnetic field can be amplified to ~MA 1 2 times the initial strength by the NRHI, where M A is the Alfvénic mach number of the shock.However, due to the limited scale range of first-principle simulations, the maximum M A studied by the PIC simulations is only M A = 100, which leads only by an order of magnitude amplification, and we need a case study of realistic M A  1000 shock wave.
In order to examine realistic growth of the NRHI, we need to know the spatial structure of the upstream CR current, which is quite hard to compute accurately because the CR current is composed of escaping CR flux whose typical energy is time dependent.Hence, many authors (sometimes implicitly) assumed the structure of the CR current and discussed the expected level of the upstream magnetic field (Zirakashvili & Ptuskin 2008;Zirakashvili et al. 2008;Bell et al. 2013;Schure & Bell 2013).
Observations of young SNRs with ages ∼1000 yr via synchrotron emission show that the magnetic field downstream of the shock can be amplified to the level of ∼100 μG or more, which is clearly stronger than that expected from a shock compression of the interstellar magnetic field, but not enough to make these SNRs PeVatrons (Vink & Laming 2003;Ballet 2006;Parizot et al. 2006).However, there is still a possibility of SNRs being PeVatrons.It has been discussed that the maximum CR energy can exceed 1 PeV when the blast wave shock propagates in a dense circumstellar medium (CSM) created by a stellar wind of a progenitor massive star (Schure & Bell 2013;Marcowith et al. 2018) although the upstream magnetic field amplification by the NRHI is still necessary.
Recently, Inoue (2019) developed a novel numerical method that enables self-consistent treatment of the CR current evolution and the resulting magnetic field amplification by the NRHI.Using such a numerical method, Inoue et al. (2021) demonstrated that the blast wave shock propagating in a dense CSM can indeed accelerate particles up to the so-called knee energy ∼3 PeV if the progenitor mass-loss rate before the supernova explosion is high enough.In their study, it is reported that the NRHI can amplify the magnetic field by an order of magnitude from the initial strength, but the final strength immediately before the shock is below the so-called saturation level expected from previous theoretical studies.They discussed that the moderate amplification by the NRHI is Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
accounted for by a small e-folding number (∼2-3) due to the limited spatial extent of the CR current.
If the moderate amplification of the NRHI is a common feature of the SNR blast wave shock even in the later interstellar medium (ISM) propagating phase, it could explain the observed magnetic field strengths of the young SNRs (∼100 μG in downstream) and explain why they are not PeVatrons.Thus, in this paper, we examine the NRHImediated CR acceleration at the SNR blast wave shock propagating in the ISM.
The paper is organized as follows: in Section 2, we provide the basic equations and numerical settings for simulations.The results of the simulations and their physical interpretation are shown in Section 3. In Section 4, we summarize the paper and discuss the implications.

Basic Equations
Similar to Inoue et al. (2021), we solve a hybrid system of the Bell magnetohydrodynamics (MHD) equations and a telegrapher-type diffusion convection equations (or Vlasov-Fokker-Planck equation for the isotropy and first-order anisotropy of CR distribution function) in the polar coordinate around θ ∼ π/2 (see also Bell et al. 2013;Inoue 2019).The system equations can be broken down into the MHD part: and the CR kinetic part: where j r ret ( ) is the return current density induced by the CR streaming current ( =j j r r ret cr ( ) ( ) ), Q inj is the injection rate, and κ(B) is the diffusion coefficient.F 0 ≡ f 0 p 3 and F 1 ≡ f 1 p 3 are, respectively, the isotropic and the first-order anisotropic (drift anisotropy) components of the CR distribution function f (r, p) = f 0 (r, p) + (p r /p)f 1 (r, p).Although they are small corrections, we adequately consider source terms ∝1/r, which are neglected in Inoue et al. (2021).
The detailed form of j r cr ( ) , Q inj , and κ(B) are given by Equations (6)-( 9) of Inoue et al. (2021).Here, we briefly state how they are described: The CR current j r cr ( ) is computed by momentum space integration of f 1 (p).The injection rate Q inj has a nonzero value only at the shock front, and the rate of injection is determined so that the fraction η of the thermal particles go into the DSA process (Blasi et al. 2005).We employ the following diffusion coefficient (Skilling 1975;Caprioli & Spitkovsky 2014b): 2 and v CR is the CR velocity at momentum p.Under small amplitude magnetic turbulence (δB B r ), Equation (11) gives the diffusion coefficient due to pitch-angle scattering, while it becomes the Bohm limit coefficient under the amplified field strength for δB > B r .This type of diffusion coefficient is supported by the results of PIC simulations (Caprioli & Spitkovsky 2014b) and test particle transport calculation in a super-Alfvénic turbulence (Roh et al. 2016).Note that, since the scale of magnetic field fluctuations induced by the NRHI is roughly 1-2 orders of magnitude smaller than the gyroscale of the maximum energy CRs, the above diffusion coefficient could be underestimated.However, it is pointed out that a filamentation instability that grows simultaneously with the NRHI helps to make the fluctuation scale larger (Reville & Bell 2012), suggesting the use of Equation (11) is reasonable.
Following the method implemented by Inoue et al. (2021), we numerically treat the momentum space from p = 100 GeV c −1 , below which we assume the standard DSA spectrum f 0 (r sh ) ∝ p −4 from the injection momentum.Thus, our method cannot describe the effect of shock structure modification by CR pressure, which is carried by CRs with p ∼ 1 GeV.In this paper, we study the cases with the injection rate η = 0.1-1.0 × 10 −4 in which the effect of CR pressure is not substantial.According to a theoretical modeling of multiwavelength emission of a historical young SNR SN1006 by Berezhko et al. (2012), η ; 10 −4 gives the best fit.Note that studying the effect of the CR pressure under a larger injection rate is still important because the injection rate can be locally enhanced, and we will further study it in our future papers.

Fluid Initial Conditions
The inner and outer boundaries of the numerical domain are located at r in = 0.5 pc and r out = 3 × 10 19 cm ;10 pc for most models.To induce a blast wave shock, we initially set an ejecta of radius R 0 = 1.0 pc, velocity v 0 = 10 4 (r/R 0 ) km s −1 , and density ρ 0 = 2.1 × 10 −23 g cm −3 .This leads the total kinetic energy of the ejecta to be E kin = 0.76 × 10 51 erg.
The initial ISM is distributed in R 0 < r < r out with uniform density ρ ISM = m n ISM and temperature T ISM = 10 6 K, where m = 1.27 m p is the mean gas particle mass.To study the impact of the ISM density on CR acceleration, we examine three types of the ISM density n ISM = 0.5, 0.1, and 0.02 cm −3 .For the lowest ISM density model (n ISM = 0.02 cm −3 ), we use r out = 4 × 10 19 cm because the shock propagates farther than the larger density models.Given that the ISM temperature T ISM = 10 6 K, the ISM sound speed is calculated to be c s ; 100 km s −1 .Thus, the Mach number of the forward shock becomes ∼v 0 /c s ; 100 for the early free expansion stage.Note that our choice of the ISM temperature is realistic for the lower ISM density models of n ISM = 0.1 and 0.02 cm −3 , but for the n ISM = 0.5 cm −3 case, the realistic temperature would be around 10 4 K, which leads to a sonic Mach number around 1000.In general, it is numerically hard to stably follow the dynamics of the shock with Mach number 1000. 4 Thus, we choose T ISM = 10 6 K even for the n ISM = 0.5 cm −3 model.Since the upstream temperature does not influence the NRHI dynamics, the result seems to be unchanged even if we set a lower ISM temperature.
We assume the constant radial component of the magnetic field B r = 3 μG for most models. 5The upstream Alfvén velocity is calculated to be For B θ and B f , we set turbulent fluctuations by superposing Alfvén waves with flat power spectrum for most runs, and we also examine the 0.01 case to study the effect of seed amplitude on the NRHI growth.We show a summary of our numerical model parameters in Table 1.

Boundary Conditions
For the boundary conditions of the MHD part of the equations, we use the free boundary conditions, except we set v = 0 and B θ = B f = 0 at the inner boundary.This set of boundary conditions reproduces the blast wave evolution starting from the free expansion phase (v sh ∝ t 0 ) and then shifts to the Sedov-Taylor phase (v sh ∝ t −3/5 ; see Figure 1  below).
For the spatial boundaries of the CR part, we assume that CRs do not penetrate into the ejecta by setting κ = 0 in the ejecta.At the outer spatial boundary, outgoing free boundary conditions are imposed: F 0 (p, r out + Δr) = F 0 (p, r out ) and . For the momentum space, as stated in Section 2.1, we assume f 0 ∝ p −4 below the boundary at p = 100 GeV c −1 (see Equation (15) in Inoue et al. 2021 for implementation), while for f 1 , we impose null values outside the numerical domain of [100 GeV c −1 , 10 PeV c −1 ].

Numerical Resolution
Since capturing the NRHI growth is the major purpose of this study, we need to resolve the most unstable scale of the instability (Bell 2004): l = c B j   To satisfy the above condition, we resolve the whole spatial numerical domain by N cell = 2 20 ; 10 6 cells, leading to the spatial resolution of Δr = (r out − r in )/N cell ; 3 × 10 13 cm.Thus, the ratio of the instability scale and the resolution is estimated to be λ B /Δr ; 30, which satisfies the resolution requirement by Inoue et al. (2021).
For the momentum space, we explicitly treat a range of .Because the CR current is mostly composed of CRs with E  10 TeV (see Figure 3 below), as long as the lower boundary of the momentum space is set below 1 TeV/c, the result will not be affected by the boundary.This resolution is tested to be enough for numerical convergence (see Section 2.6 of Inoue et al. 2021).

Advantages of Our Numerical Scheme
The detailed numerical method to solve the basic Equations (1)-( 11) is described in Section 2.5 of our previous paper (Inoue et al. 2021).Here we briefly state the advantages of our method: if we employ the usual diffusion equation and solve it explicitly, a required time step for numerical stability becomes Δt 0.5 Δr 2 /κ(B r , p U ) ∼ 10 −3 s for the initial state, indicating a 1000 yr scale integration is impossible even with the largestscale supercomputer.
On the other hand, our hyperbolic-type basic equations require s, which is still short, but the 1000 yr integration is executable with a modern supercomputer.Here, c 3 is the free-streaming velocity of CRs, and C CFL is the Courant-Friedrichs-Lewy number that is set to be 0.8.
In addition to the alleviated time step, the hyperbolic equation generally has good compatibility with parallel computing, enhancing the advantage of our method.One may think of implicit schemes or super-time-stepping methods for integration, which allow large time steps.However, the implicit scheme is incompatible with parallel computing and has insufficient accuracy to recover the DSA spectrum.Supertime-stepping can only enlarge Δt by roughly an order of magnitude if we want to keep the appropriate CR spectrum.Thus, it cannot be used under realistic situations considered in this paper.

Result of Models 1-3
Panel (a) of Figure 1 shows the density structures of the result of Model 1 (fiducial model), where we only plot outside of the contact discontinuity.We see that the density structure around the shock becomes very noisy.This is due to the back reaction of magnetic field amplification.Panel (b) shows the forward-shock velocity evolution, where the shock position is defined as the outermost position of v r > 1000 km s −1 .We can confirm that the shock evolves from the initial free expansion phase (v sh ∝ t 0 ) to the Sedov-Taylor phase (v sh ∝ t −3/5 ).As expected, the larger the initial ISM density, the faster the shock attenuates.
Figure 2 shows the structure of the magnetic field strength around the shock at t = 100, 500, and 1000 yr for Model 1 (panel (a)), Model 2 (panel (b)), and Model 3 (panel (c)).Models 1-3 have similar levels of the magnetic field on the order of tens of microgauss immediately ahead of the shock and a 100 μG behind the shock.To make a more detailed comparison, in panel (d) of Figure 2, we plot the field strength of the results of Models 1-3 when the shock passes r = 5 pc.In addition, in Table 1, we list values of so-called ò B , i.e., the ratio of downstream magnetic energy density to upstream kinetic energy, where the downstream magnetic energy is averaged over 1000 Δx region from the shock front at t = 1000 yr.We see that the degree of the amplification increases with the ISM density.This is naturally accounted for by a larger total number of CRs, which is proportional to the upstream density and thus stronger CR current.
To see how magnetic fields grow in the upstream region, we plot the structure of the CR current density in the top panel of Figure 3.We also show the upstream CR current spectrum = j dj d p ln p r CR ( ) at t = 1000 yr for Model 1 in the bottom panel.We see that CRs with energy pc ∼ 10 TeV mostly For t 0  500 yr, from the Sedov-Taylor solution, the shock velocity is given by = Using t ad , the expected e-folding number of the NRHI at distance l is computed as In Figure 4, we plot σ at t 0 = 500 and 1000 yr for Models 1, 2, and 3.An important overall feature is that, as reported by Inoue et al. (2021), the NRHI has a small e-folding number (1.5-2.5) that can amplify the seed fluctuations only by roughly a factor e afew ∼ 10.
The expected e-folding number becomes smaller as ISM density decreases (as model number increases from 1 to 3).This seems to be reasonable behavior because, at fixed time, the total number of injected particles is a decreasing function of the initial ISM density that leads to smaller escape particles for smaller-ISM-density cases.
It has been expected that, if the e-folding number is large enough, the NRHI grows until magnetic field strength reaches the so-called saturation level, which happens once the gyroradius of maximum energy CRs becomes smaller than the NRHI critical scale (Bell 2004): where p max is the maximum (escaping) CR momentum.The magnetic field strengths in the simulations do not reach the above saturation level.This result is very similar to the case of an earlier phase of SNR studied by Inoue et al. (2021).In panels (a)-(c) of Figure 2, we plot the saturation amplitude (Equation ( 18)) estimated using local j r CR ( ) and the fitted E cut for t = 1000 yr data as dashed lines.Previous simple theoretical modelings assumed a sufficient e-folding number for the NRHI to saturate.However, if we take into account the realistic spatial extent of the CR current, it does not reach saturation in the SNR environment.
Figure 5 shows the CR density spectrum (n p = f 0 p 2 ) of Model 1 immediately behind the shock front.We can reasonably fit the spectrum by a power-law function with the standard DSA spectral index and exponential cutoff: In Table 1, we list the cutoff energies E cut at t = 1000 yr for all models as a result of the fitting.The resulting cutoff energies become E cut ∼ 10 TeV at t = 1000 yr for all models, except for the no-magnetic-field-amplification model (Model 8).Since the larger-ISM-density model (Model 1) has slightly stronger magnetic field strength compared to Model 2 and 3 but has smaller shock velocities, these complementary effects seem to make similar maximum energies between Model 1 and Model 3. The numerical cutoff energy is roughly consistent with the maximum energy of the age-limited DSA acceleration theory under the amplified magnetic field strength (Blandford & Eichler 1987): 10 G 3000 km s 1000 yr . 19 We can also understand the maximum energy of Model 8, which artificially suppresses the NRHI if we substitute unamplified magnetic field strength of 3 μG with ξ B = 0.1.

Influence of Seed Amplitude and Initial Field Strength
In Model 4, we run the case that has a smaller initial magnetic field fluctuation amplitude (ξ B = 0.01) than the fiducial model.We also study the case of stronger initial field strength (B r = 10 μG) in Model 5.In the top panel of Figure 6, we compare the magnetic field structure around the shock for Models 1, 4, and 5 at t = 1000 yr.We see that the final field strength does not substantially depend on the initial seed amplitude and field strength.This is because the small upstream fluctuations allow more CR flux that results in active NRHI growth in Model 4, while the stronger initial field in Model 5 has a negative effect on the CR stream.Indeed, from the bottom panel of Figure 6, which is the same as Figure 4 but for Models 1, 4, and 5, we can confirm that the results of the smaller initial seed model show a larger e-folding number, and   The results of our self-consistent calculation of CR current structure indicate that the e-folding number depends highly on the degree of upstream magnetic field fluctuations, but the final magnetic field strength in the immediate shock upstream takes a similar value in a realistic range of parameters explored in this paper.The insensitivity of the magnetic field strength around the shock on the upstream seed fluctuations seems reasonable if we notice that total electric charge escape to the upstream is the most crucial quantity, which is constant between Models 1, 4, and 5, owing to the common injection rate and upstream density.As a consequence of the similar magnetic field strengths, the cutoff CR energies have similar values in Models 1, 4, and 5 (see Table 1).

Effect of the Injection Rate
Since the CR current has a substantial impact on the growth rate of NRHI, the injection rate η could have a strong influence on the magnetic field amplification and the resulting maximum CR energy.Figure 7 shows the structures of magnetic field strength (top) and CR current (bottom) around the shock for Models 1, 6, and 7 at t = 1000 yr.We see that, although the amplification is more effective as η increases, the influence of η variation is not so striking even though there is an order of magnitude difference of η between Model 1 and Model 7.This somewhat similar magnetic field amplification level results in similar cutoff energies (only a factor of 3 difference) between Models 1, 6, and 7 (see Table 1).
This small dependence on η can be simply explained as follows: In the very early stage, the larger η model produces a stronger CR current that induces active NRHI and boosts up E max faster than the smaller η model.However, since the CR current is composed of CRs around E max , it is weakened as E max enlarges.Because CR density at shock obeys µ - dn p dp p 2 , a factor of 2 enhancement of cutoff energy leads to a quarter of CR density around the cutoff energy and thus reduces the escaping CR current.As a consequence, in a later stage, the strength of the CR current becomes similar in level in all Models 1, 6, and 7, resulting in not substantially different magnetic field strengths.Figure 8 shows CR current spectra for Model 6 (top) and Model 7 (bottom) from which we see that the CR current in Model 6 is composed of higherenergy CRs than those of Model 7.

Summary and Discussion
have studied CR acceleration at a supernova blast wave shock propagating in the ISM under the influence of NRHI.Our method handles nonlinear interplay between CR acceleration, CR streaming, and magnetic field amplification with sufficient resolution.Our findings can be summarized as follows:  1.The NRHI is indeed effective at the young SNR forward shocks with age around 1000 yr irrespective of the ISM density.Downstream ò B (the density ratio of upstream kinetic energy and downstream magnetic energy) can be 1%-10%, while it becomes 0.002% if the effect of NRHI is switched off.2. The level of the NRHI amplification does not reach the so-called saturation level, at which the gyroradius of escaping CRs becomes smaller than the critical length of the NRHI.This is because the realistic spatial extent of the upstream CR current is finite, while previous theories assumed sufficiently large extent to achieve the saturation.3. The amplified strength of the magnetic field at t = 1000 yr does not substantially depend on the ISM density because smaller ISM density leads to faster shock propagation that compensates the CR current strength.4. The magnetic field strength at immediate shock upstream does not substantially depend on the initial field strength and the seed fluctuation amplitude.This is because smaller initial strength and/or seed fluctuations lead to more active amplification via more spatially extended CR current (due to larger CR diffusion coefficient).5.Although larger injection rate models show more effective amplification and larger resulting maximum energy of CRs, they show weaker dependence on the injection rate than the linear relation (an order of magnitude change of the injection rate results in only a factor of around 3 difference in the maximum energy).This seems to stem from the fact that the escaping CR number (current strength) is a decreasing function of the maximum energy.6.The abovementioned somewhat stable nature of the NRHI amplification against the environmental conditions always leads similar maximum CR energy of ∼30 TeV at t = 1000 yr.This is a somewhat unexpected result but consistent with recent observations of young SNRs with t age  1000 yr (Suzuki et al. 2022).
Very recently, Tibet ASγ experiment and Large High Altitude Air Shower Observatory found PeV/sub-PeV gamma rays potentially from some young and middle-aged SNRs (t age ∼ 10 3 -10 4 yr; Aharonian et al. 2021;Tibet ASγ Collaboration et al. 2021).Given the results of the present study, these ultra-high-energy gamma rays seem to be not from SNR shells but from escaped CRs that are accelerated at a much earlier phase of the SNR or from other nearby high-energy objects.However, since we have only studied the cases of moderate CR injection rate (η ∼ 10 −4 at maximum), there may be an SNR population that has a higher injection rate by a yetunknown effect.In such a higher-injection-rate model, nonnegligible CR pressure would lead to a shallower CR spectrum and thus more efficient growth of the NRHI, which will be studied in our future works.
The critical scale below which the NRHI is stabilized by magnetic tension force is given by λ B /2.The results of our fiducial model (Model 1) show that, although it depends highly on time and distance from the shock front, the CR current density roughly takes 1 cm −2 near the shock front.Using this current density, the most unstable spatial and timescales of the NRH instability can be estimated as downstream magnetic energy density to upstream kinetic energy.The downstream magnetic energy is averaged over the 1000 Δx region from the shock front.results of test simulations byInoue et al. (2021), the numerical dispersion relation can reproduce the analytic one if λ B is resolved by more than 32 numerical cells.
64 numerical cells in the logarithmic scale, i.e.,

Figure 1 .
Figure 1.Panel (a): density structure for Model 1 (fiducial model).Different colors show different snapshot times.Panel (b): forward-shock position (circles) and shock velocity divided by 1000 km s −1 (triangles) for three different models.Thin red line is a reference line to the Sedov-Taylor solution: ∝t −3/5 .

Figure 2 .
Figure 2. Panel (a): structure of the magnetic field strength around the shock at t = 100, 500, and 1000 yr for Model 1. Panel (b): same as panel (a) but for Model 2. Panel (c): same as panel (a) but for Model 3. Panel (d): magnetic field strength of the results of Models 1-3 when the shock passes r = 5 pc.Dashed lines in panels (a)-(c) represent the saturation amplitude of NRHI given by Equation (18) for t = 1000 yr.

Figure 3 .
Figure 3. Panel (a): structure of the CR current densities for Models 1, 2, and 3 at t = 500 (thin and 1000 yr (thick lines).Panel (b): upstream CR current spectrum = j r p dj d p , l n p r

Figure 5 .
Figure 5. CR density spectrum (n p = f 0 p 2 ) of Model 1 immediately behind the shock.Black line represents the DSA spectrum (∝p −2 ).We can reasonably fit the spectrum by a power-law function with the standard DSA spectral index and exponential cutoff: µ -n p p cE exp p 2 cut ( ).The result of fitting for t = 1000 yr is shown as a dotted line (E cut = 31 TeV).

Figure 8 .
Figure 8. Panel (a): upstream CR current spectrum = j r p dj d p , l n p r