© The Authors 2025

1. Introduction

For nearly a century, supernova remnants (SNRs) have been extensively studied and discussed as prime sources of Galactic cosmic rays (CRs), with energies up to the knee feature in the CR spectrum achieving a few PeV (see, e.g., Baade & Zwicky 1934; Blandford & Eichler 1987). However, despite three decades of remarkable development in gamma-ray astronomy, there is still no undisputed observational evidence that SNRs can accelerate particles above ∼100 TeV. Young SNRs such as Tycho and Casiopeia A, which had been expected to be effective particle accelerators, display even lower cutoff energies (Abeysekara et al. 2020).

Particle acceleration and high-energy radiation at shocks depend on the shock speed, the mass flux through the shock (i.e., pre-shock density), and the ambient magnetic field. Thus, we could expect that a structured circumstellar medium (CSM) would significantly modify the predicted CR population and nonthermal light curve of a supernova (SN). Murase et al. (2011) and Murase et al. (2014) were the first to investigate this, finding that SNe with a strong circumstellar interaction could accelerate particles beyond a PeV and be detectable in gamma rays out to 200 Mpc, with bright accompanying nonthermal radio emission.

It is possible that only very young SNRs evolving in dense environments may satisfy the necessary conditions to accelerate particles to peta-electronvolt (PeV) energies. Driving magnetic turbulence to the relevant scales via the excitation of the nonresonant streaming instability (Bell 2004) requires high CR currents that are possible only during the initial ∼20 years of the SNR evolution (Bell et al. 2013), while high initial shock velocity and high ambient density further boost particle acceleration (Marcowith et al. 2018; Cristofari et al. 2020b). Recent active research in this field resulted in a common conclusion that only progenitor stars with an extremely dense CSM, produced by periods of high mass-loss rates with slow wind speed, have the required shock conditions to reach PeV energies (Marcowith et al. 2018; Cristofari et al. 2020a, b; Inoue et al. 2021; Diesing 2023). If true, this means that only a relatively rare subset of SN explosions can contribute anything to the flux of Galactic CRs around the knee at 3 PeV.

In our recent paper, (Brose et al. 2022, hereafter Paper I), we investigated particle acceleration and nonthermal emission from SNRs associated with explosions of the LBV and red supergiant (RSG) progenitors during the first 20 years of their evolution, for explosions into smooth 1/r² density profiles appropriate for stellar winds without variability. We assessed the problem by numerical simulations with the Radiation Acceleration Transport Parallel Code (RATPAC) software extensively described in previous works (Telezhinsky et al. 2012, 2013; Brose et al. 2016, 2019; Sushch et al. 2018).

RATPAC uses the kinetic approach to solving the transport equation for particles in time and space simultaneously with the solution of the hydrodynamic equations of the thermal plasma. It also solves the transport equation for the energy density in magnetic turbulence that is self-generated by accelerated particles. This approach allows us to self-consistently calculate the spectrum of accelerated particles and estimate the maximum energy achievable in the acceleration process, naturally accounting for the time required to build up the turbulence. We showed that for the most favorable values of the stellar wind parameters, the maximum energy saturates at about 600 TeV for the LBV case. More common densities of the CSM facilitate acceleration only up to 200 TeV and 70 TeV for a LBV and RSG progenitor, respectively. These results agree with other works that predict acceleration up to at most a few PeV for the most optimistic scenarios (see, e.g., Marcowith et al. 2018). Differences in the estimated maximum energy could be attributed to assumptions implemented for the microphysics of the turbulence generation that vary from one study to another. The feasibility of these assumptions are discussed further in this paper.

So far all the studies, including Paper I, have considered the expansion of the SNR into a smooth featureless wind while reality could strongly differ from this assumption. While many aspects of mass loss are poorly understood during late evolutionary stages of massive stars, stars may undergo rapid changes in radius with associated change in mass-loss rate and wind speed (e.g., Langer 2012). Some of the most extreme examples are the LBVs that undergo cyclic variations in radius of a factor of up to 10 on a timescale of a few years (Humphreys & Davidson 1994; Grassitelli et al. 2021), resulting in repeated episodes of mass loss through fast and then slow winds. Some LBVs also undergo giant eruptions where a few solar masses of material are ejected on a timescale of a year, expanding with ∼100 km s⁻¹; the most famous example of this being η Carinae (Smith 2014). At least some of these giant eruptions occur shortly before core-collapse, as in the case of SN 2009ip (Mauerhan et al. 2013).

Binary interaction has emerged as a key aspect of the evolution of massive stars (e.g., de Mink et al. 2014) with strong consequences for the CSM due to rapid and asymmetric mass loss through a common-envelope or Roche-lobe-overflow phase. Common envelope followed by a stellar merger is a possible explanation for the dramatic three-ring structure revealed in the CSM of SN1987A (Morris & Podsiadlowski 2006; Menon et al. 2019). Stellar mergers appear to generate strong magnetic fields in massive stars (Schneider et al. 2019), which will, in turn, produce a strongly magnetised CSM. Through these processes, stellar evolution from single and binary stars may produce a highly structured CSM around a star at the time of explosion, with sufficiently high density and inertia to affect the hydrodynamic expansion of the SN shock in the first months and years of expansion.

