# Origin of Spectral Hardening of Secondary Cosmic-Ray Nuclei

Norita Kawanaka Department of Astronomy, Graduate School of Science, Kyoto
University,
Kitashirakawa Oiwake-cho, Sakyo-ku, Kyoto, 606-8502, Japan  Hakubi Center,
Kyoto University, Yoshida-honmachi, Sakyo-ku, Kyoto, 606-8501, Japan  Shiu-
Hang Lee Department of Astronomy, Graduate School of Science, Kyoto
University,
Kitashirakawa Oiwake-cho, Sakyo-ku, Kyoto, 606-8502, Japan

###### Abstract

We discuss the acceleration and escape of secondary cosmic-ray (CR) nuclei,
such as lithium, beryllium and boron, produced by spallation of primary CR
nuclei like carbon, nitrogen, and oxygen accelerated at the shock in supernova
remnants (SNRs) surrounded by the interstellar medium (ISM) or a circumstellar
medium (CSM). We take into account the energy-dependent escape of CR particles
from the SNR shocks, which is supported by gamma-ray observations of SNRs, to
calculate the spectra of primary and secondary CR nuclei running away into the
ambient medium. We find that if the SNR is surrounded by a CSM with a wind-
like density distribution (i.e., $n_{\rm CSM}\propto r^{-2}$), the spectra of
the escaping secondary nuclei are harder than those of the escaping primary
nuclei, while if the SNR is surrounded by a uniform ISM, the spectra of the
escaping secondaries are always softer than those of the escaping primaries.
Using this result, we show that if there was a past supernova surrounded by a
dense wind-like CSM ($\sim 2.5\times 10^{-3}M_{\odot}~{}{\rm yr}^{-1}$) which
happened $\sim 1.6\times 10^{5}~{}{\rm yr}$ ago at a distance of $\sim
1.6~{}{\rm kpc}$, we can simultaneously reproduce the spectral hardening of
primary and secondary CRs above $\sim 200~{}{\rm GV}$ that have recently been
reported by AMS-02.

cosmic rays – shock waves – supernova remnants

##  1 Introduction

Cosmic-ray (CR) secondary nuclei such as lithium, beryllium, and boron are
considered to be produced via spallation of heavier nuclei such as carbon,
nitrogen, and oxygen, which are mainly produced at Galactic supernova remnants
(SNRs), during their propagation in the interstellar medium (ISM). The amount
of secondary CR nuclei produced per primary nucleus is proportional to the
grammage traversed by the primaries in the ISM. Therefore, the fluxes and
rigidity dependence of secondary CRs are considered to be a probe of the
propagation of CRs in our Galaxy. The boron-to-carbon (B/C) ratio has been
measured up to $\sim{\rm TeV}$ per nucleon by AMS-02 (Aguilar et al., 2016a),
and it has been shown that above the rigidity $R$ of $65~{}{\rm GV}$ the B/C
ratio can be well fitted with a power law, $\propto R^{\Delta}$ with the index
$\Delta\simeq-0.333$. This agrees well with what is expected from the
Kolmogorov turbulence theory, in which $\Delta$ is asymptotically equal to
$-1/3$ (Kolmogorov, 1941), and therefore it is suggested that the spatial
diffusion coefficient $D$ as a function of rigidity is proportional to
$R^{1/3}$. However, AMS-02 later reported that the spectra of CR lithium,
beryllium, and boron deviate from a single power law above $\rm 200~{}{\rm
GV}$ in an identical way and that they harden even more than the primary CRs
(Aguilar et al., 2018) which have already been known to harden at higher
rigidities from measurements by PAMELA (Adriani et al., 2011), CREAM (Ahn et
al., 2010), and AMS-02 (Aguilar et al., 2015a, b, 2017).

Several scenarios have been discussed to account for the deviation of the
index of the secondary-to-primary ratio at high rigidities from what is
expected from the diffusion coefficient in the bulk of the ISM, including the
effect of propagation in the Galaxy (Cowsik & Burch, 2010; Tomassetti, 2012;
Cowsik et al., 2014; Cowsik & Madziwa-Nussinov, 2016; Tomassetti, 2015), re-
acceleration of secondary CRs (Bresci et al., 2019; Yuan et al., 2020),
contribution from different kinds of sources (Kawanaka & Yanagita, 2018;
Boschini et al., 2020; Niu, 2021; Malkov & Moskalenko, 2021) and the effect of
acceleration of secondary CRs in supernova remnants (SNRs) (Mertsch & Sarkar,
2009; Ahlers et al., 2009; Kachelrieß & Ostapchenko, 2013; Mertsch & Sarkar,
2014; Cholis & Hooper, 2014). These models are proposed mainly to explain the
positron fraction measured by PAMELA and AMS-02 (Adriani et al., 2009; Aguilar
et al., 2013), the antiproton-to-proton ratio measured by AMS-02 (Aguilar et
al., 2016b), or the spectral hardening of proton and helium (Adriani et al.,
2011; Aguilar et al., 2015a, b). Especially, in the last scenario, secondary
CR nuclei are produced inside the SNRs via the spallation of primary CR nuclei
being accelerated at the SNR shock, and are accelerated at the same shock.
Since the spectra of secondary CR nuclei produced in this way are harder than
that of primary CR nuclei, if their contributions are comparable to the
background CR fluxes above a certain rigidity, the rigidity dependence of the
secondary-to-primary ratio would deviate from what is expected from simple ISM
diffusion. Very recently, Mertsch et al. (2020) performed a comprehensive
study on this scenario by searching for the best fit parameters to account for
the recent AMS-02 measurements of CR protons, helium, positrons, anti-protons,
and primary/secondary nuclei. However, they discuss the spectra of secondary
CRs that are shock-accelerated and advected to the downstream (i.e., inside
the SNRs), and they do not consider the spectra of CRs escaping from the SNRs,
which are what we observe.

In this paper we discuss the energy-dependent escape of primary and secondary
CRs that are accelerated in the SNR shock, and predict their energy spectrum
assuming that they come from a local SNR. The escape of CR particles from an
SNR has been intensively investigated to interpret the observed CR spectra and
the gamma-ray emission feature of SNRs interacting with molecular clouds
(Ptuskin & Zirakashvili, 2005; Reville et al., 2009; Gabici et al., 2009;
Caprioli et al., 2010; Ohira et al., 2010, 2011, 2012, 2016; Ohira & Ioka,
2011; Drury, 2011; Kawanaka et al., 2011; Bell et al., 2013). If CRs escape
the SNR in an energy-dependent way, the spectrum of escaping CRs would be
steeper than that of CRs trapped inside the SNR (Ohira et al., 2010; Caprioli
et al., 2010), which may account for the steepness of the observed CR spectral
index. Therefore, it is worth investigating theoretically the production,
acceleration, and escape of secondary CRs in an SNR, and comparing them with
observational data. In this work, we consider CR accelaration in SNRs not only
within a uniform ISM but also those surrounded by a dense circumstellar medium
(CSM). Recent intense optical and near-infrared transient searches have
revealed that many SNe show signs of strong interaction between their ejecta
and the CSM surrounding them. There are various classes of CSM-interacting
SNe: Type IIn, Ibn, Ia-CSM, type I superluminous SN, etc (for a comprehensive
review, see Smith, 2017). The origin of such dense CSM is still uncertain, but
many observations indicate that their progenitors had expelled the stellar
material shortly before their explosions, and the mass-loss rates estimated
from their observational properties are up to $\sim 10^{-3}M_{\odot}~{}{\rm
yr}^{-1}$ (assuming their wind velocities are $\sim 100~{}{\rm km}~{}{\rm
s}^{-1}$), continuing for decades (Moriya et al., 2014). The non-thermal
emissions of photons or neutrinos due to the strong interactions between an SN
ejecta and a dense CSM have been investigated in some studies (Zirakashvili &
Ptuskin, 2016; Murase et al., 2019; Wang et al., 2019; Matsuoka et al., 2019;
Cristofari et al., 2020; Matsuoka & Maeda, 2020), and the nature of such a
system as the origin of CRs has also been investigated (Murase et al., 2011,
2014). However, the production and acceleration of secondary CRs in such a
system, as well as their escape into the ISM, have not been discussed yet.

