Cosmic-ray Acceleration in Core-Collapse Supernova Remnants with the Wind Termination Shock

In this work, we consider propagation of a core-collapse SNR shock in the CSM with the Parker-spiral magnetic field and current sheet until the SNR shock collides with the WTS. Figure 1 shows the schematic picture of the case of aligned rotators. The inner and outer black circles are the SNR shock and WTS, respectively. The region between the SNR shock and WTS is the free wind region. The SNR shock with the shock velocity, $u_{\rm SNR}$, expands in the free wind region until the SNB shock collides with the WTS. $B_{{\rm w},\phi}$ is the toroidal component of the Parker-spiral magnetic field in the free wind region. The black solid arrow is the rotation axis of the progenitors. $\theta$ and $\phi$ are the polar and azimuthal angles, respectively. The regions at $\theta=0,\pi$ and at $\theta=\pi/2$ are the poles and equator, respectively. The gray line at the equator is the current sheet. The quantities in the wind region are assumed to be spherically symmetric except for the Parker-spiral magnetic field in the free wind region, $\vec{B}_{\rm w}$. We perform test particle simulations to reveal the particle acceleration and attainable maximum energy in this system. The numerical setups in this work are almost the same as that in our previous work (kamijima22, ). The important difference of numerical setups from our previous work is that the structure of the WTS is included although our previous work did not consider the WTS because we focused on the escape process of accelerated particles from the SNR shock. In our test particle simulations, protons with energy much larger than the thermal energy are considered. We focused on WR stars, which are thought to be progenitors of type Ib/Ic supernovae. The WTS of red supergiants (RSGs) is not considered because the WTS of RSGs does not play an important role in this work [see Sec. III.3 for details]. Particles with $100~{}{\rm GeV}$ are uniformly injected on the whole SNR shock surface only at $100~{}{\rm yr}$ after the supernova explosion, $t_{\rm inj}=100~{}{\rm yr}$. To improve the statistics of high-energy particles, the particle splitting method is used in our simulations.

In this work, the magnetic field fluctuation in the free wind region (shock upstream region) is assumed to be zero. Then, the magnetic field in the free wind region consists only of the Parker-spiral magnetic field in our simulations [see Sec. II.3]. In the free wind region, we solve the equation of motion in the explosion center rest frame, $d\vec{u}/dt=(e/m_{\rm p})[\vec{E}_{\rm w}+\vec{u}/(\gamma c)\times\vec{B}_{\rm w}]$, to calculate the particle trajectory, where $\vec{u}$ and $\gamma=\sqrt{1+(u/c)^{2}}$ is the spatial three components of four velocities and Lorentz factor of particles, respectively. $\vec{E}_{\rm w}$ and $\vec{B}_{\rm w}$ are the electric and magnetic fields in the free wind region measured in the explosion center rest frame. Contrary to the free wind region, the magnetic field fluctuation in the downstream region of both the SNR shock and WTS is assumed to be highly turbulent. In the downstream region of both the SNR shock and WTS, we consider only the magnetic field strength although the magnetic field configuration is not considered. The magnetic field strength in the downstream region of both the SNR shock and WTS is given by the condition that a fraction of the upstream kinetic energy flux is converted to the downstream magnetic field energy flux in the shock rest frame [see Sec. II.4]. In the downstream region of both the SNR shock and WTS, the particle trajectory is solved by the Monte-Carlo method. The downstream particle motion is assumed to be the Bohm diffusion under the highly turbulent magnetic field and downstream particles are isotropically scattered in the downstream rest frame. The scattering angle in the downstream region is randomly chosen between 0 and 4$\pi$. The mean free path of downstream particles is given by the downstream gyroradius because the particle motion in the downstream region is assumed to be the Bohm diffusion.

Both the SNR shock and WTS are assumed to be spherical discontinuities. This is because the gyroradius of high-energy protons focused in our simulations is assumed to be larger than the width of the shock. The SNR shock velocity, $u_{\rm SNR}(t)$, is given as follows (chevalier82, ):

where $E_{\rm SN}=10^{51}~{}{\rm erg}$, $M_{\rm ej}=5M_{\odot}$, $\dot{M}=10^{-4}M_{\odot}/{\rm yr}$, and $V_{\rm w}=3000~{}{\rm km/s}$ are the explosion energy, ejecta mass, mass loss rate, wind velocity of WR stars, respectively (crowther07, ; niedzielski02, ). $t$ is the elapsed time from the supernova explosion. The density in the free wind region and ejecta profile are assumed to be