In this paper, using the same approach as in Paper I, we explore such SN-CSM interactions with spherically symmetric calculations and their impact on the maximum energy achievable in the acceleration process. The resulting emission both thermal and nonthermal will be examined in a follow-up publication. Section 2 describes the model setup and introduces our assumptions for the density and magnetic field of a pre-SN CSM. We also further discuss implications of certain aspects of our numeric simulations and assumptions that they rely upon. In Section 3 we describe the time evolution of the accelerated particle spectra. Section 4 presents our conclusions.

2. Model setup

The numerical setup of simulations as well as the basic assumptions and equations that are solved are the same as in Paper I. We use the RATPAC code to combine a kinetic treatment of the CRs with a thermal leakage injection model, a fully time-dependent treatment of the magnetic turbulence, and a PLUTO-based hydrodynamic calculation. The essential ansatz to describe the evolution of the CR-distribution is solving the kinetic diffusion-advection equation for CR-transport (Skilling 1975) as follows:

$$\begin{array}{cc}\hfill \frac{\partial N}{\partial t}=& \nabla (D_r\nabla N-\text{u}N)\hfill \\ \hfill & -\frac{\partial }{\partial p}\left((N\dot{p})-\frac{\nabla ·\text{u}}{3}\mathit{Np}\right)+Q,\hspace{0.5em}\hfill \end{array}$$(1)

where D_r denotes the spatial diffusion coefficient, u the advective velocity, ṗ the energy losses, and Q is the source of thermal particles.

For detailed description of the setup we refer the reader to Paper I, while here we focus on the only difference in the setup compared to our previous work which is the structure of the CSM into which the SNR is expanding. Additionally we further discuss and analyse the assumptions we adopt for the microphysics of the self-generated turbulence and their impact on the results.

2.1. Circumstellar medium

In this work, we focus on the most extreme cases of possible CSM shells. We have identified LBV eruptions, where the mass-loss rate can reach values as high as 1 M_⊙ yr⁻¹ for periods of a few years and where we would expect the densest CSM (Smith 2014). However, Type IIn SNe that are observationally associated with LBVs make up only ∼5% of the CCSN rate (Cold & Hjorth 2023).

The simplest model for an LBV eruption is to impose a dramatically increased mass-loss rate for a short period of time. We modelled an evolving stellar wind with the astrophysical fluid-dynamics code PION (Mackey et al. 2021) in spherical symmetry (Mackey et al. 2014) by simulating the expansion of a stellar wind in three phases. The stellar wind is injected at its terminal velocity, v_∞, with the density set by the mass-loss rate, Ṁ, and the continuity equation. Wind pressure is set by assuming adiabatic expansion from the stellar surface, but its exact value is unimportant for large Mach number winds. The wind injection region is typically the first two grid cells closest to the origin. The first phase has mass-loss rate Ṁ = 10⁻⁴ M_⊙ yr⁻¹ and wind terminal velocity v_∞ = 100 km s⁻¹ and lasts for 1800 years; the second (outburst) phase has Ṁ = 1 M_⊙ yr⁻¹ and v_∞ = 100 km s⁻¹ and lasts for two years; the third is the same as the first, but we integrated it for 10 000 years, so that the shell would expand to a radius of about 1 pc. The density structure of the shell is then used to construct the circumstellar shell in the RATPAC simulations as described below.

In Paper I, we modeled the CSM density profile following a ∝1/r² distribution applicable for the free wind. Here, we additionally introduce shells of enhanced mass following a Gaussian distribution,

$$\rho (r)=\frac{M_{\mathrm{shell}}}{4\pi r^2\sqrt{2\pi d_{\mathrm{shell}}^2}}exp\left(-\frac{(r-R_{\mathrm{shell}}{)}^2}{2d_{\mathrm{shell}}^2}\right),$$(2)

where M_shell is the shell mass, d_shell the shell thickness, and R_shell the radial position of the shell. The parameters for the free winds and the corresponding shells are given in Table 1.

We obtained the shell parameters by fitting the results of the stellar-wind simulations described before with the Gaussian profiles described by Equation (2). An example for the CSM structure as used in the RATPAC simulations is shown in Figure 1.

Fig. 1. Upper panel: Density profile of the CSM for the two of our LBV progenitor models. The shell parameters were obained from PION simulations. The parameters of the shells used in RATPAC are listed in Table 1. Lower panel: Shock velocity of the models as a function of the shock-position. Upward steps in the shock velocity are from interactions of the forward shock with reflected shocks catching-up from behind.