In Section 2, we describe our model of CR production inside a SNR surrounded
by a CSM and energy-dependent escape of primary and secondary CRs. In Section
3 we solve the transport equations for primary and secondary CR nuclei to
evaluate the spectra of the escaping CRs, compare them with observations, and
discuss our results in light of stellar evolution scenarios. We conclude our
work in Section 4.

##  2 Model

###  2.1 Energy-dependent Cosmic-ray Escape Scenario

Here we describe how CR particles accelerated at the SNR shock escape the
source into the ISM depending on their energy, following Ohira et al. (2010).
In the context of the diffusive shock acceleration (DSA) theory (Blandford &
Eichler, 1987), particles can go back and forth across the shock front and
gain kinetic energy because they are scattered by the turbulent magnetic field
around the shock. Especially, the turbulence in the upstream (i.e. the outside
of the SNR shock) may be amplified by the streaming instability caused by the
accelerating particles themselves (Bell, 1978; Lucek & Bell, 2000). Since the
turbulence is generated only in the vicinity of the shock front, if
accelerated particles have sufficiently high energy, they can escape the SNR
into the far upstream region without being trapped by the turbulent magnetic
field. Using the diffusion coefficient of accelerated particles as a function
of their momentum $D(p)$, the shock velocity $u_{\rm sh}$, and the distance
from the shock beyond which the turbulence is negligible $l$, the escape
condition for a particle accelerated at the shock can be described as

$\displaystyle\frac{D(p)}{u_{\rm sh}}\gtrsim l,$ (1)

where the left hand side represents the diffusion length of a particle with
momentum $p$. Zirakashvili & Ptuskin (2008) investigated the generation of
magnetohydrodynamic (MHD) turbulence driven by the non-resonant streaming
instability using numerical MHD calculations, and show that $l$ should be the
same order as the radius of the SNR, $R_{\rm sh}$.

As a SNR evolves with time, the magnetic field around the shock decays and the
speed of the SNR shell slows down. As a result, the escape condition Eq.(1)
also evolves with time. Especially, if the diffusion length of a particle with
a specific momentum increases faster than the SNR radius ($\sim l$), the
required momentum of particles for escaping the SNR decreases with time. In
other words, CR particles accelerated at the SNR shock can escape the SNR to
far-upstream in an energy-dependent way. In the following subsections, we will
evaluate the energy distribution of CR particles produced in a SNR taking into
account the spallation of primary CRs, the production of secondary CRs, and
their energy-dependent escape.

###  2.2 Primary and Secondary Cosmic-rays in the Supernova Remnant Shock

Here we present our formalism to derive the energy distributions of primary
and secondary CRs produced in the SNR and their escaping fluxes as functions
of time. In the following discussion, we focus on the CR nuclei of lithium,
beryllium, boron, carbon, nitrogen, and oxygen. Among them, only oxygen is
regarded as a pure primary element, while other lighter elements are partly or
entirely produced via spallation of heavier nuclei during their propagation.
Such an approximation is often used in calculating the fluxes of light nuclei
(Mertsch & Sarkar, 2009; Niu, 2021).

We assume that the CR particles can be regarded as test particles during DSA
in a SNR. Letting the shock front be at $x=0$, the stationary transport
equation for the distribution functions of CR nuclei $f_{i}(x,p)$ ($i$
represents the type of nuclei) in the shock rest frame is

$\displaystyle u(x)\frac{\partial f_{i}}{\partial x}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial x}\left[D_{i}(p)\frac{\partial
f_{i}}{\partial x}\right]+\frac{p}{3}\frac{du}{dx}\frac{\partial
f_{i}}{\partial p}-\Gamma_{i}f_{i}+q_{i}$ (2)
$\displaystyle+u_{-}Q_{i}\delta(x)\delta(p-p_{0}),$

where $u(x)$ is the fluid velocity, $D_{i}(p)$ is the diffusion coefficient
for a nuclei of $i$ with momentum $p$, $\Gamma_{i}$ is the total spallation
rate of a nuclei $i$ (i.e. $\Gamma_{i}=\sum_{i>j}\Gamma_{i\rightarrow j}$),
$q_{i}$ is the source term due to the spallation of parent particles, and
$Q_{i}$ is the injection rate of a nuclei $i$ at the shock front (the
injection momentum is $p_{0}$). Considering that the kinetic energy per
nucleon of a nucleus is conserved before and after the spallation, $q_{i}$ is
given by

$\displaystyle 4\pi
p^{2}\frac{dp}{d\varepsilon_{k}}q_{i}(p)=\sum_{i<j}\Gamma_{j\rightarrow
i}N_{j}(\varepsilon_{k}),$ (3)

where $N_{j}(\varepsilon_{k})=4\pi p^{2}f_{j}(p)(dp_{j}/d\varepsilon_{k})$
($p_{j}$ is the momentum of a nucleus $j$ with kinetic energy of
$\varepsilon_{k}$) is the kinetic energy distribution function of nuclei $j$,
and $\Gamma_{j\rightarrow i}$ is the rate at which a nucleus $i$ is produced
via spallation of a heavier nucleus $j$. Here we adopt the table of the
spallation cross sections in Berezinskii et al. (1990). The fluid velocity is
given by

$\displaystyle u(x)=\begin{cases}u_{-}&(x<0)\\\ u_{+}&(x>0),\end{cases}$ (4)

where $u_{-}=u_{\rm sh}$ and $u_{+}=u_{\rm sh}/r$ are constants, and $r$ is
the shock compression ratio, which is assumed to be equal to 4 (the strong
shock limit, ignoring non-linear effects from CR feedback) hereafter.

We then solve the transport equation (2) by imposing the following boundary
conditions:

$\displaystyle{\rm(i)}$
$\displaystyle\lim_{x\rightarrow-0}f_{i}=\lim_{x\rightarrow+0}f_{i},$ (5)
$\displaystyle{\rm(ii)}$ $\displaystyle\lim_{x\rightarrow-l}f_{i}=0,$ (6)
$\displaystyle{\rm(iii)}$
$\displaystyle\left|\lim_{x\rightarrow+\infty}f_{i}\right|<\infty,$
$\displaystyle{\rm(iv)}$ $\displaystyle\left[D_{i}(p)\frac{\partial
f_{i}}{\partial x}\right]^{x=-0}_{x=+0}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(u_{+}-u_{-})p\frac{\partial f_{i,0}}{\partial
p}+u_{-}Q_{i}\delta(p-p_{0}),$ (8)

where condition (i) means that the distribution functions should be continuous
across the shock, and condition (ii) means the free escape of CR particles
from the outer boundary. Condition (iv) comes from the integration of Eq.(2)
across the shock front ($x=0$), and it yields the differential equation for
the distribution function at $x=0$ with respect to $p$, $f_{i,0}(p)$.
Hereafter we solve Eq.(2) for the relativistic regime, i.e., the kinetic
energy per nucleon $\varepsilon_{k}$ is greater than a few ${\rm GeV}/{\rm
n}$, so that we can approximate $p\approx A\varepsilon_{k}/c$.

Following Mertsch & Sarkar (2009) we can solve for the energy distribution
function of nuclei $i$, $N_{i}(\varepsilon_{k})=4\pi
p_{i}^{2}f_{i}(p)(dp_{i}/d\varepsilon_{k})$ separately in the downstream
($x>0$) and the upstream ($x<0$), where we can neglect the second and fifth
terms on the right hand side of Eq.(2). In the downstream and upstream, the
solutions of Eq.(2) are described respectively as

$\displaystyle N_{i}^{+}$ $\displaystyle=$ $\displaystyle\sum_{j\geq
i}E_{ji}e^{\lambda_{j}x/2},$ (9) $\displaystyle N_{i}^{-}$ $\displaystyle=$
$\displaystyle\sum_{j\geq i}F_{ji}e^{\kappa_{j}x/2}+G_{i},$ (10)

where