where $n$ is set to be 10 (chevalier82, ). $r$ and $t_{\rm t}$ are the distance from the explosion center and the time when the reverse shock reaches the inner ejecta.

$R_{\rm SNR}=\int^{t}u_{\rm SNR}(t^{\prime})dt^{\prime}$ is the SNR shock radius. As for the velocity profile of the downstream region of the SNR, $u_{\rm d,SNR}(r,t)$, we use the following approximate formula:

where $u_{\rm d,SNR}(R_{\rm SNR},t)=(3u_{\rm SNR}(t)+V_{\rm w})/4$ is derived from the Rankin-Hugoniot relation at the strong shock limit. Particles in the downstream region of the SNR shock lose their energy by the adiabatic cooling because of the expansion of the downstream region of the SNR shock (${\rm div}\vec{u}_{\rm d,SNR}>0$).

The self-similar solution is used as the radius of the WTS (weaver77, ). The radius of the WTS, $R_{\rm WTS}$, is determined by the condition that the ram pressure in the free wind region at the WTS, $\rho_{\rm w}(R_{\rm WTS})V_{\rm w}^{2}$, is equal to the pressure of the shocked wind region, $P(t+t_{\rm life})\approx P(t_{\rm life})$. $t_{\rm life}\sim 10^{5}~{}{\rm yr}$ is the lifetime of WR stars. $t+t_{\rm life}$ is almost the same as $t_{\rm life}$ because $t_{\rm life}$ is much larger than the SNR age. Hence, $R_{\rm WTS}(t_{\rm life})=\sqrt{\dot{M}V_{\rm w}/(4\pi P(t_{\rm life}))}$ is calculated as follows (weaver77, ; dwarkadas07, ):

The time evolution of the WTS can be negligible between the time of the supernova explosion and the time when the SNR collides with the WTS because the dynamical timescale is much shorter than the SNR age. Therefore, $R_{\rm WTS}$ is fixed to be $8~{}{\rm pc}$ in our simulations. The density in the shocked wind region is constant because the pressure in the shocked wind region is constant and the shocked wind region is adiabatic. The following velocity profile in the shocked wind region, $u_{\rm d,WTS}(r)$, is given by the mass flux conservation between the upstream and downstream regions of the WTS (weaver77, ):

$u_{\rm d,WTS}(R_{\rm WTS})=V_{\rm w}/4$ is derived from the Rankin-Hugoniot relation at the strong shock limit. In the shocked wind region, the adiabatic cooling does not work because ${\rm div}\vec{u}_{\rm d,WTS}=0$.

The electromagnetic field in the free wind region (shock upstream region) is shown in this section. As with our previous work kamijima22 , for simplicity, we consider only the Parker-spiral magnetic field as an unperturbed magnetic field. Thus, we do not consider the magnetic field fluctuation in the free wind region. The rotation axis of progenitors is set to be the polar axis of the spherical coordinate. $\theta$ and $\phi$ are the polar and azimuthal angles. $\phi$ is the same direction as the rotation of progenitors. The regions at $\theta=0,\pi$ and at $\theta=\pi/2$ are the polar and equatorial regions, respectively. The magnetic field in the free wind region, $\vec{B}_{\rm w}=B_{{\rm w},r}\vec{e}_{r}+B_{{\rm w},\phi}\vec{e}_{\phi}$, is as follows (weber67, ):

where $H(\theta)$, $\Omega_{*}=2\pi/P_{*}$, $R_{\rm A}$, and $B_{\rm A}$ are the Heaviside step function, angular frequency of progenitors, Alfvén radius, and the magnetic field strength at the Alfvén radius, respectively. $\vec{e}_{r}$ and $\vec{e}_{\phi}$ are unit vectors of $r$ and $\phi$ directions. The rotation period of progenitors, $P_{*}$, is set to be $10~{}{\rm days}$ (chene08, ). $R_{\rm A}$ is approximately given by

where $R_{*}$ is the radius of progenitors (ud-doula08, ). $\eta_{*}=B_{*}^{2}R_{*}^{2}/(\dot{M}V_{\rm w})$ is called as the magnetic confinement parameter (ud-doula08, ). $B_{*}$ is the surface magnetic field strength of progenitors. $q$ is the index of the radial dependency of the magnetic field inside the Alfvén radius. In this work, we assume that the magnetic field configuration inside the Alfvén radius is the dipole magnetic field ($q=3$). The Alfvén radius, $R_{\rm A}$, is

