Pulsar Wind Nebulae inside Supernova Remnants as Cosmic-Ray PeVatrons
aa r X i v : . [ a s t r o - ph . H E ] J un MNRAS , 1–6 (0000) Preprint 8 June 2018 Compiled using MNRAS L A TEX style file v3.0
Pulsar Wind Nebulae inside Supernova Remnants asCosmic-Ray PeVatrons
Yutaka Ohira , , ⋆ Shota Kisaka , and Ryo Yamazaki Department of Earth and Planetary Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan Department of Physics and Mathematics, Aoyama Gakuin University, 5-10-1 Fuchinobe, Sagamihara 252-5258, Japan
ABSTRACT
We propose that cosmic-ray PeVatrons are pulsar wind nebulae (PWNe) inside super-nova remnants (SNRs). The PWN initially expands into the freely expanding stellarejecta. Then, the PWN catches up with the shocked region of the SNR, where particlescan be slightly accelerated by the back and forth motion between the PWN and theSNR, and some particles diffuse into the PWN. Afterwards the PWN is compressedby the SNR, where the particles in the PWN are accelerated by the adiabatic com-pression. Using a Monte Carlo simulation, we show that particles accelerated by theSNR to . can be reaccelerated to until the end of the PWN compression. Key words: cosmic rays – acceleration of particles – shock waves – pulsars – ISM:supernova remnants.
The origin of cosmic-ray (CR) PeVatrons is a long standingproblem in the astrophysics. The CR spectrum has a spectralbreak at ∼ eV = PeV (so called, the knee energy). Thediffusive shock acceleration (DSA) at supernova remnants(SNRs) is believed to be the acceleration mechanism of CRsup to the knee energy (Axford et al. 1977; Krymskii 1977;Bell 1978; Blandford & Ostriker 1978). Although recentgamma-ray observations support the idea (Ohira et al. 2011;Ackermann et al. 2013), there are still many problems. Oneis the knee problem. It was estimated that SNRs cannot ac-celerate CRs to the knee energy for a parallel shock withoutstrong magnetic field (Lagage & Cesarsky 1983). In order toaccelerate CRs to the knee energy, magnetic fields must beamplified in the shock upstream region. Several mechanismsof the magnetic field amplification in the shock upstreamregion have been proposed (Bell 2004; Malkov et al. 2010;Ohira & Takahara 2010; Ohira 2012). However, no simula-tions have demonstrated that the upstream magnetic field issufficiently amplified to accelerate CRs to the knee energy.In contrast to the shock upstream region, magnetic fieldsare expected to be easily amplified to the equipartition levelin the shock downstream region (Giacalone & Jokipii 2007;Inoue et al. 2009; Guo et al. 2012; Caprioli & Spitkovsky2013; Ohira 2016a). Super-Alfv´enic turbulence amplifies thefield by stretching the field line. The downstream turbulenceis generated by interactions between upstream density fluc-tuations and the shock front. In addition, the downstreamturbulence is generated by the Rayleigh-Taylor instability ⋆ E-mail: [email protected] (YO) at the contact discontinuity. As a result, the magnetic fieldin the shocked region is amplified by the turbulence. If onlythe downstream magnetic field is amplified, the accelera-tion time scale of DSA is predominantly determined by theupstream residence time of accelerated particles, which de-pends on the shock velocity, u sh , and the diffusion coefficient, D , (Ohira & Yamazaki 2017; Drury 1983). Then the accel-eration time scale is given by t acc ≈ Du = cE eu B up ≈ yr (cid:18) E (cid:19) (cid:18) u sh × km s − (cid:19) − (cid:18) B up µ G (cid:19) − , (1)where we assume the shock compression ratio of 4, the Bohmdiffusion coefficient, D B = cE / eB up , and c , e , E and B up arethe speed of light, elementary charge, particle energy, andupstream magnetic field strength, respectively. After the freeexpansion phase ( t > t S ), the velocity of the forward shockdecreases with time. The Sedov time, t S , is given by t S ≈ yr (cid:18) E SN erg (cid:19) − (cid:18) M ej ⊙ (cid:19) (cid:18) n . − (cid:19) − , (2)where E SN , M ej and n are the explosion energy, ejectamass, and the ambient number density, respectively(McKee & Truelove 1995). From the condition, t acc = t S ,the maximum energy of particles accelerated at the forwardshock is given by E max ≈ . (cid:18) E SN erg (cid:19) (cid:18) M ej ⊙ (cid:19) − (cid:18) n . − (cid:19) − (cid:18) B up µ G (cid:19) , (3) © Y. Ohira et al. where u sh = ( E SN / M ej ) / is assumed(McKee & Truelove 1995). This is about 10 times smallerthan the knee energy. The maximum energy weakly dependson parameters of supernovae and their environment exceptfor the upstream magnetic field. Therefore, in oder forSNRs to accelerate CRs to the knee energy, the upstreammagnetic field needs to be amplified to about µ G.The other possible solution for the knee problem is DSAat the perpendicular shocks (Jokipii 1987). Since acceleratedparticles cannot propagate to the far upstream region, theacceleration time scale becomes small for the perpendicu-lar shock. Although the injection to DSA was thought tobe difficult for the perpendicular shock, it was shown bythree-dimensional hybrid simulations that particles are in-jected to DSA at the perpendicular shock in a partially ion-ized plasma, so that particles are rapidly accelerated there(Ohira 2016b). However, DSA at the perpendicular shockshas another problem. For such cases, the maximum energyis limited by the size of acceleration region, R . The avail-able potential drop is ∆ φ = RB up u sh / c , so that the maximumenergy of accelerated protons is given by E max = e ∆ φ ≈ . (cid:18) R
10 pc (cid:19) (cid:18) B up µ G (cid:19) (cid:18) u sh × km s − (cid:19) , (4)which is again 10 times smaller than the knee energy. Since u sh = ( E SN / M ej ) / and R = u sh t S , equation (4) becomesidentical to equation (3). In order to accelerate CRs to theknee energy at the perpendicular shock, we need an excep-tional condition (Takamoto & Kirk 2015).In this paper, we propose a reacceleration mechanismfrom . PeV to PeV by pulsar wind nebulae (PWNe)inside SNRs. Recent observations of young γ -ray pulsarssuggest that most core-collapse supernovae generate pul-sars and their spindown luminosity is typically L sd ∼ × erg s − (Watters & Romani 2011). As mentioned above,the magnetic field in the shocked region of SNRs is strongenough to scatter high-energy particles. The magnetic fieldin young PWNe is also strong compared with that inthe interstellar medium, which is about B PWN ∼ µ G (Tanaka & Takahara 2010; Torres et al. 2014). The PWNinitially expands into the freely expanding stellar ejecta to-ward the shocked region of the SNR. Since this system canbe interpreted as two walls approaching each other, parti-cles are accelerated, shuttling between the PWN and theshocked region of the SNR. After the PWN reaches the re-verse shock of the SNR, the PWN is compressed and parti-cles inside the PWN are accelerated by the adiabatic com-pression (Blondin et al. 2001; van der Swaluw et al. 2001).In the next section, using Monte Carlo simulation, we showthat the PWN-SNR system actually accelerates particlesfrom . PeV to PeV.
In order to investigate the particle acceleration by the PWN-SNR system, we first provide evolution of an SNR and aPWN inside the SNR. As a first step, we consider a spheri-cally symmetric structure. For constant ejecta and ambientdensity profiles, the approximate time evolution of the for-ward and reverse shock radii, R SNR , fs and R SNR , rs , are given by McKee & Truelove (1995). They got R SNR , fs R S = . tt S ( + . (cid:18) tt S (cid:19) / ) − / , (5) R SNR , rs R S = . tt S ( + . (cid:18) tt S (cid:19) / ) − / , (6)for the free expansion phase ( t < t S ) , and R SNR , fs R S = (cid:18) . tt S − . (cid:19) / , (7) R SNR , rs R S = tt S (cid:26) . − . tt S − .
37 ln (cid:18) tt S (cid:19)(cid:27) , (8)for the Sedov phase ( t ≥ t S ) , where R S is given by(McKee & Truelove 1995), R S = . t S (cid:18) E SN M ej (cid:19) / ≈ (cid:18) M ej M ⊙ (cid:19) / (cid:18) n . − (cid:19) − / . (9)For a constant spindown luminosity of a pulsar andan uniform ejecta profile, the analytical solution for thetime evolution of the PWN radius, R PWN , is given by(van der Swaluw et al. 2001), R PWN R S = . (cid:18) L sd t S E SN (cid:19) / (cid:18) tt S (cid:19) / , (10)where the solution can be applied until the PWN interactswith the reverse shock of the SNR. It should be noted thatthe PWN radius, R PWN , does not significantly depend on thespindown luminosity, L sd .Fig. 1 shows the time evolutions of the forward andreverse shock radii of the SNR and the radius of PWN,where we assume the uniform ejecta profile with the ejectamass of M ej = ⊙ , the explosion energy of E SN = erg,the uniform ambient matter profile with the density of n = . − , and the constant pulsar spindown luminosity of L sd = × erg s − . For these parameters, the PWN catchesup with the reverse shock of the SNR at t c ≈ × yr. Af-terwards, the PWN is compressed by the larger pressure ofthe shocked region of the SNR. In this paper, the constantspindown luminosity is assumed for t ≤ t c . Then we simplyassume that the velocity of the PWN during the compres-sion is v PWN = − v PWN ( t c )/ and the final size of the PWNis R PWN ( t end ) = R PWN ( t c )/ . Then, the PWN size becomes ≈ . pc at t end ≈ × yr. These assumptions are reasonableto simulate evolution of a spherical PWN (e.g. Gelfand et al.2009).We next give the velocity field in the PWN-SNR sys-tem. Before the PWN reaches the reverse shock of the SNR,the expansion velocity of the shocked ejecta just outside thereverse shock in the observer frame is given by v SNR , shocked = R SNR , rs / t − v SNR , rs + v SNR , rs , (11)where we assume the compression ratio at the SNR reverseshock is , and v SNR , rs = dR SNR , rs / dt is the propagation ve-locity of the reverse shock in the observer frame. In thispaper, the velocity field between the forward and reverse MNRAS000
37 ln (cid:18) tt S (cid:19)(cid:27) , (8)for the Sedov phase ( t ≥ t S ) , where R S is given by(McKee & Truelove 1995), R S = . t S (cid:18) E SN M ej (cid:19) / ≈ (cid:18) M ej M ⊙ (cid:19) / (cid:18) n . − (cid:19) − / . (9)For a constant spindown luminosity of a pulsar andan uniform ejecta profile, the analytical solution for thetime evolution of the PWN radius, R PWN , is given by(van der Swaluw et al. 2001), R PWN R S = . (cid:18) L sd t S E SN (cid:19) / (cid:18) tt S (cid:19) / , (10)where the solution can be applied until the PWN interactswith the reverse shock of the SNR. It should be noted thatthe PWN radius, R PWN , does not significantly depend on thespindown luminosity, L sd .Fig. 1 shows the time evolutions of the forward andreverse shock radii of the SNR and the radius of PWN,where we assume the uniform ejecta profile with the ejectamass of M ej = ⊙ , the explosion energy of E SN = erg,the uniform ambient matter profile with the density of n = . − , and the constant pulsar spindown luminosity of L sd = × erg s − . For these parameters, the PWN catchesup with the reverse shock of the SNR at t c ≈ × yr. Af-terwards, the PWN is compressed by the larger pressure ofthe shocked region of the SNR. In this paper, the constantspindown luminosity is assumed for t ≤ t c . Then we simplyassume that the velocity of the PWN during the compres-sion is v PWN = − v PWN ( t c )/ and the final size of the PWNis R PWN ( t end ) = R PWN ( t c )/ . Then, the PWN size becomes ≈ . pc at t end ≈ × yr. These assumptions are reasonableto simulate evolution of a spherical PWN (e.g. Gelfand et al.2009).We next give the velocity field in the PWN-SNR sys-tem. Before the PWN reaches the reverse shock of the SNR,the expansion velocity of the shocked ejecta just outside thereverse shock in the observer frame is given by v SNR , shocked = R SNR , rs / t − v SNR , rs + v SNR , rs , (11)where we assume the compression ratio at the SNR reverseshock is , and v SNR , rs = dR SNR , rs / dt is the propagation ve-locity of the reverse shock in the observer frame. In thispaper, the velocity field between the forward and reverse MNRAS000 , 1–6 (0000)
WNe inside SNRs as CR PeVatrons R ad i u s ( p c ) Time ( yr )t c t end PWNShocked SNRR
SNR,fs R SNR,rs R PWN
Figure 1.
Time evolution of the PWN radius (red solid) witha spindown luminosity of L sd = × erg s − , and radii of theforward (blue dashed) and the reverse shocks (magenta dotted)of the SNR with an explosion energy of E SN = erg, an ejectamass of M ej = M ⊙ , and an ambient number density of n = . cm − . -3 -2 -1 E d N / d E E ( PeV )
Figure 2.
Energy spectra of reaccelerated particles for Model A.The black dashed and red solid histograms are energy spectra at t = t c and t end , respectively. The initial energy is . PeV. shocks is approximated, for simplicity, as a linear interpola-tion between v SNR , shocked at r = R SNR , rs + and v SNR , fs / at r = R SNR , fs − . The expansion velocity of the PWN is givenby v PWN = dR PWN / dt and the velocity field in the PWN isassumed to be uniform.After the PWN interacts with the reverse shock of theSNR, the velocity between the forward shock and the PWNis approximately given by the linear interpolation between v PWN = − v PWN ( t c )/ and v SNR , fs / . The velocity field in thePWN is assumed as ® v PWN , in ( t , ® r ) = − v PWN ( t c ) ® rR PWN ( t ) . (12)Since the shocked region of the SNR and the PWN re-gion are expected to be highly turbulent (Porth et al. 2014),motion of high-energy particles could be approximated asthe random walk. Using the above hydrodynamical struc-ture, we perform a test-particle Monte Carlo simulation. -3 -2 -1 E d N / d E E ( PeV )
Figure 3.
Energy spectra of reaccelerated particles at the endof the PWN compression, t = t end . The red solid, green dashed,blue dotted, magenta long-dashed, cyan dot-dashed, and yellowdot-long-dashed histograms are for Model A, B, C, D, E, and F.The red histogram in this figure is the same as in Fig. 2. Theparameters for each model are tableted in Table 1.Model M ej n L sd t inj B SNR B PWN ( M ⊙ ) ( cm − ) (erg s − ) (s) ( µ G) ( µ G)A 3 0.1 × ×
300 150
B 10 0.1 × ×
300 150
C 3 1.0 × . ×
300 150
D 3 0.1 × ×
300 150
E 3 0.1 × ×
150 150
F 3 0.1 × ×
300 50
Table 1.
List of parameters for Fig. 3.
Simulation particles are isotropically scattered in the localfluid frame. The scattering time is assumed to be the Bohmscattering, t sc = Ω − ( E / m p c ) , where Ω c ≈ − s − ( B / µ G ) is the proton cyclotron frequency, E is the particle energy,and m p is the proton mass. Once particles escape from theSNR, we do not follow the particles. However, no particlesescape from the SNR in this paper.In this paper, we set the magnetic field to be B SNR = × µ G and B PWN = . × µ G in the shocked region ofthe SNR and PWN, respectively, and B ej = in the freelyexpanding ejecta. Hence, particles are not scattered in thefreely expanding ejecta which disappears after the PWN in-teracts with the reverse shock of the SNR. Since the forwardshock of the SNR can accelerate particles to about . PeV(see introduction), we set the initial energy to be . PeV.The accelerated particles are advected toward the down-stream region of the forward shock. Furthermore, they areexpected to be advected (or diffuse) to the reverse shockedregion because of the Rayleigh-Taylor instability and tur-bulence. In the next section, we will discuss on the amountof particles which the turbulence carries from the forwardshock to the revere shock front. In this paper, we impulsively
MNRAS , 1–6 (0000)
Y. Ohira et al. inject simulation particles isotropically on the reverse shocksphere, r = R SNR , rs at t inj = . × yr instead of solvingthe particle transport from the forward shock to the reverseshock, which should be addressed in the future. About half ofthe injected particles initially diffuse to the reverse shockedregion, and the rest of particles run to the freely expandingejecta.Fig. 2 shows energy spectra of the accelerated particlesfor Model A, parameters of which are listed in Table 1. Theblack dashed histogram shows the energy spectrum at thetime when the PWN reaches the reverse shock of the SNR( t = t c ). They are accelerated up to twice the initial energyby the back and forth motion between the PWN and theshocked region of the SNR. The energy gain in each cycleis ∆ E / E ∼ ∆ v / c , where ∆ v = v SNR , shocked − v PWN is the rel-ative velocity. The time scale in each cycle is ∆ t ∼ ∆ R / c ,where ∆ R = R SNR , rs − R PWN is the relative distance. Then,the acceleration time scale for the reciprocation is given by t acc = ∆ t ( E / ∆ E ) ∼ ∆ R / ∆ v that is the same as the dynamicaltime scale in which the PWN catches up with the SNR re-verse shock. Since the acceleration time scale, t acc , does notsignificantly depend on the magnetic field strength and theparticle energy as long as the scattering time is smaller thanthe time scale of the reciprocation between the SNR and thePWN, ∆ t , it takes t acc to accelerate particles to twice theinitial energy. Therefore, the maximum energy during theapproaching phase becomes twice the initial energy, whichdoes not significantly depend on parameters of the PWN-SNR system.The red solid histogram in Fig. 2 shows the energy spec-trum at t = t end . They are further accelerated to PeV by thePWN compression. The particle energy is increased by a fac-tor R PWN ( t c )/ R PWN ( t end ) = during the compression, so thatthe particles are finally accelerated to ten times the initialenergy. Hence, the maximum energy of particles acceleratedby the PWN-SNR system is given by E max ∼ (cid:18) E inj . (cid:19) (cid:18) R PWN ( t c )/ R PWN ( t end ) (cid:19) , (13)where E inj is the initial energy of particles injected to thePWN-SNR system, that corresponds to the maximum en-ergy of particles accelerated by the SNR shocks. The com-pression factor, R PWN ( t c )/ R PWN ( t end ) , is determined by thepressure balance between the PWN and the shocked regionof the SNR at t end . The rotational energy of a pulsar, E rot , isinitially stored in the PWN, but the PWN eventually losesits energy by synchrotron radiation. The synchrotron cool-ing time is given by t cool ≈ . (cid:18) B µ G (cid:19) − (cid:18) E (cid:19) − , (14)which is smaller than t end . The characteristic energy is typ-ically a few GeV − (Torres et al. 2014). Therefore,most of the rotational energy converts into synchrotron pho-tons during the compression phase of PWNe. Then, the re-maining energy in the PWN is η E rot at t end , where η is theremaining fraction of the order of 0.1 (Gelfand et al. 2009).From the equation E SN / R SNR , fs ( t end ) = η E rot / R PWN ( t end ) ,the radius of the PWN at the end of the compression isgiven by R PWN ( t end ) = . R SNR , fs ( t end ) (cid:18) η E rot / E SN − (cid:19) / , (15) where we set E SN = erg, E rot = erg and η = . .Since R SNR , fs ( t end ) ≈
13 pc and R PWN ( t c ) ≈ in this paper(see Fig. 1), the compression factor, R PWN ( t c )/ R PWN ( t end ) = , is a reasonable approximation. If the initial rotationalenergy of the pulsar is smaller, the PWN is more compressed,so that particles are accelerated to higher energies by thecompression. However, it should be noted that the maximumenergy cannot exceed the limitation by the PWN size. Oncethe gyroradius of accelerated particles becomes comparableto the PWN size, they start to escape from the PWN.In order to explore the parameter dependence of theabove results, we perform other simulations with differentparameters. The parameter sets are tableted in Table 1.As can be seen in Fig. 3, the results do not change sig-nificantly as long as R PWN ( t c )/ R PWN ( t end ) = is fixed. Thisis because dynamics of SNR and PWN do not significantlydepend on the spindown luminosity and the ambient num-ber density (see equations (5)–(10)). Although the ejectamass dependency of the Sedov time ( t S ∝ M / ) is com-paratively strong compared with other parameter depen-dences, the ejecta mass is expected not to be distributedwidely for core collapse supernovae that leave a neutronstar. The acceleration by the PWN compression does notdepend on the magnetic field strength. For t inj < t < t c ,particles diffuse into the PWN. The diffusion length scaleis given by R diff = p D ( t c − t inj ) , where D is the diffusioncoefficient. On the other hand, during the PWN compres-sion phase ( t c < t < t end ), the particles escape from thePWN by diffusion. The escape time scale, t esc is given by t esc = R / D = t c − t inj , which is independent on the diffu-sion coefficient and magnetic field strength, so that the finalspectrum does not depend on the magnetic field strength.Therefore, many PWN-SNR systems can be expected to ac-celerate particles to the knee energy. In the previous section, we assumed that . PeV CRs areinjected at the reverse shock of the SNR to reaccelerate themto PeV. In this section, we discuss some injection mecha-nisms at the reverse shock of the SNR. There are four shocksin the PWN-SNR system before the PWN interacts with thereverse shock of the SNR, the reverse and forward shocks ofthe SNR, the termination shock of the pulsar wind, and theforward shock driven by the PWN. The termination shock ofthe pulsar wind can accelerate electrons and positrons in thestandard picture, but it has not been understood whetherprotons and heavy nuclei are accelerated by the terminationshock or not.The reverse shock of the SNR and the forward shockdriven by the PWN propagate into the freely expandingejecta, where the magnetic field strength in the shock up-stream region is expected to be very weak. In this case, wenaively expect a weak CR acceleration by the shocks propa-gating into the freely expanding ejecta. However, if the mag-netic field is sufficiently amplified by some mechanisms, thereverse shock of the SNR and the forward shock driven bythe PWN can accelerate CRs to . PeV. Then, they arefurther accelerated to PeV CRs by the PWN-SNR systemas shown in the previous section.The forward shock of the SNR can easily accelerate CRs
MNRAS , 1–6 (0000)
WNe inside SNRs as CR PeVatrons to . PeV as estimated in the introduction. They are trans-ported to the reverse shock of the SNR as described in thefollowing. Since the magnetic field in the downstream regionof the SNR is amplified by turbulence, the diffusion coeffi-cient due to the particle diffusion would be in the Bohmlimit. Then, the diffusion length scale is given by R diff , p = p D B t = × cm (cid:18) E . (cid:19) (cid:18) B µ G (cid:19) − (cid:18) t (cid:19) , (16)which is comparable to, but still smaller than the distancebetween the forward and reverse shocks, ∆ R = R SNR , fs − R SNR , rs ∼ a few parsecs. Here, D B is the Bohm diffusioncoefficient of the particle diffusion. To be reaccelerated bythe back and forth motion between the PWN and the re-verse shock of the SNR, particles have to be injected withina distance of l p from the reverse shock of the SNR, wherethe diffusion length scale, l p , is given by l p = D B u = . × cm (cid:18) E . (cid:19) (cid:18) B µ G (cid:19) − (cid:18) u km s − (cid:19) , (17)and u is the downstream flow velocity in the reverse shockrest frame. If the diffusion region, l p , overlaps with R diff , p ( l p + R diff , p > ∆ R ), CRs accelerated at the forward shock can betransported to the reverse shock. Since this condition is notsatisfied as long as the magnetic field is strongly amplified,particles cannot diffuse to the freely expanding ejecta fromthe forward shock by the Bohm diffusion. If the amplifiedmagnetic field sufficiently decays in the shock downstreamregion, most of particles accelerated at the forward shockcan diffuse to the freely expanding ejecta, so that they arereaccelerated by the PWN-SNR system.Even if the amplified magnetic field does not decay suffi-ciently, the particles accelerated at the forward shock wouldbe transported to the reverse shock by the turbulent diffu-sion. There are two types of diffusion in the downstreamregion of the SNR, the particle diffusion and the turbu-lent diffusion. The particle diffusion due to magnetic tur-bulence can move particles to other fluid elements, so thatparticles can move from the downstream region to the up-stream region. On the other hand, in the context of theturbulent diffusion, particles move with a fluid element, sothat particles cannot penetrate the shock front. If there aremagnetic turbulence and large-scale fluid turbulence in thedownstream region, and if the diffusion coefficient of theparticle diffusion is smaller than that of the turbulent dif-fusion, the particle diffusion and the turbulent diffusion co-exist. In a timescale smaller than the eddy turnover time,motion of particles can be described by the particle diffu-sion, but it can be described by the turbulent diffusion ina timescale larger than the eddy turnover time. If the eddysize of turbulence is about ∆ R , particles accelerated at theforward shock can move to the vicinity of the revere shock.In this case, the injection fraction of particles at the re-verse shock is about l p / ∆ R ∼ . . If there is turbulencewith an eddy size of L edd in the vicinity of the reverse shock,particles within a distance of l turb = D turb / u from the re-verse shock can diffuse to the reverse shock by turbulence,where D turb ∼ L edd v SMNR , rs / is the diffusion coefficient due to the turbulence around the reverse shock. If the eddy sizeis L edd ∼ . , the injection fraction becomes l turb / ∆ R ∼ . .Although we have to understand turbulence in the SNR toestimate the injection fraction precisely, about − ofparticles accelerated at the forward shock of the SNR couldbe transported to the reverse shock by the turbulent diffu-sion.The turbulent diffusion around the reverse shocked re-gion can also move particles outward from the reverse shockfront. In fact, particles that stay in the reverse shocked re-gion can move to the forward shock region again by theturbulent diffusion. However, once particles escape into thefreely expanding region and go back to the reverse shockedregion, the particles can go back to the freely expand-ing region again in the timescale of D B / cu . The residencetimescale in the downstream region of accelerating particle, D B / cu ∼ r g / u , is typically much smaller than the eddyturnover time, ∼ L edd / u , where r g is the gyroradius. There-fore, once particles diffuse to the freely expanding ejecta,their turbulent diffusion can be neglected, that is, our MonteCarlo approach in section 2 is valid for particles that havealready injected at the reverse shock front. We first discuss on the energy source of our model. Therequired energy per SNR is about erg in order to supplyGalactic CRs with an energy of GeV. Since the recent CRobservations show that the source spectrum of Galactic CRsshould be dN / dE ∝ E − . , the required energy to supply PeVCRs is about × erg. In our model, the main accelerationis due to the PWN compression by the SNR. The energysource is the work done by the SNR, which is given by pdV ∼ erg (cid:18) E SN erg (cid:19) (cid:18) R PWN ( t c )/ R SNR , fs . (cid:19) , (18)where p ∼ E SN / R , fs and dV ∼ R PWN ( t c ) are the SNR pres-sure and the PWN volume compressed by the SNR. There-fore, the PWN-SNR system has enough energy to supply thePeV CRs.We next discuss on the acceleration of heavy nuclei. CRsare organized not only by protons but also by heavy nucleiwhose origin is also a long standing problem (Ohira et al.2016). Since the supernova ejecta is metal rich, the reverseshock propagating into the supernova ejecta is thought tobe the origin of heavy CR nuclei (Ptuskin et al. 2013). How-ever, the maximum energy of the accelerated particle at thereverse shock is not so large because the magnetic field inthe expanding supernova ejecta is expected to be very small.Furthermore, particles accelerated by the reverse shock suf-fer the adiabatic cooling. Our reacceleration model can boostthe maximum energy of accelerated heavy nuclei, so that thePWN-SNR system could be important for the production ofheavy CR nuclei.Next, we discuss the energy spectrum of acceleratedparticles. In this paper, to investigate whether the PWN-SNR system can accelerate CRs to the PeV scale or not,we considered only the impulsive injection of CRs with E = . PeV at t = yr. In reality, CRs accelerated by theSNR would be continuously injected with an energy spec-trum to the PWN. Furthermore, since there is the potential MNRAS , 1–6 (0000)
Y. Ohira et al. difference of about PV in the PWN, protons could be accel-erated to PeV by drifting the toroidal magnetic field (Bell1992; Bell & Lucek 1996), which was not considered in thispaper. Hence, further studies are needed to understand theenergy spectrum of particles accelerated by the PWN-SNRsystem.In addition, to understand the source spectrum ofGalactic CRs that are injected to our Galaxy, we have toconsider the spectrum of particles that have escaped fromthe SNR (Ohira et al. 2010; Ohira & Ioka 2011), and va-riety of the PWN-SNR system. The PWN could expandagain after the end of the PWN compression. We did notsolve the particle transport and the energy loss during thelater expansion phase in this paper. The time at whichthe PWN completes the re-expansion is about twice of theepoch when the PWN size becomes minimum (Gelfand et al.2009). Therefore, CRs inside the PWN-SNR system may lose − ( / ) / ∼ of their energy due to the re-expansionof the PWN and the expansion of the SNR. When andhow accelerated particles escape from the PWN-SNR sys-tem are important issues. These issues could be addressedby gamma-ray observations of the PWN-SNR system likeSNR G327.1-1.1, W44 and so on.In this paper we assumed as a first step, spherical sym-metry for the PWN-SNR system and the Bohm scatter-ing with constant magnetic field strength for the randomwalk. In reality, the magnetic field strength is not con-stant, a pulsar has a kick velocity, and the supernova ejecta,the ambient matter, the PWN have asymmetry. In addi-tion, the Rayleigh-Taylor instability amplifies the asymme-try (van der Swaluw et al. 2004), so that the PWN-SNR sys-tem is actually more complicated. In particular, the strongturbulence could play important roles, that amplifies themagnetic field, affecting the particle motion (Porth et al.2016), and accelerates particles by turbulent acceleration(Ohira 2013). The turbulence eventually decays, so that ac-celerated particles can escape from the PWN. On the otherhand, the magnetic field is compressed by the PWN com-pression. Then, the escape time scale due to diffusion be-comes long, so that more particles are accelerated by thePWN compression. In order to address above problems, weneed a more realistic magnetohydrodynamical simulation.This will be addressed in future work. We have proposed that PWNe inside SNRs are the CRPeVatron. Firstly, the SNR shock accelerates protons to ∼ . PeV. Then, the protons diffuse into the interior of theSNR and are reaccelerated to ∼ . PeV by the back andforth motion between the SNR and the PWN. Finally, theprotons diffuse into the PWN and are accelerated to ∼ PeVby the adiabatic compression while the PWN is compressedby the SNR. Our model predicts that there must be somestructures in the spectrum of CR protons around . PeV. Inaddition, we have argued that the PWN-SNR system couldbe the origin of heavy CR nuclei.
ACKNOWLEDGEMENTS
Numerical computations were carried out on the XC30 sys-tem at the Center for Computational Astrophysics (CfCA)of the National Astronomical Observatory of Japan. Thiswork was supported by JSPS KAKENHI Grant NumberJP16K17702(YO), JP16J06773(SK), JP15K05088(RY), andJP18H01232(RY). Y.O. is supported by MEXT/JSPS Lead-ing Initiative for Excellent Young Researchers.
REFERENCES
Ackermann M., et al., 2013, Science, 339, 807Axford W. I., Leer E., Skadron G., 1977, International CosmicRay Conference, 11, 132Bell A. R., 1978, MNRAS, 182, 147Bell A. R., 1992, MNRAS, 257, 493Bell A. R., 2004, MNRAS, 353, 550Bell A. R., Lucek S. G., 1996, MNRAS, 283, 1083Blandford R. D., Ostriker J. P., 1978, ApJ, 221, L29Blondin J. M., Chevalier R. A., Frierson D. M., 2001, ApJ,563, 806Caprioli D., Spitkovsky A., 2013, ApJ, 765, L20Drury L. O., 1983, Reports on Progress in Physics, 46, 973Gelfand J. D., Slane P. O., Zhang W., 2009, ApJ, 703, 2051Giacalone J., Jokipii J. R., 2007, ApJ, 663, L41Guo F., Li S., Li H., Giacalone J., Jokipii J. R., Li D., 2012, ApJ,747, 98Inoue T., Yamazaki R., Inutsuka S.-i., 2009, ApJ, 695, 825Jokipii J. R., 1987, ApJ, 313, 842Krymskii G. F., 1977, Akademiia Nauk SSSR Doklady, 234, 1306Lagage P. O., Cesarsky C. J., 1983, A&A, 125, 249Malkov M. A., Diamond P. H., Sagdeev R. Z., 2010,Plasma Physics and Controlled Fusion, 52, 124006McKee C. F., Truelove J. K., 1995, Phys. Rep., 256, 157Ohira Y., 2012, ApJ, 758, 97Ohira Y., 2013, ApJ, 767, L16Ohira Y., 2016a, ApJ, 817, 137Ohira Y., 2016b, ApJ, 827, 36Ohira Y., Ioka K., 2011, ApJ, 729, L13Ohira Y., Takahara F., 2010, ApJ, 721, L43Ohira Y., Yamazaki R., 2017, JHEAp, 13, 17Ohira Y., Murase K., Yamazaki R., 2010, A&A, 513, A17Ohira Y., Murase K., Yamazaki R., 2011, MNRAS, 410, 1577Ohira Y., Kawanaka N., Ioka K., 2016, Phys. Rev. D, 93, 083001Porth O., Komissarov S. S., Keppens R., 2014, MNRAS, 438, 278Porth O., Vorster M. J., Lyutikov M., Engelbrecht N. E., 2016,MNRAS, 460, 4135Ptuskin V., Zirakashvili V., Seo E.-S., 2013, ApJ, 763, 47Takamoto M., Kirk J. G., 2015, ApJ, 809, 29Tanaka S. J., Takahara F., 2010, ApJ, 715, 1248Torres D. F., Cillis A., Mart´ın J., de O˜na Wilhelmi E., 2014,Journal of High Energy Astrophysics, 1, 31Watters K. P., Romani R. W., 2011, ApJ, 727, 123van der Swaluw E., Achterberg A., Gallant Y. A., T´oth G., 2001,A&A, 380, 309van der Swaluw E., Downes T. P., Keegan R., 2004, A&A,420, 937This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000