$\displaystyle\lambda_{j}$ $\displaystyle=$
$\displaystyle\frac{u_{+}}{D_{j}}\left(1-\sqrt{1+4D_{j}\Gamma_{j}/u_{+}^{2}}\right),$
(11) $\displaystyle\kappa_{j}$ $\displaystyle=$
$\displaystyle\frac{u_{-}}{D_{j}}\left(1+\sqrt{1+4D_{j}\Gamma_{j}/u_{-}^{2}}\right),$
(12)

and $E_{ji}$ and $F_{ji}$ are determined recursively as

$\displaystyle E_{ji}$ $\displaystyle=$ $\displaystyle\frac{-4\sum_{m\leq
j}\Gamma_{m\rightarrow
i}E_{jm}}{D_{i}\lambda_{j}^{2}-2u_{+}\lambda_{j}-4\Gamma_{i}}~{}({\rm
for}~{}j>i),$ (13) $\displaystyle E_{ii}$ $\displaystyle=$ $\displaystyle
N_{i,0}-\sum_{j>i}E_{ji},$ (14) $\displaystyle F_{ji}$ $\displaystyle=$
$\displaystyle\frac{-4\sum_{m_{\leq}j}\Gamma_{m\rightarrow
i}F_{jm}}{D_{i}\kappa_{j}^{2}-2u_{-}\kappa_{j}-4\Gamma_{i}}~{}({\rm
for}~{}j>i),$ (15) $\displaystyle F_{ii}$ $\displaystyle=$
$\displaystyle\frac{N_{i,0}-\sum_{j>i}F_{ji}(1-e^{-\kappa_{j}l/2})}{1-e^{-\kappa_{i}l/2}},$
(16) $\displaystyle G_{i}$ $\displaystyle=$ $\displaystyle
N_{i,0}\left(1-\frac{1}{1-e^{-\kappa_{i}l/2}}\right)-\sum_{j>i}F_{ji}\frac{e^{-\kappa_{j}l/2}-e^{-\kappa_{i}l/2}}{1-e^{-\kappa_{i}l/2}},$
(17)

where $N_{i,0}=4\pi p_{i}^{2}f_{i,0}(p_{i})(dp_{i}/d\varepsilon_{k})$. The
differential equation for the distribution function at the shock front
$f_{i,0}$ can be derived from the boundary condition (iv) as

$\displaystyle p\frac{\partial f_{i,0}}{\partial p}$ $\displaystyle=$
$\displaystyle-\frac{3D_{i}}{2(u_{-}-u_{+})}\left[\left(\frac{\kappa_{i}}{1-e^{-\kappa_{i}l/2}}-\lambda_{i}\right)f_{i,0}\right.$
(18)
$\displaystyle\left.+\sum_{j>i}\left\\{\tilde{F}_{ji}\left(\kappa_{j}-\frac{1-e^{-\kappa_{j}l/2}}{1-e^{-\kappa_{i}l/2}}\kappa_{i}\right)-\tilde{E}_{ji}\left(\lambda_{j}-\lambda_{i}\right)\right\\}\right]$
$\displaystyle+\frac{3u_{-}}{u_{-}-u_{+}}Q_{i}\delta(p-p_{0}),$

where

$\displaystyle\tilde{E}_{ji}$ $\displaystyle=$
$\displaystyle\frac{E_{ji}}{4\pi p_{i}^{2}(dp_{i}/d\varepsilon_{k})},$ (19)
$\displaystyle\tilde{F}_{ji}$ $\displaystyle=$
$\displaystyle\frac{F_{ji}}{4\pi p_{i}^{2}(dp_{i}/d\varepsilon_{k})}.$ (20)

The difference between Mertsch & Sarkar (2009) and this study is the position
of the (effective) escape boundary: the former assumes that $f_{i}$ should
damp at $x=-\infty$, while we impose the outer boundary condition $f_{i}=0$ at
a finite distance, $x=-l$. When the acceleration timescale is much shorter
than the spallation timescale (i.e., $\Gamma_{i}D_{i}/u_{-}^{2}\ll 1$) and the
escape boundary is very far from the shock front (i.e., $e^{-\kappa_{i}l/2}\ll
1$), this equation is asymptotically identical to Eq. (17) in Mertsch & Sarkar
(2009).

We can evaluate the escape flux of CR nuclei from the SNR at $x=-l$ as

$\displaystyle\phi_{i}(p)$ $\displaystyle=$ $\displaystyle
u_{-}f_{i}|_{x=-l}-D_{i}(p)\left.\frac{\partial f_{i}}{\partial
x}\right|_{x=-l}$ (21) $\displaystyle=$ $\displaystyle-D_{i}(p)\sum_{j\geq
i}\frac{\kappa_{j}}{2}\tilde{F}_{ji}e^{-\kappa_{j}l/2},$

which is the function of time through the size of the escape boundary, $l$.
The energy spectrum of CR particles escaping the SNR per unit time is given by

$\displaystyle\frac{dN_{{\rm esc},i}(\varepsilon_{k})}{dt}\simeq 4\pi R_{\rm
sh}^{2}\cdot
4\pi\left(\frac{A\varepsilon_{k}}{c}\right)^{2}\left|\phi_{i}(p)\right|\frac{dp_{i}}{d\varepsilon_{k}},$
(22)

where $R_{\rm sh}$ is the radius of the SNR shock.

To understand the nature of the escape flux $\phi_{i}(p)$, let us assume that
the $i$-nuclei are purely primary (i.e., they are not produced via spallation
of heavier nuclei), and that the loss due to spallation is negligible. In this
case, the escape flux can be described as

$\displaystyle-\phi_{i}(p)=\frac{u_{-}f_{i,0}(p)}{\exp(u_{-}l/D_{i})-1},$ (23)

Here we assume Bohm-type diffusion, in which the mean free path of a charged
particle is proportional to its Larmor radius, inside the SNR:

$\displaystyle D_{i}(p)=\eta_{g}\frac{c^{2}p}{3ZeB},$ (24)

where $B$ is the magnetic field strength and $\eta_{g}$ is the gyro factor. We
can see that the momentum at which the absolute value of the escape flux of
the $i$-nuclei attains its maximum value is given by

$\displaystyle p=p_{i,m}\equiv\frac{3u_{-}lZeB}{\gamma c^{2}\eta_{g}},$ (25)

and the maximum value of the escape flux is

$\displaystyle-\phi_{i}(p_{i,m})=\frac{u_{-}f_{i,0}(p_{i,m})}{e^{\gamma}-1},$
(26)

(Caprioli et al., 2009; Ohira et al., 2010). Noting that the distribution
function at the shock front behaves as $f_{i_{0}}\propto p^{-\gamma}$ when
$p\ll p_{i,m}$ where $\gamma=3u_{-}/(u_{-}-u_{+})=3r/(r-1)$, and that
$f_{i,0}(p)\propto\exp(-p/p_{i,m})$ when $p\gg p_{i,m}$, the particles
accelerated at the SNR shock with momentum $p_{i,m}$ can escape the SNR most
efficiently, i.e., $p_{i,m}$ can be regarded as the maximum momentum of the
escaping particles.

We can also consider the particles that are accelerated at the shock but
cannot escape the shock are advected downstream. These particles are confined
in the SNR and are finally released into the ISM as CRs when the SNR is
disrupted. The energy spectrum of this CR component can be evaluated in the
following way. First, the number of the advected particles per unit momentum
per unit time is given by

$\displaystyle\frac{dN_{\rm adv}(p)}{dt}=4\pi R_{\rm sh}^{2}\cdot 4\pi
p^{2}f_{i,0}(p)\frac{u_{\rm sh}}{r},$ (27)

where $u_{\rm sh}/r$ is the fluid velocity in the downstream in the shock rest
frame. The momentum loss of particles due to the expansion of the SNR (i.e.,
adiabatic cooling; Ptuskin & Zirakashvili, 2005) is described as
$dp/dt=-(u_{\rm sh}/R_{\rm sh})p$. Assuming the shell radius expands as
$R_{\rm sh}\propto t^{q}$ where $q$ is a constant, the momentum evolves with
time as $p\propto t^{-q}$. As a result, the energy spectrum of CR particles
that were once advected downstream and are released at the end of the SNR’s
life is given by

