Multi-dimensional simulations of the expanding supernova remnant of SN 1987A
T.M. Potter, L. Staveley-Smith, B. Reville, C.-Y. Ng, G. V. Bicknell, R. S. Sutherland, A. Y. Wagner
aa r X i v : . [ a s t r o - ph . H E ] S e p Draft version July 12, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
MULTI-DIMENSIONAL SIMULATIONS OF THE EXPANDING SUPERNOVA REMNANT OF SN 1987A
T.M. Potter , L.Staveley-Smith , B. Reville , C.-Y. Ng , G. V. Bicknell , R. S. Sutherland , A. Y.Wagner Draft version July 12, 2018
ABSTRACTThe expanding remnant from SN 1987A is an excellent laboratory for investigating the physics ofsupernovae explosions. There are still a large number of outstanding questions, such the reason forthe asymmetric radio morphology, the structure of the pre-supernova environment, and the efficiencyof particle acceleration at the supernova shock. We explore these questions using three-dimensionalsimulations of the expanding remnant between days 820 and 10,000 after the supernova. We combine ahydrodynamical simulation with semi-analytic treatments of diffusive shock acceleration and magneticfield amplification to derive radio emission as part of an inverse problem. Simulations show that anasymmetric explosion, combined with magnetic field amplification at the expanding shock, is able toreplicate the persistent one-sided radio morphology of the remnant. We use an asymmetric Truelove& McKee progenitor with an envelope mass of 10 M ⊙ and an energy of 1 . × J . A terminationshock in the progenitor’s stellar wind at a distance of 0 . ′′ − . ′′
51 provides a good fit to the turn on ofradio emission around day 1200. For the H ii region, a minimum distance of 0 . ′′ ± . ′′
01 and maximumparticle number density of (7 . ± . × m − produces a good fit to the evolving average radiusand velocity of the expanding shocks from day 2000 to day 7000 after explosion. The model predictsa noticeable reduction, and possibly a temporary reversal, in the asymmetric radio morphology of theremnant after day 7000, when the forward shock left the eastern lobe of the equatorial ring. Subject headings: acceleration of particles, hydrodynamics, ISM: supernova remnants, radiation mech-anisms: non-thermal, supernovae: general, supernovae: individual: (SN 1987A) INTRODUCTION
Supernovae play an important role in the evolution ofthe Universe: providing a source of heavy elements, driv-ing winds to regulate star formation, producing cosmicrays, magnetic fields, neutron stars and black holes. Asthe brightest supernova since 1604, the type II-P super-nova SN 1987A has been the most well studied super-nova in history. It was the only supernova to be as-sociated with a neutrino detection (Hirata et al. 1987;Bionta et al. 1987; Aglietta et al. 1987; Alexeyev et al.1988), and one of the few supernovae whose progeni-tor star was observed prior to explosion. Rousseau et al.(1978) and Walborn et al. (1989), classified the progeni-tor star Sk 69 ◦
202 as a B3 I blue supergiant (BSG). TheBSG had an estimated surface temperature of 16,000K; a mass of 19 ± M ⊙ ; an envelope mass of 5 − M ⊙ (Woosley 1988); and an estimated wind velocityand mass loss rate of 450 km s − and 7 . × − M ⊙ [email protected] International Centre for Radio Astronomy Research(ICRAR) M468, The University of Western Australia, 35 Stir-ling Highway, Crawley, WA 6009, Australia School of Earth and Environment, M004, The University ofWestern Australia, 35 Stirling Highway, Crawley, WA 6009, Aus-tralia ARC Centre of Excellence for All-sky Astrophysics (CAAS-TRO) Centre for Plasma Physics, Queen’s University Belfast, Uni-versity Road, Belfast BT7 1NN, United Kingdom Department of Physics, The University of Hong Kong, Pok-fulam Road, Hong Kong Research School of Astronomy and Astrophysics, AustralianNational University, Canberra, ACT 0200, Australia. Center for Computational Sciences, Tsukuba University,Tsukuba, Ibaraki, 305-8577, Japan yr − (Chevalier & Dwarkadas 1995). This was surpris-ing, as the expected progenitors of core-collapse super-novae were red supergiants (RSG’s). Soon after core col-lapse, a UV flash from shock breakout ionised the ma-terial surrounding the BSG and revealed a central equa-torial ring, accompanied above and below by two fainterrings (Gouiffes et al. 1989; Plait et al. 1995). Subsequentanalysis of the echoes from the UV flash Sugerman et al.(2005) showed that, at an assumed distance of 50 kpc,the equatorial ring is the waist of a much larger peanutshaped structure which extends around 2 × m (6 . × m (3 . M ⊙ for the nebula. The amount of material in thecircumstellar environment suggests the progenitor pre-viously went through a phase of high mass loss as aRSG before transforming into a BSG prior to explo-sion. The cause of the transformation is uncertain. Pos-sible explanations for the transformation involve a bi-nary merger (Podsiadlowski & Joss 1989), or low metal-licity in the progenitor (Woosley et al. 1987). Generalconsensus is that the transformation took place approx-imately 20,000 years prior to the explosion, and the fastwind from the BSG interacted with the relic RSG windto form the hourglass and rings (Crotts & Heathcote1991; Blondin & Lundqvist 1993; Crotts & Heathcote2000; Podsiadlowski et al. 2007). This interactionhas been successfully modelled with hydrodynamical(Blondin & Lundqvist 1993), magnetohydrodynamical(Tanaka & Washimi 2002) and smoothed particle hydro-dynamics simulations (Podsiadlowski et al. 2007). Potter et al. 2014 Radio observations
Over the last 25 years the interaction of the expandingsupernova remnant has been monitored at wavelengthsspanning the electromagnetic spectrum. Observationsat radio frequencies ranging from 843 MHz to 92 GHzhave traced the evolution of flux density, spectral in-dex, and radio morphology of the remnant, from a fewdays after explosion to the present day (Turtle et al.1987; Staveley-Smith et al. 1992; Ball & Kirk 1992b;Gaensler et al. 1997; Ball et al. 2001a; Manchester et al.2005; Gaensler et al. 2007; Staveley-Smith et al. 2007;Ng et al. 2008; Potter et al. 2009; Zanardo et al. 2010;Ng et al. 2011; Laki´cevi´c et al. 2012). Approximatelyfour days after core collapse, radio emission peakedaround 150 mJy at 1 GHz as the supernova transi-tioned from optically thick to optically thin regimes(Turtle et al. 1987). Over the next few months theemission faded to undetectable levels as an expandingset of forward and reverse shocks propagated througha rarefied BSG wind. About 1200 days after core col-lapse, the shocks crashed into the termination shockof the pre-supernova BSG wind. Radio emission fromthis interaction became visible, and the remnant wasre-detected at 843 MHz by the Molonglo ObservatorySynthesis Telescope (MOST) (Ball et al. 2001a) and at1-8 GHz by the Australia Telescope Compact Array(ATCA) Staveley-Smith et al. (1992). The flux den-sity increased rapidly following re-detection as a newset of shocks began propagating away from the BSGwind boundary. Since day 2500, flux density hasbeen growing exponentially at all frequencies (Ball et al.2001a; Manchester et al. 2002; Staveley-Smith et al.2007; Ng et al. 2008; Zanardo et al. 2010) as the shocksencounter relics from the RSG wind in the equatorialplane, and a hot BSG wind beyond the termination re-gion at high latitudes. When the radio emission re-turned, the spectral index α ( F ( ν ) ∝ ν − α ) was between0 . .
7. Around day 2200 the spectrum had be-come its softest, with α around 1 .
05. Since day 2500the spectral index has been hardening linearly as α ( t ) =0 . − . × ( t − / t is expressed indays (Zanardo et al. 2010). Presently, (day 9200), thespectral index has returned to a value between 0 . . . ′′ ±
200 km s − ; aroundthree times faster than the 1900 ±
400 km s − obtainedfor the western lobe.In Ng et al. (2008), the topology of the radio emittingshell from SN 1987A was modelled using a shell of finitewidth and truncated to lie within a half-opening angle ofthe equatorial plane. They projected the truncated shellsto the u − v plane, and used least squares optimisation tofind shells that fitted u − v data from the 8GHz ATCAmonitoring observations. As a result, we have estimatesof the radius, opening angle, and thickness of the ex-panding shell of emitting material. The estimated shellradius from the models indicate that the emitting regionhad a minimum average expansion of 30,000 km s − from1987 to 1992 (Gaensler et al. 1997; Ng et al. 2008). Af-ter encountering the relic RSG material inside the ring(Chevalier & Dwarkadas 1995) the average speed of thesupernova shocks slowed to around 4000 ±
400 km s − and remained at that rate of expansion until day 7000when the emitting region appears to become more ring-like (Ng et al. 2013). Theoretical models of radio emission from theexpanding shocks
Radio emission from SN 1987A is thought to arise fromrelativistic electrons accelerated at the supernova shockfront. Diffusive shock acceleration (DSA) (Krymskii1977; Axford et al. 1977; Blandford & Ostriker 1978;Bell 1978) is believed to be the main source of relativis-tic electrons at such shocks (Melrose 2009). It producesa non-thermal population of energetic electrons, whoseisotropic phase-space distribution in momentum f ( p ),has a power-law form f ( p ) ∝ p − b . For a strong shockwith a compression ratio of ζ = 4, and a ratio of specificheats of γ = 5 /
3, diffusive shock acceleration predictsthe index on the distribution is b = ζζ − .Early models of the radio emission from SN 1987Awere constructed by calculating radio emission usingpower-law distributions evolving in an expanding shellof hot gas Turtle et al. (1987); Storey & Manchester(1987). The underlying hydrodynamics of the shockwere greatly simplified by the assumption that the shellof hot, radio-emitting gas underwent self-similar ex-pansion. Turtle et al. (1987) obtained b ≈ b in therange 3 . − .
33 by fitting an expanding shell model tothe same data. Their model included synchrotron self-absorption and free-free absorption. Waßmann & Kirk(1991) also proposed a model for the early rise and fallin radio emission in which shock-accelerated electrons“surf” outwards from the shock along a pre-existingspiral magnetic field line. Ball & Kirk (1992a) andKirk et al. (1994) developed a time-dependent, two-zonemodel to evolve the radio-emitting electrons in the adi-abatically expanding downstream. They were able to fitthe 4 . b ≈ .
8. Theypostulated that the softening of the electron spectrumwas due to cosmic ray feedback on the shock.The later models of Duffy et al. (1995) andBerezhko & Ksenofontov (2000), included cosmicray feedback. The resulting weakening of the shock, asit decelerated in the cosmic-ray pressure gradient, mod-ified the compression ratio at the density discontinuityto around 2.7 and softened the electron momentumspectrum index b from 4 to 4 .
8. This provided a physicalmotivation for the index observed in SN 1987A.Models of radio emission for SN 1987A up to thispoint used a pre-existing magnetic field that was com-pressed by the shock. In recent years it has beenshown that plasma-instabilities excited by cosmic rayscan amplify a background magnetic field B (Bell 2004)by up to a factor of ≈
45 (Riquelme & Spitkovsky2009, 2010). The efficiency of magnetic field ampli-fication is an active topic being studied with simula-tions. From Bell (2004), the resulting dependence ofmagnetic field B on shock velocity v s scales as B ∝ v / s .Subsequent models of radio emission from SN 1987Aincluded prescriptions for magnetic field amplification(Berezhko & Ksenofontov 2006; Berezhko et al. 2011).These models produce a strong downstream magneticfield of 2 × − T, in contrast to previous models withan assumed magnetic field of 10 − − − T (Duffy et al.1995; Berezhko & Ksenofontov 2000). The strong depen-dence of magnetic field upon shock velocity raises the in-teresting possibility that synchrotron emission is stronglydependent on the shock velocity and thus the asymmetryof the radio remnant may be a byproduct of an asym-metric explosion. The fraction of electrons injected intothe shock, χ el , is a product of the microphysics of theshock and is not well understood. Kinetic plasma sim-ulations have made progress in this direction in recentyears (McClements et al. 2001; Matsumoto et al. 2012;Caprioli & Spitkovsky 2014, e.g.). The total synchrotronluminosity from supernovae is dependent upon both thestrength of the magnetic field and χ el . With a pre-existing magnetic field model Berezhko & Ksenofontov(2000) found χ el = (1 − × − to be a good fit to radioobservations. This was later modified to χ e = 6 × − in Berezhko et al. (2011), after allowing for additionalnon-linear magnetic field amplification. The case for a multi-dimensional simulation
A feature of previous models of radio emission fromSN 1987A is their incorporation of varying degrees ofspherical-symmetry. Such models cannot account for theinteraction of the shock with the ring, nor can they repli-cate the evolving asymmetrical radio morphology of theremnant. The need for multi-dimensional simulations ofSNR 1987A has been made clear e.g Dwarkadas (2007).Until recently, modelling SNR 1987A and other super-nova remnants in multiple dimensions has been regardede.g Dewey et al. (2012), as highly complex and compu-tationally challenging, thus limiting the potential for it-erative exploration in parameter space. However, recentadvances in computing have reduced model realisationtimes, enabling more possibilities for model explorationin higher dimensions. With an aim to address the above questions and chal-lenges we present results from a new three-dimensionalsimulation of the interaction of the shock from SN 1987Awith its pre-supernova environment. This work is moti-vated by the need to overcome some of the limitationswith previous one-dimensional models, such as the inabil-ity to adequately model the evolving radio morphologyof the remnant. In a fully three dimensional simulationwe can: (1) Test a hypothesis that magnetic field ampli-fication in combination with an asymmetric explosion isthe cause of the observed persistent asymmetry in the ra-dio morphology ; (2) gain insight into the 3D structure ofthe pre-supernova material; (3) obtain an estimate of theinjection efficiency of shock acceleration by direct com-parison with observations; (4) and make a prediction onhow the remnant might evolve if the model is accurate.The quality and relative abundance of observationalmonitoring data makes SN 1987A an ideal candidate foran inverse modelling problem. SIMULATION TECHNIQUE
Non-thermal emission from the remnant is computedin two stages. First we use a hydrodynamics code to sim-ulate the expansion of the supernova shock into a modelenvironment. For this we use FLASH (Fryxell et al.2000) to propagate the fluid according to Eulerian con-servation equations of inviscid ideal gas hydrodynamicsin a Cartesian grid. The equations solved for in FLASHare as follows: ∂ρ∂t + ∇ · ρ v = 0 , (1) ∂ρ v ∂t + ∇ · ( ρ vv ) + ∇ P = 0 , (2) ∂ρE∂t + ∇ · [( ρE + P ) v ] = 0 . (3)The density and pressure is ρ and P , and v repre-sents the three components of fluid velocity. The spe-cific energy E = ǫ + v is the sum of specific internal( ǫ ) and specific kinetic (cid:0) v (cid:1) energies. We assume anideal monatomic plasma with γ = 5 / P = ( γ − ρǫ to close the set of equations. Radiativecooling is implemented using a temperature-dependentcooling function (discussed in Section 2.3.6). The sec-ond stage involves post-processing the hydrodynamicsoutput. We locate both forward and reverse shocksand apply semi-analytic models of diffusive shock ac-celeration to fix the momentum distribution f ( p ) at lo-cations in the grid immediately downstream from theshocks. The magnetic field is assumed to be amplifiedat the shock through cosmic-ray current driven instabili-ties (Bell 2004) and is evolved adiabatically downstreamof the shock. The development of both the momentumdistribution and magnetic field intensity for shocked vox-els is implemented using a simple up-winded advectionscheme, with source terms where necessary (see sections2.3.4 & 2.3.5). Synchrotron emission and absorptionis calculated from the resulting momentum distributionand magnetic field using standard analytic expressions.Comparisons of the simulated radius, flux density andmorphology with observations are used to fit parameters Potter et al. 2014in the initial environment. Initial environment
The initial conditions for the pre-supernova environ-ment are designed to be physically motivated as much aspossible, while fitting the monitoring observations. Weuse estimates from the literature to build an approximatemodel and refine it using inverse modelling by compar-ison with observations. In regions of the model wherewe do not have good data we use results from a pre-supernova environment formation simulation (EFS, seeAppendix D for details). We report the results of ourbest 3D supernova model, whose parameters have beentuned to fit radio observations such as the expandingradius, flux density and radio morphology. We do notclaim that the residual of the model has been minimised,rather we report on a model that fits the majority ofradio observations and can be used to make meaningfuldeductions about the real supernova.
Computational domain
The scope of our simulation is to model the interac-tion of the supernova shock with the inner hourglass andequatorial ring for a simulated period of 10,000 days afterexplosion. In order to simplify radiative transfer calcula-tions we used a Cartesian grid aligned with the plane ofthe sky and centred on the progenitor. The environmentwas then inclined within this grid in order to match theorientation of the equatorial ring. The positive X axiscorresponds to West on the plane of the sky, the positiveY axis is North, and the positive Z axis points to Earth.The grid is a cube with length 256 cells (3 . × m) ona side. This corresponds to an angular separation of 4 . ′′ i x = 41 ◦ , i y = − ◦ , i z = − ◦ . In prac-tice we use a series of successive X, Y, and Z rotations toachieve the observed inclination. The required rotationsare ( x rot = 41 ◦ , y rot = − ◦ , z rot = − ◦ ). Within thisinclined environment we use the Cartesian coordinates( x ′ , y ′ , z ′ ) and cylindrical coordinates ( s ′ , φ ′ , z ′ ) centredon the progenitor. If the environment were not inclined,the positive z ′ axis would point to Earth and the angle φ would be measured as a counterclockwise rotation fromthe X axis. We also use r , the radial distance from theprogenitor. Model features
Within our domain, the main components of the pre-supernova environment surrounding SN 1987A are as fol-lows. Outwards from the progenitor a supersonic, lowdensity wind extends to a termination shock located at
Figure 1.
Cartoon of the equatorial ring showing the inclinationof the environment at angles i x = 41 ◦ , i y = − ◦ , i z = − ◦ .The rotated cylindrical coordinate system radial coordinate s ′ andvertical coordinate z ′ has z ′ parallel to the plane normal. a radius approximately 3 . × m (0 . ′′ . × m (1 . ′′ ii region (Chevalier & Dwarkadas 1995). Theequatorial ring lies within the H ii region at a distance of(6 . ± . × m or (0 . ′′ ± . ′′
1) (Plait et al. 1995;Sugerman et al. 2005) and forms the waist of the hour-glass. Exterior to the hourglass the density fades to thebackground density as s ′− near the waist and as s ′− . at large | z ′ | .In Figures 2, 3 and 4 are cross sections of the model en-vironment, obtained by slicing halfway along the X axisin the log of particle number density, temperature, andvelocity. In Figure 2 we have labelled the main featuresof the model, including the progenitor, free BSG wind,shocked BSG wind, mach disk, HII region, and equato-rial ring. Details of how we arrived at the model will bediscussed in forthcoming sections. Properties of the plasma
For all simulations we assume an ideal monatomicgas with γ , the ratio of specific heats, set to .An average atomic mass per particle was derived us-ing the MAPPINGS shock and photo-ionization code(Sutherland & Dopita 1993) and the abundances of theinner ring from Table 7 of Mattila et al. (2010). The de-rived average atomic mass per particle, µ , is 0.678 amu,and the average number of particles per electron, Hydro-gen atom, ion, and nucleon is 2.62, 3.51, 3.00, and 1.62respectively. Progenitor
The progenitor star Sanduleak -69 ◦
202 was observedto be a B3 I blue supergiant (Rousseau et al. 1978; −2−1012 Distance (arcsec) −2−1012 D i s t a n c e ( a r c s e c ) −1.5 −1.0 00.5 0.0 0.5 1.0 1.5Distance (10 m)01.501.000.50.00.51.01.5 D i s t a n c e ( m ) E(uato)ial )ingS ocked BSG .indSupe)no−aen−elopeHalf-openingangle θ ho BSG freewindHII region and hourglass Mach diskTo Earth 567891011121314 L o g nu m b e r d e n s i t y ( m − ) Figure 2.
Slice through the three dimensional volume takenhalfway along the x axis. The variable shown is the base 10 log ofthe number density (m − ). Earth is to the right and north is up.Features of the plot are the central supernova envelope and BSGfree wind region, the H ii region, hourglass, and equatorial ring. −2−1012 Distance (arcsec) −2−1012 D i s t a n c e ( a r c s e c ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) L o g t e m p e r a t u r e ( K ) Figure 3.
Same slice as in Figure 2 but in the log of temperature.The highest temperature material is in the core of the progenitorand the shocked BSG wind.
Walborn et al. 1989), with an estimated surface temper-ature of 16,000 K (Arnett et al. 1989), a mass of 19 ± M ⊙ an envelope mass of 5 − M ⊙ (Woosley 1988;Nomoto et al. 1988) and an estimated wind velocity andmass loss rate of 450 km s − and 7 . × − M ⊙ yr − (Chevalier & Dwarkadas 1995). As the initial stages ofthe explosion are too small to represent at our chosengiven grid resolution, we used an analytic solution toevolve the supernova to a size large enough to repre-sent within a sphere 20.25 voxels or r sn = 2 . × m (0 . ′′
36) in radius. The self-similar analytical solutionsin Chevalier (1976) and Truelove & McKee (1999) de-scribe the propagation of a supernova into a circumstel-lar environment with density profile of ρ ( r ) ∝ r − s . Thedensity of the expanding supernova envelope varies with −2−1012 Distance (arcsec) −2−1012 D i s t a n c e ( a r c s e c ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) L o g v e l o c i t y ( m / s ) Figure 4.
Same slice as in Figure 2, but in the log of velocity.The progenitor is the central core of highest velocity in the range(10 − ) m/s. velocity as ρ ( v ) ∝ v − n . We use n = 9, s = 2 fromChevalier & Dwarkadas (1995) and modify the analyticsolution of Truelove & McKee (1999) to include inter-nal energy and introduce an asymmetric explosion in theeast-west direction. Details of these modifications aredescribed at length in Appendix A. The mechanics ofthe model supernova envelope are completely describedby: M env the mass of the supernova envelope; E tot , thetotal mechanical energy of the explosion (kinetic plusinternal), χ ke , the ratio of kinetic energy to total; and χ asym , the ratio of kinetic energy in the eastern hemi-sphere of the supernova envelope with respect to thewestern hemisphere. We adopt E tot = 1 . × J, and M env = 10 M ⊙ , consistent with the results of Woosley(1988); Arnett & Fu (1989); Bethe & Pizzochero (1990);Shigeyama & Nomoto (1990). In order to keep internalenergy low we set the initial ratio of kinetic energy tototal mechanical energy to an arbitrary value of χ = 0 . χ asym . It was fitted as 1 . ± .
05 using theradio morphology from the high-resolution October 2008observations at 36 GHz (Potter et al. 2009), however itis likely that our fitted explosion asymmetry is an upperbound. In future work the asymmetry fit needs to berefined with more high resolution radio images. Usingthis solution a supernova explosion radius of 20 .
25 vox-els corresponds to a simulated time of around 850 daysafter explosion. A cross section of the initial supernovaenvelope can be seen in Figures 2, 3, and 4.
BSG wind region
A star emitting a spherically symmetric and steadywind with a mass loss rate ˙ M and a wind velocity v w produces a wind whose density profile scales with radius(r) as ρ ( r ) = ˙ M πv w r . (4)For our model BSG wind exterior to the progenitor Potter et al. 2014we use ˙ M = 7 . × − M ⊙ yr − , v w = 450 km s − as in Chevalier & Dwarkadas (1995). This is consis-tent with the density profile for the BSG wind de-rived from our environment formation simulation. FromLundqvist & Fransson (1991) there is evidence that theBSG wind was ionized by the shock breakout and at-tained a temperature in the range (3 . − . × K.For the wind we assume an isothermal temperature of5 . × K.Previous theoretical work has shown that the BSGwind ends in a termination shock located around(2 . − . × m (0 . ′′ − . ′′
5) from the centralstar (Blondin & Lundqvist 1993; Zhekov et al. 2010;Berezhko et al. 2011), (see also Appendix D). The EFSshows that the termination shock is a prolate spheroidwhose ratio of polar to equatorial minor axes is 1 .
18 withan average radius of 4 . × m (0 . ′′
56) from the twosemi-axes. We found that the turn on in radio emis-sion around day 1200 is sensitive to the location of thetermination shock. We used the same ratio of polarto equatorial axes as for the prolate sphere and foundthat an average radius of R ts = (3 . ± . × m(0 . ′′ ± . ′′
07) provided the best fit to the return of radioemission around day 1200 as seen in Figure 11.
Shocked BSG wind
Exterior to the termination shock, but still within thehourglass and H ii region is a bubble of shocked BSGwind. Our EFS indicates that the bubble is hot, with adensity of 1 . × − kg m − a temperature of around2 . × K and a velocity of 170 km s − on average.Outwards from the termination shock the outflowing ma-terial thins slightly and becomes supersonic again at aMach disk, which we take to be at a radial distance of1 . × m (1 . ′′
8) from the progenitor. Outside the Machdisk the gas properties are approximately constant untilthe edge of the expanding bubble is reached. We alsouse EFS results to fit polynomials to the pressure anddensity profiles of the shocked BSG wind. Details of thefitted polynomial and related constants are given in Ta-ble 3. The total mass for all BSG wind structures withinthe grid is 8 . × − M ⊙ . H ii region and hourglass Measurements of radial expansion show thataround day 1800 the shock slowed significantlyfrom 30 ,
000 km s − to 3 ,
000 km s − at a radiusof 4 . × m (0 . ′′ ii gas in the vicinity ofthe equatorial ring, swept there by either the expansionof the shocked BSG wind or evaporated from thering. As neither the literature nor the EFS have anyinformation on the morphology of the H ii region wemodel the innermost edge of the H ii gas by a convexcircular profile in the toroidal ( z ′ , s ′ ) plane and set it asthe inner edge of the waist of the hourglass. The curvedescribing the inner edge of the H ii region is constructedusing an arbitrary radius of 9 . × m (1 . ′′
33) withits origin placed in the equatorial plane and beyond the equatorial ring. The height of the waist (in | z ′ | ) is2 . × m (0 . ′′ ii region. The opti-mal fit to the radius observations places the inner edgeof our model H ii region at R H ii = (4 . ± . × m (0 . ′′ ± . ′′ . × m (0 . ′′
04) provided thebest fit to turnover in shock velocity.From (Sugerman et al. 2005) the hourglass is definedbetween a cylindrical radius (1 . − . × m (1 . ′′ − . ′′
73) and a maximum height of 2 . × m (3 . ′′ | z ′ | such that it asymptotes to the outer rimof the hourglass at large | z ′ | . The parameters of the ex-ponential were chosen such that it completes 3 e-foldinglengths between corners of the inner waist and the outerrim of the hourglass.We anticipate that the H ii region gradually mergeswith the density of the hourglass at large s ′ and | z ′ | , suchthat the boundary between the H ii region and hourglassis undefined. For the density and pressure profiles ofthe H ii region we use a truncated two-dimensional raisedGaussian in the ( s ′ , z ′ ) plane. Fits to the evolving super-nova shock radius place the peak density of the HII regionat its innermost edge. The FWHM of best fit in the s ′ di-rection was s ′ FWHM = (2 . ± . × m (0 . ′′ ± . ′′ θ ho is defined as θ ho = tan − (cid:16) z ′ FWHM R H ii (cid:17) . Fits to the observed radio mor-phology of the radio emission (Potter et al. 2009) showthat the best-fit half-opening angle is 15 ◦ ± ◦ , whichyields z ′ FWHM = (2 . ± . × m (0 . ′′ ± . ′′ . ± . × m − which is consistent withthe results from Zhekov et al. (2010).For the hourglass Sugerman et al. (2005) measured aHydrogen gas density of (2 − × m − out to acylindrical radius of s ′ = 1 . × m (2 . ′′
02) . Giventhe abundances in use we set the gas number densityof the hourglass to 8 . × m − and have the H ii region gas properties asymptote to this value at large s ′ and z ′ . Outside the hourglass in the radial direc-tion we model the equatorial belt and outer walls de-scribed in Sugerman et al. (2005) using a density pro-file that scales as s ′− for | z ′ | < . × m and as s ′− . elsewhere. The particle density of all the hour-glass structures is limited to a floor value of 3 . × m − , consistent with Sugerman et al. (2005). The upperand lower boundaries of the hourglass at | z ′ | = 2 . × m (3 . ′′
16) are smoothed using a logarithmic ramp func-tions in the pressure and density profiles, and a transi-tion width of 2 . × m (0 . ′′ . × − M ⊙ .The EFS shows that the temperature of the H ii regionis about 10 K, approximately a factor of 2 higher thanthe estimated 4 ,
500 K for the H ii region prior to theUV flash, and an order of magnitude less than the 10 K for post supernova models (Lundqvist 1999). For thepressure profile of the H ii region we adopt an isother-mal temperature of 8 . × K (Lundqvist 1999) for thehourglass, H ii region, equatorial belt and outer walls. Equatorial ring
Ionisation of the equatorial ring by the supernova UVflash has enabled accurate distance measurements to thesupernova (Panagia 2005). Assuming a circular ring, theradius r and width w of the ionized ring has also beendetermined as r er = (6 . ± . × m (0 . ′′ ± . ′′ w er = (9 . ± . × m (0 . ′′ ± . ′′
02) (Plait et al.1995; Sugerman et al. 2005). While the geometry of thenon-ionized portion is unknown, model fits to HST ra-dial profiles of the ring (Plait et al. 1995) suggest a cres-cent torus geometry for the ionized region. There is alsoconvincing evidence that over-dense clumps of materialreside within the ring and form hot spots of optical andradio emission when the shocks encounter it (Pun et al.2002; Sugerman et al. 2002; Ng et al. 2011).Spectroscopic optical and u − v line emission measure-ments of the equatorial ring between days 1400 and 5000show that the characteristic atomic number density forthe ionized gas varies in the range (1 × − × )atoms m − (1 . × − − × − kg m − ), giving a to-tal ionized mass of around 5 . × − M ⊙ and a ring tem-perature of around (3 − × K (Lundqvist & Fransson1991; Mattila et al. 2010).In similar fashion to Dewey et al. (2012) we use a two-component ring model consisting of a smooth equatorialring interspersed with dense clumps. The smooth ringbegins at s ′ er , inner = 5 . × m (0 . ′′ s ′ er = 6 . × , z ′ = 0) m (0 . ′′ , . ′′ s ′ inner , is a function of the height z ′ from the equa-torial (ring) plane. The width (in the z ′ direction) ofthe Gaussian is w er = 9 . × m (0 . ′′
12) at a heightof z ′ = w er /
2, and the Gaussian asymptotes to the inneredge of the hourglass s ′ = 1 . × m (1 . ′′
3) for large z ′ .As the filling factor of the ring is uncertain, we delineatethe outer edge of the smooth ring using the same Gaus-sian profile as for the inner edge but translated outwardsby a width w er , eff = 4 . × m (0 . ′′
06) in s ′ .Within the boundaries of the inner and outer edgesof the smooth ring we specified the pressure and den-sity profiles using another raised Gaussian function as afunction of the minimum distance from the ring locus at s ′ er , z ′ = 0. The FWHM was set to w er and the floor ofthe Gaussian was set to the density and temperature ofthe hourglass. The density and pressure were truncatedat the values of the surrounding H ii material to ensure asmooth transition. The peak number density and tem-perature of the smooth ring was set to 8 . × m − and2 . × K. At the innermost edge of the ring the numberdensity and temperature are 4 . × m − and 2 . × K. At the outermost edge of the ring in the equatorialplane, the density and temperature are at peak values.The total mass of the smooth ring is 6 × − M ⊙ . Withinthe smooth ring we place 20 dense clumps of material,centred on the ring at a radius of 6 . × m and evenly Table 1
Key fixed parametersDescription ParameterLength of the grid (m) 3 . × Inclination of the environment i x = 41 ◦ i y = − ◦ i z = − ◦ Ratio of specific heats γ = 5 / µ = 0 . r sn = 2 . × Index on BSG wind density profile s = 2Index on supernova envelopedensity profile n = 9Supernova energy (J) E tot = 1 . × Supernova envelope mass (kg) M env = 1 . × Ratio of kinetic to total energy χ = 0 . − ) ˙ M = 4 . × BSG wind velocity (m/s) v w = 4 . × Ratio of polar to equatorial distancesfor BSG wind termination shock 1.18Distance to Mach disk (m) 1 . × Radius describing inner profileof HII region (m) 9 . × Height (above equatorial plane)of inner profile of HII region (m) 2 . × Temperature of the HII regionand hourglass (K) 8 . × Hourglass number density m − . × Minimum backgroundnumber density (m − ) 3 . × Equatorial ring radius (m) r er = (6 . ± . × Equatorial ring width (m) w er = (9 . ± . × Equatorial ring number density,(smooth component) (m − ) 8 . × Equatorial ring temperature (K) 2 . × Equatorial ring clumppeak number density m − . × Equatorial ring clumppeak peak temperature (K) 2 × Total mass of ring clumps (kg) 7 . × distributed in azimuth. Each clump has a diameter of4 . × m, a peak density of 3 . × m − and apeak temperature of 2 . × K. For the density andpressure profile we choose the FWHM of the Gaussiansuch that at the periphery, the density of each clump is3 . × m − and has a temperature of 2 × K. Thetotal mass of the dense clumps is 3 . × − M ⊙ . Alongwith the mass of the smooth ring, this is consistent withthe 5 . × − M ⊙ currently estimated for the ionizedmaterial in the ring (Mattila et al. 2010). Summary of parameters
In Tables 1 and 2 and 3 is a summary of fixed andfitted parameters describing the environment of the fi-nal model. Error estimates on the fitted parameters arebased on the discretisation of parameter space used inthe model search.
Modelling radio and thermal emission processes
We assume a population of ultra-relativistic particlesis produced at both forward and reverse shocks via diffu-sive shock acceleration, where particles gain energy byrepeatedly sampling the converging flows at a strongshock front. Frequent scattering on magnetic fluctua-tions maintains a near isotropic distribution, ensuring Potter et al. 2014
Table 2
Fitted parametersDescription parameterSupernova envelope asymmetry χ asym = 1 . ± . R ts = (3 . ± . × Inner boundaryof HII region (m) R HII = (4 . ± . × Peak number densityof HII region (m − ) (7 . ± . × z ′ FWHM of HII region (m) z ′ FWHM = (2 . ± . × s ′ FWHM of HII region (m) s ′ FWHM = (2 . ± . × HII region half opening angle 15 ± ◦ that, on average, a particle will cross the shock manytimes before escaping downstream. In the absence ofnon-linear effects, this results in a uniform power-lawspectrum in momentum space f ( p ) = κp − b , where κ is a normalisation term. The distribution extends overseveral decades in energy. These ultra-relativistic elec-trons cool via synchrotron radiation, and the emissioncan typically be observed in the radio band. In our sim-ulations, we determine the synchrotron radio emission,by calculating the volume emissivity J ( ν ) and absorp-tion coefficient χ ( ν ) directly from the particle momen-tum distribution f ( p ) at shocked voxels. We assume arandomly-oriented magnetic field, and inject it at theshock using the analytic estimates for cosmic-ray drivenmagnetic field amplification. Then we follow its adia-batic evolution downstream. For the emissivity and ab-sorption we use the expressions for J ( ν ) and χ ( ν ) givenin Longair (1994). The synchrotron emissivity, in unitsof Watts m − Hz − sr − is J syn ( ν ) = √ e Bκc ( b − πǫ m e (cid:18) eB πνm e c (cid:19) b − a , (5)where a is given in terms of the Gamma function Γas a ( b ) = √ π (cid:0) b − + (cid:1) Γ (cid:0) b − − (cid:1) Γ (cid:0) b +34 (cid:1) ( b − (cid:0) b +54 (cid:1) . (6)The synchrotron absorption coefficient, in units of m − is χ syn ( ν ) = √ πe κc ( b − B b πǫ m e (cid:18) e πm e c (cid:19) b − a ν − b +22 , (7)where a is a ( b ) = Γ (cid:0) b − + (cid:1) Γ (cid:0) b − + (cid:1) Γ (cid:0) b +44 (cid:1) Γ (cid:0) b +64 (cid:1) . (8)The parameters κ and b are calculated at the shocklocation at each time-step, using the dynamically deter-mined shock-jump conditions. The magnetic field inten-sity B is estimated from the local shock parameters, us-ing the saturated magnetic field amplification estimatesand evolved downstream, together with the distributionof shocked particles. Details of the model are discussed in the following subsections. In order to produce an ef-fective comparison with monitoring observations we cal-culate synchrotron emission and absorption at frequen-cies 843 MHz and 1 .
38 GHz. Synchrotron cooling can besafely neglected, as the loss timescale for microwave emit-ting electrons is on the order of 10 years, much longerthan the dynamical timescale being studied here. Shock localisation
Within the diffusion approximation, the shape of thepower-law spectrum produced from shock accelerationdepends solely on the compression ratio of the shock,which can be determined from the shock velocity and theupstream plasma conditions. The ability to accuratelylocate shock positions in the hydrodynamical simulationis clearly a necessity. In the shock rest-frame, fluid ofdensity ρ , pressure P enters from upstream with ve-locity v and exits down-stream with ρ , P and velocity v . The compression ratio of a shock, ζ = ρ /ρ , can berelated to the pressure ratio P /P using the Rankine-Hugoniot shock relations (Landau & Lifshitz 1959). ζ = ( γ −
1) + ( γ + 1) P P ( γ + 1) + ( γ − P P (9)For γ = 5 / P /P >>
1, and thecompression ratio asymptotes to 4.In hydrodynamical simulations the shock is not a thindiscontinuity, but is spread over a region several cellswide. In order to locate shocks within the simulation wehave adapted the shock locator that FLASH 3.2 uses toswitch on a hybrid Riemann scheme in the presence ofa shock (Fryxell et al. 2000). It works by finding (viathe velocity divergence) voxels where fluid is being com-pressed . If the local pressure gradient is greater than athreshold value then a voxel is deemed to be in a shock.This is not sufficient to find points outside a shock, so thepressure gradient is followed upstream and downstreamuntil the gradient relaxes at points ρ , ρ , P , and P . De-tails of this technique are in Appendix B. Once the up-stream and down-stream variables have been determined,the shock compression ratio associated with a shockedvoxel is derived from Equation 9. As ζ = v /v from theshock relations, the inbound fluid velocity in the shockframe v , is obtained in terms of the lab frame velocities v L and v L v = (cid:12)(cid:12)(cid:12)(cid:12) ( v L − v L ) ζζ − (cid:12)(cid:12)(cid:12)(cid:12) . (10) Magnetic field amplification
We assume the magnetic field energy density upstreamof the shock is amplified via the Bell instability (Bell2004), which has been shown in numerical simulations toamplify fields by more than an order of magnitude. Thisis achieved through the stretching of magnetic field lines,driven by the cosmic-ray current j cr , which acceleratesthe background plasma via the j cr × B force. Hence, thefree energy available to amplify magnetic fields, is somefraction of the cosmic-ray energy density U cr .If µ is the permeability of free space (in SI units) thenBell (2004) relates the magnetic field energy density tothe cosmic ray energy density as B sat µ ≈ v c U cr , (11)We define an efficiency factor η cr U cr v = η cr ρ v . (12)such that B sat ≈ s µ η cr ρ v v v c = r µ η cr ρ v c . (13)In practice we found that v , as calculated from Equation10, is not very stable due to the finite width of the nu-merical shock. This consequently dampens the responseto changes in shock speed from abrupt changes in theupstream environment. We therefore adopt a more con-servative approach where v is approximated from thelab-frame shock velocity v ,L and the shock normal n by v ≈ v ,L · ˆ n ζζ − . The shock normal is derived fromthe pressure gradient, and the saturated magnetic fieldis approximated by B sat ≈ s µ η cr ρ (cid:20) ( v ,L · ˆ n ) ζζ − (cid:21) v c . (14)Supernova remnants are generally thought to be theprimary source of Galactic cosmic rays, which requiresan acceleration efficiency for protons and other heavynuclei of η cr ≈ . Acceleration of electrons at the shock
A precise treatment of diffusive shock acceleration overthe entire remnant is not possible, and we are forced touse a reduced model for the acceleration of electrons atthe shock front. The standard theory of shock accelera-tion predicts an acceleration time (Drury 1983) t acc = 3 v − v (cid:18) D v + D v (cid:19) (15)where D , are the shock-normal spatial diffusion coeffi-cients in the upstream and downstream regions. Thesecoefficients are typically taken to be Bohm-like, i.e. onthe order c / Ω g , where Ω g = eB/γm is the electronrelativistic gyro frequency. Given that the peak in thesynchrotron spectrum emitted from particles at a givenLorentz factor γ is ν syn ≈ π γ Ω g , (16)it follows that the characteristic acceleration time for ra-dio emitting electrons in the GHz range is shorter thanour numerical time steps. This allows us to update theelectron spectrum at every timestep in our simulations,such that a new spectrum is deposited at the shock loca-tion at each update. In the simplest theory of DiffusiveShock Acceleration (DSA), the power-law index of thedistribution, b is related to the compression ratio by b = 3 ζζ − . (17) The index b is related to the spectral index of radioemission α ( F ( ν ) ∝ ν − α ) by α = b − . For a strongshock in our ideal monatomic gas the compression ratiois 4 and the spectral index from shock acceleration is 0 . . ii region. It reached a peak of 1 . . t as t / (Kirk et al. 1996). The resulting index on the momen-tum is modified to b = 3 ζζ − (cid:18) ζ (cid:19) . (18)A tangled magnetic field is consistent with observationsas significant polarisation is yet to be observed in SNR1987A (Potter et al. 2009).The fraction of available electrons that were injectedinto the shock, χ el , while distinct from the accelerationefficiency, can be estimated with radio observations givenassumptions about the injection momentum. Assumingthe electrons are injected into the shock at a single mo-mentum δ ( p − p ), it can be shown (e.g (Melrose 2009)) that the isotropic downstream power law distributionof electrons, f ( p ) = κp − b , (where κ is a constant) is de-fined between p and the maximum momentum, whichis taken to be indefinite. If the downstream density ofenergised electrons is a fraction χ el of the electron num-ber density n , then conservation of mass requires that χ el n = R ∞ p πp f ( p ) dp . Solving for f ( p ) shows that f ( p ) = χ el n ( b − π p b − p − b . (19)and therefore κ = χ el n ( b − π p b − . (20)We assume electrons are injected into the shock fromdownstream. The injection momentum is derived by as-suming the electrons are in thermal equilibrium with thedownstream plasma and are injected into the shock at theelectron thermal velocity. By equating thermal energy torelativistic kinetic energy, then the injection momentumis given in terms of the temperature at the downstreampoint T p ( T ) = m e c s(cid:18) k b T m e c + 1 (cid:19) − . (21)0 Potter et al. 2014 Advection of the magnetic field
Once the magnetic field has been amplified by theshock, we assume that it is frozen into the backgroundflow, satisfying ddt (cid:18) B ρ (cid:19) = (cid:18) B ρ · ∇ (cid:19) u . (22)Following Kirk (1994), we assume an homologous expan-sion inside the remnant, i.e. u ∝ r ˆ r , for which the ratio ψ = B/ρ / is constant for a given fluid element.The method implemented to track ψ is given in Ap-pendix C. As the density is calculated in the main part ofthe hydro-code, the magnetic field can be reconstructedat a later time t simply by multiplying ψ by ρ ( t ) / . Advection of the particle distribution
In order to track the evolution of the electron dis-tribution in the downstream we follow the method ofDuffy et al. (1995), where the electrons are assumedfrozen to the flow (i.e. diffusion is neglected). This allowsus to simplify the transport equation ∂f∂t + u · ∇ f −
13 ( ∇ · u ) p ∂f∂p = 0 . (23)While this equation can in principle be solved usingthe method of characteristics, with dpdt = −
13 ( ∇ · u ) p. (24)to reduce the numerical effort, we choose instead to repli-cate the approach used for the magnetic field advection.Equation (23) can be re-written in the form ∂f∂t + ∇ · ( u f ) = ( ∇ · u ) f (cid:18) ∂ ln f∂ ln p (cid:19) . (25)where it is immediately noticed that ∂ ln f∂ ln p is just the in-dex of our power law − b . We update each component, i , of the two-point power law f ( p i ) using the advectionscheme in Appendix C. To preserve conservation of par-ticle number we update the injection momentum p byevolving it along the characteristic implied by Equation24. Radiative cooling and thermal emission
Unlike synchrotron emission, radiative cooling is imple-mented by converting a small fraction of the available in-ternal energy to thermal energy as the simulation evolves.We constructed a temperature-dependent cooling func-tion using MAPPINGS (Sutherland & Dopita 1993;Sutherland et al. 2003; Sutherland & Bicknell 2007).Figure 5 shows a plot of the cooling function multipliedby the particle mass squared ( µ ).If ǫ is the specific internal energy such that P = ( γ − ρǫ , and Λ( T ) is the cooling function, then the evolutionequation for the internal energy is given by dǫdt = ǫ ( γ − dρdt − ρ Λ( T ) . (26)As the dρ/dt term is already handled within FLASHthrough operator splitting, we complete the update to Temperature (K)10 -37 -36 -35 -34 Λ m p ( W a tt s m ) Figure 5.
The cooling function Λ( T ) multiplied by the square ofthe particle mass µ . The function was derived from the abundancesin Mattila et al. (2010) using the MAPPINGS shock and photoion-ization code (Sutherland & Bicknell 2007; Sutherland et al. 2003;Sutherland & Dopita 1993) the internal energy through the first order ODE dǫdt = − ρ Λ( T ) , (27)which is solved using a fourth order Runge Kuttascheme. The thermal X-ray emissivity J therm ( T ) frommaterial is then J therm ( T ) = 14 π ρ Λ( T ) . (28) Radiative transfer
Radiative transfer is implemented using the analyticsolution of the radiative transfer equation to propagatethe brightness across the grid in the direction of the ob-server. If ∆ x is the width of each voxel then the analyticsolution gives the brightness at the edge of each voxel interms of the volume emissivity J syn ( ν ) and absorptioncoefficient χ syn ( ν ) discussed in Section 2.3 I v = J syn ( ν ) χ syn ( ν ) [1 − exp ( − χ syn ( ν )∆ x )] . (29)After propagation across the grid, the flux density isobtained by multiplying by the apparent angular size ofthe voxel face as seen from Earth. RESULTS AND DISCUSSION
The hydrodynamical simulations were evolved to day10 ,
023 after the explosion, with an average timestep of 33simulated days. In post-processing, distributions of ac-celerated electrons were placed in the downstream flow ofthe forward and reverse shocks and were advected withthe flow. Radio emissivity at each timestep was gener-ated from the electron distributions. At a resolution of256 simulations took approximately 4 hours to completewith 8 cores. In Figures 6 and 7 are slices of the log ofdensity and pressure at a number of different simulatedepochs.The plots were formed by taking a cut plane at around45% of the X axis total length. In the top row of Figure1 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 1214 L o g nu m b e r d e n s i t y ( m − ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 1214 −15−14−13−12−11−10−9−8−7−6−5 L o g p r e ss u r e ( P a ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 1998 L o g nu m b e r d e n s i t y ( m − ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 1998 −15−14−13−12−11−10−9−8−7−6−5 L o g p r e ss u r e ( P a ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 5543 L o g nu m b e r d e n s i t y ( m − ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 5543 −15−14−13−12−11−10−9−8−7−6−5 L o g p r e ss u r e ( P a ) Figure 6.
Early epochs of the evolving shockwave from SN 1987A. Shown is the log of number density and pressure for a slice of thecomputational domain. Earth is to the right and the slice has been taken at 45% along the X axis in order to intersect one of the denseblobs in the ring. In the top row is the interaction of the supernova shock with the inner edge of the hot BSG wind around day 1200.Around day 2000 (middle row) the shocks begin to interact with the inner edge of the H ii region. By day 5500 (bottom row) the supernovashocks have begun interacting with the dense blobs within the ring. −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 6797 L o g nu m b e r d e n s i t y ( m − ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 6797 −15−14−13−12−11−10−9−8−7−6−5 L o g p r e ss u r e ( P a ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 8030 L o g nu m b e r d e n s i t y ( m − ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 8030 −15−14−13−12−11−10−9−8−7−6−5 L o g p r e ss u r e ( P a ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 10023 L o g nu m b e r d e n s i t y ( m − ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Day 10023 −15−14−13−12−11−10−9−8−7−6−5 L o g p r e ss u r e ( P a ) Figure 7.
Later epochs of the evolving shockwave from SN 1987A. Around day 6800 (top row), the supernova forward shock has almostcompleted its crossing of the ring. The reverse shock continues to interact with the highest density blobs within the ring. By day 8000(middle row), the forward shock has completely left the equatorial ring. The reverse shock continues to interact with the equatorial ringuntil after the end of the simulation at day 10 ,
023 (bottom row).
36 the supernova shock reaches the hot BSG wind aroundday 1200. After this encounter the shock splits into a for-ward shock, contact discontinuity, and a reverse shock.Around day 2000 (middle row) the forward shock encoun-ters the H ii region and around day 5500 (bottom row) thesupernova forward shock begins encountering the equa-torial ring. In the top row of Figure 7 the forward shockthe forward shock has almost completed its crossing ofthe ring around day 6800. The reverse shock is beginningto encounter the highest density blobs within the ring.By day 8000 (middle row), the forward shock has com-pletely left the dense ring. The reverse shock continuesto interact with the densest part of the ring until afterthe end of the simulation at day 10 ,
023 (bottom row).
Shock radius
The expanding radius from the simulation was calcu-lated from the 3D expanding morphology by deriving theradial distribution of radio luminosity at each timestep.The expectation value of the distribution forms an esti-mate of the radius. In this way we hope to determinethe measured shock radius in the most general way pos-sible that approximates model fits to the observationaldata from Staveley-Smith et al. (1992); Ng et al. (2008,2013). The resulting radius curves are shown in Figure8. The orange background in the plot is the time-varyingradial distribution of radio luminosity. The bin width forthe distribution has been normalised by its representa-tive width of 5 . × m. A sum over all dimensionlessbins in the distribution produces the total luminosity ofthe remnant at each timestep. We used radio emissionat a simulated frequency of 1 . E ( r ) for the distri-bution, along with upper and lower bounds containing68% of the radio luminosity. The radius and boundsare plotted both with and without polar emission abovea half-opening angle of 45 ◦ from the equatorial plane.Atop this is plotted the radius from u − v domain mod-els fitted to the observations (Staveley-Smith et al. 1992;Ng et al. 2008, 2013). We have also plotted a second-degree smoothing spline fitted to the radii from the u − v domain models. The spline fit at each point was weightedby the inverse of the error of the observed radius and thesensitivity of the spline was adjusted to produce a rea-sonably smooth function for the noisy data around day2000. The radial position and width of the H ii regionand equatorial ring are also delineated by horizontal lineswithin the plot. In order to compare the simulated ra-dius with the truncated shell model of Ng et al. (2013),we also fit a truncated shell to the simulation data at thesame epochs as the observations. The midpoint radiusand accompanying errors of the shell model are overlaidas blue diamonds. From the plot it is clear that the ra-dius from the truncated shell is systematically larger thanthe expectation of radius. It appears that the truncatedshell model fit to the simulation more closely follows theforward shock of the simulated data, however caution isadvised in applying the same interpretation for the trun-cated shell fit to the observations as it is still largelyunknown how the radio emission is distributed betweenthe real forward and reverse shocks. For the simulationwe have assumed that radio emission is generated at both forward and reverse shocks. This assumption may not beaccurate.Overall, the fitted radius from observations is well ap-proximated by the expectation of radius from the simu-lation. Prior to the collision with the H ii region aroundday 2000, the shape of the luminosity distribution is abroad and steep line, a clear signature of spherical expan-sion. Beyond the H ii region the time varying distributionof radio luminosity is clearly aspherical, as indicated bythe bi-modality in the distribution after day 2250. Theencounter of the forward and reverse supernova shockswith high latitude material above the plane of the ring isresponsible for the concentration of radio luminosity atradii greater than 8 . × m (1 . ′′ . × m(1 . ′′
77) produces a dramatic reduction in the productionof radio luminosity due to a lowering of shock velocity asthe shocks restart at that interface.At lower latitudes, radio luminosity is dominated bythe interaction of the supernova shocks with the equato-rial ring and H ii region. From the plot it appears thatthe forward shock began to encounter the ring aroundday 5400 and the reverse shock began to encounter thering around day 6200. Around day 7000 there is a dis-tinct turnover in radius for all estimates. This is morelikely to be the result of a change in the distribution ofradio emitting material than a real deceleration. Theapparent deceleration might be due to the forward shockleaving the equatorial ring. This seems likely to be true,as the distribution of radio emission, and therefore theexpectation of radius, is biased toward the reverse shockafter the forward shock leaves the ring.The Drishti (Limaye 2006) volume renderings in Fig-ure 9 confirm that the forward shock does indeed leavethe equatorial ring between days 7000-8000. Shown inthe figure is a volume rendering of the shock interactionwith a cross-section of two sides of the equatorial ring atday 7000 and 8000. Plotted in greyscale is the log of en-tropy, which is particularly sensitive to the reverse andforward shocks. Contrasted with this is the equatorialring and ring blobs, rendered as red and tan features.The figure shows that the forward shock has almost leftthe eastern equatorial ring by day 7000. By day 8000, theforward shock has completely left the eastern ring, andhas almost completed its crossing of the western ring.After day 8000 the radius appears to again accelerateas the relative amount of radio luminosity in the forwardshock begins to increase relative to the luminosity in thering. Shock velocity
In Figure 10 is the average shock velocity computedas the time derivative of the smoothed radius curves inFigure 8. Also shown is the velocity derived from thespline fit to the observational data, obtained by obtain-4 Potter et al. 2014 R a d i u s ( a r c s e c ) R a d i u s ( m )
68% of emission , with polar emissionE(r) with polar emission68% of emission , without polar emissionE(r) without polar emissionSpline fit to observational dataStaveley-Smith et al. (1992)Ng et al. (2013)Shell model fit to simulation L o g l u m i n o s i t y ( W a tt s H z − B i n − ) Figure 8.
The evolving radial distribution of radio luminosity from the final simulation. In orange is the radial distribution at eachtimestep. Overlaid is the expectation of radius E ( r ), and a boundary containing 68% of the luminosity. This is plotted for distributionsboth with and without high latitude luminosity above a half-opening angle of 45 ◦ . The observed radius from Staveley-Smith et al. (1992);Ng et al. (2008, 2013) and a spline fit to the observed data is also shown for comparison. Shown in blue diamonds is the shock radiusformed by fitting the truncated shell model of Ng et al. (2013) to the simulated radio emission. The vertical lines correspond approximatelyto epochs where the forward shock encountered various hydrodynamic structures; the BSG wind at day 1200, the HII region at day 2000,the equatorial ring ingress at day 5500, egress from the eastern and western lobes of the ring around days 6800 and 8000. Figure 9.
Volumetric renderings of the shock interaction with the equatorial ring at days 7042 (top row) and 8029 (bottom row). Inthe left column is a cross section of the easternmost part of the equatorial ring, on the right is a similar cross section of the western ring.Rendered in greyscale is the log of entropy to reveal the reverse and forward shocks (in that order) from the centre of the figure. In redand tan is the equatorial ring in the log of particle density. From the Figure it is clear that the forward shock has almost left the easternequatorial ring by day 7000. By day 8000 the forward shocks on both sides of the ring have almost completed traversing the ring, asindicated in Figure 7. ing the slope of the fitted spline at the observed epochs.From the plot we see that prior to day 2000 after ex-plosion, the average supernova shock expansion velocityof around (1 . − . × km s − is due to the for-ward shock propagating through the BSG wind. Afterthe shock encounters the H ii region, the average velocityis dramatically slowed to around 2300 km s − . Follow-ing this, the apparent shock velocity climbs steadily. Forradio luminosity within 45 ◦ of the equatorial plane, thevelocity reaches a peak around 6500 km s − around day5600. As the forward shock leaves portions of the thering, both estimates of radio emission experience a se-quence of rapid drops in average shock velocity betweendays 6700 and 8000. The double dip structure in theshock velocity may arise due to the forward shock leavingthe eastern lobe first, around day 7000, then the westernlobe at day 7800. These events cause a large change inthe average position of radio emitting material, and re-sult in a perceived reversal in the shock velocity. Afterday 8000, the average shock velocity appears to graduallyaccelerate as the forward shock encounters material witha greater sound speed. At day 8000 the average shockvelocity is around 1000 km s − . By day 10 ,
000 the aver-age shock speed has accelerated to around 6 ,
000 km s − . The velocity derived from a spline fit to the observationsappears to follow the general trend from the simulations,however caution is advised in interpreting high frequencyoscillations from the general trend, as the fit to the ex-panding shock radius is determined by the sensitivity ofthe least squares spline fit. The rate of increase in theaverage shock speed for the observations appears slowerafter day 8000, this may mean that the real sound speedin the ring and/or beyond the ring may be lower thanwe expected, indicating it may have a lower temperaturethan the temperature of 8 × K we had set for the ringand H ii region. Flux density and spectral index
We calculate flux density by summing over radial binsin the time-varying radio luminosity distribution in Fig-ure 8. The free parameter of the model, the fractionof available electrons swept up by the shock χ el , wasfitted to the observations by chi-square minimisationof the simulated flux density scaled by χ e l . We per-formed this fit for 843 MHz and 1.38 GHz. In Figure 11is the fitted flux densities plotted against the observedflux densities at the two frequencies. From the plot theshape of the simulated flux density provides a good fit6 Potter et al. 2014 −15−10−50510152025300 2000 4000 6000 8000 1000021521025051015202530 V e l o c i t ( k m / s ) Witho−t (olar emissionS(line fit to Stavele1-Smith et al. (1992) and Ng et al. (2013)
Figure 10.
Time-varying shock velocity computed from thesmooth radius curves in Figure 8. Also included is the veloc-ity derived from a spline fit to the to the observations fromStaveley-Smith et al. (1992) and Ng et al. (2013). Overlaid aremarkers showing the assorted interaction events from the forwardshock and its effect on the distribution of radio emission and hencethe derived shock velocity. The vertical lines are at the same epochsfor the hydrodynamical events as discussed in Figure 8. -1 F l u x d e n s i t y ( m J y )
843 MHz, with polar emission, η e =0.041.38 GHz, with polar emission, η e =0.04843 MHz, without polar emission, η e =0.041.38 GHz, without polar emission, η e =0.04Observed 843 MHz, Ball et al. (2001)Observed 1.38 GHz, Zanardo et al. (2010)Observed 1.38 GHz, Zanardo et al. (2013) Figure 11.
Simulated flux density plotted against observations at843 MHz (green) and 1.38 GHz (blue). The flux density both withand without polar emission is plotted for comparison. In black isthe observed fluxes at 843 MHz from Ball et al. (2001a) and at 1.4GHz from Zanardo et al. (2010) and Zanardo et al. (2013). Thevertical lines correspond to the same hydrodynamical epochs as inFigure 10. to the observational data at both frequencies betweendays 1200 and 5400. The scaling factor, χ el = 4%, pro-vides an optimal fit to the flux density, assuming theelectrons are in thermal equilibrium with the ions andare injected into the shock from the downstream regionat a momentum consistent with their thermal velocity.This fraction is consistent with the range of (1 − χ el = 6 × − obtained in Berezhko et al. (2011).The scaling factor is derived using the additional assump-tion that a constant fraction of the electrons are injectedinto the shock at all times. Since magnetic field ampli-fication is highly non-linear and is still an area of ac- tive research, these assumptions may not be correct andwe consider our derived value of χ el a preliminary re-sult. The sudden increase in flux density around day1200 is particularly sensitive to the location of the ter-mination shock in the relic BSG wind. In our simula-tion the termination shock was placed at distances inthe range (3 . − . × m (0 . ′′ − . ′′
51) from theprogenitor, with corresponding gas densities in the range(7 . − . × − kg m − . The flux densities around day1900 are discrepant with the observational data becausethe shock velocity (and hence the magnetic field) slowsconsiderably at the H ii region prior to restarting. This isprobably an artefact of a comparatively large numericalshock width, and might be resolved with an increase ingrid resolution in future studies. Another possibility isthat radius as reported from the truncated shell modelfits may be overestimated, as suggested by the truncatedshell model fits to the simulated data in Figure 8. Wenote that the spurious dip in flux density disappears insome of our models if we move the HII region and BSGtermination shock closer to the progenitor. Around day5500 there is an even greater discrepancy between the ob-served and simulated fluxes. It is interesting that the ob-served flux density does not also display a similar markedjump as the shock encounters the ring. There are manypossible reasons for this. The mass or filling factor of thesimulated ring may be overestimated. Alternatively, thevelocity of the simulated shock or the injection efficiencymay be overestimated during the crossing of the ring,thus more radio emission is produced than is observed.It may also be that the ambient magnetic field within theequatorial ring is lower than expected, hence the shockencounter with the ring is not producing as much syn-chrotron emission as the simulations predict. A spectralindex was calculated from the two frequencies. Howeverwe see little deviation from α = 0 .
75, which is expectedfor a strong shock and sub-diffusive shock accelerationwithout cosmic ray feedback.
Morphology and opening angle
The exact reason for the persistent asymmetry in theradio morphology of the remnant has been a longstand-ing mystery. Magnetic field amplification may providea solution to the problem by explaining the asymmet-ric radio morphology as a consequence of an asymmetricexplosion. From Equation 5, we see that synchrontronemissivity is a nonlinear function of b , B and κ . If weemploy shock acceleration to generate the particle dis-tribution f ( p ) and magnetic field amplification to obtain B , then radio emissivity should scale with shock velocity v s as J ν ∝ ρ ( b +3) / v b − / s ν − ( b − / . (30)From Landau & Lifshitz (1959) the shock velocity ofa strong forward shock propagating into a stationarymedium is proportional to downstream pressure P andupstream density ρ as ( P /ρ ) / . Radio emissivity thenscales as J ν ∝ P b − / ρ (9 − b ) / ν − ( b − / . (31)For a strong shock, b is in the range 4 − .
5, and radioemission is more sensitive to shock strength than den-sity. Conversely, thermal X-ray emission is more sen-sitive to density. Observations of thermal X-rays from7SN 1987A on day 7736 (Ng et al. 2009) show that theeast-west asymmetry is around 3 − Morphological comparisons with the observations
In order to test the magnetic field amplification hy-pothesis we derived 36 GHz synthetic images of the radiomorphology at day 7900 and compared them with obser-vations at the same epoch (Potter et al. 2009). Figure 12contains the result. At the top left is the observed image.At top right is the imaged model at day 7900. The fluxdensity of the model has been scaled to match that of theobservations. At bottom left is the imaged model wherethe model has been transformed to the u − v domain andthe u − v data of the transformed model is used to re-place corresponding u − v data of the observations. Theresult has been convolved with the (0 . ′′ × . ′′
2) beam fromthe October 2008 (day 7900) observation in Potter et al.(2009). At the lower right is the imaged residual, wherethe u − v data of the model has been subtracted from theobservation prior to imaging.From the plot it is clear that the simulated image andmodel show similar morphologies. Due to the faster east-ern shock, radio emission in the eastern lobe of the modelhas more radio emission than the western lobe. Theresidual images shows that extra radio emission from themodel occurs at the eastern lobe outside the position ofthe ring. This is unlikely to be due to image registrationas the difference in position between the dark patcheson the eastern lobe is greater than the registration errorof 0 . ′′
03. It may be that the speed of the eastern shockis faster than expected, or that the real forward shockis interacting with more high-latitude material than thesimulation indicates.
The evolving asymmetry
In order to track the evolution of asymmetry in thesimulation we obtained the ratio of the total integratedflux density either side of the origin in the rotated X ′ co-ordinate. The resulting evolution in asymmetry for thesimulation is shown as blue diamonds in Figure 13. Insimilar fashion we integrated the flux density over trun-cated shell models that were fitted to both the observedand simulated data. The results are shown as black andblue points.From the figure we see that the evolution in asymmetryprovides a reasonably good fit between the asymmetrymeasured from the simulation and the asymmetry ob-tained via a truncated shell model fit to the observationsin Ng et al. (2013). The truncated shell model fit to thesimulation appears to have a very high level of asymme-try. We suspect the truncated shell model is biased byhigh-latitude components of radio emission. Overall, theeastern lobe in both simulated and observed remnantshas consistently more flux than the western lobe fromday 2000 to day 7000. This shows that an asymmetric explosion combined with magnetic field amplification atthe shock is a viable physical model for reproducing theasymmetry in the remnant. The sudden positive jumpsin the simulated asymmetry in Figure 13 appear to becorrelated with hydrodynamical events such as the in-teraction with the H ii region around day 2000 and theencounter with the equatorial ring around days (5000-6000). Such behaviour indicates that the 3D model maynot be smooth enough.Both observed and simulated asymmetries experiencea decline around day 7000. This suggests that eithereastern lobe of the remnant loses a large portion of itsflux density relative to the western lobe at that time.Such a decline may be due to the forward shock exitingthe equatorial ring first, as is expected for an asymmetricexplosion. It may also indicate that the real shock hasencountered a significant overdensity in the western lobeof the ring. However the consequent X-ray emission froma shock encounter with such an overdensity has not beenobserved in X-ray images taken around the same timeNg et al. (2009).The rapid decline in the simulated asymmetry aroundday 7000, in contrast to that derived from observation, islikely to be the result of placing the ring at points equidis-tant from the progenitor. The decline in asymmetry fromfits to the observations is more gradual. This might indi-cate the ring is more broadly distributed in radius thanwe have simulated. The plot shows that the timing forthe simulated events is sooner than the observed eventsand that we may have overestimated the asymmetry inthe simulated explosion. Morphological predictions
An interesting prediction from the simulation is thatthe asymmetry will at least temporarily reverse directionin coming years, as evidenced by the asymmetry of thesimulation dipping below parity after day 8000. In Figure14 is synthetic images of the radio morphology betweendays 8700 and 9900. The images show that the westernlobe of the ring will dim more slowly than its easterncounterpart due to the lower shock speed.Measurements of the radius obtained through trun-cated shell modelling have shown that the radius curveshows an apparent deceleration at that time, (Ng et al.2013), suggesting that the shock has slowed down or therelative contribution of radio emission from the forwardshock has decreased. Fits to the radius obtained withring and torus models suggest that the radio emission isbecoming more ringlike with age. These observations areconsistent with the hypothesis put forward in Ng et al.(2013) - that day 7000 corresponds to the time the for-ward shock left the ring, leaving the reverse shock buriedwithin the ring.
Radio luminosity at different half-opening angles
In Figure 15 is the evolving distribution of radio emis-sion as a function of half-opening angle from the equa-torial plane. Overlaid on the distribution is the expec-tation of half-opening angle from the distribution, alongwith contours containing 68 and 95% of the radio emis-sion. For comparison, we have overlaid the half-openingangle derived from truncated shell model fits to boththe simulation (in blue) and the observations (in black)from Ng et al. (2013). At early times the supernova8 Potter et al. 2014
Figure 12.
Comparison of the real 36 GHz image with a corresponding synthetic model at day 7900. At top left is the 36 GHz imageat day 7900 from Potter et al. (2009). At top right is the model image formed from radio emission in the simulation. At lower left is themodel imaged with the same beam as the observations. At lower right is the residual image formed by subtracting synthetic u − v data ofthe model from the u − v data of the observation and imaging the result. Units are in Jy/Beam for the images and Jy/pixel for the model. F l u x d e n s i t y r a t i o Simulated asymmetry, X ′ axisTorus model, Ng et al. (2013) Figure 13.
The evolving asymmetry of the supernova model. Inblack is the asymmetry from the model fit in Ng et al. (2013).In blue is the the asymmetry obtained by summing flux densitieseither side of the origin in the rotated X ′ coordinate. shock is spherical, as evidenced by a half-opening an-gle around 90 ◦ seen before day 2000. After day 2000,the half-opening angle from the observations appears toseparate into two distributions representing componentsfrom high latitude material and the shock interaction with the H ii region. The expectation value of the simu-lated half-opening angle follows the radio emission nearthe H ii region and drops sharply, reaching a half-openingangle of 8 ◦ by day 4000. Between days 4000 and 6000,both simulated and fitted models show a flat slope for theevolving half-opening angle. The 95% boundary of thesimulated distribution appears to diverge from the fit-ted model after day 4000. This is due to radio emissionfrom high latitude material between days 4000 − − ◦ between days 8000 and 10000. Conversely,points from the truncated shell model fits to the simu-lation appear to be scattered around the expectation ofhalf-opening angle at early times, however they soon di-verge from the expected half-opening angle around day2000 and appear to follow the 95% confidence contourfrom the distribution, presumably as a result of high-latitude radio emission. The truncated shell failed toconverge to a solution after day 8000. It is suspectedthis is caused by hotspots in the in the simulated west-9 Figure 14.
Simulated 8.7 GHz radio images of SNR 1987A, between days 8000-9900. The model images (left column) were convolvedwith 0 . ′′ H a l f o p e n i n g a n g l e ( d e g r ee s )
68% of emission 95% of emission E(r)Torus model, Ng et al. (2013)Shell model fit to simulation L o g l u m i n o s i t y ( W a tt s H z − B i n − ) Figure 15.
The evolving distribution of radio emission as a function of half opening angle. Overlaid is the expectation of half openingangle and the contours containing 68 and 95% of the radio luminosity. Also plotted is the half-opening angle from the truncated shellmodel fit to synthetic images made from the simulation (in blue), and truncated shell model fits to observations (in black) from Ng et al.(2013). The vertical lines at days 1200, 2000, 5500, 6800 and 8000 delineate the shock interaction events discussed in Figure 10.
Injection parameters
As an independent consistency check to the semi-analytic injection physics of Section 2.3 we obtained theinjection parameters set for newly shocked cells in thesimulation. The injection parameters at each timestepwere obtained from a radio-weighted average of pa-rameters from cells shocked during the previous threetimesteps. We used the evolving radio emissivity at 843MHz as the weight for the average at each timestep. InFigure 16 is the result. At top left is the average num-ber density of injected electrons, scaled by the fractionof injected electrons required to reproduce the observa-tions. It is interesting to note two main jumps in theinjected electron density. The first is from the shock en-counter with the H ii region at day 2000 and second isthe encounter with the equatorial ring around day 5500.The injected electron density as the shock crosses theH ii region is around 10 m − . This is consistent withthe pre-supernova electron density of the H ii region, af-ter scaling by the injection efficiency and compressionratio. The higher electron injection density of 10 m − obtained at late times is consistent with the maximumscaled electron density of the ring. This suggests mostof the radio emission at these times is arising from com-paratively dense regions of the equatorial ring.At top right of Figure 16 is the average injection mo-mentum, in units of m e c . As seen in Equation 21 we de-rived the injection momentum from the gas temperatureat the downstream point. The maximum normalised in-jection momentum of 1 . ii region is equivalent to a shock temperatureof 4 × K. During propagation through the H ii regionthe injection momentum of 0.8 is equivalent to a temper-ature of 1 × K. During the shock crossing of the ring,the injection momentum drops to 0.4, or a temperatureof 4 × K.The injected magnetic field in the lower left panel ofFigure 16 shows that the amplified magnetic field is inthe range 8 × − − × − T. This is within an or-der of magnitude of the amplified magnetic field esti-mates in Duffy et al. (1995) and Berezhko et al. (2011),but is an order of magnitude larger than the estimate inBerezhko & Ksenofontov (2000).The magnetic fields form these other works are also in-cluded in the figure for comparison. We believe the dipof the injected magnetic field around day 2000 is due tothe lowering of shock velocity as the shocks crashed intothe H ii region. This is probably an effect of the low res-olution or an overestimated distance for the H II region.As the magnetic field is the only injection parameter toexperience a dip at day 2000, we are confident this is re-sponsible for the anomalous dip in seen around day 2000in the flux density plots of Figure 11. For completeness, the compression ratio ζ and the in-dex b obtained at the shocks is plotted in the bottomright panel of Figure 16. Overall the compression ratioreturned is fairly stable at the expected compression ra-tio of ζ = 4 for a strong shock, and an index of b = 4 . α = 0 .
75. A brief period of instability isobserved at early times when the shock was establishedfrom the initial conditions of the simulation.
Energy density at newly shocked cells
We also looked at the balance of energy density be-tween kinetic, thermal, and magnetic processes at newlyshocked points. Shown in Figure 17 is the evolving en-ergy density obtained by averaging in the same fashionas for the injection parameters. The magnetic energydensity is nearly two orders of magnitude below the ki-netic energy. This suggests that magnetic field amplifi-cation will have a negligible effect on shock evolution ifthe energy expended in amplifying a magnetic field is in-cluded in the energy budget. At early times prior to theencounter of the shock with the H ii region, the kineticenergy is clearly dominant. However this picture reversessoon after the encounter with the H ii region and thermalenergy at the shock front becomes dominant for the restof the simulation. CONCLUSIONS
In this work we have sought to: (1) Test magnetic fieldamplification and an asymmetric explosion as the causefor the long term asymmetry in the radio remnant; (2)Refine the structure of the pre-supernova environment,(3) Obtain an estimate of the injection efficiency at thesupernova shock, and (4) Provide a model that predictsfuture behaviour of the expanding radio remnant. Wehave addressed these questions by using a hydrodynami-cal simulation and a semi-analytic method incorporatingDiffusive Shock Acceleration and magnetic field ampli-fication to estimate power law distributions of electronmomenta and the magnetic field in the downstream re-gion of the shock. The distributions and magnetic fieldwere evolved with the downstream flow. Morphologicalcomparisons of the simulated radio emission with real ob-servations shows that magnetic field amplification com-bined with an asymmetric explosion is able to reproducethe persistent asymmetry seen in radio observations ofSN 1987A. The asymmetry in radio emission is primar-ily the result of non-linear dependence of the amplifiedmagnetic field on the shock velocity. The evolving radioemission from the simulation was compared to a num-ber of time-varying observations from SN 1987A, suchas: radius, flux density, morphology, opening angle, andspectral index. These comparisons formed the objectivefunctions for an inverse problem, and allowed us to refinethe model of the initial supernova environment.Essential features of the model are an asymmetric ex-plosion, a blue supergiant wind, an H ii region and anequatorial ring at the waist of an hourglass. We fixedthe energy and mass of the explosion at 1 . × Jand 10 solar masses. From the radius and flux densitycomparisons we find that a termination shock distance of(3 . − . × m (0 . ′′ − . ′′
51) provides a good fit forthe turn on of radio emission around day 1200. An H ii region with an innermost radius of (4 . ± . × χ e l n e m − Number density of energetic electrons ( p / m e c ) Injection momentum -9 -8 -7 -6 T e s l a Injected magnetic field
This workDuffy et al. (1995)Berezhko and Ksenofontov (2000)Berezhko et al. (2011) D i m e n s i o n l e ss r a t i o Compression ratio and momentum index
Figure 16.
Injection parameters from the final model. In the top left is the injected density of energetic electrons χ el n e ; top right is the in-jected momentum p ; bottom left is the injected magnetic field B with the magnetic fields from Duffy et al. (1995), Berezhko & Ksenofontov(2000) and Berezhko et al. (2011). At bottom right is the compression ratio ζ and index b on the isotropic momentum distribution. -11 -10 -9 -8 -7 -6 -5 J / m − Energy density
Downstream kineticDownstream magneticDownstream thermal
Figure 17.
Kinetic, thermal, electron, and magnetic field energydensity at the shock. The kinetic energy density is dominant beforethe interaction of the shock with the H ii region. m (0 . ′′ ± . ′′
01) and a maximum gas number densityof (7 . ± . × m − provides a good fit to theshock deceleration around day 2000 and subsequent ra-dius and flux density evolution to day 5500. The additionof clouds within the ring with a radius of 2 . × m,a peak number density of 3 . × m − and a totalmass 3 . × − M ⊙ results in an abrupt increase in theflux density around day 5500, given a constant injectionefficiency. It also results in a rapid reduction in open-ing angle and beading in the radio morphology after day7000.Three dimensional renderings of the computational do-main show that the period of apparent deceleration inshock velocity between days 7000 and 8000 may be theresult of the forward shock leaving the equatorial ring.The forward shock emerged from the eastern lobe of thering first around day 7000. It then emerged from thewestern lobe around day 8000. Following day 7000, theexit of the forward shock from the eastern lobe of equato-rial ring leaves strong radio-emitting components in thewestern lobe.The shock radii returned by truncated shell model fitsto the simulation appear to be substantially larger thanthe expectation of radius from the simulation, and ap-pears to follow the forward shock of the simulation. Aswe do not know how the radio emission of the real rem-nant is distributed between forward and reverse shockswe are unable to determine if truncated shell modellingis also similarly biased toward the forward shock of thereal remnant.Comparisons between simulated and observed flux den-sity show that during the supernova shock traversal of theHII region, the flux densities of simulation and observa-tion are in agreement if the fraction of electrons injectedinto the shock is around 4%. We arrive at this figure bymaking the somewhat speculative assumption that theelectrons are in thermal equilibrium with the ions andare injected into the shock from the downstream regionat a momentum consistent with their thermal velocity.There is a discrepancy between simulated and observedflux density around day 2000 due to a reduction in theamplified magnetic field caused by a stalled shock veloc-ity as the shock encountered the H ii region. This prob- lem might be rectified with higher resolution simulations.It may also mean that the radial distance of of the H ii region has been overestimated. The flux density is alsonot in agreement with observations from day 5500 on-wards as the shock encounters the thickest parts of thering. This may be due to the reasonably coarse resolu-tion of the simulation or lower χ el arising from yet to beunderstood microphysics at the shock as it collides withthe ring.As a result of the absence of cosmic ray feedback, thecompression ratio, and hence the index on the inversepower law for the electron distribution remains constant.This produces a spectral index for radio emission whichis inconsistent with the large dip seen in spectral index offrom the real remnant (Zanardo et al. 2010). We expectthat future models of the remnant that incorporate non-linear feedback (Lee et al. 2014; Ferrand et al. 2014) ormagnetic field topology (Bell et al. 2011) will be able toaddress this discrepancy.By capturing the injection parameters at the shockand performing a radio emission weighted average, wealso obtained estimates of the number density, momen-tum, magnetic field, compression ratio, and energy den-sity of the supernova shock. This permitted a consis-tency test of the semi-analytic method in use. The den-sity of electrons injected into the shock is consistent withthe upstream density (scaled by χ el ) of the medium intowhich the shock propagates. The injection momentumis consistent with a shock that has a temperature in therange 1 − × K. The injected magnetic field is in therange 8 × − − × − T, which is broadly consis-tent with Duffy et al. (1995) and Berezhko et al. (2011)but an order of magnitude higher than the estimate inBerezhko & Ksenofontov (2000).The ratio of energy densities at the shock clearly showsthat kinetic and thermal energy are approximately twoorders of magnitude stronger than magnetic energy den-sity. It is interesting to note that the downstream ther-mal energy occupies the largest fraction of the availableshock energy after the shock encounters the H ii regionaround day 2000.In terms of future predictions, the model indicates thatthe asymmetry in radio morphology may temporarilyreverse in coming years as radio emitting spots in thewestern lobe of the ring decrease in brightness moreslowly than their eastern counterparts. This is becausethe shock leaves the eastern ring more quickly. Syntheticimages of the future radio morphology indicate thatradio emission is concentrated in hotspots centred onoverdense blobs within the equatorial ring.We look forward to how this amazing young supernovaremnant evolves in years to come.We are grateful to John Kirk for his insightful and ex-tremely helpful input in all stages of this project. Weappreciate the support of staff at iVEC and ICRAR forproviding supercomputing resources for our use. In addi-tion we thank the Centre for Petroleum Geoscience andCO2 Sequestration for providing visualisation infrastruc-ture to produce volumetric renderings of our simulations.The software (FLASH) used in this work was in part de-veloped by the DOE - supported ASC/Alliances Center4 Potter et al. 2014for Astrophysical Thermonuclear Flashes at the Univer-sity of Chicago. REFERENCESAglietta, M., Badino, G., Bologna, G., Castagnoli, C., Castellina,A., Dadykin, V., Fulgione, W., Galeotti, P., Kalchukov, F., &Kortchaguin, B. 1987, Europhysics Letters, 3, 1315Alexeyev, E., Alexeyeva, L., Krivosheina, I., & Volchenko, V.1988, Physics Letters B, 205, 209Arnett, W. D., Bahcall, J. N., Kirshner, R. P., & Woosley, S. E.1989, IN: ARA&A, 27, 629Arnett, W. D. & Fu, A. 1989, ApJ, 340, 396Axford, W. I., Leer, E., & Skadron, G. 1977, 15th InternationalCosmic Ray Conference, 11, 132Ball, Lewis, Crawford, F, D., Hunstead, W, R., Klamer, I,McIntyre, & J, V. 2001a, ApJ, 549, 599Ball, L. & Kirk, J. 1992a, ApJ, Part 2 - Letters, 396, 1, 39Ball, L. & Kirk, J. G. 1992b, Astronomical Society of Australia,Proceedings, 10, 1, 38Bell, A. R. 1978, MNRAS, 182, 147Bell, A. R. 2004, MNRAS, 353, 550Bell, A. R., Schure, K. M., Reville, B. 2011, MNRAS, 418, 2,1208-1216Berezhko, E. G. & Ksenofontov, L. T. 2000, Astronomy Letters,26, 639Berezhko, E. G. & Ksenofontov, L. T.. 2006, ApJ, 650, 59Berezhko, E. G., Ksenofontov, L. T., & V¨olk, H. J. 2011, ApJ,732, 58Bethe, H. A. & Pizzochero, P. 1990, ApJ, 350, L33Bionta, R. M., Blewitt, G., Bratton, C. B., Caspere, D., & Ciocio,A. 1987, Physical Review Letters, 58, 1494Blandford, R. D. & Ostriker, J. P. 1978, ApJ, 221, 29Blondin, J. M. & Lundqvist, P. 1993, ApJ, 405, 337Caprioli, D. & Spitkovsky, A. 2014, ApJL, 765:L20Chevalier, R. A. 1976, ApJ, 207, 872Chevalier, R. A. 1982, ApJ, 259, 302Chevalier, R. A. & Dwarkadas, V. V. 1995, ApJL, 452, 45Crotts, A. & Heathcote, S. 1991, Nature, 350, 683Crotts, A. & Heathcote, S. 2000, ApJ, 528, 426Dewey, D., Dwarkadas, V. V., Haberl, F., Sturm, R., &Canizares, C. R. 2012, ApJ, 752, 103Dwarkadas, Vikram V., 2007, AIP Conference Proceedings, 937,120-124Drury, L. O. 1983, Reports on Progress in Physics, 46, 973Duffy, P., Ball, L., & Kirk, J. G. 1995, ApJ, 447, 364Ferrand, G. Danos, R. J., Shalchi, A., Safi-Harb, S., Edmon, P.,Mendygral, P. 2014, accepted for publication in ApJFryxell, B., Olson, K., Ricker, P., Timmes, F. X., Zingale, M.,Lamb, D. Q., MacNeice, P., Rosner, R., Truran, J. W., & Tufo,H. 2000, ApJS, 131, 273Gaensler, B. M., Manchester, R. N., Staveley-Smith, L., Tzioumis,A. K., Reynolds, J. E., & Kesteven, M. J. 1997, ApJ, 479, 845Gaensler, B. M., Staveley-Smith, L., Manchester, R. N.,Kesteven, M. J., Ball, L., & Tzioumis, A. K. 2007,SUPERNOVA 1987A: 20 YEARS AFTER: Supernovae andGamma-Ray Bursters. AIP Conference Proceedings, 937, 86Gouiffes, G., Wampler, E. J., Baade, D., & Wang, L.-F. 1989,The Messenger, 58, 11Hirata, K., Kajita, T., Koshiba, M., Nakahata, M., & Oyama, Y.1987, Physical Review Letters, 58, 1490Kirk, J. G. 1994 in Saas-Fee Advanced Course 24, PlasmaAstrophysics, Lecture Notes of the Swiss Society for Astronomyand Astrophysics (SSAA), ed. A. O. Benz and T. J.-L.Courvoisier, (Berlin, New York: Springer) 225Kirk, J. G., Duffy, P., & Ball, L. 1994, ApJS, 90, 807Kirk, J. G., Duffy, P., & Gallant, Y. A. 1996, A&A, 314, 1010Krymskii, G. F. 1977, Akademiia Nauk SSSR, 234, 1306Laki´cevi´c, M., Zanardo, G., van Loon, J. T., Staveley-Smith, L.,Potter, T., Ng, C.-Y., & Gaensler, B. M. 2012, A&A, 541, L2Landau, L. D. & Lifshitz, E. M. 1959, Fluid mechanics, 2nd ed.;Course of theoretical physics / by L. D. Landau and E. M.Lifshitz, Vol. 6 (Pergamon Press, Oxford)Lee, Shiu-Hang., Patnaude, D. J., Ellison, D. C., Nagtaki, S.,Slane, P. O. 2014, ApJ, 791, 2, 97 Limaye, A,
Drishti - Volume Exploration and Presentation Tool ,Poster presentation, Vis 2006, Baltimore.Longair, M. 1994, High Energy Astrophysics, 2nd ed., Vol. 2 (ThePitt Building, Trumpington Street, Cambridge, UnitedKingdom: The Press Syndicate of the University of Cambridge)Lundqvist, P. 1999, ApJ, 511, 389Lundqvist, P. & Fransson, C. 1991, ApJ, 380, 575Manchester, R. N., Gaensler, B. M., Wheaton, V. C.,Staveley-Smith, L., Tzioumis, A. K., Bizunok, N. S., Kesteven,M. J., Reynolds, J. E. 2002, PASA, 19, 207-221Manchester, R. N., Gaensler, B. M., Staveley-Smith, L.,Kesteven, M. J., & Tzioumis, A. K. 2005, ApJ, 628, 131Matsumoto, Y., Amano, T., & Hoshino, M. 2012, ApJ, 755, 109Mattila, S., Lundqvist, P., Gr¨oningsson, P., Meikle, P., Stathakis,R., Fransson, C., & Cannon, R. 2010, ApJ, 717, 1140McClements, K. G., Dieckmann, M. E., Ynnerman, A., Chapman,S. C., & Dendy, R. O. 2001, Physical Review Letters, 87, 25Melrose, D. B. 2009, eprint arXiv, 0902, 1803Ng, C.-Y., Gaensler, B. M., Murray, S. S., Slane, P. O., Park, S.,Staveley-Smith, L., Manchester, R. N., & Burrows, D. N. 2009,ApJ, 706, 100Ng, C.-Y., Gaensler, B. M., Staveley-Smith, L., Manchester,R. N., Kesteven, M. J., Ball, L., & Tzioumis, A. K. 2008, ApJ,684, 481Ng, C.-Y., Potter, T. M., Staveley-Smith, L., Tingay, S.,Gaensler, B. M., Phillips, C., Tzioumis, A. K., & Zanardo, G.2011, ApJL, 728, 15Ng, C.-Y., Zanardo, G., Potter, T. M., Staveley-Smith, L.,Gaensler, B. M., Manchester, R. N., & Tzioumis, A. K. 2013,ApJ, 777, 131Nomoto, K., Shigeyama, T., Kumaga, S., & Hashimoto, M.-A.1988, in Elizabeth and Frederick White Research Conference onSupernova 1987A, (Astronomical Society of Australia), 7, 4, 490Panagia, N. 2005, Cosmic Explosions. Springer Proceedings inPhysics, 585-592Plait, P. C., Lundqvist, P., Chevalier, R. A., & Kirshner, R. P.1995, ApJ, 439, 730Podsiadlowski, P. & Joss, P. C. 1989, Nature, 338, 401Podsiadlowski, P., Morris, T. S., & Ivanova, N. 2007,SUPERNOVA 1987A: 20 YEARS AFTER: Supernovae andGamma-Ray Bursters. AIP Conference Proceedings, 937, 125Potter, T. M., Staveley-Smith, L., Ng, C.-Y., Ball, L., Gaensler,B. M., Kesteven, M. J., Manchester, R. N., Tzioumis, A. K., &Zanardo, G. 2009, ApJ, 705, 261Pun, C. S. J., Michael, E., Zhekov, S. A., McCray, R., Garnavich,P. M., Challis, P. M., Kirshner, R. P., Baron, E., Branch, D.,Chevalier, R. A., Filippenko, A. V., Fransson, C., Leibundgut,B., Lundqvist, P., Panagia, N., Phillips, M. M., Schmidt, B.,Sonneborn, G., Suntzeff, N. B., Wang, L., & Wheeler, J. C.2002, ApJ, 572, 906Reville, B. & Bell, A. R. 2013, MNRAS, 430, 2873Reynolds, J. E., Jauncey, D. L., Staveley-Smith, L., Tzioumis,A. K., de Vegt, C., Zacharias, N., Perryman, M. A. C., vanLeeuwen, F., King, E. A., McCulloch, P. M., Russell, J. L.,Johnston, K. J., Hindsley, R., Malin, D. F., Argue, A. N.,Manchester, R. N., Kesteven, M. J., White, G. L., & Jones,P. A. 1995, A&A, 304, 116Riquelme, M. A. & Spitkovsky, A. 2009, ApJ, 694, 626Riquelme, M. A. & Spitkovsky, A. 2010, ApJ, 717, 1054Rousseau, J., Martin, N., Pr´evot, L., Rebeirot, E., Robin, A., &Brunet, J. P. 1978, A&A, 31, 243Shigeyama, T. & Nomoto, K. 1990, ApJ, 360, 242Soker, N. 1999, MNRAS, 303, 611Staveley-Smith, L., Briggs, D. S., Rowe, A. C. H., Manchester,R. N., Reynolds, J. E., Tzioumis, A. K., & Kesteven, M. J.1993, Nature, 366, 136Staveley-Smith, L., Gaensler, B. M., Manchester, R. N., Ball, L.,Kesteven, M. J., & Tzioumis, A. K. 2007, SUPERNOVA1987A: 20 YEARS AFTER: Supernovae and Gamma-RayBursters. AIP Conference Proceedings, 937, 96Got to hereStaveley-Smith, L., Manchester, R. N., Kesteven, M. J.,Reynolds, J. E., Tzioumis, A. K., Killeen, N. E. B., Jauncey,D. L., Campbell-Wilson, D., Crawford, D. F., & Turtle, A. J.1992, Nature, 355, 147Storey, M. C. & Manchester, R. N. 1987, Nature, 329, 421 Sugerman, B. E. K., Crotts, A. P. S., Kunkel, W. E., Heathcote,S. R., & Lawrence, S. S. 2005, ApJS, 159, 60Sugerman, B. E. K., Lawrence, S. S., Crotts, A. P. S., Bouchet,P., & Heathcote, S. R. 2002, ApJ, 572, 209Sutherland, R. S. & Bicknell, G. V. 2007, ApJS, 173, 37Sutherland, R. S., Bicknell, G. V., & Dopita, M. A. 2003, ApJ,591, 238Sutherland, R. S. & Dopita, M. A. 1993, ApJS, 88, 253Tanaka, T. & Washimi, H. 2002, Science, 296, 5566, 321Toro, E. F. 2009, Riemann Solvers and Numerical Methods forFluid Dynamics (Springer Berlin Heidelberg)Truelove, J. K. & McKee, C. F. 1999, ApJS, 120, 299Turtle, A. J., Campbell-Wilson, D., Bunton, J. D., Jauncey,D. L., & Kesteven, M. J. 1987, Nature, 327, 38V¨olk, H. J., Berezhko, E. G., & Ksenofontov, L. T. 2005, A&A,433, 229 Walborn, N. R., Prevot, M. L., Prevot, L., Wamsteker, W.,Gonzalez, R., Gilmozzi, R., & Fitzpatrick, E. L. 1989, A&A,219, 229Waßmann, M. & Kirk, J. G. 1991, Astronomische GesellschaftAbstract Series, 6, 86Woosley, S. E. 1988, ApJ, 330, 218Woosley, S. E., Pinto, P. A., Martin, P. G., & Weaver, T. A.1987, ApJ, 318, 664Zanardo, G., Staveley-Smith, L., Ball, L., Gaensler, B. M.,Kesteven, M. J., Manchester, R. N., Ng, C.-Y., Tzioumis,A. K., & Potter, T. M. 2010, ApJ, 710, 1515Zanardo, G., Staveley-Smith, L., Ng, C.-Y., Gaensler, B. M.,Potter, T. M., Manchester, R. N., & Tzioumis, A. K. 2013,ApJ, 767, 98Zhekov, S. A., Park, S., McCray, R., Racusin, J. L., & Burrows,D. N. 2010, MNRAS, 407, 2, 1157APPENDIX
PROGENITOR
From Truelove & McKee (1999) if an explosion were to propagate into a vacuum, it would expand with the velocity v ej and have radius R ej = v ej t . The velocity as a function of r and t is given by v ( r, t ) = (cid:26) rt , r < R ej , r > R ej . (A1)The density profile of the supernova exploding into a power law environment with density ρ ( r ) = ρ s r − s is given interms of a structure function f ( v/v ej ) = f ( w, n ), and the ejecta mass M env . We introduce asymmetry in the progenitorby multiplying the density in Truelove & McKee (1999) by the asymmetry function (1 + k sin θ cos φ ). ρ ( r, t ) = ( M env ( v ej t ) f (cid:16) vv ej (cid:17) (1 + k sin θ cos φ ) , r < R ej ρ s r − s , r > R ej . (A2)Where the constant k controls the degree of asymmetry in the progenitor. We chose s = 2 for the environmentsurrounding the progenitor as we assume a constant velocity wind. The density scaling ρ s is determined by the windvelocity v wind and progenitor mass loss rate ˙ M ρ s = ˙ M πv wind . (A3)The structure function f ( w, n ) specifies the shape of the solution given the exponent n f ( w, n ) = (cid:26) f n w − n core , ≤ w ≤ w core f n w − n , w core ≤ w ≤ . (A4)Chevalier & Dwarkadas (1995) used n = 9 as a best fit to the supernova. The purpose of a core is to avoid asingularity when n is large with w core = v core v ej as a free parameter. Truelove and McKee recommend small values( w core =0.001-0.1) thus setting a small core velocity. We have adopted w core = 0 . M env reveals f n as f n = 34 π (cid:20) − n − nw − n core (cid:21) . (A5)The kinetic energy of the explosion is determined by integration: E kin = 12 M ej v Z π dθ Z π dφ Z dww f ( w, n )(1 + k sin θ cos φ ) sin θ (A6)If the ratio of kinetic energies in the Eastern hemisphere to the Western hemisphere is χ ke then the constant k maybe obtained by taking the ratio of kinetic energy across the two hemispheres k = − (cid:18) χ ke − χ ke + 1 (cid:19) . (A7)Since E kin has been specified as a fraction χ of the total explosion energy, we can solve for v ej to obtain6 Potter et al. 2014 v ej = s E kin M ej (cid:18) − n − n (cid:19) (cid:18) − nw − n core − nw − n core (cid:19) . (A8)In order to incorporate pressure we assume it is related to density via an adiabatic process P = k P ρ ( r, t ) γ . Internalenergy is derived from pressure through the ideal gas equation of state: E int = (cid:18) M ej ( v ej t ) (cid:19) γ ( v ej t ) k P ( γ − Z π dθ Z π dφ Z dww f ( w, n ) γ (1 + k sin θ cos φ ) γ sin θ. (A9)Since the internal energy is (1 − χ ) E tot , the constant k P is given by k P = E int ( γ − f γn ( v ej t ) M ej t v ! − γ − nγ )3 − nγw − nγ core R π dθ R π dφ (1 + k sin θ cos φ ) γ sin θ . (A10)The position of the forward shock ( w b ) as a function of time is calculated from the differential equation of the shockmotion, dw b dt = − w b t (cid:18) M env v − s ej ρ s (cid:19) / l ed h f ( w b /l ed ) σ ed i / w s/ b t ( s − / . (A11)The constants σ ed = 0 .
212 and l ed = 1 .
19 are adopted for n = 9 from Table 6 of Truelove & McKee (1999).Regarding the position of the blast wave as a function of time, Truelove and McKee adopt the following characteristicvalues for position, time and mass R ch = M / (3 − s )env ρ − / (3 − s ) s (A12) t ch = E − / M [(5 − s ) / (2(3 − s ))]env ρ − / (3 − s ) s (A13) M ch = M env . (A14)Assuming an initial condition w b (0) = l ed , Truelove & McKee (1999) integrate Equation A11 to find the normalisedforward blast position ( R ⋆b ) as a function of time R ⋆b = R b R ch = (3 − s )2 (cid:18) l ed σ ed (cid:19) Z w b /l ed [ wf ( w )] dw ! s − . (A15)They then invert this solution for the normalised time taken to get to R b t ⋆ = E kin M env v ! R ⋆b l ed " − (cid:18) − n − s (cid:19) (cid:18) σ ed l ed f n (cid:19) R ⋆ − s b − − n . (A16)Therefore, given the energies E kin and E int , a progenitor mass M env , n, w core , and s , we can completely describe theearly expansion of the supernova. LOCALISING A SHOCK AND DETERMINING UPSTREAM AND DOWNSTREAM FLUID VARIABLES
The scheme FLASH 3.2 uses to locate voxels undergoing a shock is to look for compression as well as a significantpressure gradient. We look for compression by finding velocity divergence ( ∇ · v ) using central differencing. A negativevelocity divergence indicates compression. If the voxel is indeed in a shock then it should also have a pressure gradient ||∇ P || defined over the width of a numerical shock s w . If P is the upstream pressure and P is the downstreampressure then the pressure gradient ||∇ P || is ||∇ P || = (cid:18) P (cid:18) P P − (cid:19) / ( s w ) (cid:19) . (B1)Supposing we are looking for shocks with a compression ratio of at least r min . We can derive a minimum pressureratio r p,min from Equation 9 r p,min = ( γ + 1) r min − ( γ − γ + 1) − ( γ − r min . (B2)7Further supposing that the pressure P in a voxel is at least the upstream pressure of a shock, P , then the pressuregradient in a voxel should be larger than ||∇ P || > P ( r p,min − / ( s w ) . (B3)A shock is crossing a voxel if B3 is true and ( ∇ · v ) <
0. Using s w = 10 cells at maximum mesh refinement and r min = 2 . P and P by following the pressure gradient in the upstream and downstream directions until the slopefails to fulfil equation B3 or the pressure gradient turns back on itself.
310 315 320 325 330 335 3400.10.150.20.250.3 Voxel N o r m a li s ed p r e ss u r e Downstream point (P ) Upstream point (P ) −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5Distance (10 m)−1.5−1.0−0.50.00.51.01.5 D i s t a n c e ( m ) Log pressure at da : 8030 −15−14−13−12−11−10−9−8−7−6−5 L o g p r e ss u r e ( P a ) Figure 18. (Left) Points localised at the shock of a Sod shock tube problem. Crosses represent the localised points and the stars representthe upstream and downstream points P and P located by following the pressure gradient. (Right) The shock location algorithm appliedto the supernova simulation at day 8030. Shown in blue is a slice across the computational domain in the log of pressure variable. Overlaidis a contour plot (with 1 level) of voxels that were determined to be within a shock. Outwards from the centre, the inner and outer contoursmark the reverse and forward shock. ADVECTING A SCALAR VARIABLE
Given the advection equation for a scalar variable Ψ, a velocity field u , and constant κ∂ Ψ ∂t + ∇ · ( u Ψ) = ( ∇ · u )Ψ κ. (C1)The method of solution is similar to that (Toro 2009), pages 533-535. We first solve the associated homogeneousequation ∂ Ψ ∂t + ∇ · ( u Ψ) = 0 . (C2)Solving this equation gives a temporary solution Ψ ∗ which is the solution to d Ψ dt = 0. The full solution of equationC1 is completed by solving the following ODE in way following the prescription of (Toro 2009), pages 533-535 d Ψ dt = ( ∇ · u )Ψ ∗ κ. (C3)In order to solve equation C2 we integrate over the spatial-temporal cell dimensions to obtain the exact solution interms of fluxes entering each interface. Given at timestep ∆ t and grid spacing ∆ x the exact solution isΨ ∗ ( n +1) = Ψ n − ∆ t ∆ x (cid:16) F n +1 / x +1 / − F n +1 / x − / + F n +1 / y +1 / − F n +1 / y − / + F n +1 / z +1 / − F n +1 / z − / (cid:17) . (C4)To determine the fluxes we use the framework from page 457-461 of (Toro 2009). Using the velocity field v =[ v nx , v ny , v nz ] from the hydro simulation, the flux entering the cell from the left x direction is given byF n +1 / x − / = 12 (1 + σ x − / ) v nx,x − Ψ nx − + 12 (1 − σ x − / ) v nx,x Ψ nx , (C5)where σ x − / is a flux limiter function. We use the simple up-wind flux limiter. Given the average velocity acrossthe cell interface v n +1 / x,x − / = ( v nx,x − + v nx,x ) the flux limiter is defined as8 Potter et al. 2014 σ x − / = ( v n +1 / x,x − / > − v n +1 / x,x − / < = 0 . (C6)Once the solution to C2 has been approximated the full solution is obtained through the analytic solution to equationC3 Ψ n +1 = exp (cid:16) ( ∇ · u )Ψ ∗ ( n +1) κ ∆ t (cid:17) . (C7) PRE-SUPERNOVA ENVIRONMENT FORMATION SIMULATION
The beautiful hourglass structure of SN 1987A is thought to arise as a blue supergiant wind interacts with materialfrom past evolutionary phases of the progenitor. Simulations of remnant formation have more or less been able toreplicate the beautiful hourglass surrounding SN 1987A by placing a spherically-symmetric blue supergiant wind insidean asymmetric environment (Blondin & Lundqvist 1993; Soker 1999; Tanaka & Washimi 2002; Podsiadlowski et al.2007).Previous simulations of remnant formation of SN 1987A Blondin & Lundqvist (1993) have shown that a supersonicblue supergiant (BSG) wind extends radially outwards from the progenitor. The density profile of the wind scaleswith radius as r − since the flow is essentially a free-flowing wind. The free wind ends in a termination shock around(2 . − . × m (0 . ′′ − . ′′
5) from the progenitor. Material downstream from the termination shock is hot dueto adiabatic compression and forms the bubble responsible for inflating the hourglass.In order to obtain the density and temperature profiles of material in the free-wind and shocked-wind regions priorto the explosion, we simulated the formation of the pre-supernova environment in three dimensions. As with thesupernova simulation we used the standard hydrodynamics solver in FLASH with radiative cooling as discussed inSection 2.3.6. The computational domain was constructed as a rectangular grid of dimensions 256 × ×
640 at thefinest level of mesh refinement. This corresponds to a box of dimensions (2 . × . × . × m or (3 . ′′ × . ′′ × . ′′ M and radial wind velocity v w , the radial density profile of the free-wind is ρ ( r ) = ˙ M πv w r . (D1)Under the assumption of adiabatic flow, pressure is expected to scale with density as P ( r ) ∝ ρ ( r ) γ . In the centre of thegrid and at the origin we fixed a ”star” - a spherical region of radius r s = 9 . r s = 8 . × m) at the highestlevel of refinement. Within the star we set constant boundary conditions using equation D1 and ˙ M = 7 . × − M ⊙ yr − , 450 km s − from Chevalier & Dwarkadas (1995). Everywhere within the star we set a constant radial velocityof v w = 450 km s − .The pressure profile within the “star” was generated assuming adiabatic flow and a wind temperature of 16 ,
000 Kat the stellar surface where r = 3 . × m (Woosley 1988).For the the initial environment of the asymmetric RSG wind we used the wind profile from Blondin & Lundqvist(1993). If θ = sin − ( z ′ /r ) is the angle from the z ′ axis (see section 2.1 where these axes are defined), ˙ M RSG is the massloss rate of the red supergiant (RSG), v w, RSG is the RSG wind velocity, the density of the environment is described interms of the asymmetry parameter A: ρ ( r, θ ) = 3 ˙ M RSG v w, RSG πr (3 − A ) (1 − A cos θ ) . (D2)We used ˙ M RSG = 2 . × − M ⊙ yr − , v w, RSG = 5 km s − and A = 0 .
95 from the best-fit model ofBlondin & Lundqvist (1993). For A = 0 .
95, half of the RSG mass was lost within a half-opening angle of 21 ◦ from theequatorial plane and the equatorial-to-polar density ratio is 20:1. The pressure profile of the relic RSG wind was setby keeping the wind temperature constant at 500 K.The simulation was evolved until the waist of the bipolar inflated bubble matched the radius of the equatorial ringfrom the observations. For A = 0 .
95 this occurred around 19,725 simulated years from the initial conditions. Thisis consistent with other estimates of around 20,000 years for the time taken for the BSG to inflate the hourglass(Podsiadlowski et al. 2007). Figure 19 shows slices of the computational domain formed by cuts halfway along the x axis. Shown are the slices in different hydrodynamical variables overlaid by one-dimensional profiles. The horizontaland vertical profiles are represented by solid and dashed lines.Outwards from the star, a rarefied and fast, blue supergiant (BSG) wind extends to a termination shock located ata radial distance of (3 . − . × m (0 . ′′ − . ′′
55) along the polar axis, and (4 . − . × m (0 . ′′ − . ′′ − . Given the density and pressures of the environment this corresponds Mach9 −0.50.00.5 D i s t a n c e ( m ) N u m b e r d e s i t y ( m − ) (0.50.00.5 D i s t a c e ( m ) -21 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 P r e ss u r e ( P a ) −0.50.00.5 D i s t a n c e ( m ) -2 -1 T e m p e r a t u r e ( K ) −0.50.00.5 D i s t a n c e ( m ) V e l o c i t y ( m / s ) (2.5 (2.0 (1.5 (1.0 (0.5 0.0 0.5 1.0 1.5 2.0 2.5Dista ce (m) 1e16(0.50.00.5 D i s t a c e ( m ) M a c h N u m b e r Figure 19.
Slices of the formation simulation around 19,725 years after the initial conditions when the distance to waist of the hourglassapproximates the radius of the observed equatorial ring. The slices are formed by cuts halfway along the x axis. Shown are log-scaledimages in different hydrodynamical variables overlaid by one-dimensional profiles. The horizontal (polar) and vertical (equatorial) profilesare represented by solid and dashed lines. Note the highly supersonic BSG wind bubble in the centre, surrounded by a hot bubble ofshocked BSG wind at a temperature of around 10 K. Table 3
Parameters of the fit to hydrodynamical variables in the environment formation simulationRadial DistanceRegion (m) Parameter a b c
BSG wind 7 × − . × Density − . − . − . − . − . − . − . . . . − . × Density 0 . − . − . . − . − . . . − . . − . × Density 0 . . − . . . − . . . . numbers ranging from 19 at the BSG surface, 23 ,
000 at the inner boundary conditions of the “star” and 72 ,
000 justinside the polar termination shock, where the Mach number crosses unity.Exterior to the termination shock is a hot bubble comprising shocked BSG wind. Overall, the shocked BSG windhas approximately constant gas properties with a particle density of 1 . × m − and a temperature of 2 . × K.Due to the bipolar nature of the outflow, another shock known as a Mach disk forms at a distance (1 . − . × m (1 . ′′ − . ′′
81) along the polar axes.At at the expanding edge of the bubble, the hot BSG wind interacts with the relic RSG wind in two places. Theinner edge of the bubble is the interface between BSG and RSG winds. From the inner edge a forward disturbancepropagates outwards to become the outer edge of the expanding bubble. The “shocked” RSG material between theinner and outer edges of the bubble is associated with the H ii region from (Chevalier & Dwarkadas 1995). The H ii region in this simulation has a particle density in the range 10 − m − and a temperature around 10 K.In order to derive profiles for use in the supernova simulations polynomials were fitted to the density, pressure andvelocity profiles within the expanding bubble. If r is the radius of the BSG at 3 . × m and r ′ = rr then thepolynomial to be fitted is y = 10 a (log r ′ ) + b (log r ′ )+ c with the resulting units are in SI (units of density are in kgm − ). The coefficients of the fit are listed in Table 3.We compared these fits of the BSG wind profiles to theoretical estimates of the density, pressure and velocity profileson the assumption of ballistic flow. We find that the properties of the free BSG wind is in agreement with the adiabaticapproximation. The average deviation of the fits to the theoretical profiles are 6% in density, 9% in pressure and 0 ..