Numerical Simulations of Cosmic-Ray Acceleration at Core-Collapse Supernovae
Gwenael Giacinti, Vikram Dwarkadas, Alexandre Marcowith, Andrea Chiavassa
NNumerical Simulations of Cosmic-Ray Accelerationat Core-Collapse Supernovae
Gwenael Giacinti ∗ Max-Planck-Institut für Kernphysik, Postfach 103980, 69029 Heidelberg, GermanyE-mail: [email protected]
Vikram Dwarkadas
Department of Astronomy and Astrophysics, University of Chicago, 5640 S Ellis Ave, Chicago,IL 60637, USAE-mail: [email protected]
Alexandre Marcowith
Laboratoire Univers et Particules de Montpellier (LUPM) Université Montpellier, CNRS/IN2P3,CC72, place Eugène Bataillon, 34095, Montpellier Cedex 5, FranceE-mail:
Andrea Chiavassa
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Lagrange, CS 34229, Nice,FranceE-mail: [email protected]
Core-collapse supernovae exploding in dense winds are favorable sites for cosmic-ray (CR) ac-celeration to very high energies. We present our CR-radiation-hydrodynamics simulations of theexplosion of a red supergiant. We study the evolution of the shock wave during the first day fol-lowing core collapse, and estimate the time at which CR acceleration can start. We then calculatethe maximum CR energy at the forward shock as a function of time, and show that it may alreadyexceed 100 TeV only a few hours after shock breakout from the surface of the star. ∗ Speaker. c (cid:13) Copyright owned by the author(s) under the terms of the Creative CommonsAttribution-NonCommercial-NoDerivatives 4.0 International License (CC BY-NC-ND 4.0). http://pos.sissa.it/ a r X i v : . [ a s t r o - ph . H E ] A ug R Acceleration at Core-Collapse Supernovae
Gwenael Giacinti
1. Introduction
Several studies have suggested that core-collapse supernovae (SN) exploding in dense windscould accelerate cosmic-rays (CR) to PeV energies for a few decades [1, 2, 3, 4]. The situationduring the first few days of the SN is somewhat more complicated. After core collapse, a radiation-mediated shock travels through the progenitor. Such a shock is not expected to accelerate CRsbecause its width is larger than the gyroradius of suprathermal particles. Later, at shock breakout(SB), it stalls in the outer layers of the star at an optical depth τ ∼ c / U s , where U s denotes theshock velocity. At this point, the radiation in the immediate shock downstream escapes, producinga flash of photons [5], which accelerates the circumstellar wind. A collisionless shock (CS) laterforms [6]: The dilution of photons as 1 / R , where R is the stellar radius ( R =
2. Numerical simulations
Stellar convection modelling is carried out using three-dimensional (3D) radiative hydrody-namical (RHD) code CO BOLD [7]. It solves the coupled non-linear equations of compressiblehydrodynamics and non-local radiative energy transfer in the presence of a fixed external gravita-tional field. The computational domain is a cubic grid equidistant in all directions, and the sameopen boundary condition is employed for all sides of the computational box (see, e.g., Ref. [8]).The code employs realistic input physics either for the equation of state [7] or for the radiativetransfer opacities [9]. A detailed and precise solution of the radiative transfer is essential for arealistic treatment of convection because it is the radiative losses in the surface layers of the starthat drain the convective movements and thus influence the whole simulation domain.The fundamental stellar parameters of the RHD simulation used in this work are reported inTable 1 of Ref. [10]. We plot in Fig. 1 a snapshot of the resulting density profile. In the left panel,we plot the logarithm of the density in a plane cut of the 3D simulated domain. The plane containsthe centre of the star, R =
0, in the centre of the panel. The very dense central region of the staris not represented. The dense, optically thick layers of the progenitor are in yellow/orange, whilethe surrounding dark regions correspond to the optically thin circumstellar material. Deviationsfrom spherical symmetry are clearly visible. In the right panel, we plot with the thick red line theradial density profile of this star, averaged over all radial directions. The thin lines, denoted “D1”to “D4”, represent the density profiles along four given radial directions.The simulation of the SN is carried out using our Eulerian 1D-spherical radiation-hydro-dynamics code, presented in Ref. [11]. We use the average radial density and temperature profilesof the above star as an input. CRs have been added as a pressure term in our code, see Ref. [12] formore details. The code is two-temperature, i.e. electron and proton temperatures are assumed to1
R Acceleration at Core-Collapse Supernovae
Gwenael Giacinti -20-18-16-14-12-10-8-6 ρ [ g / c m ] R [m] AverageD1D2D3D4
Figure 1: (Initial) density profile of the progenitor red supergiant star considered in this study, and of itssurroundings. Left panel: log (cid:0) ρ / ( g · cm − ) (cid:1) in a plane containing the centre of the star R =
0, located inthe centre of the panel. See the colour code in the bar on the right-hand side. Right panel: Correspondingradial density profile averaged over all radial directions (thick red line), and density profiles along four givenradial directions, “D1” to “D4” (thin lines, see the key). be equal. For the radiation, we use a gray frequency average, and represent it by its internal energy E rad with characteristic temperature T rad = ( cE rad / σ ) / , where σ denotes the Stefan-Boltzmannconstant. The radiation transport is solved using a “square-root” flux-limited diffusion approxima-tion. We take into account Compton cooling and bremsstrahlung for the transfer of energy betweenfluid and radiation, using the formulae from Ref. [6]. At t = erg at the centre of the star as a spike in temperature. We follow the evolutionof the physical conditions around the shock in the simulation, and determine the formation time ofthe CS, using the method of Ref. [11]. Once the CS is formed, we calculate, at each time in thesimulation, the typical Coulomb loss time of suprathermal particles at the CS, τ Coul , T e , u , and theirtypical acceleration time to 1 GeV, τ acc , . For the latter, we assume Bohm diffusion of the CRs,in an initial circumstellar magnetic field of 1 mG strength. For stronger magnetic fields, particleacceleration would proceed faster. We also calculate the characteristic times for CR energy lossesdue to inelastic pp collisions, τ pp , and adiabatic losses, τ adiab . For the latter, we take the conserva-tive lower value τ adiab = R s / U s , where R s denotes the radius of the CS. See Refs. [11, 12] for moredetails. We assume that particle acceleration starts in our simulation once the acceleration time to1 GeV is shorter than all loss times.We assume an E − spectrum for the CRs accelerated at the CS. We further assume that mag-netic field amplification at the shock is due to Bell’s instability [13], which is driven by the escapeof the highest energy CRs in the upstream of the CS [2]. We calculate the maximum CR energy atthe CS, E max , using the results of Ref. [2]. At each time step in the simulation, we calculate the CRcurrent escaping ahead of the CS, the growth rate of the instability in the upstream, and E max . SeeRef. [12] for technical details.Finally, we calculate the CR acceleration time to E max in the amplified magnetic field fromBell’s instability, τ acc , E max , and verify that it does not exceed the characteristic CR loss times τ adiab ,2 R Acceleration at Core-Collapse Supernovae
Gwenael Giacinti τ pp , and τ p γ , where τ p γ is the loss time due to pion production in inelastic p γ collisions on thehigh-energy photons emitted by the CS. For τ p γ , we calculate here a conservative (and pessimistic)lower value, τ p γ , min , which assumes that the CS radiates all the energy it processes in photons withenergies set to the threshold for pion production.
3. Results
We present in Fig. 2 results from our CR-radiation-hydrodynamics simulation of the explosionof the star, at three different times: Shortly before SB in the upper row, during SB in the middlerow, and soon after SB in the lower row. In the left column, we show U / c , the velocity of thefluid normalized to c , as a function of the radius R , counted from the centre of the star. In the rightcolumn, we plot the electron temperature T e (black lines), the fluid density ρ (red), the radiationenergy density ε rad (orange), and the CR energy density ε CR (green), as functions of R . See thekeys in the Figure for the units and normalizations. In the initial profile, the optical depth is equalto 1 at R ∼ . · m. In the first row, the shock is still inside the optically thick layers of the star,and is radiation-mediated. It is located at R (cid:39) . · m, see the vertical jump in the blue curve inthe left panel. The shock thickness is several photon mean free paths, and the post-shock radiationis well confined, see the orange line in the right panel. In the middle row, the shock reaches theouter layers of the star. The radiation in the region immediately behind the shock starts to escape,and creates a flash of photons: The inflection around R (cid:39) · m in the curve for the radiationenergy density (orange line in the right panel) corresponds to the front of SB photons which escapeat the speed of light towards R → ∞ . The escaping radiation accelerates the circumstellar medium,which is visible in the left panel: The discontinuity in U / c is now significantly wider, and spansfrom R (cid:39) · m to R (cid:39) m. In the lower row, the CS has formed: In the left panel, it appearsas the abrupt discontinuity in U / c at R (cid:39) . · m, embedded within a significantly broadertransition region smoothed by the radiation from SB. One can see in the right panel a spike intemperature ( ∼
100 keV), which corresponds to the shock-heated region in the downstream of thenewly formed CS. CR acceleration has started at the CS, cf. the green curve.We show in Fig. 3 our calculations of the typical timescales for particle acceleration and energylosses, as functions of time since core collapse. In the left panel, the solid red line correspondsto τ acc , , the CR acceleration time to 1 GeV in an initial circumstellar magnetic field set to1 mG. It is plotted against the loss time for inelastic pp collisions ( τ pp , solid blue line), as wellas the adiabatic ( τ adiab , dotted black), and Coulomb loss times using the electron temperature inthe CS upstream from the simulation ( τ Coul , T e , u , solid orange). For information, we also show thelimiting Coulomb loss time one would have assuming that the electron temperature is zero ( τ Coul , ,dashed orange). This plot shows that Coulomb losses are the limiting factor for the onset of CRacceleration here. In this simulation, the CS forms around t (cid:39) t (cid:39) ∼ E max ( τ acc , E max , solid red line) and to 1 TeV ( τ acc , , dashed red) in the amplified magnetic fieldand assuming Bohm diffusion. We plot it against the pp and adiabatic loss times, as well as the3 R Acceleration at Core-Collapse Supernovae
Gwenael Giacinti U / c R [m] e /keV ρ /(10 -10 g/cm ) ε rad /(10 J/m ) ε CR /(J/m ) U / c R [m] e /keV ρ /(10 -10 g/cm ) ε rad /(10 J/m ) ε CR /(J/m ) U / c R [m] e /keV ρ /(10 -10 g/cm ) ε rad /(10 J/m ) ε CR /(J/m ) Figure 2:
Simulation of the explosion of a red supergiant. Velocity of the fluid normalized to c (left column),and electron temperature T e , fluid density ρ , radiation and CR energy densities ε rad , CR (right column) at threedifferent times: shortly before (upper row), during (middle row) and soon after (lower row) SB —the opticaldepth is equal to 1 at R ∼ . · m in the initial profile. See the keys in the right column for the linecolours and normalizations. R Acceleration at Core-Collapse Supernovae
Gwenael Giacinti c ha r a c t e r i s t i c t i m e sc a l e s [ s ] time since core collapse [s] τ acc,1GeV τ Coul,0 τ Coul,T e,u τ pp τ adiab c ha r a c t e r i s t i c t i m e sc a l e s [ s ] time since core collapse [s] τ acc,E max τ acc,1TeV τ pp τ p γ ,min τ adiab Figure 3:
Characteristic timescales for particle acceleration and energy losses, versus time since corecollapse. See the keys for the line types and colours, and see the text for more details. minimum p γ loss time for CRs with energy E max ( τ p γ , min , dashed blue line). Fig. 3 (right panel)demonstrates that CR acceleration to E max is not hindered by any of these losses. E m a x [ e V ] time since core collapse [s]Average E m a x [ e V ] time since core collapse [s]D1D2D3D4D5D6 Figure 4:
Maximum CR energy, E max , versus time since core collapse, for the initial (circum)stellar densityprofile averaged over all radial directions (left panel) and for the initial density profiles along six given radialdirections, “D1” to “D6” (right panel). In the left panel, the solid line corresponds to negligible Coulomblosses, and the dashed line assumes that CR acceleration starts when τ acc , ≤ τ Coul , T e , u . We plot in Fig. 4 our calculations of the maximum CR energy at the forward shock, E max ,versus time since core collapse. Results in the left panel are for a simulation in the average densityprofile of the SN progenitor and its surroundings. The solid red line is for the optimistic case ofCR acceleration starting as soon as the CS is formed, i.e. for negligible Coulomb losses. This maybe possible if the magnetic field strength around the progenitor star is (cid:29) τ acc , ≤ τ Coul , T e , u . These results show that CRsreach ≥
100 TeV energies typically ∼ E max startsto plateau to values of a few hundreds of TeV by the end of the first day following core collapse.The fact that both lines reach almost the same limiting value shows that Coulomb losses have littleimpact on E max , and only slightly delay the beginning of particle acceleration. The results in the5 R Acceleration at Core-Collapse Supernovae
Gwenael Giacinti right panel are for six density profiles along six given radial directions, denoted as “D1”, “D2”,...,“D6”. The profiles D1 to D4 correspond to those plotted in Fig. 1 (right panel) —the profiles D5and D6 are not shown there due to limited space in the plot. The diversity of the curves in Fig. 4(right panel) provides an estimate of how the 3D geometry of the progenitor affects the results for E max versus time. First, the times at which the CS forms and CR acceleration starts depend on theconsidered direction around the star. Indeed, the progenitor is not perfectly spherically symmetric,which implies that SB and thence the onset of CR acceleration occurs earlier in some directionsthan in others. However, we caution that a 3D simulation would be needed for a more watertightestimate of the spread in SB times: In 3D, the forward shock can get around dense clumps, whichcannot be taken into account in our simulation. Second, the limiting value of E max at which eachcurve plateaus depends on the direction. This means that the limiting E max is not exactly the same ineach region of the forward shock, at least during the first day of the SN. For example, the magentacurve reaches a value that is a few times larger than that of the black curve. Looking at Fig. 1 (rightpanel), one can see that the circumstellar material at R ≥ · m is about an order of magnitudedenser in the direction D3 (magenta curve) than in the direction D1 (black one). These results arein line with the calculations of Ref. [2], where E max is found to be larger in denser regions —at U s fixed. All six curves nonetheless tend towards the same order-of-magnitude values of E max , whichshows that taking the average density profile provides a reasonable first estimate.
4. Conclusions
We study here the formation of a CS and the beginning of CR acceleration at a Type II SN, dur-ing the first day that follows core collapse. We present the CR-radiation-hydrodynamics simulationof the explosion of a red supergiant with a realistic density profile. We find that a CS soon formsafter SB, and that CR acceleration starts soon after. For the progenitor considered here, Coulomblosses do not delay CR acceleration by more than an hour. The maximum CR energy at the forwardshock quickly increases, and already reaches a few hundreds of TeV only hours after SB. It startsto plateau towards its maximum value by the end of the first day of the SN. The density of thecircumstellar medium in the immediate surroundings of the progenitor depends on the consideredradial direction around the star, and we find that CRs can be accelerated to higher energies in denserregions. Our results suggest that core-collapse SNe are plausible sources of very high-energy CRs,even in their earliest phases.
Acknowledgments
We thank Matthieu Renaud, Vincent Tatischeff, and Pierre Cristofari for useful discussions.This research collaboration is supported by a grant from the FACCTS program to the University ofChicago (PI: VVD). This work is supported by the ANR-14-CE33-0019 MACH project.
References [1] V. Tatischeff,
Radio emission and nonlinear diffusive shock acceleration of cosmic rays in thesupernova SN 1993J , Astron. Astrophys. (2009) 191 [ arXiv:0903.2944 ]. R Acceleration at Core-Collapse Supernovae
Gwenael Giacinti[2] A. R. Bell, K. Schure, B. Reville, G. Giacinti,
Cosmic ray acceleration and escape from supernovaremnants , MNRAS (2013) 415 [ arXiv:1301.7264 ].[3] A. Marcowith, M. Renaud, V. Dwarkadas, V. Tatischeff,
Cosmic-ray acceleration and gamma-raysignals from radio supernovæ , Nucl. Phys. B Proc. Suppl. (2014) 94 [ arXiv:1409.3670 ].[4] A. Marcowith, V. Dwarkadas, M. Renaud, V. Tatischeff, G. Giacinti,
Core collapse supernovae asCosmic Ray sources , MNRAS (2018) 4470 [ arXiv:1806.09700 ].[5] S. A. Colgate,
ApJ (1974) 321.[6] R. A. Chevalier, R. I. Klein,
Nonequilibrium processes in the evolution of type II supernovae , ApJ (1979) 597.[7] B. Freytag, M. Steffen, H.-G. Ludwig, S. Wedemeyer-Böhm, W. Schaffenberger, O. Steiner,
Simulations of stellar convection with CO5BOLD , J. Comput. Phys. (2012) 919[ arXiv:1110.6844 ].[8] A. Chiavassa, B. Freytag, T. Masseron, B. Plez,
Radiative hydrodynamics simulations of redsupergiant stars. IV. Gray versus non-gray opacities , Astron. Astrophys. (2011) A22[ arXiv:1109.3619 ].[9] B. Gustafsson, B. Edvardsson, K. Eriksson, U. G. Jørgensen, Å Nordlund, B. Plez,
A grid of MARCSmodel atmospheres for late-type stars. I. Methods and general properties , Astron. Astrophys. (2008) 951 [ arXiv:0805.0554 ].[10] K. Kravchenko, S. Van Eck, A. Chiavassa, A. Jorissen, B. Freytag, B. Plez,
Tomography of cool giantand supergiant star atmospheres. I. Validation of the method , Astron. Astrophys. (2018) A29[ arXiv:1711.08327 ].[11] G. Giacinti, A. R. Bell,
Collisionless shocks and TeV neutrinos before supernova shock breakout froman optically thick wind , MNRAS (2015) 3693 [ arXiv:1503.04170 ].[12] G. Giacinti et al. , in preparation.[13] A. R. Bell,
Turbulent amplification of magnetic field and diffusive shock acceleration of cosmic rays , MNRAS (2004) 550.(2004) 550.