$\displaystyle\int_{t_{\rm ini}}^{t_{\rm fin}}dt^{\prime}\frac{dN_{{\rm
adv},i}(\varepsilon_{k})}{dt^{\prime}}$ (28) $\displaystyle\simeq$
$\displaystyle\int_{t_{\rm ini}}^{t_{\rm fin}}dt^{\prime}4\pi R_{\rm
sh}^{2}\cdot
4\pi\left(\frac{A\varepsilon_{k}}{c}\right)^{2}f_{i,0}(p^{\prime})\frac{u_{\rm
sh}}{r}\frac{p^{\prime 2}dp^{\prime}}{p^{2}dp},$

where $p^{\prime}$ is the momentum at the time of injection $t^{\prime}$ of a
particle whose momentum is $p$ at the time of release into the ISM, $t_{\rm
fin}$ (i.e., $p^{\prime}=p(t^{\prime}/t_{\rm fin})^{-q}$).

To evaluate the resulting CR spectra, we need a model for the spatial
diffusion coefficient, $D_{i}(p)$, which in general depends on time because of
the evolution of the magnetic field around the shock, and the primary CR
injection rate at the shock, $Q_{i}$. In the next section, we describe a
phenomenological model to determine these quantities based on analytic
formulae of SNR evolution and the observed spectra of Galactic CRs.

###  2.3 Evolution of SNRs and CR fluxes

Neglecting the back reaction of particle acceleration, the dynamics of the SNR
shock can be well described by analytic solutions. When an SNR is surrounded
by the ISM or a CSM, the SNR expands freely until it sweeps up a mass
comparable to its own mass from the surrounding medium, and afterward the SNR
shell is decelerated and expands in a self-similar way. This is so called the
Sedov-Taylor phase, and the evolution of the SNR shell radius is determined by
the explosion energy of a supernova $E_{\rm SN}$, the ejecta mass $M_{\rm
ej}$, and the number density of the ambient medium $n$. In the case with the
uniform ISM, the SNR shell radius is given by

$\displaystyle R_{\rm sh}=R_{\rm S}\left(\frac{t}{t_{\rm S}}\right)^{2/5}$
(29)

where $R_{\rm S}$ and $t_{\rm S}$ are the Sedov radius and Sedov time,
respectively, and they are given by

$\displaystyle R_{\rm S}$ $\displaystyle=$ $\displaystyle 4.59~{}{\rm
pc}\left(\frac{M_{\rm ej}}{1M_{\odot}}\right)^{1/3}\left(\frac{n}{0.1~{}{\rm
cm}^{-3}}\right)^{-1/3},$ (30) $\displaystyle t_{\rm S}$ $\displaystyle=$
$\displaystyle 450~{}{\rm yr}\left(\frac{E_{\rm SN}}{10^{51}~{}{\rm
erg}}\right)^{-1/2}\left(\frac{M_{\rm
ej}}{1M_{\odot}}\right)^{5/6}\left(\frac{n}{0.1~{}{\rm
cm}^{-3}}\right)^{-1/3},$ (31)

while in the case that the ambient density profile is wind-like (i.e.
$n\propto r^{-2}$), the SNR shell radius is given by

$\displaystyle R_{\rm sh}=R_{\rm S}\left(\frac{t}{t_{\rm S}}\right)^{2/3}$
(32)

where

$\displaystyle R_{\rm S}$ $\displaystyle=$ $\displaystyle 1.26~{}{\rm
pc}\left(\frac{M_{\rm
ej}}{1M_{\odot}}\right)\left(\frac{\dot{M}}{10^{-3}M_{\odot}~{}{\rm
yr}^{-1}}\right)^{-1}\left(\frac{v_{w}}{100~{}{\rm km}~{}{\rm
s}^{-1}}\right),$ (33) $\displaystyle t_{\rm S}$ $\displaystyle=$
$\displaystyle 178~{}{\rm yr}\left(\frac{E_{\rm SN}}{10^{51}~{}{\rm
erg}}\right)^{-1/2}\left(\frac{M_{\rm
ej}}{1M_{\odot}}\right)^{3/2}\left(\frac{\dot{M}}{10^{-3}M_{\odot}~{}{\rm
yr}^{-1}}\right)^{-1}$ (34) $\displaystyle\times\left(\frac{v_{w}}{100~{}{\rm
km}~{}{\rm s}^{-1}}\right),$

where $\dot{M}$ and $v_{w}$ are the mass loss rate and wind velocity,
respectively. Here we use the density profile for a steady wind, i.e.,
$n=\dot{M}/(4\pi m_{p}v_{w}r^{2})$.

We here adopt a phenomenological model proposed by Gabici et al. (2009) (see
also Ohira et al., 2010), which is based on the assumption that Galactic SNRs
are responsible for CRs with energy below the knee ($\sim 10^{15.5}~{}{\rm
eV}$). In this model it is assumed that the maximum momentum of a CR particle
$p_{i,m}$ accelerated at the SNR is limited by its escape during the Sedov-
Taylor phase, and that the maximum energy of an escaping particle, $\sim
cp_{i,m}$, at the beginning of the Sedov phase is equal to the knee energy
($\sim 10^{15.5}~{}{\rm eV}$). Using a variable $\chi$ to describe the
evolution of an SNR (e.g., the SNR’s age, the shock radius of an SNR, etc.),
we assume that the maximum momentum of an escaping CR particle decreases with
time as $p_{i,m}\propto\chi^{-\alpha}$, and that the spectrum of CRs
accelerated at the SNR is proportional to $\chi^{\beta}p^{-s}$ (i.e., CRs
inside the SNR have a power-law spectrum with index of $s$ and their total
number increases with time). One can then see that the spectrum of CRs
escaping the SNR is proportional to $p^{-s+\beta/\alpha}$. This means that if
$\alpha$ and $\beta$ are positive (i.e., the maximum momentum decreases and
the total number of accelerated CRs increases with time), the spectrum of
escaping CRs become steeper than that inside the SNR. Here $\alpha$ and
$\beta$ are the phenomenological parameters that are determined so that the
resulting CR spectra are consistent with observations. Following Ohira et al.
(2010), we hereafter adopt $R_{\rm sh}$ as $\chi$, and fix the parameters as
$\alpha=6.5$ and $\beta=2.1$. Actually, with these parameters, the maximum
energy of escaping CRs evolves from the knee energy at the beginning of the
self-similar phase to $1~{}{\rm GeV}$ at the end of the Sedov phase, when the
radius of the SNR shell becomes 10 times larger than $R_{S}$, and the spectral
index of escaping primary CRs would be similar to what is inferred from
observations, assuming that the diffusion coefficient in the ISM, $D_{\rm
ISM}(R)$, is proportional to a power of rigidity, $R$.

One can determine the diffusion coefficient and the normalization factor of
the escaping CR flux in the following way. First, the time dependence of
$p_{i,m}$ ($\propto\chi^{-\alpha}$) is related to the evolution of the
diffusion coefficient, which depends on the amplification and decay of
magnetic field and the turbulence around the shock (Ptuskin & Zirakashvili,
2003, 2005; Yan et al., 2012). Assuming the Bohm diffusion inside the SNR, one
can determine the evolution of the diffusion coefficient in the SNR using the
relation in Eq.(25). Assuming $l\propto R_{\rm sh}$, we can describe the
evolution of the diffusion coefficient at a specific momentum as
$\propto\chi^{\alpha-1/2}$ in the case with uniform ISM, while
$\propto\chi^{\alpha+1/2}$ in the case with the wind-like ambient density
profile. Second, the normalization factor of the spectrum of CRs accelerated
at the SNR ($\propto\chi^{\beta}$) is related to the CR injection rate at the
shock, $Q_{i}$. If we assume that $Q_{i}$ evolves with time as
$\propto\chi^{\beta^{\prime}}$ where $\beta^{\prime}$ is a constant, since the
escaping CR spectrum per unit time is expressed as Eq.(22), where $\phi(p)$ is
proportional to $u_{\rm sh}f_{i,0}$, one can see that the time-integrated
spectrum of escaping CRs would be