The magnetic field strength at the Alfvén radius, $|B_{\rm A}|$, is

In this work, $B_{*}$ and $R_{*}$ are set to be $100~{}{\rm G}$ and $10R_{\odot}$, respectively (hubrig20, ; crowther07, ; niedzielski02, ). The sign of $B_{\rm A}=\pm B_{*}\left(R_{*}/R_{\rm A}\right)^{q}$ is given by whether the angle between the rotation and magnetic axes, $\alpha_{\rm inc}$, is larger than $\pi/2$ radian. When $\alpha_{\rm inc}\leq(\geq)\pi/2$, the sign of $B_{\rm A}$ is positive (negative). In this work, $\alpha_{\rm inc}$ is set to be $0,\pi/6,5\pi/6$, and $\pi$. Aligned rotators are the case of $\alpha_{\rm inc}=0$ or $\pi$, and oblique rotators are the case of $\alpha_{\rm inc}\neq 0$ or $\pi$. The position of the current sheet, $\theta_{\rm CS}$, is determined as follows (alanko-huotari07, ):

The width of the current sheet is assumed to be infinitely thin because the gyroradius of high-energy protons is assumed to be much larger than the width of the current sheet. The wind velocity is assumed to have only the radial component ($\vec{V}_{\rm w}=V_{\rm w}\vec{e}_{r}$). Hence, in the simulation frame (explosion center rest frame), the electric field, $\vec{E}_{\rm w}=-(\vec{V}_{\rm w}/c)\times\vec{B}_{\rm w}$, emerges.

As we mentioned above, the magnetic field in the downstream region of both the SNR shock and WTS is assumed to be highly turbulent, which is suggested by recent observations and simulations (prajapati19, ; berezhko03, ; bamp, ). The downstream magnetic field strength is determined by the condition that a fraction, $\epsilon_{B}$, of the upstream kinetic energy flux measured in the shock rest frame is converted to the downstream magnetic field energy flux. In this work, $\epsilon_{B}$ is set to be $0.1$ for both the SNR shock and WTS. As for the SNR shock, the upstream kinetic energy flux measured in the SNR shock rest frame is $(1/2)\rho_{\rm w}(R_{\rm SNR})(u_{\rm SNR}-V_{\rm w})^{3}$ and the energy flux of the downstream magnetic field, $B_{\rm d,SNR}$, is $(B_{\rm d,SNR}^{2}/(8\pi))\times(u_{\rm SNR}-V_{\rm w})/4$. Thus, $B_{\rm d,SNR}$ is

As for the WTS, the upstream kinetic energy flux measured in the WTS rest frame is $(1/2)\rho_{\rm w}(R_{\rm WTS})V_{\rm w}^{3}$ and the energy flux of the downstream magnetic field, $B_{\rm d,WTS}$, is $(B_{\rm d,WTS}^{2}/(8\pi))\times V_{\rm w}/4$. Here, the velocity of the WTS is negligible because the velocity of the WTS is much smaller than the wind velocity. Then, $B_{\rm d,WTS}$ is

