XXV European Cosmic Ray Symposium, Turin, Sept. 4-9 2016 1 SNRs W28 and W44: old cosmic ray accelerators in molecular clouds V.N.Zirakashvili, V.S.Ptuskin Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radiowave Propagation, 108840 Moscow Troitsk, Russia Nonlinearmodelofdiffusiveshockaccelerationisusedforinvestigationoftheparticleacceleration in old supernova remnants W28 and W44. We modeled the hydrodynamical evolution of the rem- nants, shockmodification and streaming instability produced byaccelerated particles. Comparison with available radio and gamma-ray data is given. 7 I. INTRODUCTION 1 0 ∂u ∂u 1 ∂Pg ∂Pc ∂Pm = u + + (2) 2 The diffusive shock acceleration (DSA) process [1– ∂t − ∂r − ρ(cid:18) ∂r ∂r ∂r (cid:19) 4]isconsideredastheprincipalmechanismforproduc- n a tion of galactic cosmic rays (CR) in supernova rem- 2 1 ∂P ∂P γ P ∂r u J nants (SNRs). During the last decade the excellent g g g g +u + = 5 results of X-ray and gamma-ray astronomy supplied γg 1(cid:18) ∂t ∂r r2 ∂r (cid:19) − 1 the observational evidence of the presence of multi- P ∂P 2 m c TeV energetic particles in these objects (see e.g. [5]). Λ(Te)n +Hc+2Γn ξAVAr(1 hm) (3) − γ 1 − − ∂r ] m− E MostofexistingDSAmodelswereappliedtoyoung H SNRs(seehowever[6]). Thisisprobablybecauseitis 2 ∂P ∂P P ∂r (γ u+ξ V ) h. ethxepreec.teHdotwheavteCrlRowsweritehnehriggyhepsatretnicelregsiaesrearperopdruodceudceind ∂tm +(u+ξAVAr) ∂rm + rm2 m ∂r A Ar = p old SNRs either. So the investigationof CR accelera- ∂Pc h (γ 1)ξ V 2Γ P (4) - m m A Ar n m o tioninoldSNRsisimportantforcalculationofoverall − − ∂r − r CR spectra produced by SNRs. t s In this paper we describe the modifications of our ∂N 1 ∂ ∂N ∂N ∂N p ∂r2w a non-linearDSAmodel[7]designedforinvestigationof = r2D(p,r,t) w + [ ∂t r2∂r ∂r − ∂r ∂p 3r2 ∂r DSAinoldSNRs. WeapplyitfortwoGeVbrightold 1 ∂ 3 SNRs W28 and W44. + p2b(p)N + v p2∂p 4 η δ(p p ) 4 f − f ρ(R +0,t)(R˙ u(R +0,t))δ(r R (t)) 8 II. NONLINEAR DIFFUSIVE SHOCK 4πp2fm f f − f − f 0 0 ACCELERATION MODEL +ηbδ(p−pb)ρ(R 0,t)(u(R 0,t) R˙ )δ(r R (t))(5) 4πp2m b− b− − b − b . b 1 0 DetailsofourmodelofnonlinearDSAcanbefound Here P = 4π dpp3vN/3 is the cosmic ray pressure, c 7 in [7]. The model contains coupled spherically sym- w(r,t) is the aRdvection velocity of cosmic rays, Te, γg 1 metric hydrodynamic equations and the transport and n are the gas temperature, adiabatic index and : equations for energetic protons, ions and electrons. v numberdensity respectively,γ isthe waveadiabatic m i The forward and reverse shocks are included in the index, D(r,t,p) is the cosmic raydiffusion coefficient. X consideration. The radiative cooling of gas is described by the cool- r a Damping of magnetohydrodynamic (MHD) waves ing function Λ(Te). The function b(p) describes the due to the presence of neutral atoms is important for energy losses of particles. In particular the Coulomb oldSNRs. To takethis effect intoaccountwe addthe losses of sub GeV ions and the radiative cooling are equation that describes the transport and generation important in old SNRs. The energy of sub GeV ions of MHD waves (see Eq.(4) below). goes to the gas heating described by the term H in c The hydrodynamical equations for the gas den- Eq. (3). sity ρ(r,t), gas velocity u(r,t), gas pressure P (r,t), Cosmicraydiffusionisdeterminedbyparticlescat- g wave pressure P (r,t), and the equationfor isotropic tering on magnetic inhomogeneities. The cosmic ray m partofthecosmicrayprotonmomentumdistribution streaming instability increases the level of MHD tur- N(r,t,p)inthesphericallysymmetricalcasearegiven bulenceintheshockvicinity[2]andevensignificantly by amplifiestheabsolutevalueofmagneticfieldinyoung SNRs [8, 9]. It decreases the diffusion coefficient and increases the maximum energy of accelerated parti- ∂ρ 1 ∂ cles. The results of continuing theoretical study of 2 ∂t =−r2∂rr uρ (1) this effect can be found in review papers [10, 11]. eConf C16-09-04.3 2 XXV European Cosmic Ray Symposium, Turin, Sept. 4-9 2016 Cosmicrayparticlesarescatteredbymovingwaves TABLE I: Physical parameters of SNRsW28 and W44 anditiswhythecosmicrayadvectionvelocityw may differ from the gas velocity u by the value of the ra- dial component of the Alfv´en velocity V = V /√3 d Rf ESN Mej nH nn Kep T Vf Bf calculated in the isotropic random maAgrneticAfield: kpc pc 1051erg M⊙ cm−3 cm−3 kyr km/s µG w = u+ξ V . The factor ξ describes the possi- W28 1.9 13.4 1.3 6.8 4.0 0.2 0.008 37 121 79 A Ar A ble deviationof the cosmic raydrift velocityfrom the W44 2.8 12.4 1.6 7.1 6.0 0.3 0.006 33 130 102 gas velocity. We use values ξ = 1 and ξ = 1 up- A A − streamoftheforwardandreverseshocksrespectively, where Alfv´en waves are generated by the cosmic ray streaminginstabilityandpropagateinthecorrespond- The injection efficiency is taken to be independent ing directions. oftime η =0.001,andthe particleinjectionmomen- f The pressure of generated waves Pm determines tumispf =2m(R˙f u(Rf+0,t)). Protonsofmassm − thescatteringanddiffusionofenergeticparticleswith are injected at the forward shock. The high injection charge q, momentum p and speed v efficiency results in the significant shock modification alreadyat earlystage ofSNR expansion. Since inthe B2 cpv old remnants the reverse shock absent we put ηb = 0 D =D , D = , B = B2+8πP (6) in the modeling below. B B 0 m 8πPm 3qB q We neglect the pressure of energetic electrons and where B is the total magnetic field strength, while treatthemastestparticles. Theevolutionoftheelec- B0 is the strength of the mean field. At high wave tron distribution is described by equation analogous amplitudes the diffusioncoefficient coincides with the to Eq. (5) with function b(p) describing Coulomb, Bohm diffusion coefficient D . synchrotronand inverse Compton(IC) losses andad- B The parameter h in Eqs. (3,4) describes the frac- ditionalterms describingthe productionof secondary m tionofmagnetic energyproducedby streaminginsta- leptons by energetic protons and nuclei. bility. We use the following dependence h (B) m III. MODELING OF DIFFUSIVE SHOCK B B ACCELERATION IN THE OLD SNRS h =1, <3; h =0.5, >3. (7) m m B0 B0 At high amplitudes the waves are damped and the A significant part of core collapse supernova explo- fraction 1 h of energy goes into the gas heating sion occurs in molecular gas. The stars with initial m − upstream of the shocks [12] that is described by the massesbelow 12 M⊙ haveno power stellarwinds and lastterminEq. (3). Theheatingandwavegeneration thereforedonotproduceastrongmodificationoftheir limits the total compression ratio of cosmic ray mod- circumstellarmedium. Themolecularcloudhavebeen ified shocks. Inthe downstreamregionofthe forward totally destroyed by stellar winds and supernova ex- and reverse shock at R <r <R we put ξ =0 and plosionsofmoremassivestarsatthe instantofexplo- b f A therefore w=u. sion. As a result the star explodes in the inter-clump Intheshocktransitionregionthemagneticpressure mediumwiththedensity5 25cm−3 [13]. Manysuch − is increased by a factor of σγm, where σ is the shock SNRs are observed in gamma rays now. compression ratio. Its impact on the shock dynamics Bright in GeV gamma rays SNRs W44 and W28 is taken into account via the Hugoniot conditions. at distance d 2 3 kpc from the Earth are at the ∼ − Below we use the adiabatic index of Alfv´en waves radiativephasenowandshowsignsofinteractionwith γ = 3/2. For this value of the adiabatic index, the molecular gas [14]. It is believed that circumstellar m wave pressure P = (δB)2/8π equals to the wave mediumisalmostfullyionizedbyultravioletradiation m magnetic energy density. from the remnant interior at the radiative stage [13]. The rate of the neutral damping Γ = 0.5σ v n The same is true for young SNRs because the gas is n ex T n depends on the number density of neutrals n , ther- ionized by the radiationfrom the shock breakoutand n malspeedv 10kms−1 andchargeexchangecross- in the hot shock precursor produced by accelerated T section σ ∼10−14 cm2. We take the damping into particles. Probably there exist an intermediate phase ex ∼ account in the upstream region only and put Γ = 0 when the shock propagates in the neutral medium. n downstream of the shock. We leave the detail description of the gas ionization Two last terms in Eq. (5) correspond to the injec- to the future work. Here we shall use a simplified tion of thermal protons with momenta p=p , p=p approachandconsiderthenumberdensityofneutrals f b andmassmattheforwardandreverseshockslocated n as a free parameter. n at r =R (t) and r =R (t) respectively . The dimen- The parametersofsupernovamodeling aregivenin f b sionlessparametersη andη determinetheefficiency Table I. The explosion energy E was taken within f b SN of injection. the range (1 2) 1051 erg. The ambient number − · eConf C16-09-04.3 XXV European Cosmic Ray Symposium, Turin, Sept. 4-9 2016 3 2 1 luo,g1(P02c/kρm0R˙sf2−)1,,lologg((Pρg//ρρ00)R,˙lf2o)g,(lToeg/(kPemV/)ρ0R˙f2) 0.1 4πNint(p)p4c/ESN u ρ 0 0.01 P c 0.001 −2 Pg 10−4 p/mc Pm Te 0.01 1 100 104 106 r,pc 4 − 4 6 8 10 12 14 FIG. 3: Spectra of particles produced in the SNR W28 during 37 kyr after explosion. The spatially integrated spectrum of protons (thick solid line), the spectrum of FIG. 1: Radial dependencies of the gas density (thick protonsescapedfrom theremnant(dashedline),theinte- solid line), the gas velocity (dotted line), CR pressure gratedspectrumofelectrons(thinsolidline)andthespec- (thick dashed line), the wave magnetic energy density trum of electrons escaped from the remnant (thin dashed (δB)2/8π = Pm (thin dotted line), the gas temperature line) are shown. Te (thin solid line) and the gas pressure Pg (dashed line) at T =37 kyrin SNR W28. logF(E)E2,ergcm−2s−1 9 TB 10 − Fermi 4πN(p)p4c/ρ0R˙f2 −10 pp −−−− 0.11 −−1121◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦ synch NTB−−−− −−−−−−−−H−ESS 13 − IC 0.01 logE,eV 14 − 5 0 5 10 15 − 0.001 FIG. 4: The results of modeling of electromagnetic radia- p/mc tion of W28. The following radiation processes are taken into account: synchrotron radiation of accelerated elec- 0.01 1 100 104 trons (solid curve on the left), IC emission (dashed line), gamma-ray emission from pion decay (solid line on the FIG. 2: Spectra of accelerated particles at the forward right), thermal bremsstrahlung (dotted line on the left), shock in SNR W28 at T =37 kyr. The spectrum of pro- nonthermalbremsstrahlung(dottedlineontheright). Ex- tons (thick solid line) and electrons (thin solid line) are perimentaldataingamma-rayFermiLAT[16];HESS[17] shown. (data with error-bars) and radio-bands [18] (circles) are also shown. density n was adjusted to reproduce the observ- H able gamma ray fluxes. The number density of neu- contactdiscontinuitybetweentheejectaandtheinter- trals was adjusted to reproduce the spectral shape of stellargasisatr =Rc =6.6pc. Thegastemperature gamma emission. The remnant evolution was calcu- drops sharply downstream of the forward shock due lated up to the instant of time T when the radius of to the radiative cooling. However a thin dense shell the forwardshockequalsto the observableradiusR . is not formed because of CR pressure (cf. [6]). The f The electron to proton ratio K was adjusted to re- centralpartoftheremnantisfilledbythehotrarefied ep produce the observable radio-flux. The numbers in gas with temperature 106 107K. − three last columns of Table I that is the age T, shock SpectraofacceleratedinW28protonsandelectrons speed V and magnetic field strength B just down- at T = 37 kyr are shown in Fig.2. At this point the f f stream of the shock were obtained in the modeling. maximum energy of accelerated protons is about 10 Figures (1)-(5) illustrate the results of our numeri- GeV. This value is in accordance with with the esti- cal calculations. mate of the maximum energy in the partially ionized Radial dependencies of physical quantities in SNR medium Emax = u38n1H/2n−n1 TeV [15]. For shock ve- W28atpresent(T =37kyr)areshownin Fig.1. The locityVf 100kms−1 thatisu8 =0.1andusingthe ∼ eConf C16-09-04.3 4 XXV European Cosmic Ray Symposium, Turin, Sept. 4-9 2016 logF(E)E2,ergcm−2s−1 W44 are shown in Figures 4,5. Thermal emission 9 TB shows two components. One is produced by the hot − Fermi gasintheremnantinteriorwhilelowerenergycompo- −−−111210◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦◦ synch pp−−−−N−−−−T−−B−−−−−−−−−−−− nanendtriescopmrobdinuecdedbbehyinthdethdeenssheogckasfrroandti.atively cooled 13 IV. DISCUSSION − IC logE,eV 14 − 5 0 5 10 15 − High acceleration efficiency (85%) in our modeling FIG. 5: The results of modeling of electromagnetic radia- seems at odds with the lower energetics of 10 20% − tion of W44. The following radiation processes are taken expected for Galactic SNRs. Higher density of cir- into account: synchrotron radiation of accelerated elec- cumstellar medium might permit the lower acceler- trons (solid curve on the left), IC emission (dashed line), ation efficiency. However this will result in the over- gamma-ray emission from pion decay (solid line on the productionofthethermalradioemissioninSNRW44. right), thermal bremsstrahlung (dotted line on the left), Indeed the thermal radio-flux is only slightly below nonthermal bremsstrahlung (dotted line on the right). than the Planck data [21] at high radio-frequencies Experimental data in gamma-ray Fermi LAT [19], (data (see Fig. 5). In this regard the direct DSA scenario witherror-bars) andradio-bands[20,21](circles) arealso shown. seems more probable in comparison with the reaccel- eration scenario in SNR W44 [22]. ThemainpartofgammaemissioninW28andW44 parametersfromthe TableIweindeedgetEmax 10 isproducedbyenergeticprotonsviappcollisions. The ∼ GeV. gamma emission at GeV energies is produced by pro- The spectra of particles Nint produced during 37 tons recently accelerated at the forward shock. Their kyr after supernova explosion are shown in Fig.3. energy is regulated by the neutral damping of MHD Theyarecalculatedviatheintegrationthroughoutthe waves upstream of the shock. Higher energy gamma simulation domain and via the integration on time of emission is produced by particles accelerated earlier the outwarddiffusive flux atthe simulationboundary whentheshockspeedwashigher. Theconfinementof atr =2Rf. About85%ofthekineticenergyofexplo- these particles is very efficient because of the Bohm sionistransferredtocosmicrays. Mostofthisenergy like diffusion in the ionized remnant interior. isgonebyescapedparticles. Themaximumenergyof So we conclude that the neutral damping naturally escaped particles is 100 TeV for this SNR. results in appearance of high energy tails of gamma While themaximumenergyofprotonsattheshock emission in old SNRs. is about10 GeV the interior ofSNR containsprotons The work was supported by Russian Foundation of with energies up to several TeV. They were acceler- Fundamental Research grant 16-02-00255. ated earlier when the shock speed was higher. Results of multi-band modeling of SNRs W28 and [1] Krymsky,G.F.1977, SovietPhysics-Doklady,22,327 [13] Chevalier R.1999, Astrophys.J. 511, 798. [2] Bell, A.R.,1978, MNRAS,182, 147 [14] Reach, W.T., Rho, J., & Jarrett, T.H., 2005, ApJ [3] Axford, W.I., Leer, E. & Skadron, G., 1977, Proc. 618,297 15th ICRC, Plovdiv,90, 937 [15] Drury,L.O.,Duffy,P.,&Kirk,J.G.1996,A&A,309, [4] Blandford,R.D.,&Ostriker,J.P.1978,ApJ,221,L29 1002 [5] Lemoine-Goumard M. Proceedings of IAU Sympo- [16] Abdo, A.A., Ackermann, M., Ajello, M., et al. 2010, sium 2014, 296, 287 ApJ, 718, 348 [6] Lee S.H., Patnaude D.J., Raymond J.C. et al. 2015, [17] Aharonian, F., et al. 2008, A&A,481, 401 ApJ806, 71 [18] Dubner,G.M.,Velazquez,P.F.,Goss,W.M.,&Hold- [7] Zirakashvili, V.N. & Ptuskin V.S. 2012, Astropart. away, M.A., 2000, ApJ, 120, 1933 Phys.,39, 12 [19] Ackermann, M., Ajello, M., Allafort, A., et al. 2013, [8] Bell, A.R.,2004, MNRAS,353, 550 Sci, 339, 807 [9] Zirakashvili, V.N., & Ptuskin, V.S., 2008, ApJ, 678, [20] Castelletti1, G., Dubner, G., Brogan, C., & Kassim, 939 N.E., 2007, A&A 471, 537 [10] Bell, A.R.,2014, Astropart. Phys., 43, 56 [21] Arnaud M., Ashdown M., Atrio-Barandela F. et al., [11] Caprioli, D.,2014, Nuclear Physics B (Proc. Suppl.), 2016, A&A586, 134 256, 48 [22] Cardillo, M.,Amato,E.,&Blasi, P.,2016,A&A595, [12] McKenzie, J.F., & V¨olk, H.J., 1982, A&A,116, 191 58 eConf C16-09-04.3