$\displaystyle\int dt\frac{dN_{{\rm esc},i}}{dt}$ $\displaystyle\propto$
$\displaystyle\chi^{\beta^{\prime}+3}p_{i,m}^{\gamma-2}$ (35)
$\displaystyle\propto$ $\displaystyle
p^{-(\gamma-2)-\frac{\beta^{\prime}+3}{\alpha}},$

where we use $p_{i,m}\propto\chi^{-\alpha}$. Considering that the power-law
index of the on-site CR spectrum is $\sim\gamma-2$, we can see how $Q_{i}$
should increase with time as

$\displaystyle\beta^{\prime}=\beta-3.$ (36)

###  2.4 Observed CR Spectra from a Local Single SNR

In this study, we try to explain the observed CR spectral hardening above
$\sim 200~{}{\rm GV}$ by introducing a single local SNR surrounded by the ISM
or a dense CSM. The propagation of CR particles in the ISM can be described by
the diffusion equation,

$\displaystyle\frac{\partial}{\partial t}f_{i,{\rm ISM}}(r,p,t)=D_{\rm
ISM}(p)\nabla^{2}f_{i,{\rm ISM}}(r,p,t)+q_{i}(r,p,t),$ (37)

where $f_{i,{\rm ISM}}(r,p,t)$ is the distribution function of CR particles in
the ISM at a distance $r$ from the source and the time $t$, with momentum $p$,
and $q_{i}(r,p,t)$ is the CR injection rate from the source, which is located
at $r=0$ (i.e., $q_{i}\propto\delta(r)$). In our case, this diffusion equation
can be solved as

$\displaystyle f_{i,{\rm ISM}}(r,p,t)=\int_{0}^{t}dt^{\prime}\frac{4\pi R_{\rm
sh}(t^{\prime})^{2}p^{2}\left|\phi_{i}(p)\right|}{\left(4\pi D_{\rm
ISM}t\right)^{3/2}}\exp\left(-\frac{r^{2}}{4D_{\rm ISM}t}\right).$ (38)

. The observed spectra can then be obtained by summing this distribution
function and the background flux due to the myriad of sources in the Galaxy.

##  3 Results and Discussion

###  3.1 Time integrated energy spectra of CR nuclei escaping the SNR

Figure 1: Time-integrated energy spectra of CR nuclei escaping the SNR
surrounded by a uniform interstellar medium (solid lines) and those once
advected downstream and confined in the SNR (dashed lines). The ambient
density is assumed to be $0.1~{}{\rm cm}^{-1}$. Figure 2: Time-integrated
energy spectra of CR nuclei escaping the SNR surrounded by a wind-like
circumstellar medium (solid lines) and those once advected downstream and
confined in the SNR (dashed lines). The mass loss rate and wind velocity are
assumed to be $3\times 10^{-3}M_{\odot}~{}{\rm yr}^{-1}$ and $100~{}{\rm
km}~{}{\rm s}^{-1}$, respectively.

Fig. 1 depicts the time-integrated spectra of CRs that have escaped the SNR
surrounded by a uniform ISM during the Sedov phase and those that had been
advected downstream and confined in the SNR until it was disrupted at the end
of the Sedov phase111Here we assume that the advected CRs are released
instantaneously when the SNR is disrupted, which may overestimate the effect
of adiabatic cooling.. Here the ISM density is assumed as $0.1~{}{\rm
cm}^{-3}$). Fig. 2 depicts the same spectra but in the case that the SNR is
surrounded by a wind-like medium whose mass loss rate and velocity are
$3\times 10^{-3}M_{\odot}~{}{\rm yr}^{-1}$ and $100~{}{\rm km}~{}{\rm
s}^{-1}$, respectively. In both cases, the explosion energy and ejecta mass of
the supernova are assumed to be $10^{51}~{}{\rm erg}$ and 3 $M_{\odot}$
respectively, and the elemental abundance in CRs are assumed to be identical
to that of the background CRs. First of all, we can see that in both cases the
contributions of CRs that had once been advected into the SNR are much lower
than those escaping the SNR during the Sedov phase. This means that we can
consider only the CR particles that escape the SNR during the Sedov phase to
discuss the observed spectral feature of CRs. Second, in the former case, the
spectra of escaping secondary nuclei (lithium, beryllium, and boron) are
softer than those of the primary nuclei (carbon and oxygen), while in the
latter case the spectra of secondary nuclei are a bit harder than those of
primary nuclei. These spectral features and the difference between two cases
can be interpreted in the following way. Since the secondary CR nuclei
escaping the SNR are produced by the spallation of primary CR nuclei being
accelerated at the shock, the number of CR particles escaping the SNR per unit
time is proportional to $u_{-}f_{\rm pr}n_{\rm amb}t_{\rm int}$, where $f_{\rm
pr}$ is the distribution function of their parent nuclei accelerated at the
shock ($\chi^{\beta}$), $n_{\rm amb}$ is the density of the ambient gas
(constant, or $\propto\chi^{-2}$), and $t_{\rm int}$ is the time during which
their parent nuclei interact with the ambient gas before escape. Especially,
in the escape-limited regime, $t_{\rm int}$ is equal to their acceleration
time in the SNR, $\sim D_{i}/u_{-}^{2}$. Since we are assuming Bohm diffusion
in the SNR, $D_{i}\propto p$, the intrinsic spectral index of secondary nuclei
inside the SNR would be $s-1$. On the other hand, the normalization factor of
the secondary CR flux evolves with time as $\propto
D_{i}/u_{-}^{2}\cdot\chi^{\beta}$ in the case of a uniform ISM and $\propto
D_{i}/u_{-}^{2}\cdot\chi^{-2}\cdot\chi^{\beta}$ in the case of a wind-like
ambient medium. Taking into account the time-dependence of the shell expansion
velocity ($u_{-}\propto\chi^{-3/2}$ in the case of uniform ISM and
$u_{-}\propto\chi^{-1/2}$ in the case with the wind-like ambient medium) and
the diffusion coefficient (see the previous section), we can evaluate the
spectral index of the distribution function of secondary CR nuclei escaping
the SNR per unit time as

$\displaystyle(s-1)+\frac{-3/2+\beta+(\alpha+5/2)}{\alpha}=s+\frac{\beta+1}{\alpha},$
(39)

in the case of uniform ISM, and

$\displaystyle(s-1)+\frac{-1/2+\beta-2+(\alpha+3/2)}{\alpha}=s+\frac{\beta-1}{\alpha},$
(40)

in the case of wind-like ambient medium. Now we can see that in the former
case the spectrum of secondary nuclei is softer, while in the latter case it
is harder compared to that of the primary nuclei, regardless of the value of
$\alpha$ or $\beta$. This difference is caused mainly by the difference in the
density profile of the ambient medium. Generally, the CR particles with higher
energy would escape the SNR earlier, so the time for primary CRs with higher
energy to produce secondary CRs due to the interaction with ambient matter
would be shorter. This makes the number of escaping secondary CRs with higher
energy smaller when the ambient matter distribution is uniform. By contrast,
when the SNR is surrounded by a wind-like CSM, the primary CRs escaping into
the ISM earlier can interact with larger amount of matter than those escaping
later, and then the number of secondary CRs produced by them is enhanced. As a
result, the energy spectrum of secondary CRs escaping the SNR would be harder
than that of primary CRs. Note that this spectral hardening of secondary
nuclei is essentially different from what have been studied in Mertsch &
Sarkar (2009) and Mertsch & Sarkar (2014), in that they discussed the
acceleration of secondary CRs produced inside the SNR assuming that the
maximum energy of accelerated CRs is limited by the age of the SNR. In this
scenario, the primary CRs with higher energy can interact with ambient medium
longer and produce more secondary CRs than those with lower energy. However,
they did not take into account the energy-dependent escape of primary and
secondary CRs, or the non-uniform distribution of ambient medium, both of
which have been implied by recent observations of SNRs and SNe.