Particles with the charge, $Ze$, are considered in this section. This section presents the global particle motion between the SNR shock and WTS, the condition for realizing this global particle motion, and the maximum attainable energy. First, we consider the expected global particle motion between the SNR shock and WTS. Figure 2 shows the schematic picture of the expected global particle motion between the SNR shock and WTS for the case of aligned rotators ($\alpha_{\rm inc}=0$). The inner and outer circles are the SNR shock and WTS, respectively. The red arrows are the expected motion of accelerating particles. The black arrows are the mean velocities of accelerating particles along the SNR shock, $v_{\theta,\rm SNR}$, current sheet, $v_{\rm CS}$, WTS, $v_{\theta,\rm WTS}$, and pole, $v_{\rm pl}$. $B_{{\rm w},r}$ and $B_{{\rm w},\phi}$ are the radial and toroidal components of the Parker-spiral magnetic field in the free wind region. When $\alpha_{\rm inc}\leq(\geq)\pi/2$, particles injected on the SNR shock surface are accelerated in perpendicular shocks, and some particles are advected to the far downstream region of the SNR shock with certain probabilities. Accelerating particles that are not advected to the far downstream region of the SNR shock move to the equator (poles) along the SNR shock. For $\alpha_{\rm inc}\leq\pi/2$, particles around the equator move to the WTS along the current sheet while performing the meandering motion (levy75, ; kamijima22, ). For $\alpha_{\rm inc}\geq\pi/2$, particles around the poles move to the WTS along $B_{{\rm w},r}$ because $B_{{\rm w},r}$ is larger than $B_{{\rm w},\phi}$ around the poles. Hence, for both $\alpha_{\rm inc}\leq(\geq)\pi/2$, particles accelerated at the SNR shock escape from the SNR shock towards the WTS (kamijima22, ). Escaped particles eventually reach the WTS. These particles are accelerated at the WTS and some of these particles are advected to the far downstream region of the WTS similar to the SNR shock. Accelerating particles that are not advected to the far downstream region move to the poles (equator) along the WTS. Here, the direction of the particle motion along the SNR shock is opposite to the direction of the particle motion along the WTS. For $\alpha_{\rm inc}\leq\pi/2$, particles that reach the poles return to the SNR shock along $B_{{\rm w},r}$. This is because $B_{{\rm w},r}$ is larger than $B_{{\rm w},\phi}$ around the poles. For $\alpha_{\rm inc}\geq\pi/2$, particles around the equator of the WTS are expected to return the SNR shock while moving along the current sheet (meandering motion). Particles that return from the poles (equator) to the SNR shock can be accelerated at the SNR shock again. These particles move to the equator (poles) and escape from the equator (poles) to the WTS again. Thus, the cyclic particle motion between the SNR shock and WTS is expected.

We estimate the attainable energy of particles accelerated by the cyclic motion between the SNR shock and WTS, $\varepsilon_{\rm cyc}$. $\varepsilon_{\rm cyc}$ is given by

$N_{\rm cyc}$ is the maximum number of cycles until the SNR shock collides with the WTS. $\varepsilon_{\rm SNR}$ ($\varepsilon_{\rm WTS}$) is the energy gain of particles accelerated in the region between the pole and equator of the SNR shock (WTS). $t_{\rm i}$ is the time when the i-th cycle starts.

First, we estimate $\varepsilon_{\rm SNR}$. In the SNR shock rest frame, the flow velocity in the free wind region (shock upstream region) at the time, $t$, is $u_{\rm SNR}(t)-V_{\rm w}$. Then, the motional electric field measured in the SNR shock rest frame, $\vec{E}_{\rm w,SNR}$, emerges in the free wind region:

where $\vec{e}_{\theta}$ and $\langle B_{{\rm w},\phi}\rangle$ are the unit vector of the $\theta$ direction and mean magnetic field strength that particles in the free wind region feel. For the oblique rotator case ($\alpha_{\rm inc}\neq 0,\pi$), $\langle B_{{\rm w},\phi}\rangle$ inside the wavy current sheet region ($\pi/2-\alpha_{\rm inc}\leq\theta\leq\pi/2+\alpha_{\rm inc}$) is

which is derived in our previous paper (kamijima22, ). Outside the wavy current sheet region ($0\leq\theta\leq\pi/2-\alpha_{\rm inc}$ and $\pi/2+\alpha_{\rm inc}\leq\theta\leq\pi$), $\langle B_{{\rm w},\phi}\rangle$ is $B_{{\rm w},\phi}$ in Eq. (14). For the aligned rotator case ($\alpha_{\rm inc}=0,\pi$), $\langle B_{{\rm w},\phi}\rangle=B_{{\rm w},\phi}$. In the SNR shock rest frame, particles that move along the SNR shock are accelerated by $\vec{E}_{\rm w,SNR}$ (kamijima22, ). Therefore, $\varepsilon_{\rm SNR}$ is given by the potential difference between the pole and equator of the SNR shock in the SNR shock rest frame. Thus, $\varepsilon_{\rm SNR}$ is

Here, we ignored the time evolution of the SNR shock until particles injected around the poles reach the equator. This is because the time when particles injected around the pole of the SNR shock reach the equator, $R_{\rm SNR}/v_{\theta,{\rm SNR}}\approx R_{\rm SNR}/(c/2)$, is much shorter than the dynamical timescale of the SNR shock, $R_{\rm SNR}/u_{\rm SNR}$. $v_{\theta,{\rm SNR}}$ is the mean velocity in the $\theta$ direction of accelerating particles [see Eq. (59) in Appendix A].