The total mass of material in the CSM is important not only for the shock dynamics but also for absorption in the ambient medium, for instance of X-rays and radio waves. It has to be noted that we ignore the effects of radiative cooling in our calculations. Strictly speaking, especially around the reverse-shock, the shock will be radiative initially but the intense photon-fields would require the use of a different equation of state for the plasma. The interactions with the dense shells will push the shock-dynamics into the radiative regime for a brief period even for the forward-shock for our LBV cases with the earliest interactions. However, resolving the changed dynamics of radiative shocks would be very challenging and the particle injection algorithms would need significant modification if the post-shock cooling length could not be spatially resolved.

We modeled the large-scale magnetic field in the wind in an identical way to the process detailed in Paper I. Hence, the field is described as

$$B(r)=B_{*}\frac{R_{*}}{r},$$(3)

where B_* is the surface magnetic field and R_* the radius of the progenitor. We assumed a constant ratio of B_*(R_*/R_⊙ = 1000 G for both progenitor cases. In this work, we do not account for any potential compression of the magnetic field in the shells, as studying this aspect would require a more profound understanding of the variation of the surface field of LBVs and RSGs throughout their lifetime. Additional structure in the magnetic field (especially compression in the RSG case) potentially introduces effects such as those described in Sushch et al. (2022), whereas this work aims to understand the impact of the density structure of the CSM.

2.2. Self-generation of magnetic turbulence

Any treatment of diffusive shock acceleration (DSA) needs a prescription for the scattering of CRs, as they need to be returned to the shock from the downstream region in order to undergo repeated shock-crossings. The various models on the market follow different approaches in terms of describing the scattering or diffusion of CRs and we briefly review the current status of the field, to place our ansatz in the appropriate context below.

The only numerical method that is able to capture all the microphysics involved is particle-in-cell simulation. However, to this point the computational costs are too large to simulate timescales beyond a few ion-cyclotron times, ω_i, and usually need to use reduced ion-to-electron mass ratios (e.g., Amano & Hoshino 2009; Bohdan et al. 2019a, b; Amano et al. 2022). Hybrid simulations, which treat electrons as a fluid are able to cover longer time periods; however, the simulated timescales are still on the order of a few hours and it has at least be questioned that the observed acceleration is yet fully DSA-governed (Caprioli & Spitkovsky 2014a, b, c; Caprioli et al. 2020). All other approaches to tackling the shock acceleration problem rely on simplifications of the underlying microphysics in one or multiple aspects.

2.2.1. Scattering problem

One central aspect is what we call the “scattering problem” in the following. CRs have to be returned to the shock in order to participate in the acceleration process. Fermi's initial description of DSA just relies on the assumption that magnetic structures (or mirrors) exist, which scatter the particles isotropically in the rest-frame of the local plasma (Fermi 1949). By now, it has been established that CRs are scattered by (magnetic) turbulence in the plasma, although the question of which wave-modes exactly are responsible for the scattering is still subject to ongoing research. Among the types of turbulence that are discussed are Alfven-waves (e.g., Skilling 1975), Bell's mode waves (Lucek & Bell 2000) or fast-mode waves (Yan & Lazarian 2004). Each of these has a different driving mechanism and different resonance kernels when it comes to CR scattering.

In the early calculations of DSA, it was believed that CRs stream along the large-scale magnetic field and were scattered back towards the shock by resonant interactions with Alfven-waves (e.g., Skilling 1975). With the realization that the magnetic field in SNRs needs to be amplified beyond the value of the large-scale field, the resonant growth of Alfven waves (Skilling 1975; Bell 1978) was deemed not to be efficient enough (Lagage & Cesarsky 1983). An apparent solution is the nonresonant streaming instability (Lucek & Bell 2000; Bell 2004). However, while in this case the turbulent field can grow to values of δB≈10B₀, the instability initially grows on very small scales and the nature of the wave-particle interaction changes. The central assumption to many studies of CR acceleration is that particles diffuse at the Bohm limit as inferred from quasilinear theory for the ∝E⁻² particle spectrum (e.g., Vladimirov et al. 2006; Caprioli et al. 2009b; Kobashi et al. 2022; Diesing 2023, and associated works). However, the relation between a certain energy-spectrum (and type) of magnetic turbulence and CR scattering is neither straightforward nor is it well understood at the moment. Despite the existence of some evidence of Bohm-like diffusion in the presence of nonresonant magnetic field amplification as inferred from numeric simulations, this evidence is far from being conclusive.

High-resolution MHD simulations of the nonresonant instability performed by Reville et al. (2008) indicate that diffusion in amplified field approaches the Bohm limit¹ for particles at lower energies than those driving the turbulence growth. Further studies utilizing three-dimensional (3D) hybrid MHD–kinetic simulations implied a diffusion coefficient that is about four times larger than Bohm for particles driving the field growth; however, the simulations end before the saturation is reached (Reville & Bell 2013). Kinetic-protons-fluid-electrons hybrid simulations imply that after a sufficient amount of time diffusion reaches the Bohm limit in the total amplified field at all energies except for the highest (Caprioli & Spitkovsky 2014c), but it should be understood that these simulations operate on much smaller scales and it is unclear how we should apply these results to length-scales of order the remnant size.

Finally, it should be noted that X-ray observations of young SNRs suggest that diffusion of electrons falls short of the Bohm limit in most cases (Tsuji et al. 2021; Sapienza et al. 2024). While approximations are needed to tackle the DSA-problem at the scales of real astrophysical sources, it must be clear that the assumption of Bohm-like scaling of diffusion and a Bohm-like turbulence spectrum constitutes the most optimistic case possible.

2.2.2. Steady-state assumptions

A second, a common simplification is the assumption that one or more of the involved processes are operating fast enough to reach a steady state. Commonly, this is assumed to be the case for the turbulent magnetic field by assuming that the saturation level of the nonresonant instability can be reached (e.g., Bykov et al. 2014; Marcowith et al. 2018; Diesing 2023).

A few fully time-dependent simulations in the literature (Brose et al. 2016, 2022; Inoue et al. 2021) demonstrated that the shock-capture time is sometimes not sufficient to allow enough growth-times for the nonresonant instability to saturate, which is in accordance with earlier analytical estimates by Niemiec et al. (2008). This can easily be understood as a “chicken-and-egg-problem”. Simulations of the nonresonant instability usually assume a beam of monoenergetic high-energy particles to trigger the growth of the instability. For example, recent hybrid simulations considering a current of ions with a monochromatic initial momentum show that the saturation of the Bell instability takes place after ∼10 growth times when the length scale of perturbations becomes comparable to the ion gyroradius (Zacharegkas et al. 2024). In reality, the particles of energy E have to provide turbulence scattering particles with E+dE. However, the scattering for particles of E+dE is less efficient unless enough turbulence was build up at the right scales. This interdependence of turbulence and maximum particle energy slows down the acceleration process and can limit E_max (Brose et al. 2016).

The fact that the nonresonant instability might not actually reach saturation has implications again for the so-called scattering problem. The instability saturates when the scale of the driven turbulence becomes comparable to the gyroradius of the most-energetic particles in the system. It has not been studied so far, whether the nonsaturated and nonresonant instability is as efficient in scattering particles as the saturated case studied by Reville et al. (2008). Steady-state approaches assume the Bohm-diffusion regime is applicable over the entire growth time of the instability.

2.2.3. Turbulence-driving CR flux

Another difference in the literature (especially concerning the maximum energy) is which part of the CRs is driving the nonlinear instability. The two possibilities at the extreme ends are that all CRs contribute, which requires that the field-amplification caused by escaping CRs rebuilds a quasi-planar shock geometry far upstream of the shock. This assumption is used in the approaches described in the literature (Bell & Lucek 2001; Bykov et al. 2014; Cristofari et al. 2020b), while the field-growth of Paper I and Inoue et al. (2021) is driven by particles in the precursor only; hence, it is lower than in the other models. Diesing (2023) calculated the maximum energy for both assumptions and finds a difference that is roughly one order of magnitude in E_max.

2.2.4. Energy fraction of injected CRs

Besides the particularities of turbulence growth, the actual energy carried by CRs plays the single most important role for determining the achievable maximum energy. Many analytic or semi-analytic approaches for the DSA problem rely on assuming the actual fraction of the explosion energy that converted to CRs. The usual assumption is that 10% of the explosion energy get converted to CRs. However, some works assume that this fraction of energy is already converted very early on (e.g., Murase et al. 2011). However, conceptually assuming a certain energy fraction a priori is unphysical, but it is needed to make the problem mathematically tractable (Blasi 2002; Caprioli et al. 2009a; Murase et al. 2011; Kobashi et al. 2022). Recent Fermi-LAT observations of SN 2023ixf indicate, that the energy fraction of CRs right after explosion ≤1% (Martí-Devesa et al. 2024). Hybrid simulations on the other hand suggest injection efficiencies on the order of 10–15% of the energy flux passing through the shock (Caprioli et al. 2020). Our simulations result in a comparable but slightly lower injection fraction, reaching the order of 5–10% of the explosion energy in CRs after 50 kyrs (Brose et al. 2020).

The actual fraction of thermal particles that participate in the DSA process is not clear and subject to ongoing studies. Our injection fraction is tailored to observations of historical SNRs, assuming that the injection fraction is roughly constant over time and shock conditions (Brose et al. 2020).

The injection efficiency that we have applied here automatically ensures that the CR pressure at the shock stays below 10% of the shock's ram pressure. Thus, the flow structure ahead of the shock is not modified by the CRs themselves and our simulation operated within the so-called test-particle limit. While we need to make this simplification due to computational constraints, we note that modifications to DSA are discussed in the community at the moment Caprioli et al. (2020, and references therein). They usually result in a softening of the spectral index at low energies and, hence, in a lower CR flux at the highest energies compared to our model. In this respect, any modification to our ansatz softening the spectrum will likely also reduce the maximum energy that is achieved.

2.2.5. Our ansatz

In this work, we assess the magnetic turbulence by solving the transport equation for the turbulence spectrum in parallel to the transport equation of CRs (Brose et al. 2016). While this method was initially developed to describe Alfvenic turbulence, we include parts that are related to the nonresonant modes without adding an additional turbulence component. In any case, the diffusion coefficient varies strongly in space and time and is coupled to the spectral energy's density per unit logarithmic bandwidth, E_w(r,k,t). The evolution of E_w is described by

$$\begin{array}{cc}\hfill \frac{\partial E_w}{\partial t}+\nabla ·(\text{u}{E}_w)+k\frac{\partial }{\partial k}\left(k^2{D}_k\frac{\partial }{\partial k}\frac{E_w}{k^3}\right)=& \hfill \\ \hfill =2({\mathrm{\Gamma }}_g-{\mathrm{\Gamma }}_d)E_w+Q.\hspace{0.5em}& \hfill \end{array}$$(4)

Here, u denotes the advection velocity, k the wavenumber, D_k the diffusion coefficient in wavenumber space, and Γ_g and Γ_d the growth and damping terms, respectively (Brose et al. 2016).

The diffusion coefficient of CRs is coupled to E_w by

$$D_r=\frac{4v}{3\pi }{r}_g\frac{U_m}{E_w},\_$$(5)

where U_m denotes the energy density of the large-scale magnetic field, v is the particle velocity, and r_g the gyro-radius of the particle.

As initial condition, we used a turbulence spectrum derived from the diffusion coefficient, as suggested by Galactic propagation modeling (Trotta et al. 2011), reduced by a factor of 10 on account of numerical constraints:

$$D_0={10}^{28}{\left(\frac{\mathit{pc}}{10\text{GeV}}\right)}^{1/3}{\left(\frac{B_0}{3\mu \text{G}}\right)}^{-1/3}\frac{{\text{cm}}^2}{\text{s}}.$$(6)

We used a growth rate based on the resonant streaming instability (Skilling 1975; Bell 1978) as follows:

$${\mathrm{\Gamma }}_g=A·\frac{v_{\mathrm{A}}{p}^{2}v}{3E_{\mathrm{w}}}\left|\frac{\partial N}{\partial r}\right|,\hspace{0.5em}$$(7)

where v_A is the Alfvén velocity. We introduced a linear scaling factor, A, to artificially enhance the amplification. We used A = 10 throughout this paper to mimic the more efficient amplification due to the nonresonant streaming instability (Lucek & Bell 2000; Bell 2004). Using A as a scaling parameter is of course only a crude approximation of the process. However, as outlined before, there is at the moment no fully self-consistent prescription of the nonlinear instability over the relevant timescales as described below.

We found in earlier works that a value of 10 seems to be consistent with observations of historical SNRs (Brose et al. 2021) and estimates of the growth rates operating at the early stages of CR acceleration (Marcowith et al. 2018). Since we have exceeded the growth rate of the resonant streaming instability (Bell 1978) by a factor of 10, the turbulent field was amplified to δB≫B₀.

Given this state of CR-acceleration theory in SNRs, we have to emphasize that our approach is not free of simplifying assumptions. We can only mimic the growth rate of the nonresonant instability in our fully time-dependent calculation, where we resolve the energy spectrum of the magnetic turbulence at any location (noting that we are not numerically able to track the return current or, thus, the onset of the saturation. Additionally, we assume resonant scattering of CRs on the turbulence that we are aiming to drive

$$D(p)\propto r_g(p)\frac{B_0^2}{\delta B(p{)}^2}.$$(8)

Our approach gives a more realistic description of D(p), especially at larger distances from the shock compared to the usual assumption of Bohm-like diffusion everywhere. In general, our estimates for E_max are within a factor of 2, compared to Inoue et al. (2021) when similar ambient parameters are concerned. Furthermore, we include cascading, which acts as a damping mechanism.

The spectral energy transfer by cascading balances the growth of the magnetic turbulence and hence the magnetic field level in our simulations. This process is described as a diffusion process in wavenumber space and the diffusion coefficient is given by (Zhou & Matthaeus 1990; Schlickeiser 2002) as:

$$D_{\mathrm{k}}=k^3{v}_{\mathrm{A}}\sqrt{\frac{E_{\mathrm{w}}}{2B_0^2}}.$$(9)

This phenomenological treatment will result in a Kolmogorov-like spectrum, if cascading is dominant. The cascading rate will depend on the level of turbulence in two regimes, since v_A∝B_tot

$$D_k\propto \{\begin{array}{cc}\sqrt{E_{\mathrm{w}}}\hfill & \text{for}\hspace{0.5em}{E}_{\mathrm{w}}\ll B_0^2/8\pi ,\hfill \\ E_{\mathrm{w}}\hfill & \text{for}\hspace{0.5em}{E}_{\mathrm{w}}\gg B_0^2/8\pi .\hfill \end{array}$$(10)

When the turbulent field starts to dominate over the background field, the cascading rate depends more sensitively on the energy density of magnetic turbulence.

Again, this mechanism is efficiently shaping the diffusion coefficient at larger distances from the shock by redistributing scattering turbulence from scales resonant with high-energy to low-energy particles. We did not include any other damping mechanism in this work as we assume that the dense shells get photoionized by the SN flash. However, in comparison to other models, we are able to already single out a less optimistic scenario and additional damping can only bring the maximum energy further down.

3. Results and discussion

Using the wind as well as the circumstellar shell and SN parameters in Table 1, we simulated the evolution of the remnants for 20 years for each of the six cases. In this paper, we focus on the evolution of the particle distribution and maximum energy (Sections 3.1 and 3.2) and the amplified magnetic field (Section 3.3). A follow-up publication will address the emission signatures from Radio, optical, and (thermal) X-rays to gamma rays.

3.1. Maximum particle energies

To obtain the maximum energy we fitted the simulated proton spectrum below the cutoff to a power law with a fixed spectral index of s=−2 to obtain the spectral normalization N₀². We then define E_max as the lowest energy where

$$N(E_{\max })E_{\max }^2/N_0=1/e.$$(11)

The time evolution of the maximum energy for our LBV simulations is shown in Figure 2. The shock-shell interaction decreases E_max initially. For the readability of the plots, we employed E_max(t) = max([E_max(0),E_max(t)]). During the interaction, substructure arises in the spectra due to injection of many particles at low energies. The cutoff shape changes over time and is initially super exponential, with its shape following $exp\left(-{\left(\frac{E}{E_{\max }}\right)}^a\right)$ with a>20. The definition of E_max according to Equation (11) avoids the influence of a on E_max obtained by fitting. However, we note that this treatment of the spectrum is an assumption we make to extract the maximum energy in a consistent way across all simulations. The actual particle (and, subsequently, the emission) spectra can feature substructures and spectral breaks and strongly deviate from s=−2 spectra at energies below the cutoff. Observations indeed suggest very soft radio-spectra for Type-Ib and Type-Ic SNe, indicating electron-spectra as soft as s=−3 (Chevalier & Fransson 2006). However, Type-IIP explosions which are more closely related to the cases that we focus on based on the properties of the progenitors, have been modeled with spectra closer to s=−2.2 and it was pointed out that cooling can significantly soften the spectra in these cases (Chevalier et al. 2006). We discuss the emission spectra in detail in a follow-up publication.

Fig. 2. Top: Maximum energy of protons for the simulated LBV progenitors plotted over the simulation time. Bottom: Maximum energy obtained in our simulation as a function of the onset time of shell interaction for the LBV case (points) and the fit using the function defined in Equation (12).

The interaction of the SN with the dense shell introduces additional structure in evolution of the maximum energy for the LBV progenitors, if the interaction takes place in the first two years after the explosion. In all these cases, the maximum energy is enhanced after the shock passes through the shell, although the jump in the maximum energy depends on the time of the interaction. In general, the earlier the interaction, the larger E_max becomes. For an early interaction at ≈0.1 yrs, E_max surpasses 10 PeV. Fitting the maximum energy dependent on the interaction time (see Figure 2) by

$$E_{\max }=E_0{\left(\frac{t}{36\text{days}}\right)}^{-\beta }+E_1,$$(12)

yields β=−1.80±0.03, E₀ = 10.8±1.1 PeV and E₁ = 90 TeV. Hence, we predict that early interactions before ∼140 days lead to a maximum energy above 1 PeV.

This is not affected by the fact that the case of t_shell = 0.1yr shows more substructure in the upper panel of Figure 2 between 100-1000d. We can see our ansatz using Equation (11) tracks a spectral feature below E_max in the CR spectrum for an intermediate time period, which fulfills the same condition. After 1000d, our approach finds the right cutoff energy.

Simulations presented in Paper I suggest that the nonlinear instability only undergoes between three and five growth cycles and, thus, it is not able to reach its saturation level. This finding is supported by the work of Inoue et al. (2021). The reason here is that particles reside in the upstream in a region called the precursor up to a distance of L=D(E)/v_shock from the shock. The time available to grow turbulence is limited to $\tau =L/v_{\mathrm{shock}}=D(E)/v_{\mathrm{shock}}^2$. We found that E_max was thus limited to sub-PeV energies, whereas many other works found higher energies. In the case of interaction with dense circumstellar shells, E_max is now boosted by three mechanisms:

1. The shock-shell interaction slows down the shock considerably and suddenly enhances the precursor-scale D(E)/v_sh. The time available to grow turbulence in the precursor is enhanced. After the shock passes through the shell, the precursor scale decreases again and the shock runs through a medium with a now pre-amplified field, boosting E_max.

2. The shell interaction boosts the CR current and, hence, the magnetic field amplification during and after the shell interaction. The densities are higher for earlier interactions and, hence, the E_max achieved is also higher.

3. The collision of the forward-shock with a dense shell creates reflected shocks, which can be rereflected at the contact discontinuity and catch up with the forward shock from behind. Such interactions with the subsequent sharp increases in the shock velocity are shown in the lower panel of Figure 1. These interactions enhance the forward shock speed and slightly boosts E_max, as seen around a few hundred days for the scenario with the earliest shell interaction. Similar effects have been described earlier for reflected shocks produced from the interaction of the SNR shock with the termination shocks of the progenitor star's wind bubble (Sushch et al. 2022).

In our case, the shell interaction enhances the available number of growth cycles for a brief period of time, while the shock accelerates after passing through the shell. This brings our estimates of E_max closer to the values obtained by approaches that assume a steady state or saturation of the nonlinear instability.

3.2. Energy in comic rays

We made no a priori assumption on the fraction of the explosion energy that gets converted into CRs and fixed only the fraction of the thermal plasma particles that get injected as CRs at the shock. The total energy in CRs is an outcome of our simulation and depends on the energy flux through the shock or, alternatively, how the shock converts kinetic to thermal energy. We can easily find that the energy passing through the shock is given by³

$$E_{\mathrm{acc}}=\int \frac{1}{2}{\rho }_u{v}_{\mathrm{shock}}^{2}4\pi R_{\mathrm{shock}}^2{v}_{\mathrm{shock}}\text{d}t\propto t^{3a-2},$$(13)

where ρ_u is the upstream density, v_shock the shock velocity, R_shock the shock radius, and a the expansion parameter, where R_shock∝t^a. The point where E_acc roughly equals the explosion energy marks the transition to the Sedov stage of the remnant's evolution. In the first few years, E_acc is only small fraction of the 10⁵¹ erg of explosion energy.

Figure 3 shows the evolution of the total energy in CRs. Within the first day, the energy in CRs rises fast and then continues to increase roughly as ∝t, following from Equation (13) as our initial expansion parameters is a∼1. However, the interaction with the shells breaks this power-law dependence, greatly increasing E_acc and, hence, the energy in CRs.

Fig. 3. Total energy in CRs in units of 1051 erg for the LBV interaction simulations with different shell radii, as a function of time.

The early interaction case reaches a peak energy-fraction of ∼2%, marking the case with the highest conversion efficiency among our Type IIn cases. The increase of available energy and the increased total energy in CRs supports the strong increase of E_max. We note, however, that the yield in reality might be lower than our predictions and our predictions have to be considered as upper limits. The medium around the SN can feature strong inhomogeneities itself or the explosions can be asymmetric as observed for many core collapse explosions. This reduces and/or smears the effects resulting from the shell interactions that we discussed in this paper over longer time periods.

3.3. Magnetic field evolution

The most crucial aspect with regard to the acceleration of particles is the evolution of the turbulent field. We calculated the component of the turbulent field from an integral over the turbulence spectrum and plotted the time evolution in Figure 4.

Fig. 4. Log ratio of the amplified magnetic field in the immediate upstream relative to the field of a simulation without any shell interactions, shown as the solid lines show. The vertical dashed red lines mark from left to right the times of the first shock-shell interaction, the time of when the shock reached the peak density of the shell, and the time when the shock is propagating in the free wind again for the simulation with the earliest shell interaction (strong line). The pale lines show the behavior for later shell interactins.

Figure 4 shows the amplified field relative to the amplified field of a simulation without a shell interaction. As shell interaction begins (the first red line from the left marks the point where the preshock density increased by 3×), there is an initial surge in the field immediately downstream of the shock. This is partially due to the increasing CR current, but it also partially arises because particles of the highest energies are no longer confined and can escape the remnant from deep downstream. This escape creates a limited spike and the field drops again as the importance of the particles escaping from deep down fades compared to particles that get injected freshly into the upstream.

After this first spike, the field increases again till after the peak of the downstream density (marked by the middle red line in Figure 4) is reached. The period of rising field is caused by the steady increase of the CR current due to the increasing density and the fact that the surge of particles escaping during the onset of the shell-interaction pre-amplified the upstream field.

The upstream field strength starts to fall shortly before the shock completely passed the dense shell, when the shock starts to accelerate again. At its peak, the field is roughly a factor of 100 higher than compared to the case without a shell interaction. Many CRs that were injected during the shell interaction are still present and can drive field growth, but the field decays over time along with the density of injected particles. At the time when E_max is reached after roughly 100 days, the upstream field is still a factor of 5 higher compared to the case without a shell interaction. The same structure of field growth can be seen in the scenarios with a later shell-interaction, although the upstream field strength at the end of the simulation and during the peak is lower in all other scenarios.

3.4. Comparison to other models

As described in Section 2.2, there is no model available at the moment that captures all the underlying physics of CR acceleration in a satisfactory manner. As a consequence, our predictions reflect the assumptions and shortcomings of our ansatz.

In general, we have obtained maximum energies that are at the low end of the models in the literature, with values of ≈100 TeV before the shell-wind interactions. This is considerably lower than the values obtained⁴ by Bell et al. (2013, 21 PeV), Bykov et al. (2018, 8 PeV), Marcowith et al. (2018, 60 PeV), Inoue et al. (2021, 7.5 PeV), and Diesing (2023, 2.4–60 PeV). The disparity is due to the fact that we use an injection fraction of CRs that is on the low end of the range assumed in the above-cited works. Before the shell interaction, only 0.01% of the explosion energy is converted into CRs. Additionally, the fact that we are not considering the saturation of the nonresonant instability and that the number of turbulence growth cycles is limited further limits E_max. Our model is most comparable with the work of Inoue et al. (2021); however, on account of turbulence cascading which acts as a damping mechanism, we achieved an E_max value that is about a factor of ≈2 below the values obtained by Inoue et al. (2021). Both works also only consider the driving of turbulence by CRs in the precursor.

The early shock-shell interactions alter a few of the points above. Most importantly, the deceleration of the shock, which enhances CR escape – followed by an acceleration that lets the shock run through a medium where escaping CRs have driven the turbulence – greatly enhances E_max by basically allowing a larger number of growth cycles. As a consequence, our maximum energies become comparable to more optimistic models such as those of Bell et al. (2013) and Bykov et al. (2018). Despite the fact that we are not switching off the driving of turbulence when the (analytical) saturation limit is reached, efficient cascading continues to limit E_max, as the damping rate of turbulence starts to scale with E_w when $E_w\gg B_0^2/8\pi$ (see also Equation (10)). Furthermore, the interactions boost the number of CRs and enhance the energy converted to CRs greatly, reaching up to ∼2% after the interaction. This certainly provides more energy that can eventually be converted to turbulence at larger scales, essential for the enhancement of E_max.

4. Conclusions

We performed numerical simulations of particle acceleration in very young SNRs expanding in a dense CSM, featuring dense shells created by LBV progenitors. We solved time-dependent transport equations of CRs and magnetic turbulence in the test-particle limit alongside the standard gas dynamical equations for CC-SNRs. We derived the CR diffusion coefficient from the spectrum of magnetic turbulence that evolves as a result of driving by the CR pressure gradient, as well as cascading and wave damping.

We found that the maximum proton energy that we observe in our simulations exceeds PeV energies when the shock interaction with a shell of about 2M_⊙ takes place prior to ≈140 days. Later interactions still boost E_max, but not beyond the PeV frontier. The maximum energy as function of the interaction time with the shell follows a power law of the form E_max∝t^β with β=−1.80±0.03. The increase of the maximum energy is driven by a combination of (i) the increased CR current due to the increasing ambient density during the interaction and (ii) the increase and then decrease in the precursor scale caused by the temporary drop in the shock velocity during the shell interaction. It should be emphasized that such a scenario which reaches PeV energies is limited to very young SNRs and is plausible only for a small subset of SNRs. In this sense, our results agree with previous works arguing that SNRs are unlikely to be responsible for most of the CRs around the knee.

The upstream magnetic field is largely boosted when the shock passes through dense material and, simultaneously, the fraction of the explosion energy that is processed through the shock increases. Hence, a much larger fraction of the explosion energy is available for particle acceleration, which will consequently also boost the gamma-ray emission. These effects will be addressed in a follow-up publication.

Data availability

The data underlying this article will be shared on reasonable request to the corresponding author.

Acknowledgments

RB and JM acknowledge funding from an Irish Research Council Starting Laureate Award. JM acknowledges funding from a Royal Society-Science Foundation Ireland University Research Fellowship. This publication results from research conducted with the financial support of Taighde Éireann – Research Ireland under Grant numbers 20/RS-URF-R/3712, 22/RS-EA/3810, IRCLA\2017\83. We acknowledge the SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support (project dsast026c). IS acknowledges funding from Comunidad de Madrid through the Atracción de Talento “César Nombela” grant with reference number 2023-T1/TEC-29126. And finally we thank Martin Pohl and Pasquale Blasi for the useful discussions and suggestions.

References

All Tables

All Figures

Fig. 1. Upper panel: Density profile of the CSM for the two of our LBV progenitor models. The shell parameters were obained from PION simulations. The parameters of the shells used in RATPAC are listed in Table 1. Lower panel: Shock velocity of the models as a function of the shock-position. Upward steps in the shock velocity are from interactions of the forward shock with reflected shocks catching-up from behind.

In the text

	Fig. 2. Top: Maximum energy of protons for the simulated LBV progenitors plotted over the simulation time. Bottom: Maximum energy obtained in our simulation as a function of the onset time of shell interaction for the LBV case (points) and the fit using the function defined in Equation (12).
In the text

	Fig. 3. Total energy in CRs in units of 1051 erg for the LBV interaction simulations with different shell radii, as a function of time.
In the text

Fig. 4. Log ratio of the amplified magnetic field in the immediate upstream relative to the field of a simulation without any shell interactions, shown as the solid lines show. The vertical dashed red lines mark from left to right the times of the first shock-shell interaction, the time of when the shock reached the peak density of the shell, and the time when the shock is propagating in the free wind again for the simulation with the earliest shell interaction (strong line). The pale lines show the behavior for later shell interactins.

In the text