According to the recent measurements of secondary CR nuclei by AMS-02 (Aguilar
et al., 2018), the spectra of lithium, beryllium, and boron are hardened at
$\sim 200~{}{\rm GV}$ as has been observed for primary CR nuclei, but they
harden more than the primaries. To account for such a feature, we will
introduce the CR contribution from a past local SN. In this context, the case
with a wind-like CSM, which can make the spectra of secondary CRs harder, is
more relevant. In the following discussions, therefore, we will show results
for the wind-like CSM case only. In addition to the elements shown in Figure
2, we also show the spectra of protons and helium nuclei expected from the SNR
in the next section.

###  3.2 Observed spectra of CR nuclei

In our scenario, we assume that the flux of CR nuclei observed by AMS-02 is a
superposition of two components: one is from a local SNR with a wind-like CSM
that is supposed to reproduce the hardenings observed in the CR nuclei spectra
at $\sim 200~{}{\rm GV}$, and the other is the background flux due to average
SNRs without dense CSM (i.e., these SNRs do not produce hard secondary CRs
efficiently inside themselves). As for the first component, we can calculate
its flux using Eq.(38), choosing the parameters such as the distance and age
of the SNR, total CR energy, CSM properties, abundance ratio in CRs, etc. On
the other hand, the background flux of $i$-th nuclei ${\mathcal{N}}_{i}$ is
given recursively by

$\displaystyle{\mathcal{N}}_{i}(\varepsilon_{k})=\frac{\sum_{i<j}\Gamma_{j\rightarrow
i}{\mathcal{N}}_{j}(\varepsilon_{k})+\mathcal{R}N_{{\rm
esc},i}(\varepsilon_{k})}{1/\tau_{{\rm esc},i}(\varepsilon_{k})+\Gamma_{i}},$
(41)

where ${\mathcal{R}}\simeq 0.03~{}{\rm yr}^{-1}$ is the Galactic SN rate and
$\tau_{{\rm esc},i}(\varepsilon_{k})$ is the timescale for an $i$-th nucleus
with kinetic energy per nucleon of $\varepsilon_{k}$ to escape from the
Galaxy, which is modeled using an expression of $D_{\rm ISM}$,

$\displaystyle\tau_{{\rm esc},i}\simeq\frac{H^{2}}{D_{\rm
ISM}(A_{i}\varepsilon_{k}/Z_{i})},$ (42)

where $H\simeq 3~{}{\rm kpc}$ is half the thickness of the Galactic halo. In
our calculation, we assume the diffusion coefficient in the ISM as $D_{\rm
ISM}(R)=D_{0}(R/1~{}{\rm GV})^{0.45}$ where $D_{0}=2\times 10^{28}~{}{\rm
cm}^{2}~{}{\rm s}^{-1}$ (see e.g., Boschini et al., 2018).

Figure 3: Comparisons of our model spectra of CR oxygen, nitrogen, and carbon
with AMS-02 data. The background fluxes (dashed lines) and the contributions
from our hypothetical past CSM-interacting SN (dotted lines) are also shown.
Figure 4: Comparisons of our model spectra of CR boron, beryllium, and lithium
with AMS-02 data. The background fluxes (dashed lines) and the contributions
from our hypothetical past CSM-interacting SN (dotted lines) are also shown.
---
Figure 5: Li/C (left), Be/C (center), and B/C (right) ratios in the observed CRs reported by AMS-02 along with our model predictions. age | $1.6\times 10^{5}~{}{\rm year}$
---|---
distance | $1.6~{}{\rm kpc}$
total energy | $10^{51}~{}{\rm erg}$
ejecta mass | $3~{}M_{\odot}$
mass loss rate | $2.5\times 10^{-3}M_{\odot}~{}{\rm yr}^{-1}$
wind velocity | $100~{}{\rm km}~{}{\rm s}^{-1}$
Table 1: Properties of the past local SN we introduced to explain the AMS-02
data.

Figure 3 depicts the observed spectra of proton, helium, carbon, nitrogen, and
oxygen nuclei, and Figure 4 depicts the observed spectra of lithium,
beryllium, and boron nuclei as functions of rigidity reported by AMS-02,
fitted by our model. Figure 5 depicts the predicted ratios of lithium to
carbon, beryllium to carbon, and boron to carbon in the observed CRs as
functions of rigidity compared with the ratios reported by AMS-02. In making
these plots, we assume a local SNR surrounded by dense wind-like CSM with
parameters summarized in Table 1. As for the elemental abundance in CRs, we
assume that the fractions of carbon, nitrogen, and oxygen are three times as
much as those in the background CRs. We can see that the hard component
appearing above $\sim 200~{}{\rm GV}$ in each spectrum is well explained by
introducing this local SNR. The important point here is that we can fit not
only the spectra of primary CRs (i.e., oxygen and carbon) but also that of
secondary CRs (i.e., lithium, beryllium and boron) by taking into account the
production, acceleration, and escape of secondary CRs inside a local SNR
surrounded by a dense CSM. Another remarkable thing is that our model can be
tested by the secondary to primary ratios at rigidity of $\gtrsim 1~{}{\rm
TeV}$: within our model, they are predicted to show rising or flattening with
rigidity. Such a feature is never expected as long as we consider the ISM
propagation effects as the origin of spectral hardening of CR nuclei. As shown
in Section 3.1, the energy spectra of secondary CR nuclei escaping a SNR
surrounded by wind-like medium would be harder than those of the primary CR
nuclei. Since the observed CR spectra above $\sim 200~{}{\rm GV}$ are
dominated by the contribution from a hypothetical SNR in our model, the
predicted secondary-to-primary ratios in this rigidity range would reflect the
hardening of secondary CRs escaping the SNR. The secondary-to-primary ratios
above $\sim~{}{\rm TeV}$ will be explored by future experiments such as
AMS-100 (Schael et al., 2019), which will provide a critical test for the
existence of such a local SNR with dense wind-like CSM.

###  3.3 Origin of the local SNR with a dense CSM

While it is beyond the scope of this study to discuss in depth about the
progenitor nature and detailed mass loss mechanism of the proposed local SNR,
it is instructive to briefly account on its possible origin, and the
feasibility of invoking such a local source without violating current
observations. There exists a few possible scenarios for an enhanced mass loss
rate during a certain period prior to core-collapse, such as giant eruptions
at the envelope, envelope stripping in a binary system by Roche-lobe overflow
(RLOF) to a companion star and a possible common envelope phase, and so on
(see, e.g., Smith, 2014). The hypothetical close-by SNR invoked here in
particular can possibly be the result of a past SN of type Ib/c, originating
from the explosion of a Wolf-Rayet (WR) star in a close binary system (e.g.,
Yoon et al., 2010; Dessart et al., 2012, and reference therein). For example,
a RLOF phase that lasts for $\sim 10^{4}$ yrs with a mass loss rate of $\sim
10^{-3}$ M$\mathrm{sun}$/yr and a velocity $\sim 100$ km/s can result in a
dense CSM extending up to a radius $R_{w}\sim 1$ pc. Assuming an ejecta mass
of $3M_{\odot}$ typical for SN Ib/c, and for simplicity a $r^{-2}$ density
profile for the CSM structure and the mass loss parameters in Table 1, the
Sedov radius $R_{S}$ of the SNR is about $0.1$ pc.222We are ignoring any
anisotropy and episodic history of the mass loss from the progenitor here.
This is much smaller than $R_{w}$, and the transition to Sedov phase happened
at an age of roughly 10 yrs old. Therefore, the escape of the highest energy
CRs must have occurred when the SNR blastwave was still in the midst of
interacting with the dense CSM material, which is consistent with our picture
above.