Next, we estimate $\varepsilon_{\rm WTS}$. Here, we ignored the WTS velocity because the WTS velocity is much smaller than the wind velocity. Then, in the WTS rest frame, the motional electric field in the free wind region, $\vec{E}_{\rm w,WTS}$, is

Similar to $\varepsilon_{\rm SNR}$, $\varepsilon_{\rm WTS}$ is given by the potential difference between the equator and pole of the WTS as follows:

For $\alpha_{\rm inc}\geq\pi/2$, $\varepsilon_{\rm SNR}$ and $\varepsilon_{\rm WTS}$ are given as follows:

$\langle B_{{\rm w},\phi}\rangle$ for $\alpha_{\rm inc}\geq\pi/2$ is the opposite sign of $\langle B_{{\rm w},\phi}\rangle$ for $\alpha_{\rm inc}\leq\pi/2$ (kamijima22, ).

The maximum number of cycles between the SNR shock and WTS, $N_{\rm cyc}$, is estimated in this section. We consider the situation that particles injected on the SNR shock at $t=t_{\rm inj}$ continue to be accelerated until the SNR shock collides with the WTS. Hence, we do not consider particles advected to the far downstream region of the SNR shock or WTS until the SNR collides with the WTS. The SNR shock collides with the WTS at $t=t_{\rm col}$. $N_{\rm cyc}$ is

where $\Delta t_{\rm cyc}$ is one cycle time between the SNR shock and WTS (SNR $\rightarrow$ WTS $\rightarrow$ SNR). $\Delta t_{\rm cyc}$ is given as follows:

where $v_{\theta,\rm SNR},v_{\rm eq},v_{\theta,\rm WTS}$, and $v_{\rm pl}$ are the mean velocities of accelerating particles moving along the SNR shock, equator, WTS, and poles, respectively [see Fig. 2]. As shown in Appendix A, $v_{\theta,\rm SNR},v_{\rm eq},v_{\theta,\rm WTS}$, and $v_{\rm pl}$ are approximately $c/2$ if the residence time in the downstream region of the SNR shock and WTS is much smaller than that in the free wind region (shock upstream region) [see Eqs. (59),(60), and (62) in Appendix A]. In our model, the downstream residence time can be negligible if the downstream magnetic field strength is much larger than the magnetic field strength in the free wind region (kamijima20, ). The first term in Eq. (41) means the elapsed time when accelerating particles move along the SNR shock from the pole to the equator. The second term means the elapsed time when particles move along the equator between the SNR shock and WTS. The third term means the elapsed time when accelerating particles move along the WTS from the equator to the pole. The fourth term means the elapsed time when particles move along the pole between the SNR shock and WTS. We ignored the first term in Eq. (42) because $R_{\rm SNR}<R_{\rm WTS}$. Therefore, $N_{\rm cyc}$ is calculated to be

where $\Delta t_{\rm col}=t_{\rm col}-t_{\rm inj}$. In this work, $t_{\rm col}$ is approximately $1500~{}{\rm yr}$. Here, we assume that $t_{\rm inj}$ is much earlier than $t_{\rm col}$. Thus, $\Delta t_{\rm col}=t_{\rm col}-t_{\rm inj}\approx t_{\rm col}$. Furthermore, $N_{\rm cyc}$ is determined by $u_{\rm SNR}$ because $R_{\rm WTS}\approx u_{\rm SNR}t_{\rm col}$. Then, thanks to the WTS, the attainable energy of particles accelerated in the SNR-WTS system, $\varepsilon_{\rm cyc}$, can be about 8 times the maximum energy of particles accelerated in SNR shock, $\varepsilon_{\rm SNR}$. If particles are injected at a time close to $t_{\rm col}$, these particles cannot experience the cyclic motion until the SNR shock collides with the WTS. $N_{\rm cyc}$ is larger than unity if $t_{\rm inj}<1300~{}{\rm yr}$. Therefore, particles injected at $t_{\rm inj}<1300~{}{\rm yr}$ could experience one or more cycles between the SNR shock and WTS. From Eqs. (30), (LABEL:eq:emax_snr), (37), and (44), $\varepsilon_{\rm cyc}$ is

