<figure><img src="image_1.png" /></figure>

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

Tsuyoshi Inoue<fnref n="1" />, Alexandre Marcowith<fnref n="2" />, and Gwenael Giacinti<fnref n="3" />

## Abstract

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.

Unified Astronomy Thesaurus concepts: Supernova remnants (1667); Cosmic rays (329)

## 1. 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 ~$M_A^{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 \gtrsim 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

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Physics, Konan University, Okamoto 8-9-1, Higashinada-ku, Kobe 658-8501, Japan; tsuyoshi.inoue@konan-u.ac.jp</footnote>

<footnote n="2">²Laboratoire Universe et Particules de Montpellier (LUPM) Université Montpellier, CNRS/IN2P3, CC72, place Eugène Bataillon, 34095, Montpellier Cedex 5, France</footnote>

<footnote n="3">³Tsung-Dao Lee Institute and School of Physics and Astronomy, Shanghai 200240, People's Republic of China
Received 2023 November 4; revised 2024 February 22; accepted 2024 March 5; published 2024 April 12</footnote>
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 NRHI-mediated 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.

## 2. Basic Equations and Numerical Setups

### 2.1. 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 $\theta \sim \pi/2$ (see also Bell et al. 2013; Inoue 2019). The system equations can be broken down into the MHD part:

$$ \frac{\partial \rho}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} r^2 (\rho v_r) = 0, \quad (1) $$

$$ \begin{aligned} \frac{\partial}{\partial t}(\rho v_r) + \frac{1}{r^2} \frac{\partial}{\partial r} r^2 \left( \rho v_r^2 + p + \frac{B_{\theta}^2 + B_{\phi}^2}{8 \pi} \right) \\ = \frac{1}{r} \left\{ \rho(v_{\theta}^2 + v_{\phi}^2) + 2p + \frac{B_r^2}{4\pi} \right\}, \end{aligned} \quad (2) $$

$$ \begin{aligned} \frac{\partial}{\partial t}(\rho v_{\theta}) + \frac{1}{r^2} \frac{\partial}{\partial r} r^2 \left( \rho v_r v_{\theta} - \frac{B_r B_{\theta}}{4\pi} \right) \\ = -\frac{1}{c} j_r^{(\text{ret})} B_{\phi} + \frac{1}{r} \left( \frac{B_r B_{\theta}}{4\pi} - \rho v_r v_{\theta} \right), \end{aligned} \quad (3) $$

$$ \begin{aligned} \frac{\partial}{\partial t}(\rho v_{\phi}) + \frac{1}{r^2} \frac{\partial}{\partial r} r^2 \left( \rho v_r v_{\phi} - \frac{B_r B_{\phi}}{4\pi} \right) \\ = \frac{1}{c} j_r^{(\text{ret})} B_{\theta} + \frac{1}{r} \left( \frac{B_r B_{\phi}}{4\pi} - \rho v_r v_{\phi} \right), \end{aligned} \quad (4) $$

$$ \begin{aligned} \frac{\partial \epsilon}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} r^2 \left\{ v_r \left( \epsilon + p + \frac{B^2}{8\pi} \right) - B_r \frac{\mathbf{B} \cdot \mathbf{v}}{4\pi} \right\} \\ = -\frac{j_r^{(\text{ret})}}{c} (\mathbf{v} \times \mathbf{B})_r, \end{aligned} \quad (5) $$

$$ \epsilon = \frac{p}{\gamma - 1} + \frac{1}{2} \rho v^2 + \frac{B^2}{8\pi}, \quad (6) $$

$$ \frac{\partial B_{\theta}}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r} r^2 (B_r v_{\theta} - B_{\theta} v_r) + \frac{1}{r} (v_r B_{\theta} - v_{\theta} B_r), \quad (7) $$

$$ \frac{\partial B_{\phi}}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r} r^2 (B_r v_{\phi} - B_{\phi} v_r) + \frac{1}{r} (v_r B_{\phi} - v_{\phi} B_r) \quad (8) $$

and the CR kinetic part:

$$ \begin{aligned} & \frac{\partial F_0(r, p)}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} \{r^2 v_r F_0(r, p)\} - \frac{1}{3} \frac{\partial v_r}{\partial r} \frac{\partial F_0(r, p)}{\partial \ln p} \\ & = -\frac{c}{3} \frac{1}{r^2} \frac{\partial}{\partial r} \{r^2 F_1(r, p)\} + Q_{\text{inj}}(\eta, r, p) p^3, \end{aligned} \quad (9) $$