On the other hand, the required chemical enhancement of heavy metal abundances
in the CRs produced by this local accelerator is also consistent with a
stripped envelope SN origin, in that the CSM in the vicinity of the ejecta of
a type Ib/c progenitor is expected to be rich in metals in comparison to that
around a SN IIP from the explosion of a red super-giant star (e.g., Yoon,
2015). The blast wave colliding with this metal-rich CSM shell in the early
phase after SN should then be able to accelerate these heavy ions to high
energies333We note that ejecta of SN Ib/c can reach trans-relativistic
velocities right after explosions, for which particle acceleration at the
shock can be modelled self-consistently using methods such as Monte-Carlo
simulations (e.g., Ellison et al., 2013; Warren et al., 2015). A more accurate
approach with treatment of DSA at trans-relativistic shocks will be done in a
future work.. This interaction between the blastwave and a chemically-enriched
CSM shell has been suggested by recent observations of type Ibn SNe (e.g., see
review in Smith, 2017), or even the hypothetical type Icn that may be detected
in the future. Meanwhile, at a current age of $1.6\times 10^{5}$ yrs, the
shock should have long died out already due to radiative loss, and the SNR
should have merged with and become indistinguishable from the surrounding ISM,
rendering it invisible in various wavebands at a distance $>1$ kpc. A more
detailed model for the evolution and broadband emission of such a stripped-
envelope SNR interacting with a dense CSM using a method similar to Yasuda &
Lee (2019) will be presented in a follow-up paper (Matsuoka et al., in
preparation).

###  3.4 About the nonlinear effects

In this study, we assume that the CR particles in the SNR can be treated as
test particles, and neglect the nonlinear effects from the back reaction of
CRs on the incoming flow that is known to modify the shock structure and
amplify the magnetic fields. In general, we should take into account these
nonlinear effects in the calculation of CR production by DSA coupled to the
hydrodynamics, in a way similar to some existing frameworks in the literature
(e.g., Lee et al., 2012). First, the acceleration efficiency and maximum
energies of CRs are expected to increase due to the nonlinear effects. Second,
shock modification by CR feedback can harden the CR spectra around the maximum
energies to some extent (for reviews see Malkov & Drury, 2001; Schure et al.,
2012). The inclusion of nonlinear effects in our model can therefore
necessitate some adjustments of our parameters to maintain a good fit to the
data. However, even in that case, our main conclusion that shock interaction
of SNR with a dense wind in the CSM can explain the enhanced hardening of the
secondary CR nuclei will not change. In this work, we decide to focus on
demonstrating this new possibility using a simpler test-particle picture, and
reserve a nonlinear description for a follow-up paper.

##  4 Summary

In this paper, we investigated the effect of the production of secondary CR
nuclei at supernova remnants by the diffusive shock acceleration mechanism on
the CR population measured near the Earth using an analytic approach. By
including nuclear spallation effects and acceleration locally at the
acceleration sites, we predicted the spectra of the accelerated ions escaping
from the SNR shocks. Our results show that the secondary CR nuclei escaping
the SNR have a softer spectral index compared to the primary CRs when the SNR
is surrounded by uniform ISM, while they have a harder spectral index when
they are surrounded by the CSM with a wind-like density profile. We show that,
by introducing a past and relatively close-by SN event surrounded by a wind-
like CSM, the current CR measurements including the spectral hardening
recently discovered above a few $100$ GV in ion species up to oxygen can be
successfully reproduced. We suggest that this past local accelerator can be a
SN of type Ib/c with its progenitor enclosed by a CSM from pre-SN mass loss
highly enriched in heavy metals. Our model also predicts a characteristic
spectral flattening of CR secondary-to-primary ratios, such as Li/C, Be/C and
so on, above a rigidity $\sim 1$ TV. Future CR measurements by next-generation
experiments such as AMS-100 with an extended high energy range will put this
to the test.

## Acknowledgement

We are grateful to Kohta Murase for his comments. N.K. acknowledges support by
the Hakubi project at Kyoto University. S.H.L. acknowledges support by JSPS
Grant No. JP19K03913 and the World Premier International Research Center
Initiative (WPI), MEXT, Japan.