where $u_{\rm SNR}(t_{\rm inj})\Delta t_{\rm col}\approx u_{\rm SNR}(t_{\rm inj})t_{\rm col}\approx R_{\rm WTS}$ because $t_{\rm inj}\ll t_{\rm col}$ is assumed. Deceleration of the SNR shock velocity can be negligible under the parameters we use in our simulations. This is because the mass of the circumstellar matter that the SNR shock sweeps up until the SNR shock collides with the WTS is $0.1-1M_{\odot}$, which is much smaller than the ejecta mass, $M_{\rm ej}\sim 10M_{\odot}$. The SNR shock velocity at the i-th cycle, $u_{\rm SNR}(t_{\rm i})$, is almost same as the SNR shock velocity at $t_{\rm inj}$, $u_{\rm SNR}(t_{\rm inj})$. Hence, $\sum_{\rm i=1}^{N_{\rm cyc}}u_{\rm SNR}(t_{\rm i})\approx N_{\rm cyc}u_{\rm SNR}(t_{\rm inj})$. As one can see Eq. (50), $\varepsilon_{\rm cyc}$ depends on not SNR parameters but WR star parameters. From Eqs. (21) and (24), $\varepsilon_{\rm cyc}$ is

In Eq. (50), the attainable energy becomes large when the magnetic field strength in the free wind region of WR stars is strong ($B_{{\rm w},\phi}\propto B_{\rm A}R_{\rm A}^{2}\Omega_{*}/V_{\rm w}$). However, the residence time in the free wind region (upstream region) becomes comparable to that in the downstream region when the magnetic field strength in the free wind region is strong. $v_{\theta,\rm SNR}$ and $v_{\theta,\rm WTS}$ can be smaller than $c/2$ when the downstream residence time is not negligible. This is because the particle motion in the downstream region is assumed to be diffusion and the downstream diffusion velocity is much smaller than the drift velocity in the free wind region. If the downstream residence time is larger than the residence time in the free wind region, particles are almost resident in the downstream region. Downstream diffusing particles are hard to spread along the $\theta$ direction compared with the case that particles are almost resident in the free wind region. Therefore, $v_{\theta,\rm SNR}$ and $v_{\theta,\rm WTS}$ can be smaller than $c/2$. In our acceleration model, the downstream residence time can be larger than the residence time in the free wind region when the downstream magnetic field strength is not much larger than the magnetic field strength in the free wind region (kamijima20, ). $v_{\theta,\rm WTS}$ rather than $v_{\theta,\rm SNR}$ can be the bottleneck for the cyclic motion between the SNR shock and WTS. From Eq. (59), $v_{\theta,\rm WTS}$ is

When $B_{\rm d,WTS}/B_{{\rm w},\phi}(R_{\rm WTS})\geq 16c/(3\pi V_{\rm w})$, the residence time in the downstream region of the WTS is smaller than the residence time in the free wind region. Thus, $v_{\theta,\rm WTS}$ approximately becomes $c/2$. This condition, $B_{\rm d,WTS}/B_{\rm w}(R_{\rm WTS})\geq 16c/(3\pi V_{\rm w})$, can be rewritten as the condition of the rotation period of progenitors, $P_{*}$, as follows:

where $P_{*,\rm cr}$ is the critical rotation period, which is given by the condition that the residence time in the downstream region of the WTS is equal to the residence time in the free wind region. The shorter rotation period leads to the larger toroidal component of the Parker-spiral magnetic field in the free wind region ($B_{{\rm w},\phi}\propto P_{*}^{-1}$). However, the downstream magnetic field strength in our model does not depend on $P_{*}$. Then, the shorter $P_{*}$ leads to the smaller $B_{\rm d,WTS}/B_{{\rm w},\phi}(R_{\rm WTS})$. The larger $B_{\rm d,WTS}/B_{{\rm w},\phi}(R_{\rm WTS})$ is favorable for the larger $v_{\theta,\rm WTS}$, which means that the longer $P_{*}$ is favorable. For the case of the longer $P_{*}$, the energy gain in the SNR shock becomes small due to the smaller magnetic field strength in the free wind region although $v_{\theta,\rm WTS}$ approaches $c/2$. The attainable energy, $\varepsilon_{\rm cyc}$, becomes the maximum value, $\varepsilon_{\rm cyc,max}$, when $P_{*}=P_{\rm*,cr}$. From Eqs. (50), (55), and $P_{*}=2\pi/\Omega_{*}$, $\varepsilon_{\rm cyc,max}$ is

For protons ($Z=1$), unrealistic mass-loss rate and wind velocity are required to accelerate protons to the PeV scale without the upstream magnetic field amplification.