$$ \begin{aligned} & \frac{\partial F_1(r, p)}{\partial t} + \frac{1}{r^2} \frac{\partial}{\partial r} \{r^2 v_r F_1(r, p)\} \\ & = -\frac{c}{r} \frac{\partial}{\partial r} \{r F_0(r, p)\} + \frac{c F_0}{r} - \frac{c^2}{3 \kappa(p, \mathbf{B})} F_1(r, p), \end{aligned} \quad (10) $$

where $j_r^{(\text{ret})}$ is the return current density induced by the CR streaming current ($j_r^{(\text{ret})} = -j_r^{(\text{cr})}$), $Q_{\text{inj}}$ is the injection rate, and $\kappa(\mathbf{B})$ is the diffusion coefficient. $F_0 \equiv f_0 p^3$ and $F_1 \equiv f_1 p^3$ are, respectively, the isotropic and the first-order anisotropic (drift anisotropy) components of the CR distribution function $f(r, \mathbf{p}) = f_0(r, p) + (p_r/p) f_1(r, p)$. Although they are small corrections, we adequately consider source terms $\propto 1/r$, which are neglected in Inoue et al. (2021).

The detailed form of $j_r^{(\text{cr})}$, $Q_{\text{inj}}$, and $\kappa(\mathbf{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^{(\text{cr})}$ is computed by momentum space integration of $f_1(p)$. The injection rate $Q_{\text{inj}}$ has a nonzero value only at the shock front, and the rate of injection is determined so that the fraction $\eta$ of the thermal particles go into the DSA process (Blasi et al. 2005). We employ the following diffusion coefficient (Skilling 1975; Caprioli & Spitkovsky 2014b):

$$ \kappa(p, \mathbf{B}) = \frac{4}{3\pi} \frac{\max(B_r^2, \delta B^2)}{\delta B^2} \frac{v_{\text{CR}} p c}{e \max(|B_r|, \delta B)}, \quad (11) $$

where $\delta B^2 = B_{\theta}^2 + B_{\phi}^2$ and $v_{\text{CR}}$ is the CR velocity at momentum $p$. Under small amplitude magnetic turbulence ($\delta B \ll 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 $\delta 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_{\text{sh}}) \propto 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 \sim 1$ GeV. In this paper, we study the cases with the injection rate $\eta = 0.1–1.0 \times 10^{-4}$ in which the effect of CR pressure is not substantial. According to a theoretical modeling of multi-wavelength emission of a historical young SNR SN1006 by Berezhko et al. (2012), $\eta \simeq 10^{-4}$ gives the best fit. Note that
<table><caption>**Table 1**
Model Parameters</caption><thead><tr><th>Model ID</th><th>n<sub>ISM</sub></th><th>B<sub>r</sub><br>(μG)</th><th>ξ<sub>B,ini</sub></th><th>η</th><th>Bell term<sup>a</sup></th><th>E<sub>cut</sub><br>(TeV)</th><th>ε<sub>B</sub><sup>b</sup></th></tr></thead><tbody><tr><td>1</td><td>0.5</td><td>3.0</td><td>0.1</td><td>1 × 10<sup>-4</sup></td><td>yes</td><td>31 ± 1.2</td><td>4.0 × 10<sup>-2</sup></td></tr><tr><td>2</td><td>0.1</td><td>3.0</td><td>0.1</td><td>1 × 10<sup>-4</sup></td><td>yes</td><td>41 ± 0.7</td><td>3.6 × 10<sup>-2</sup></td></tr><tr><td>3</td><td>0.02</td><td>3.0</td><td>0.1</td><td>1 × 10<sup>-4</sup></td><td>yes</td><td>29 ± 0.6</td><td>1.3 × 10<sup>-2</sup></td></tr><tr><td>4</td><td>0.5</td><td>3.0</td><td>0.01</td><td>1 × 10<sup>-4</sup></td><td>yes</td><td>29 ± 0.9</td><td>5.0 × 10<sup>-2</sup></td></tr><tr><td>5</td><td>0.5</td><td>10.0</td><td>0.1</td><td>1 × 10<sup>-4</sup></td><td>yes</td><td>37 ± 1.3</td><td>12.6 × 10<sup>-2</sup></td></tr><tr><td>6</td><td>0.5</td><td>3.0</td><td>0.1</td><td>3 × 10<sup>-5</sup></td><td>yes</td><td>19 ± 0.5</td><td>1.5 × 10<sup>-2</sup></td></tr><tr><td>7</td><td>0.5</td><td>3.0</td><td>0.1</td><td>1 × 10<sup>-5</sup></td><td>yes</td><td>9.3 ± 0.4</td><td>0.24 × 10<sup>-2</sup></td></tr><tr><td>8</td><td>0.5</td><td>3.0</td><td>0.1</td><td>1 × 10<sup>-4</sup></td><td>no</td><td>0.82 ± 0.02</td><td>2.4 × 10<sup>-5</sup></td></tr></tbody></table>

**Notes.**

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.

## 2.2. Fluid Initial Conditions

The inner and outer boundaries of the numerical domain are located at $r_{in} = 0.5$ pc and $r_{out} = 3 \times 10^{19}$ cm $\simeq 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 $\rho_0 = 2.1 \times 10^{-23}$ g cm$^{-3}$. This leads the total kinetic energy of the ejecta to be $E_{kin} = 0.76 \times 10^{51}$ erg.

The initial ISM is distributed in $R_0 < r < r_{out}$ with uniform density $\rho_{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 \times 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 \simeq 100$ km s$^{-1}$. Thus, the Mach number of the forward shock becomes $\sim v_0/c_s \simeq 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.<fnref n="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.<fnref n="5" /> The upstream Alfvén velocity is calculated to be $v_A = B_r / \sqrt{4\pi\rho_{ISM}} = 8.2$ km s$^{-1}(B_r/3 \mu\text{G})(n_{ISM}/0.5 \text{ cm}^{-3})^{-1/2}$. For $B_\theta$ and $B_\phi$, we

set turbulent fluctuations by superposing Alfvén waves with flat power spectrum ($B_r^2 = \text{const.}$). The fluctuation dispersion is set as $\xi_B \equiv \langle(B_\theta^2 + B_\phi^2)/B_r^2\rangle = 0.1$ 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.

## 2.3. 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_\theta = B_\phi = 0$ at the inner boundary. This set of boundary conditions reproduces the blast wave evolution starting from the free expansion phase ($v_{sh} \propto t^0$) and then shifts to the Sedov–Taylor phase ($v_{sh} \propto 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 $\kappa = 0$ in the ejecta. At the outer spatial boundary, outgoing free boundary conditions are imposed: $F_0(p, r_{out} + \Delta r) = F_0(p, r_{out})$ and $F_1(p, r_{out} + \Delta r) = \max[F_1(p, r_{out}), 0]$. For the momentum space, as stated in Section 2.1, we assume $f_0 \propto 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 \text{ GeV } c^{-1}, 10 \text{ PeV } c^{-1}]$.

## 2.4. 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): $\lambda_B = c B_r / |j_r^{(CR)}|$. The critical scale below which the NRHI is stabilized by magnetic tension force is given by $\lambda_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 a value of $j_r^{(cr)} \sim 10^{-10}$ esu s$^{-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

$$ \lambda_B \simeq 0.9 \times 10^{15} \text{ cm} \\ \times \left( \frac{j^{(cr)}}{10^{-10} \text{ esu s}^{-1} \text{ cm}^{-2}} \right)^{-1} \left( \frac{B_r}{3 \mu\text{G}} \right), \quad (12) $$

<footnote>a If no, we always set null return current $j^{(ret)} = 0$.</footnote>

<footnote>b Ratio of downstream magnetic energy density to upstream kinetic energy. The downstream magnetic energy is averaged over the 1000 $\Delta x$ region from the shock front.</footnote>

<footnote n="4">⁴ We employ the HLLD Riemann solver (Miyoshi & Kusano 2005) in the integration of the MHD part, which is known to be a very robust scheme.</footnote>

<footnote n="5">⁵ This is because we set inconsistent radial dependence of $B_r$ to the divergence-free condition of the magnetic field, which brings erroneous force term $\propto 1/r$ to the Lorentz force. However, this erroneous force is orders of magnitude smaller than other physical terms whose ratio is estimated to be $\lesssim \lambda_{NRH}/r_{in} \sim 10^{-5}$. Thus, our initial magnetic field condition does not raise any noticeable issues.</footnote>
<figure><img src="image_3.png" /><figcaption>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⁻¹ (triangles) for three different models. Thin red line is a reference line to the Sedov–Taylor solution: ∝t⁻³/⁵.</figcaption></figure>

$$ \omega_B^{-1} = \frac{\lambda_B}{2\pi v_A} = \frac{c \rho^{1/2} |j_r^{(cr)}|}{\pi^{1/2} |j_r^{(cr)}|} \approx 5.5\text{yr} \\ \times \left( \frac{j^{(cr)}}{10^{-10} \text{esu s}^{-1} \text{cm}^{-2}} \right)^{-1} \left( \frac{n}{0.5 \text{cm}^{-3}} \right)^{1/2} . \quad (13) $$

According to the results of test simulations by Inoue et al. (2021), the numerical dispersion relation can reproduce the analytic one if $\lambda_B$ is resolved by more than 32 numerical cells. To satisfy the above condition, we resolve the whole spatial numerical domain by $N_{cell} = 2^{20} \approx 10^6$ cells, leading to the spatial resolution of $\Delta r = (r_{out} - r_{in})/N_{cell} \approx 3 \times 10^{13}$ cm. Thus, the ratio of the instability scale and the resolution is estimated to be $\lambda_B/\Delta r \approx 30$, which satisfies the resolution requirement by Inoue et al. (2021).

For the momentum space, we explicitly treat a range of $[p_L = 100 \text{ GeV c}^{-1}, p_U = 1 \text{ PeV c}^{-1}]$, which is divided into uniform $N_{cell,p} = 64$ numerical cells in the logarithmic scale, i.e., $\Delta \ln p = \ln(p_U/p_L)/N_{cell,p}$. Because the CR current is mostly composed of CRs with $E \gtrsim 10 \text{ TeV}$ (see Figure 3 below), as long as the lower boundary of the momentum space is set below $1 \text{ 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).

## 2.5. 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 $\Delta t \le 0.5 \Delta r^2/\kappa(B_r, p_U) \sim 10^{-3}$ s for the initial state, indicating a 1000 yr scale integration is impossible even with the largest-scale supercomputer.

On the other hand, our hyperbolic-type basic equations require $\Delta t \le C_{CFL} \Delta r/(c/\sqrt{3}) \sim 10^3$ s, which is still short, but the 1000 yr integration is executable with a modern supercomputer. Here, $c/\sqrt{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. Super-time-stepping can only enlarge $\Delta 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.

# 3. Results

## 3.1. 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 \text{ km s}^{-1}$. We can confirm that the shock evolves from the initial free expansion phase ($v_{sh} \propto t^0$) to the Sedov–Taylor phase ($v_{sh} \propto 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, \text{ and } 1000 \text{ 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 $\epsilon_B$, i.e., the ratio of downstream magnetic energy density to upstream kinetic energy, where the downstream magnetic energy is averaged over 1000 $\Delta x$ region from the shock front at $t = 1000 \text{ 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_p = d j_r^{(CR)}/d \ln p$ at $t = 1000 \text{ yr}$ for Model 1 in the bottom panel. We see that CRs with energy $pc \sim 10 \text{ TeV}$ mostly
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

compose the upstream CR current. The spatial extent of the CR current is an order of magnitude larger than the diffusion length of the 10 TeV CRs:

$$l_{\text{diff}} = \kappa(p)/v_{\text{sh}} \\ \approx 10^{17} \text{cm} \ \xi_B^{-1} \left( \frac{p_{\text{max}}}{10 \text{ TeV} \ c^{-1}} \right) \left( \frac{B}{10 \ \mu\text{G}} \right)^{-1} \\ \times \left( \frac{v_{\text{sh}}}{3 \times 10^3 \ \text{km} \ \text{s}^{-1}} \right)^{-1} \quad (14)$$

At time $t_0$, advection time $t_{\text{ad}}$ for the gas at distance $l = r - r_{\text{sh}}$ from the shock front is obtained by solving

$$l = \int_{t_0}^{t_{\text{ad}}} v_{\text{sh}} dt, \quad (15)$$

For $t_0 \gtrsim 500$ yr, from the Sedov–Taylor solution, the shock velocity is given by $v_{\text{sh}} = v_{\text{sh}}(t_0)(t/t_0)^{-3/5}$, which leads to

$$t_{\text{ad}} = \left( \frac{2}{5} \frac{l}{v_{\text{sh}}(t_0)t_0} + 1 \right)^{5/2} t_0. \quad (16)$$

Using $t_{\text{ad}}$, the expected e-folding number of the NRHI at distance $l$ is computed as

$$\sigma(t_0, r) = (t_{\text{ad}} - t_0)\omega_B. \quad (17)$$

In Figure 4, we plot $\sigma$ 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^{\text{afew}} \sim 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):

$$B_{\text{sat}} = \left( \frac{4\pi j_r^{(\text{CR})} p_{\text{max}}}{e} \right)^{1/2} \\ \approx 64 \ \mu\text{Gauss} \left( \frac{j_r^{(\text{CR})}}{10^{-10} \text{ esu s}^{-1} \text{ cm}^{-2}} \right)^{1/2} \left( \frac{p_{\text{max}}}{30 \text{ TeV c}^{-1}} \right)^{1/2}, \quad (18)$$

where $p_{\text{max}}$ is the maximum (escaping) CR momentum. The magnetic field strengths in the simulations do not reach the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Panel (a): structure of the CR current densities for Models 1, 2, and 3 at $t = 500$ (thin lines) and 1000 yr (thick lines). Panel (b): upstream CR current spectrum $j_p(r, p) = d j_r^{(CR)}/d \ln p$ at $t = 1000$ yr for Model 1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Expected e-folding number of the NRHI at distance $r - r_{sh}$ (see Equation (17)) for Model 1 (magenta), Model 2 (green), and Model 3 (light blue) at $t_0 = 500$ yr (circle) and 1000 yr (triangle).</figcaption></figure>

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><img src="image_7.png" /><figcaption>**Figure 5.** CR density spectrum ($n_p = f_0 p^2$) of Model 1 immediately behind the shock. Black line represents the DSA spectrum ($\propto p^{-2}$). We can reasonably fit the spectrum by a power-law function with the standard DSA spectral index and exponential cutoff: $n_p \propto p^{-2} \exp(-pc/E_{cut})$. The result of fitting for $t = 1000$ yr is shown as a dotted line ($E_{cut} = 31$ TeV).</figcaption></figure>

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: $n_p \propto p^{-2} \exp(-pc/E_{cut})$. 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} \sim 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):

$$E_{\max} \approx 43\text{TeV} \xi_B \left( \frac{B}{10 \ \mu\text{G}} \right) \left( \frac{v_{\text{sh}}}{3000 \ \text{km s}^{-1}} \right)^2 \left( \frac{t}{1000 \ \text{yr}} \right). \quad (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 $\xi_B = 0.1$.

## 3.2. Influence of Seed Amplitude and Initial Field Strength

In Model 4, we run the case that has a smaller initial magnetic field fluctuation amplitude ($\xi_B = 0.01$) than the fiducial model. We also study the case of stronger initial field strength ($B_r = 10 \ \mu\text{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
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Panel (a): magnetic field structure around the shock for Model 1 ($B_r = 3$ µG, $\xi_B = 0.1$), Model 4 ($B_r = 3$ µG, $\xi_B = 0.01$), and Model 5 ($B_r = 10$ µG, $\xi_B = 0.1$) at $t = 1000$ yr. Panel (b): expected e-folding number of the NRHI at distance $r - r_{sh}$ (see Equation (17)) for Model 1 (magenta), Model 4 (green), and Model 5 (light blue) at $t_0 = 500$ yr (circle) and 1000 yr (triangle).</figcaption></figure>

larger initial strength model shows a more compact NRHI active region.

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

### 3.3. Effect of the Injection Rate

Since the CR current has a substantial impact on the growth rate of NRHI, the injection rate $\eta$ 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 $\eta$ increases, the influence of $\eta$ variation is not so striking even though there is an order of

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Panel (a): magnetic field structure around the shock for Model 1 ($\eta = 1 \times 10^{-4}$), Model 6 ($\eta = 3 \times 10^{-5}$), and Model 7 ($\eta = 1 \times 10^{-5}$) at $t = 1000$ yr. Panel (b): structure of the CR current densities for Models 1, 6, and 7 at $t = 1000$ yr.</figcaption></figure>

magnitude difference of $\eta$ 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 $\eta$ can be simply explained as follows: In the very early stage, the larger $\eta$ model produces a stronger CR current that induces active NRHI and boosts up $E_{max}$ faster than the smaller $\eta$ 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 \propto p^{-2} dp$, 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 higher-energy CRs than those of Model 7.

## 4. Summary and Discussion

We 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:
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Panel (a): upstream CR current spectrum $j_p(r, p) = d_j^{(\text{CR})}/d \ln p$ at $t = 1000$ yr for Model 6. Panel (b): same as panel (a) but for Model 7.</figcaption></figure>

1. The NRHI is indeed effective at the young SNR forward shocks with age around 1000 yr irrespective of the ISM density. Downstream $\epsilon_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_{\text{age}} \lesssim 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_{\text{age}} \sim 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 ($\eta \sim 10^{-4}$ at maximum), there may be an SNR population that has a higher injection rate by a yet-unknown 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.

## Acknowledgments

The numerical computations were carried out using the supercomputer “Flow” at Information Technology Center, Nagoya University and XC50 system at the Center for Computational Astrophysics (CfCA) of National Astronomical Observatory of Japan. This work is supported by grants-in-aid from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan (20H01944 and 23H01211).

## ORCID iDs

Tsuyoshi Inoue https://orcid.org/0000-0002-7935-8771

Alexandre Marcowith https://orcid.org/0000-0002-3971-0910

Gwenael Giacinti https://orcid.org/0000-0001-9745-5738

## References

Aharonian, F., An, Q., Axikegu, B. L. X., et al. 2021, Natur, 594, 33

Aharonian, F. A., Buckley, J., Kifune, T., & Sinnis, G. 2008, RPPh, 71, 096901

Ballet, J. 2006, AdSpR, 37, 1902

Bell, A. R. 2004, MNRAS, 353, 550

Bell, A. R., Schure, B. R., & Giacinti, G. 2013, MNRAS, 431, 415

Berezhko, E. G., Ksenofontov, L. T., & Völk, H. J. 2012, ApJ, 759, 12

Blandford, R. D., & Eichler, D. 1987, PhR, 154, 1

Blasi, P., Gabici, S., & Vannoni, G. 2005, MNRAS, 361, 907

Caprioli, D., & Spitkovsky, A. 2014a, ApJ, 794, 46

Caprioli, D., & Spitkovsky, A. 2014b, ApJ, 794, 47

Drury, L. O'. C. 1983, RPPh, 46, 973

Inoue, T. 2019, ApJ, 872, 46

Inoue, T., Marcowith, A., Giacinti, G., van Marle, A. J., & Nishino, S. 2021, ApJ, 922, 7

Koyama, K., Petre, R., Gotthelf, E. V., et al. 1995, Natur, 378, 255

Marcowith, A., Bret, A., Bykov, A., et al. 2016, RPPh, 79, 046901

Marcowith, A., Dwarkadas, V. V., Renaud, M., et al. 2018, MNRAS, 479, 4470
Miyoshi, T., & Kusano, K. 2005, JCoPh, 208, 315

Parizot, E., Marcowith, A., Ballet, J., & Gallant, Y. A. 2006, A&A, 453, 387

Reville, B., & Bell, A. R. 2012, MNRAS, 419, 2433

Roh, S., Inutsuka, S., & Inoue, T. 2016, APh, 73, 1

Schure, K. M., & Bell, A. R. 2013, MNRAS, 435, 1174

Skilling, F. 1975, MNRAS, 172, 557

Suzuki, H., Bamba, A., Yamazaki, R., & Ohira, Y. 2022, ApJ, 924, 45

Tibet ASγ Collaboration, Amemori, M., Bao, Y. W., et al. 2021, NatAs, 5, 460

Vink, J., & Laming, J. M. 2003, ApJ, 584, 758

Zirakashvili, V. N., & Ptuskin, V. P. 2008, ApJ, 678, 939

Zirakashvili, V. N., Ptuskin, V. P., & Volk, H. J. 2008, ApJ, 678, 255