## References

  * Adriani et al. (2009) Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2009, Nature, 458, 607, doi: 10.1038/nature07942
  * Adriani et al. (2011) —. 2011, Science, 332, 69, doi: 10.1126/science.1199172
  * Aguilar et al. (2013) Aguilar, M., Alberti, G., Alpat, B., et al. 2013, Physical Review Letters, 110, 141102, doi: 10.1103/PhysRevLett.110.141102
  * Aguilar et al. (2015a) Aguilar, M., Aisa, D., Alpat, B., et al. 2015a, Physical Review Letters, 114, 171103, doi: 10.1103/PhysRevLett.114.171103
  * Aguilar et al. (2015b) —. 2015b, Physical Review Letters, 115, 211101, doi: 10.1103/PhysRevLett.115.211101
  * Aguilar et al. (2016a) Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2016a, Physical Review Letters, 117, 231102, doi: 10.1103/PhysRevLett.117.231102
  * Aguilar et al. (2016b) Aguilar, M., Ali Cavasonza, L., Alpat, B., et al. 2016b, Physical Review Letters, 117, 091103, doi: 10.1103/PhysRevLett.117.091103
  * Aguilar et al. (2017) Aguilar, M., Ali Cavasonza, L., Alpat, B., et al. 2017, Phys. Rev. Lett., 119, 251101, doi: 10.1103/PhysRevLett.119.251101
  * Aguilar et al. (2018) Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2018, Phys. Rev. Lett., 120, 021101, doi: 10.1103/PhysRevLett.120.021101
  * Ahlers et al. (2009) Ahlers, M., Mertsch, P., & Sarkar, S. 2009, Phys. Rev. D, 80, 123017, doi: 10.1103/PhysRevD.80.123017
  * Ahn et al. (2010) Ahn, H. S., Allison, P., Bagliesi, M. G., et al. 2010, ApJ, 714, L89, doi: 10.1088/2041-8205/714/1/L89
  * Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147, doi: 10.1093/mnras/182.2.147
  * Bell et al. (2013) Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415, doi: 10.1093/mnras/stt179
  * Berezinskii et al. (1990) Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., & Ptuskin, V. S. 1990, Astrophysics of cosmic rays
  * Blandford & Eichler (1987) Blandford, R., & Eichler, D. 1987, Phys. Rep., 154, 1, doi: 10.1016/0370-1573(87)90134-7
  * Boschini et al. (2018) Boschini, M. J., Della Torre, S., Gervasi, M., et al. 2018, ApJ, 858, 61, doi: 10.3847/1538-4357/aabc54
  * Boschini et al. (2020) —. 2020, ApJ, 889, 167, doi: 10.3847/1538-4357/ab64f1
  * Bresci et al. (2019) Bresci, V., Amato, E., Blasi, P., & Morlino, G. 2019, MNRAS, 488, 2068, doi: 10.1093/mnras/stz1806
  * Caprioli et al. (2009) Caprioli, D., Blasi, P., & Amato, E. 2009, MNRAS, 396, 2065, doi: 10.1111/j.1365-2966.2008.14298.x
  * Caprioli et al. (2010) Caprioli, D., Kang, H., Vladimirov, A. E., & Jones, T. W. 2010, MNRAS, 407, 1773, doi: 10.1111/j.1365-2966.2010.17013.x
  * Cholis & Hooper (2014) Cholis, I., & Hooper, D. 2014, Phys. Rev. D, 89, 043013, doi: 10.1103/PhysRevD.89.043013
  * Cowsik & Burch (2010) Cowsik, R., & Burch, B. 2010, Phys. Rev. D, 82, 023009, doi: 10.1103/PhysRevD.82.023009
  * Cowsik et al. (2014) Cowsik, R., Burch, B., & Madziwa-Nussinov, T. 2014, ApJ, 786, 124, doi: 10.1088/0004-637X/786/2/124
  * Cowsik & Madziwa-Nussinov (2016) Cowsik, R., & Madziwa-Nussinov, T. 2016, ApJ, 827, 119, doi: 10.3847/0004-637X/827/2/119
  * Cristofari et al. (2020) Cristofari, P., Renaud, M., Marcowith, A., Dwarkadas, V. V., & Tatischeff, V. 2020, MNRAS, 494, 2760, doi: 10.1093/mnras/staa984
  * Dessart et al. (2012) Dessart, L., Hillier, D. J., Li, C., & Woosley, S. 2012, MNRAS, 424, 2139, doi: 10.1111/j.1365-2966.2012.21374.x
  * Drury (2011) Drury, L. O. 2011, MNRAS, 415, 1807, doi: 10.1111/j.1365-2966.2011.18824.x
  * Ellison et al. (2013) Ellison, D. C., Warren, D. C., & Bykov, A. M. 2013, The Astrophysical Journal, 776, 46, doi: 10.1088/0004-637x/776/1/46
  * Gabici et al. (2009) Gabici, S., Aharonian, F. A., & Casanova, S. 2009, MNRAS, 396, 1629, doi: 10.1111/j.1365-2966.2009.14832.x
  * Kachelrieß & Ostapchenko (2013) Kachelrieß, M., & Ostapchenko, S. 2013, Phys. Rev. D, 87, 047301, doi: 10.1103/PhysRevD.87.047301
  * Kawanaka et al. (2011) Kawanaka, N., Ioka, K., Ohira, Y., & Kashiyama, K. 2011, ApJ, 729, 93, doi: 10.1088/0004-637X/729/2/93
  * Kawanaka & Yanagita (2018) Kawanaka, N., & Yanagita, S. 2018, Phys. Rev. Lett., 120, 041103, doi: 10.1103/PhysRevLett.120.041103
  * Kolmogorov (1941) Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301
  * Lee et al. (2012) Lee, S.-H., Ellison, D. C., & Nagataki, S. 2012, ApJ, 750, 156.  https://arxiv.org/abs/1203.3614
  * Lucek & Bell (2000) Lucek, S. G., & Bell, A. R. 2000, MNRAS, 314, 65, doi: 10.1046/j.1365-8711.2000.03363.x
  * Malkov & Drury (2001) Malkov, M. A., & Drury, L. O. 2001, Reports on Progress in Physics, 64, 429, doi: 10.1088/0034-4885/64/4/201
  * Malkov & Moskalenko (2021) Malkov, M. A., & Moskalenko, I. V. 2021, ApJ, 911, 151, doi: 10.3847/1538-4357/abe855
  * Matsuoka & Maeda (2020) Matsuoka, T., & Maeda, K. 2020, ApJ, 898, 158, doi: 10.3847/1538-4357/ab9c1b
  * Matsuoka et al. (2019) Matsuoka, T., Maeda, K., Lee, S.-H., & Yasuda, H. 2019, ApJ, 885, 41, doi: 10.3847/1538-4357/ab4421
  * Mertsch & Sarkar (2009) Mertsch, P., & Sarkar, S. 2009, Physical Review Letters, 103, 081104, doi: 10.1103/PhysRevLett.103.081104
  * Mertsch & Sarkar (2014) —. 2014, Phys. Rev. D, 90, 061301, doi: 10.1103/PhysRevD.90.061301
  * Mertsch et al. (2020) Mertsch, P., Vittino, A., & Sarkar, S. 2020, arXiv e-prints, arXiv:2012.12853.  https://arxiv.org/abs/2012.12853
  * Moriya et al. (2014) Moriya, T. J., Maeda, K., Taddia, F., et al. 2014, MNRAS, 439, 2917, doi: 10.1093/mnras/stu163
  * Murase et al. (2019) Murase, K., Franckowiak, A., Maeda, K., Margutti, R., & Beacom, J. F. 2019, ApJ, 874, 80, doi: 10.3847/1538-4357/ab0422
  * Murase et al. (2011) Murase, K., Thompson, T. A., Lacki, B. C., & Beacom, J. F. 2011, Phys. Rev. D, 84, 043003, doi: 10.1103/PhysRevD.84.043003
  * Murase et al. (2014) Murase, K., Thompson, T. A., & Ofek, E. O. 2014, MNRAS, 440, 2528, doi: 10.1093/mnras/stu384
  * Niu (2021) Niu, J.-S. 2021, Chinese Physics C, 45, 041004, doi: 10.1088/1674-1137/abe03d
  * Ohira & Ioka (2011) Ohira, Y., & Ioka, K. 2011, ApJ, 729, L13, doi: 10.1088/2041-8205/729/1/L13
  * Ohira et al. (2016) Ohira, Y., Kawanaka, N., & Ioka, K. 2016, Phys. Rev. D, 93, 083001, doi: 10.1103/PhysRevD.93.083001
  * Ohira et al. (2010) Ohira, Y., Murase, K., & Yamazaki, R. 2010, A&A, 513, A17, doi: 10.1051/0004-6361/200913495
  * Ohira et al. (2011) —. 2011, MNRAS, 410, 1577, doi: 10.1111/j.1365-2966.2010.17539.x
  * Ohira et al. (2012) Ohira, Y., Yamazaki, R., Kawanaka, N., & Ioka, K. 2012, MNRAS, 427, 91, doi: 10.1111/j.1365-2966.2012.21908.x
  * Ptuskin & Zirakashvili (2003) Ptuskin, V. S., & Zirakashvili, V. N. 2003, A&A, 403, 1, doi: 10.1051/0004-6361:20030323
  * Ptuskin & Zirakashvili (2005) —. 2005, A&A, 429, 755, doi: 10.1051/0004-6361:20041517
  * Reville et al. (2009) Reville, B., Kirk, J. G., & Duffy, P. 2009, ApJ, 694, 951, doi: 10.1088/0004-637X/694/2/951
  * Schael et al. (2019) Schael, S., Atanasyan, A., Berdugo, J., et al. 2019, Nuclear Instruments and Methods in Physics Research A, 944, 162561, doi: 10.1016/j.nima.2019.162561
  * Schure et al. (2012) Schure, K. M., Bell, A. R., O’C Drury, L., & Bykov, A. M. 2012, Space Sci. Rev., 173, 491, doi: 10.1007/s11214-012-9871-7
  * Smith (2014) Smith, N. 2014, Annual Review of Astronomy and Astrophysics, 52, 487, doi: 10.1146/annurev-astro-081913-040025
  * Smith (2017) Smith, N. 2017, Interacting Supernovae: Types IIn and Ibn, ed. A. W. Alsabti & P. Murdin, 403, doi: 10.1007/978-3-319-21846-5_38
  * Tomassetti (2012) Tomassetti, N. 2012, ApJ, 752, L13, doi: 10.1088/2041-8205/752/1/L13
  * Tomassetti (2015) —. 2015, Phys. Rev. D, 92, 081301, doi: 10.1103/PhysRevD.92.081301
  * Wang et al. (2019) Wang, K., Huang, T.-Q., & Li, Z. 2019, ApJ, 872, 157, doi: 10.3847/1538-4357/aaffd9
  * Warren et al. (2015) Warren, D. C., Ellison, D. C., Bykov, A. M., & Lee, S.-H. 2015, MNRAS, 452, 431, doi: 10.1093/mnras/stv1304
  * Yan et al. (2012) Yan, H., Lazarian, A., & Schlickeiser, R. 2012, ApJ, 745, 140, doi: 10.1088/0004-637X/745/2/140
  * Yasuda & Lee (2019) Yasuda, H., & Lee, S.-H. 2019, ApJ, 876, 27, doi: 10.3847/1538-4357/ab13ab
  * Yoon (2015) Yoon, S.-C. 2015, Publications of the Astronomical Society of Australia, 32, e015, doi: 10.1017/pasa.2015.16
  * Yoon et al. (2010) Yoon, S. C., Woosley, S. E., & Langer, N. 2010, ApJ, 725, 940, doi: 10.1088/0004-637X/725/1/940
  * Yuan et al. (2020) Yuan, Q., Zhu, C.-R., Bi, X.-J., & Wei, D.-M. 2020, J. Cosmology Astropart. Phys, 2020, 027, doi: 10.1088/1475-7516/2020/11/027
  * Zirakashvili & Ptuskin (2008) Zirakashvili, V. N., & Ptuskin, V. S. 2008, ApJ, 678, 939, doi: 10.1086/529580
  * Zirakashvili & Ptuskin (2016) —. 2016, Astroparticle Physics, 78, 28, doi: 10.1016/j.astropartphys.2016.02.004