This work focuses on the WTS of WR stars. RSGs are also expected to create the WTS similar to WR stars. However, the cyclic motion between the SNR shock and WTS of RSGs could be hard to occur compared with the WTS of WR stars. This is because, in the explosion center rest frame, the velocity of the WTS of RSGs is almost the same as the wind velocity of RSGs (dwarkadas05, ). The upstream kinetic energy flux measured in the WTS rest frame is much smaller than that for WR stars. This leads to a quite small magnetic field in the downstream region of the WTS of RSGs. Then, the residence time in the downstream region of the WTS is much larger than that in the free wind region. Therefore, the WTS of RSGs could not play an important role in the cyclic motion between the SNR shock and WTS.

First, we show the simulation results for the case of aligned rotators. The top four panels in Fig. 3 are the time evolution of the particle distribution for the case of $\alpha_{\rm inc}=0$. The vertical axis is the $z$ component of the spatial coordinate. The $z$ direction is parallel to the rotation axis of progenitors. The horizontal axis shows the distance from the rotation axis of progenitors, $\sqrt{x^{2}+y^{2}}$. Both the horizontal and vertical axes are normalized by the WTS radius, $R_{\rm WTS}$. The inner and outer black semicircles are the SNR shock and WTS, respectively. The points and their color show the positions of simulation particles and their particle energies, respectively. The gray line at the equator ($z=0$) is the current sheet. Particles injected on the SNR shock surface are accelerated at the SNR shock while moving to the equator, reaching the equator eventually. Once particles interact with the current sheet, particles escape from the SNR shock while moving along the current sheet (meandering motion). These acceleration and escape processes are shown in our previous work (kamijima22, ). Particles escaped from the SNR shock move towards the WTS ($t=120.7~{}{\rm yr}$ in Fig. 3). Particles that reach the WTS are accelerated at the WTS while moving to the poles ($t=193.2~{}{\rm yr}$ in Fig. 3). The toroidal component of the Parker-spiral magnetic field, $B_{{\rm w},\phi}$, becomes small towards the poles. $B_{{\rm w},\phi}$ around the poles is smaller than the radial component of the Parker-spiral magnetic field, $B_{{\rm w},r}$. Then, the shock around the poles is the parallel shock. Particles around the poles move along $B_{{\rm w},r}$ to the SNR shock ($t=227.0~{}{\rm yr}$ in Fig. 3). These returned particles are accelerated at the SNR shock again while moving to the equator. These particles eventually interact with the current sheet at the equator and escape from the SNR shock again ($t=273.9~{}{\rm yr}$ in Fig. 3). This cycle process between the SNR shock and WTS lasts until the SNR shock collides with the WTS ($t=1449.4~{}{\rm yr}$) and accelerates particles between the SNR shock and WTS.

The top two panels in Fig. 4 show the energy spectra of all particles for the case of $\alpha_{\rm inc}=0$. The horizontal and vertical axes are the particle energy, $\varepsilon$, and $\varepsilon^{2}dN/d\varepsilon$, respectively. The vertical red and blue lines are the energy limited by the potential difference between the pole and equator of the SNR shock (the energy gain at the SNR shock) and energy limited by the cyclic motion between the SNR shock and WTS. In the top two panels, the value of the vertical red and blue lines are given by substituting $\alpha_{\rm inc}=0$ into Eq. (LABEL:eq:emax_snr) and Eq. (50). The top left panel shows the energy spectrum at the time when particles escaped from the SNR shock still do not reach the WTS ($t=120.7~{}{\rm yr}$). The top right panel shows the energy spectrum at the time when the SNR shock collides with the WTS ($t=1449.4~{}{\rm yr}$). As shown in Ref. kamijima22 , the maximum energy of particles escaped from the SNR shock is limited by the potential difference between the pole and equator of the SNR shock ($t=120.7~{}{\rm yr}$). The reason why the cutoff energy of the energy spectrum at $t=120.7~{}{\rm yr}$ is slightly smaller than the theoretical estimate of the energy gain at the SNR shock (red line) is because the number of particles injected around the poles is small due to the small solid angle around the poles although particles injected around the poles are considered in the theoretical estimate of the energy gain at the SNR shock. Escaped particles are accelerated at the WTS while moving to the poles and the particle energy can exceed the energy gain at the SNR shock (red line). Particles around the pole return to the SNR shock while moving along $B_{{\rm w},r}$ and are accelerated at the SNR shock again. Thanks to this cyclic motion between the SNR shock and WTS, the maximum energy continues to increase until the SNR shock collides with the WTS ($t=1449.4~{}{\rm yr}$). The cutoff energy of the energy spectrum at $t=1449.4~{}{\rm yr}$ is almost in good agreement with the theoretical estimate of the maximum energy limited by the cyclic motion (blue line) within a factor of two. The energy spectrum of accelerated particles is the same as the energy spectrum of the standard DSA prediction ($dN/d\varepsilon\propto\varepsilon^{-2}$) because particles are accelerated at the SNR shock and WTS by the DSA in our simulations. Many of injected particles are advected to the far downstream of the SNR shock and lose their energies by the adiabatic cooling. Thus, there are particles with energy below the injected particle energy in Fig. 4.

Next, we show the results for the case of $\alpha_{\rm inc}=\pi$. The bottom four panels in Fig. 3 show the time evolution of the particle distribution for the case of $\alpha_{\rm inc}=\pi$. The bottom two panels in Fig. 4 show the energy spectra of all particles for the case of $\alpha_{\rm inc}=\pi$. The sign of the magnetic field in the free wind region is opposite to that in the case of $\alpha_{\rm inc}=0$. Hence, the direction of the particle motion between the SNR shock and WTS is opposite to that in the case of $\alpha_{\rm inc}=0$. The results for $\alpha_{\rm inc}=\pi$ are the same as the results for $\alpha_{\rm inc}=0$ except for the direction of the particle motion.

The schematic picture of the oblique rotator case ($\alpha_{\rm inc}\leq\pi/2$) is shown in Fig. 5. The differences from aligned rotators are that the magnetic axis of progenitors is tilted by the angle, $\alpha_{\rm inc}$, from the rotation axis of progenitors and the current sheet has a wavy structure.

Next, we show the results for $\alpha_{\rm inc}=\pi/6$. The top four panels in Fig. 6 are the time evolution of the particle distribution for $\alpha_{\rm inc}=\pi/6$. The format of Fig. 6 is the same as that of Fig. 3. The gray region ($\pi/3\leq\theta\leq 2\pi/3$) is the wavy current sheet region. Particles accelerated at the SNR shock escape from the SNR shock ($t=120.7~{}{\rm yr}$ in Fig. 6). Particles can escape from the SNR shock even though the current sheet is wavy. Almost similar to the aligned rotator case, the cyclic motion between the SNR shock and WTS occurs. However, contract to the aligned rotator case, particles feel the magnetic field averaged over the reversal magnetic field structure inside the wavy current sheet region, which is smaller than the magnetic field for the aligned rotator case [see Eq. (LABEL:eq:mean_bphi), and Eq. (17) and Fig. (10) in Ref. (kamijima22, )].

The top two panels in Fig. 7 show the energy spectra of all particles for $\alpha_{\rm inc}=\pi/6$. Particles inside the wavy current sheet region can spread around the equator due to the weak mean magnetic field, which prevents the ideal cyclic motion. Therefore, the deviation between the theoretical estimate of the cycle-limited maximum energy (blue line) and the cutoff energy of the energy spectrum for $\alpha_{\rm inc}=\pi/6$ is slightly larger than that for $\alpha_{\rm inc}=0$. This deviation at $t=1449.4~{}{\rm yr}$ is within a factor of three. As with the case of $\alpha_{\rm inc}=0$, the energy spectrum of particles accelerated by the cyclic motion is the same as the standard DSA prediction ($dN/d\varepsilon\propto\varepsilon^{-2}$) because particles are accelerated at both the SNR shock and WTS by the DSA.

Next, we show the results for the case of $\alpha_{\rm inc}=5\pi/6$. The bottom four panels in Fig. 6 show the time evolution of the particle distribution for the case of $\alpha_{\rm inc}=5\pi/6$. The bottom two panels in Fig. 7 show the energy spectra of all particles for the case of $\alpha_{\rm inc}=5\pi/6$. The sign of the magnetic field in the free wind region is opposite to that in the case of $\alpha_{\rm inc}=\pi/6$. Hence, the direction of the particle motion between the SNR shock and WTS is opposite to that in the case of $\alpha_{\rm inc}=\pi/6$. The results for $\alpha_{\rm inc}=5\pi/6$ are the same as the results for $\alpha_{\rm inc}=\pi/6$ except for the direction of the particle motion.