The Los Alamos Supernova Light Curve Project: Computational Methods

3 MB·

by  Lucille H. Frey

Draftversion January28,2013
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
THE LOS ALAMOS SUPERNOVA LIGHT-CURVE PROJECT: COMPUTATIONAL METHODS
Lucille H. Frey1,2, Wesley Even2, Daniel J. Whalen3, Chris L. Fryer4,5,6, Aimee L. Hungerford2, Christopher
J. Fontes7, and James Colgan8
Draft version January 28, 2013
ABSTRACT
We haveenteredthe eraof explosivetransientastronomy,in whichcurrentand upcoming real-time
3 surveys such as the LargeSynoptic Survey Telescope (LSST), the PalomarTransient Factory (PTF)
1 andPanoramicSurveyTelescopeandRapidResponseSystem(Pan-STARRS)willdetectsupernovae
0
in unprecedented numbers. Future telescopes such as the James Webb Space Telescope may discover
2
supernovaefromtheearlieststarsintheuniverseandrevealtheirmasses. Theobservationalsignatures

n of these astrophysical transients are the key to unveiling their central engines, the environments in
a whichtheyoccur,andtowhatprecisiontheywillpinpointcosmicaccelerationandthenatureofdark
J energy. We presentanew methodfor modeling supernovalightcurvesandspectrawith the radiation

5 hydrodynamicscodeRAGEcoupledwithdetailedmonochromaticopacitiesintheSPECTRUMcode.
2 Weincludeasuiteofteststhatdemonstratehowtheimprovedphysicsandopacitiesareindispensable
to modeling shock breakout and light curves when radiation and matter are tightly coupled.

] Subject headings: atomic data - methods: numerical - radiative transfer - supernovae: general
E
H
. 1. INTRODUCTION curves and spectra, with varying degrees of physics.
h
They include codes with radiation transport linked
p Wearenowintheeraofexplosivetransientastronomy,
to simple hydrodynamics such as EDDINGTON
- in which real-time all-sky surveys with the Large Syn-
o optic Survey Telescope (Ivezic et al. 2008), the Palomar (Eastman & Pinto 1993) and others (Hoeflich et al.
r TransientFactory(Law et al.2009) andPanoramicSur- 1993), one-temperature radiation diffusion codes such
t as Kepler (Weaver et al. 1978; Woosley et al. 2002),
s veyTelescopeandRapidResponseSystem(Kaiser et al.
a and radiation hydrodynamics (RHD) codes that
2002) are or soon will be detecting supernovae (SNe) in
[ evolve radiation and material temperatures sepa-
unprecedented numbers. Gamma-ray bursts now probe
rately such as VISPHOT (Ensman & Burrows 1992),
3 the universe out to z & 8 and upcoming telescopes
TITAN (Gehmeyr & Mihalas 1994) and STELLA
v suchastheJamesWebbSpaceTelescope(Gardner et al.
(Blinnikov et al. 1998). There are many codes which
2 2006)coulddiscoverSNe fromthe earlieststars(PopIII
have been developed to post-process profiles from
3 SNe) and reveal their masses. The observational signa-
hydrodynamical codes with radiation transport, in-
8 turesofthesecosmicexplosionsarekeytounveilingtheir
5 central engines, the environments in which they occur, cluding SEDONA (Kasen & Woosley 2009), SYNOW,
3. andto whatprecisionthey canbe usedas standardcan- (Branch et al. 1985) and several Monte Carlo codes
(Mazzali & Lucy1993, Kromer & Sim2009andothers).
0 dlestopinpointcosmicexpansionratesandconstrainthe
RHD simulations have also been post-processed to ob-
2 natureofdarkenergy. Lightcurvesonlyprovideindirect
tainnon-LTE(localthermodynamicequilibrium)spectra
1 measures of the physics at work in transients but are by
(Hauschildt & Ensman 1994).
: far the most abundant type of data. New advances in
v In this paper we present a new method developed at
numerical codes can now produce SN light curves and
i Los Alamos National Laboratory (LANL) in which we
X spectra with greater accuracy and time resolution than
start from a full RHD simulation and post-process the
in previous studies, mandating a new generation of sim-
r resulting structure with detailed monochromatic opac-
a ulations to better understand these cosmic explosions.
ities to calculate SN light curves and spectra with an
Numerous codes have been developed in the past
unprecedented level of accuracy. This level of detail is
three decades to simulate SNe and their light
necessary for times when radiation-matter coupling is
important, from shock breakout until the start of the
1Department of Computer Science, University of New Mex- nebularphasewhentheSNejectabecomeopticallythin.
ico,Albuquerque,NM87131,USA ModelswithfullRHDarerequiredwhereradiationplays
2XTD-6,Los Alamos National Laboratory, Los Alamos, NM
an important role in the hydrodynamics, such as where
87545, USA
3Department of Physics, Carnegie Mellon University, Pitts- shock heating is important in determining the tempera-
burgh,PA15213, USA tureprofile. Thecodespresentedheredonotincorporate
4CCS-2, Los Alamos National Laboratory, Los Alamos, NM non-LTE effects, which become more significant at late
87545, USA
5Physics Department, University of Arizona, Tucson, AZ times in SNe.
85721, USA Inthe initialstagesofsupernovae,reboundfromgrav-
6Physics and Astronomy Department, University of New itational collapse or thermonuclear burning drives a
Me7xXicCo,PA-5l,bLuqouseArqlaume,oNsMNa8t7io1n3a1l,LUaSbAoratory, Los Alamos, NM highly radiative shock into the upper layers of the star.
The shock is not visible until it breaks through the sur-
87545, USA
8T-1, Los Alamos National Laboratory, Los Alamos, NM faceofthestarbecauseopacityduetoelectronscattering
87545, USA in the intervening layers traps photons in the shock and
2
they are advected outwardby the flow. When it reaches core-collapse and Type Ia explosions and for 2 – 3 years
the surface of the star, the shock abruptly accelerates in pair-instability explosions. The fireball cools as it ex-
in the steep density gradient there, becoming even hot- pands and the peak of its spectrum evolves to longer
ter and releasing a sharp burst of photons into the sur- wavelengths. At the same time, the envelope surround-
rounding medium. This results in an extremely bright ing the star that was ionized by the breakout pulse be-
breakout transient with peak bolometric luminosities & gins to recombine and blanket the spectrum with lines,
1045 erg s−1, chiefly in the form of x-rays, hard UV and particularly at shorter wavelengths. In the frame of the
occasionallygammarays. Theseradiativelossesalterthe shock, the photosphere from which photons escape de-
hydrodynamics of the system and so cannot be ignored. scends deeper into the ejecta as it expands. The heavy
The ambient medium into which the shock breaks elementsitexposestothesurroundingmediumovertime
out can strongly affect light curves and spectra in ways depend on how Rayleigh-Taylor instabilities and asym-
that may have been underestimated in past studies. metries mix the interior of the star as the shock prop-
The luminosity and hardness of the breakout transient agates to its surface. When there is little mixing, lines
are governed in part by how much the shock acceler- duetoheavierelementsdonotappearinthespectraun-
ates as it exits the star: the steeper the gradient, the til later times, but when mixing is extensive such lines
brighter and harder the pulse. Radiation breakout can can become manifest at early times. The order in which
occur after shock breakout in dense circumstellar en- metallinesappearinSNspectraovertimecanbeapow-
velopes because photons can be trapped by the wind af- erful probe of mixing during the explosion, and of the
ter the shock exits the surface of the star (Ofek et al. central engine itself.
2010; Balberg & Loeb 2011; Chevalier & Irwin 2011, Bolometric luminosities can surge at late times if
L.H. Frey et al., in preparation). Shock breakout the ejecta is dense and delays the escape of IR pho-
has been the subject of numerous analytical studies tons or if the photosphere moves into a warm re-
(Colgate 1974; Matzner & McKee 1999; Nakar & Sari gion such as a 56Ni layer. Dense winds can also im-
2010; Piro et al. 2010; Katz et al. 2012) and numerical pact spectra at late times by more thoroughly blan-
simulations (Ensman & Burrows 1992; Blinnikov et al. keting the spectrum with lines, particularly at short
2000; Tominaga et al. 2009; Fryer et al. 2009; Tolstov wavelengths. The light curve can also brighten if the
2010; Fryer et al. 2010; Kasen et al. 2011) in the three ejecta collides with a dense shell from a violent lumi-
decades since it was first predicted (Colgate 1974; nous blue variable-type mass ejection (as in some Type
Klein & Chevalier 1978). However, there have only IIn SNe). Such collisions can change even dim SNe
recently been observations which are interpreted by intoextremelybrightevents(e.g.Smith & McCray2007;
some as shock breakout. These include direct ob- van Marle et al. 2010; Moriya et al. 2010; Roming et al.
servations of X-ray pulses in GRB 060218/SN2006aj 2012). After breakout, a radiative precursor, which
(Campana et al. 2006) and XRO 080109/SN2008D would not be seen in hydro-only codes, propagates in
(Soderberg et al. 2008), early-time peaks in UV light front of the shock and interacts with any surrounding
curves (Schawinski et al. 2008; Gezari et al. 2008), and material. If a shell was previously ejected, this precur-
indirectobservationsthroughIRechoesintheCasASN sorimpacts the shellbeforethe mainshock andcausesa
remnant (Dwek & Arendt 2008). bump in the light curve (D.J. Whalen et al., in prepara-
Thelocationandtime atwhichradiationbreakoutoc- tion). During the nebular phase, the influence of shock
curs depend on the radius at which the star or its envi- heating and radiative losses on the hydrodynamics are
ronmentbecomesopticallythinandisafunctionofwave- minimal and hydro-only codes can be used for accurate
length. Standard blackbody and light-crossing time ar- simulations. As departures from LTE become more sig-
guments for peak breakout luminosities and pulse width nificant during the nebular phase, the accuracy of the
assumeagrayopacityandnoattenuationbeyondtheτ = codes and opacities presented here will decrease.
1surface. Inreality,theopacitycanbeastrongfunction WedescribetheLANLRHDcodeRAGEandhowitis
of wavelength in SN ejecta and is poorly approximated usedto model supernovaexplosionsin Section 2. We re-
by a gray opacity. Analyses of shock breakout obser- viewtheLANLSPECTRUMcode,whichpost-processes
vations (i.e. Waxman et al. 2007; Chevalier & Fransson RAGE profiles to compute time-dependent spectra and
2008;Gezari et al.2008;Soderberg et al.2008)havecon- light curves in Section 3, and discuss how atomic opac-
sidereda variety of progenitorsandstellar environments ities derived from the LANL OPLIB database are im-
but treat breakout as a wavelength-independent event plemented in both codes in Section 4. In Section 5 we
thatoccursatasingleradius,temperatureandtime. De- presentthesuiteoftestsusedtoverifySPECTRUM.We
pendingonthewavelengthband,thisassumptioncanin- analyzeaseriesoftestsinwhichmorephysicalprocesses
troduce significanterrorsto estimates of explosionprop- are progressively activated in RAGE in order to assess
erties derived from observations. theireffects onspectraandlightcurvesinSection6. We
The breakout transient is generally followed by a rise show why full RHD and accurate opacities are so im-
and then gradual decline in bolometric luminosity that portanttorealisticsupernovalightcurvesandspectrain
is powered by both the conversionof kinetic energy into Section7,andinSection8wesummarizeourconclusions
thermal energy and the decay of 56Ni. Shock heating and discuss future applications.
fromthesupernovaejectainteractingwiththesurround-
ingwindsalsocontributestothepeakluminosityinType 2. RAGE
II SNe, and to a lesser extent in Ib/c SNe. Shocks in RAGE (Radiation Adaptive Grid Eulerian) is a mul-
TypeIaSNsurroundingscanalsocontributetothelight tidimensional, multispecies Eulerian adaptive mesh re-
curves near peak luminosity. This luminosity powered finement RHD code developed at LANL (Gittings et al.
by shocks and 56Ni decay persists for several months in 2008). Itcouplessecond-orderconservativeGodunovhy-
3
drodynamics to either gray or multi-group flux-limited thermal and radiation equilibrium with their surround-
diffusion radiation transport. RAGE models radi- ings. When this happens, the true opacity can deviate
ation hydroflows on one-dimensional (1D) Cartesian, fromthoseinourtables. ThetestswedescribeinSection
cylindrical and spherical polar coordinate meshes, two- 6 show that LTE is a good approximation during shock
dimensional (2D) Cartesian and cylindrical grids, and breakout. We alsoinclude anadvectionterminthe gray
three-dimensional (3D) Cartesianboxes. The mass frac- diffusion calculation. This advection term conserves en-
tion of each species is evolved with its own continuity ergy(consistentlyincorporatestheworkdefinedbypres-
equation under the assumption that all species have the sure and change in volume, P*dV) and includes mate-
same velocity at a given mesh point. rial motion effects (assuming the radiation flows with
the matter). Although this is not correct in the free-
2.1. Microphysics
streaminglimit,thiserrorislikelytobesmallerthanthe
AlthoughRAGE hasthree-temperature(electron,ion, fact that we are using flux-limited diffusion to approxi-
and radiation) capability, we use two-temperature (2T) mate the transition from diffusion to free-streaming.
physics in our simulations, in which radiation and mat-
ter, while coupled, can be at different temperatures. 2.3. Gravity
This can occur during shock breakout, when matter
Gravitycanplayanimportantroleincore-collapsesu-
is known to be out of thermal equilibrium with pho- pernovae at very early times, causing ejecta to fall back
tons (Nakar & Sari 2010). This 2T physics is an im- onto the newly formed compact remnant. In most su-
portant improvement over the pure hydrodynamics and
pernovae,thebulkofthefallbackoccursinthefirst1000
one-temperature radiation models used in previous SN
s of the explosion (Fryer et al. 1999, 2006a; Zhang et al.
light-curve calculations. However, deviations between 2008; Fryer 2009). In our calculations, we first evolve
the matter and radiation temperatures can imply de-
SNe until after fallback and nuclear burning is complete
partures from local thermodynamic equilibrium (LTE)
incodessuchasKepler,CASTRO(Almgren et al.2010),
and lead to inaccuracies in our opacities during break- SNSPH (Fryer et al. 2006b) and others (Fryer 1999).
out, as we discuss below. We do not explicitly transport
Thisisdoneinordertodeterminenucleosyntheticyields
gamma rays from the radioactive decay of 56Ni in the
andhowmuchmaterialfallsbackontothecompactrem-
ejecta. Instead, their energy is deposited locally in the nant so this material can be removed from the profiles
gas according to
prior to mapping them into RAGE. After fallback and
burning is complete, gravity no longer plays a large role
dE E E
= Nie−t/τNi+ Co (e−t/τCo e−t/τNi), (1) in the explosionand we neglect it in RAGE.As a conse-
dt τNi τCo τNi − quence,thestarcanslightlyexpandintothesurrounding
−
medium in the time it takes the shock to reach the sur-
where τ = 7.6 105 s, τ = 9.6 106 s and the mean
Ni Co face,butthis motionisminor. Ifweregridthe explosion
energiesreleased×peratomforthed×ecayof56Niand56Co
prior to breakout, the original profile of the star is re-
are E = 1.7 MeV and E = 3.67 MeV (Fryer et al.
Ni Co stored to the grid.
2009). This approximation is valid at early times when
theejectaisdenseandopticallythickbutcanbreakdown 2.4. Supernova Profiles
at later times when gamma-rays are no longer trapped.
Gas densities, velocities, specific internal energies (erg
This occurs after about50 days inType Ia SNe and,de-
pending on the degree of mixing, after 150 days in core- g−1), radiation energy densities (erg cm−3) and species
massfractionsarerequiredtoinitializeanexplosionpro-
collapse SNe. The energy deposition described above
file in RAGE. Gas densities, mass fractions and veloc-
has been tested by comparing it to non-thermal γ-ray
transport,whichproducesthe sameresultatearlytimes ities are directedly imported from codes such as those
listed above. However, rather than converting the gas
(Fryer et al. 2009). We apply the same constant specific
energies in those codes into the specific internal energies
heat C to each constituent material, so only ions con-
v
tribute to the total pressure. We also use an ideal gas required by RAGE, we calculate gas and radiation en-
ergies from the temperature of the gas, assumed to be
equationofstate (EOS)because the densities commonly
ideal and monatomic, which is also usually supplied by
foundinourastrophysicalapplicationsareoffofexisting
EOS tables in RAGE by many orders of magnitude. the codes:
3
2.2. Radiation Transport Egas = RT (2)
2
Both gray and multi-group flux-limited diffusion are
implemented in RAGE with opacities derived from the and
LANL OPLIB database using the TOPS code, as ex- E =aT4, (3)
rad
plained in Section 4. Rosseland mean opacities are used
toevaluatethediffusioncoefficientinthetransportequa- where R is the ideal gas constant and a is the radiation
tionsandPlanckmeanopacitiesareusedtocalculatethe constant. This avoids subtleties in how gas energies are
emissionandabsorptionterms. Allthetestsinourstudy defined in various codes, some of which lump contribu-
weredonewithgrayflux-limiteddiffusion. OPLIBopac- tions due to the ionizationstates ofatoms and radiation
ities are calculated under the assumption that the gas is in with the internal energy. In contrast to internal en-
inLTEandarethus functions offrequency,temperature ergy, the temperature is unambiguously defined in most
and density only. At late times, when the SN becomes codes.
transparent, LTE may no longer be valid because indi- We map the profile, which consists of the explosion,
vidual fluid elements become too diffuse to remain in the surrounding star and the circumstellar environment,
4
onto a 1D spherical grid in RAGE. We choose the grid The 1D RAGE data, sampled as described above, is
resolution so that the blast and its photosphere are well mapped onto a 2D grid of radial and angular bins in
resolved, and we permit up to five levels of refinement SPECTRUM as shown in Figure 1, where the angular
in both the initial interpolation of the profile onto the bins are uniform in µ = cosθ. SPECTRUM takes the
grid and during the simulation. We set reflecting and fluid at each mesh point to be in local thermodynamic
outflowboundaryconditionsonfluidandradiationflows equilibrium(LTE)withitssurroundingsanditsradiation
at the inner and outer boundaries of the mesh, respec- spectrum to be Planckian. In this approximation, each
tively. Whenthecalculationisbegun,Couranttimesare fluid element emits radiationat the same rate it absorbs
initially small due to high temperatures, large velocities itfromits environment,sothe Kirchhoff-Planckrelation
and small cell sizes. To reduce execution times and ac- holds:
commodate the expansion of the ejecta, we periodically
η = κ I =κ B , (4)
regrid the profiles onto a larger mesh as the explosion abs abs ν
grows. Thissignificantlyincreasesthetimesteponwhich
where η is the emissivity of the fluid parcel, κ is its
abs
the model evolves. We remap just the explosion itself,
monochromatic absorption opacity at frequency ν and
ignoring any medium beyond the shock or its radiation
front, and then graft the original environment lying be- 2hν3 1
yondthis radiusbackontothe shockorfrontonthe new B = . (5)
grid. We apply the same criteria in choosing a new grid ν c2 ehν/kTgas 1
−
as in the original setup. The entire calculation typically
From Equations 4 and 5 we construct the luminosity of
requires 20,000 CPU hours on LANL platforms.
each fluid element into a given line of sight:
3. SPECTRUM
2hν3 1 (1 vµ/c)2
We post-process RAGE hydrodynamical profiles with Lr,µ =κ m − . (6)
the LANL SPECTRUM code to construct light curves ν abs c2 ehν/kTgas 1p(1 (v/c)2)
− −
and spectra for our simulations. Stand-alone post-
processingallowsus toexploitthe comprehensiveLANL Here, r and µ define the radial and angular position of
OPLIBopacitydatabasetocomputedetailedabsorption
the fluid parcel, m,v and T are its mass, velocity and
gas
and emission features in our spectra without resorting
temperature, ν is the frequency at which it radiates in
to thousands of energy groups in the transport calcu- its own frame, and h,c and k are Planck’s constant, the
lations. These monochromatic opacities are calculated speed of light and Boltzmann’s constant, respectively.
assuming LTE conditions for a range of temperatures, Notethatthesphericalsymmetryoftheproblemimplies
densities andmaterials; this database is describedin de- thattheluminositydefinedinEquation6isfortheentire
tail in Section 4. SPECTRUM is general enough that it
ringofmaterialpassingthroughthe givenmeshpoint in
canaccommodate hydrodynamic profiles fromany code. Figure 1 and centered on the µ = 1 axis. Consequently,
The current version of SPECTRUM is parallelized and m is the mass of the ring:
reads in pre-computed opacity tables, derived from the
LANL OPLIB database, unlike the earlier version sum- 4 π
marized in Fryer et al. 2009. This previous version was mi = 3n (ri3+1−ri3)ρi, (7)
serial,didnotefficientlyreferencethefrequencygridand A
required an additional external calculation to determine where ρ is the ejecta density at r and n is number of
i i A
the opacities for each radial grid point. angular bins.
The calculations described below, with the monochro- Discretizing the spherically symmetric blast profile
matic OPLIB opacities in the spectrum calculations, in- uniformly in µ guarantees that every fluid element at
clude the effects of line and continuum opacities, fluo- a given radius has the same volume and so the same lu-
rescence, Doppler shifting, time dilation and limb dark- minosity. The v/c term in the denominator adjusts the
ening. Other ways of calculating these effects, such as luminositytoaccountfortherelativistictimedilationas-
the Sobolev approximation, are not necessary since we sociatedwiththemotionofthefluidandthev/ctermin
include them explicitly. the numerator accounts for Doppler shifts in photon en-
For 1D RAGE simulations, densities, velocities, mass ergyduetomotionrelativetothelineofsight. Asshown
fractionsandradiationtemperaturesfromthefinestlevel in Figure 1, if the ejecta is expanding photons from µ<
ofrefinementoftheadaptivemeshgridareextractedand
0 are redshifted and from µ> 0 are blueshifted.
ordered by radius. These profiles typically contain more
than 50,000 radial data points, but constraints on ma- 3.2. Escape Luminosity
chinememoryandtimepreventusfromusingallofthem
Tocalculatethe luminosityfromeachfluidparcelthat
inaSPECTRUMcalculation. Asubsetofthedata from
actually escapes into the interstellar medium we attenu-
each RAGE profile must be carefully chosen to fully re-
solveallradiatingregionsoftheflowfromwhichphotons ateLrνi,µj bye−τri,µj, where τri,µj is the integratedopti-
can escape, such as the photosphere of an SN shock or cal depth from the mesh point out of the ejecta to pho-
thecollisionofashockwithadensestructure. Theradial tons at frequency ν. It is constructed from consecutive
grid used in SPECTRUM for one set of SN simulations line segments as shown in Figure 2:
is described in more detail in Section 5.2.1.
3.1. Luminosity of a Fluid Element τri,µj =Xκabs,iρi|ri+1−ri|, (8)
i
5
Fig.1.— The SPECTRUM mesh. The radius of a given zone is defined at its lower face but µ = cosθ is cell-centered. The grid is
uniformlypartitionedinµ.
where ρ is the density along the segment (assumed to
i
be that of the mesh point), and the sum is over all line
segments extending from the source to the outer edge of
the grid along the line of sight. As shown in Figure 2,
the r and µ of the fluid parcel uniquely determine the
i j
lengths of all segments from its position out to the edge |r - r|
i+1 i
of the grid because all the r are known and the ends of
each segment are at the same height above the axis of
symmetry. Each zone along the line of sight from the
radiating cell intercepts photons at an energy different
from that at the source that is determined by the total
relative Doppler shift between them. Thus, the κ for
abs
each line segment is computed at the energy into which
source photons are Doppler shifted upon reaching that r r
zone. i i+1
Since Equation 8 is used to calculate the contribution
to the escaping luminosity, the use of κ here is appro-
abs
priate as only the photons that are actually destroyed
do not escape. However, in situations where there is a
strongly scattering dominated atmosphere with little or
no absorption, the thermalization depth of the photons µ
i
willbesignificantlydifferentfromthedepthoflastinter-
action. Whilethesescatteredphotonsarenotdestroyed,
their properties may change due to a scattering event
and their emergence from the ejecta will be undoubt- Fig.2.— The paths over which the τri,µj are calculated. The
ri and µj of a given fluid element set the lengths of all segments
edly delayed in time. The approximation in Equation fromitspositionouttotheedgeofthegridalongthelineofsight
8 will remain adequate for modeling the emergent spec- since the |r| aresimplythe radii of adjacent mesh points and the
trum so long as two conditions are met: (1) the added ends of the segments all lie at the same height above the axis of
symmetry.
delay from scattering processes near the thermalization
surface (τ 1) is short compared to the ejecta evo- shown in Figure 3.
abs
∼
lution timescale, and (2) the spectral information in the Duetothehighvelocitiesseenattheshockfront,even
emerging photons is not strongly altered by the scatter- with well-resolved data there can be a significant veloc-
ing interactions. For typical conditions in low energy X- ity gradientacrossjust one zone. Consequently,photons
ray emission from supernova ejecta, Thompson electron from one end of the zone can be removed at a differ-
scattering is sufficient and its coherent nature assures ent energy at the other end that is determined by the
thatthesecondcriteriaismet. InrealityThompsonscat- difference in velocity between the two ends. Depending
tering limit is not perfectly coherentand attention must on the source frequency ν, severalenergy bins may span
be paid to the number of scatters to ensure that the in- theDopplershiftofthevelocitygradientacrossthezone.
tegrated energy shift is not significant. The first criteria Wethereforefrequencyaveragetheopacityofeveryzone
is one that must be checkedfor eachnew ejecta scenario by first computing its total opacity at each energy bin
andisanimportantcaveatintheuseoftheSPECTRUM bracketed by the Doppler shift. We then sum them and
code. Scatteringdominatedatmospherescanoccurinex- divide by the number of bins across the shift to obtain
tremely energetic explosions, such as the pair-instability theaverageopacityofthezone,effectivelygivingussub-
SNe of blue compact giant stars (Whalen et al. 2012a). grid resolution for our opacities. The total luminosity
In these cases, the flux calculated in SPECTRUM with escaping the explosion along the line of sight at source
κ inEquation8insteadofκ canbe10ordersofmag- frequency ν is then just the sum of the luminosities of
tot abs
nitude lower,withapeakwavelength10times longer,as eachmeshpoint attenuatedby the matter alongthe line
6
TABLE 1
OPLIB OpacityFrequency Grid
∆u urange #ofgridpoints
0.00125 0.00125–12 9600
0.005 12.005–20 1600
0.01 20.01–30 1000
0.1 30.1–100 700
1 101–1000 900
10 1010–10000 900
100 10100–30000 200
3.3. Light Curves
BolometriclightcurvescanbeconstructedwithSPEC-
TRUM by summing over the full energy range of the
spectra at eachtime step. Light curves for a wavelength
band from a specific instrument can be created by con-
volving the spectra and a filter response function. Each
Fig.3.— Spectra froma225 M⊙ Pop III pair-instabilitySN at
1.3hours,calculatedwitheither κtot orκabs inEquation 8. Note SPECTRUMcalculationaverages200CPUhourssothe
thatluminosityispreferentiallylostatshortwavelengths. total CPU time required to compute a 250-point light
curve is one to two times that of the RAGE run. How-
ever, because many SPECTRUM jobs can run in simul-
of sight
taneouslyandeachoneiscompletedquickly,itisusually
i=nRj=nA possibletocalculatealightcurveperdayonLANLplat-
Ltνot = X X Lrνi,µje−τri,µj, (9) forms.
r=1 j=1
4. OPACITIES
where n and n are the number of radial and angu-
R A
lar bins, respectively. Equations 6 and 9 together with The Rosseland and Planck mean opacities in RAGE
the symmetry of the grid ensure that limb darkening, andthemonochromaticopacitiesinSPECTRUMarecal-
Doppler shifts and time dilation due to relativistic ex- culated with the LANL TOPS code9 from cross sections
pansionoftheejectaareincludedinourspectrumcalcu- compiled in the OPLIB database (Magee et al. 1995).
lations. We also account for atomic emission lines with An advantage to using this database is that all of the
theκ terminEquation6andabsorptionlineswiththe opacities are calculated in a consistent manner, i.e. the
abs
e−τri,µj term in Equation 9. Note that departures from results have not been culled from a variety of sources.
A disadvantage is that only electric dipole transitions
LTE in the ejecta can lead to errorsin the SPECTRUM
are considered in the bound-bound contributions, while
opacities, and hence luminosities. Furthermore, devia-
other databases include forbidden lines. However, for
tions from LTE may also partially invalidate Equation 4
LTE conditions, E1-allowed transitions should be suf-
and cause additional errors in our spectra, especially at
ficient for capturing the relevant physics. We also note
late times.
thattheOPLIBdatabaseaveragestogethersomenearby
The14,900energybinsforwhichthereareopacitiesin
spectral lines, but this level of accuracy is sufficient for
theLANLOPLIBdatabase,describedinSection4.1,can
the light-curve modeling considered here as well as for
changefrompointtopointonthegridbecausetheprofile
identifying general spectral features.
temperature varies with radius. These bins range from
We now describe these cross sections and opacities in
0.8˚A–2.0 107 ˚Aat0.5eVto4.0 10−6 ˚A–100˚Aat
105 eV. To×calculate the Lri,µj one×must choose a fixed detail and how mixed opacities are calculated for multi-
ν species flows in RAGE and SPECTRUM.
setof14,900ν forallthespectra. Thebinsweadoptare
those for OPLIB opacities at T = 1 eV. These energies
correspondtowavelengthsof0.415to9.1 106˚A,which
× 4.1. The LANL OPLIB Database
capture all spectral features from shock breakout to the
end of the light curve. We calculate Lrνi,µj using the TOPS can calculate external tables of absorption and
energy grid for the temperature in that cell and then scatteringopacitiesκforasetofdensityandtemperature
map each Lν into the standard T = 1 eV energy grid. points that the user chooses. The κ are derived from
The SPECTRUM code can be extended to post- cross sections σ in the OPLIB database that have been
process data from 2D RHD simulations. The luminosity compiledforpurechemicalelementsthatareassumedto
and the attenuation along a line of sight are calculated be in LTE, according to
with the same methods described by Equations 6 and
8. The 2Ddiscretizationis done by taking parallelslices N
0
κ (ρ,T)= σ (η,T), (10)
throughthegridthatalignwiththedesiredline ofsight. ν M ν
Each slice can then be solved as separate 1D spectra
calculation and the total luminosity can be obtained by whereN0 isAvogadro’snumber, M is the atomicweight
summing the spectra of individual slices. How we ex- ofthe elementandρis the density. The σ arecalculated
tend SPECTRUM to general 2D data sets is described
in detail in (W. Even et al., in preparation). 9 http://aphysics2/www.t4.lanl.gov/cgi-bin/opacity/tops.pl/
7
exist, TOPS will interpolate the cross section and com-
pute an opacity for that point. If the density point falls
above or below the corresponding range of ρ for which
σ are tabulated, TOPS returns the opacity of the clos-
est corresponding on-table ρ and does not extrapolate.
This procedure can lead to serious errors in opacity at
very high or low densities that are beyond the OPLIB
limits. The astrophysical simulations for which SPEC-
TRUM was developed frequently have densities that are
off-table on the low end of the original OPLIB grid. We
therefore extended the database downward by 10 orders
ofmagnitudeindensityforH,He,C,NandO,asshown
inFigure4forH.TheregioninρandT forwhichopaci-
tiesofalltheelementscanbecalculatedisverysimilarto
the redregionshownfor H; the extended regionsfor He,
C, N and O are very similar to the black region show-
ing our extension to OPLIB for H. The TOPS opacity
tables are output to files that are read in by RAGE and
Fig. 4.— The densities and temperatures for which hydrogen
SPECTRUM at execution time.
crosssectionsexistintheLANLOPLIBdatabase. Ingeneral,the
densityrangevarieswithbothtemperatureandmaterial. Thered We showhowopacitiesforhydrogenvarywithdensity
hatched lines indicate the original density range of the hydrogen and temperature in Figure 5. Absorption opacities at
cross sections and the black lines delineate our extension to this any given density generally fall as the temperature rises
dataset.
and the atoms become increasingly ionized (note that
the total opacity falls to the Thomson scattering floor
from first principles based on quantum mechanics, and
at sufficiently short wavelengths). Also, the absorption
are tabulated on a grid of temperatures T, electron de-
linesthataresoprominentatlowertemperaturesfadeat
generaciesη, and dimensionless reducedphoton energies
higheronesasthenumberofboundelectronsdecreasesin
u given by
the increasingly ionized gas. These trends demonstrate
hν thatthe atomicphysicsincludedinthe OPLIBdatabase
u . (11)
≡ kT capturestheeffectoftemperatureanddensityonopacity.
In certain applications it is advantageousto create ta-
This (T,η,u) grid is fixed for all elements. For a given η blesofopacitiesforamulti-materialmixturewithTOPS.
and T, there are 14,900 u that range from 1.25 10−3 To calculate mixed opacities, TOPS assumes that all
to 3 104, which we summarize in Table 1. × constituents are homogeneously mixed in each cell (i.e.
TO×PS can calculate gray, multi-group, or monochro- atomically mixed) and imposes the EOS condition of
matic opacities for pure elements or mixtures of ele- temperature-η (T-η) equilibrium on the constituents of
mentsandchemicalcompounds. OpacitiesinRAGEand the mixture. By virtue of Equation 13, imposing T-η
SPECTRUM are referenced by ρ and T but, as men- equilibrium is equivalent to imposing T-Ne equilibrium,
tioned above,cross sections are tabulated by η and T in whichisconsistentwiththephysicalpictureofeachcon-
OPLIB. TOPS determines the ρ and T associated with stituentexperiencingthesameglobalfree-electrondistri-
an η and T according to the standard relation bution. Furthermore,sinceanEOSprovidestherelation
between three physical quantities, imposing T-N equi-
e
libriumamongtheconstituentsisequivalenttoimposing
MN (η,T) MN (η,T)
ρ(η,T)= I = e , (12) T-P equilibrium,i.e. thetemperatureandelectronpres-
N Z¯(η,T)N e
0 0 suresassociatedwitheachconstituentareequaltothose
ofthe mixture. Undertheseassumptions,TOPSobtains
where ρ is the density of the pure element or mixture,
the opacity for a mixture at conditions ρ and T by sum-
N is its ionnumber density andZ¯ is its averagecharge.
I ming the cross sections for pure elements according to
The relation between N and η is assumed to be that of
e
a Fermi degenerate ideal gas:
N
4 (2πm kT)3/2 κν(ρ,T)= M0 Xfiσνi(η,T), (14)
e i
N = F (η), (13)
e √π h3 1/2
whereN isAvogadro’snumber, M is the atomicweight
0
of the mix, f is the atomic fraction of the ith element,
whereF (η) is the Fermi-Diracintegraloforder1/2,k i
1/2 andη is chosento satisfy Equation12. TOPScalculates
istheBoltzmannconstant,hisPlanck’sconstantandm
e
grayor multi-groupmixedopacities by integratingthese
is the mass of an electron. When constructing an opac-
monochromaticopacitieswiththe appropriateweighting
ity table, TOPS will produce data only at temperatures
function.
that existin the OPLIB database,i.e. neither interpola-
tionnorextrapolationareperformedforrequestedTthat
4.2. RAGE Opacities
are not in the database. For this reason, the T points
we adopt in our opacity tables are the same as those in We use TOPS to precompile Rosseland and Planck
OPLIB. On the other hand, if a requested density point meanopacities for eachpure element of the SN ejecta in
falls between two corresponding ρ for which OPLIB σ a single external file for RAGE. However, since a mesh
8
Fig. 5.—Totalopacity(upperrow)andabsorptionopacity(lowerrow)forhydrogenasafunctionoftemperatureandwavelength. Left
toright: T =0.5eV,2.0eV,and100eV.
point on the RAGE grid might be populated by several more accurate than those obtained under the assump-
elements, RAGE must construct a single opacity for the tion that electrons dominate the pressure, as in TOPS
mixture. RAGE does not construct mixed opacities in (Zeeb & Fontes 2006). However, because we are using
thesamewayasTOPSbecauseitimposesdifferentEOS an ideal gas EOS rather than the RAGE tabular EOS,
mix conditions on the constituents. Instead of taking we ignore electron pressure when calculating the ρ .
i
themtobeatomicallymixedineachcell,RAGEassumes In principle, one should first construct the sum in
that they remain physically segregated (i.e. consistent Equation 15 at each frequency and then average over
with a heterogeneous or chunk mix), each with its own frequency to compute the mean mixed opacity rather
volume fraction and partial density (for additional de- than sum opacities that have already been frequency-
tails concerning the atomic and chunk mix assumptions, averaged. In the Planck mean, the two approaches yield
see Smith 2003; Pomraning 1991). Like TOPS, RAGE thesameopacitybecausethemeanisasimplearithmetic
imposes T-P equilibrium onthe components of the mix- averagefor whichsumming andaveragingare commuta-
ture, but considers both the ion and electron pressure. tive. ThetwooperationsarenotcommutativeforRosse-
Morespecifically, inRAGE the ionpressureofeachcon- land opacities because the mean is harmonic, not arith-
stituent is equal to the ion pressure of the mixture and metic. Unlike TOPS, which frequency averages mixed
the electron pressure of each component is the same as monochromaticopacitiesto compute meanmixedopaci-
theelectronpressureofthe mixture. After obtainingthe ties,RAGEsumsmeanopacities,eventhoughitissome-
partial densities from the above procedure, the opacity what less accurate, because computational constraints
of the mixture is obtained from preclude a RAGE calculation for all 14,900 u. Finally,
to construct the mixed opacity at a given mesh point at
κ¯(ρ,T)=XfiMκ¯i(ρi,T), (15) some ρ and T, RAGE performs a bi-linear interpolation
i of log κ¯ from Equation 15 over log ρ and log T.
whereκ¯representsamean(grayorgroup)opacity,fM is
i 4.3. SPECTRUM Opacities
themassfractionandρ thepartialdensityoftheithcon-
i
stituent,whichiscalculatedfromitsvolumefractionand WeuseTOPStocreatemonochromaticabsorptionand
mass fraction. The RAGE simulations presented in this total opacity tables for SPECTRUM, one for each ele-
paper use gray opacities, so κ¯ represents the Rosseland ment. We calculated tables for the first 30 elements at
or Planck mean opacity. In regions where ion pressure the 50 OPLIB temperatures from 5 10−4 keV to 100
is non-negligible and material is not atomically mixed, keV and a grid of 81 densities from 1×0−15 to 105 g/cm3.
Equation15isexpectedtoyieldmixedopacitiesthatare Opacities for H, He, C, N, and O were calculated for an
9
additional20densitiesfrom10−25 to10−15 g/cm3,using
data from the extended tables described earlier. Opaci-
ties for all 14,900u are tabulated for each ρ and T. The
OPLIBugridstructure (showninTable 1)simplifies re-
trieval of the desired opacity from the tables because af-
tertheuintervalinwhichagivenfrequencyν resideshas
been determined, a simple index that directly references
that interval is easily computed. This approach results
inspectrumcalculationtimesthataremuchshorterthan
those done with bisection searches.
To calculate the mixed opacity for a cell at a given ρ
andT,SPECTRUMfirstcomputestheopacityκ ofeach
i
constituent of the mixture at its partial density ρ . At
i
thegivenfrequency,logκi isextractedfromthe tableat
ν
thetwodensitiesthatbracketρ atthetabletemperature
i
closestto T. SPECTRUMthenperformsa simple linear
interpolationofthelogarithmofthesetwoopacitieswith
respect to log ρ while holding T constant to obtain the
i
desired κi, i.e. as in TOPS, SPECTRUM does not in- Fig. 6.—Analytical(solidlines)andSPECTRUM(dottedlines)
ν blackbodyspectraat1eVand100eV.
terpolate κi in T. Like RAGE,the mixed opacity of the
ν
cell is then taken to be the sum of the opacity of each
constituent weighted by its mass fraction, according to withamasslossrateof10−5M⊙/yrandavelocityof108
cm/s. Thisstarwasoriginallyinvestigatedasapotential
κν(ρ,T)=XfiMκiν(ρi,T), (16) CthaeseAffepcrtosgoefnpitroer-caonrdeicsonlloawpsbeeminagssstluodsisedontoshuoncdkerbsrteaankd-
i outandSNlightcurves. Theconvergencetestspresented
which is the monochromatic version of Equation 15. hereareforasingletimeaftershockbreakout,400s,but
we have performed identical tests at a variety of times
5. SPECTRUMTESTS
beforeandafterbreakoutandfindsimilarresults. These
WeranaseriesofteststoverifytheSPECTRUMcode radialandangularbinconvergencetests wereconducted
andtodeterminethenumberandstructureofradialbins in tandem over a 2D parameter space, 500-9000 radial
and the number of angular bins required to converge to bins and 2-400 angular bins, to find the best combina-
accuratespectragivenmemoryandtimeconstraints. We tion of radial and angular resolution.
find that convergence is problem specific, and should be
investigated for each new type of simulation. 5.2.1. Radial Convergence
As noted in Section 3, the RAGE simulations used
5.1. Blackbody Tests
as input to SPECTRUM contain at least 50,000 radial
One basic test every spectra code should be able to datapointsbutcurrentmemoryconstraintslimitSPEC-
pass is to reproduce an analytic blackbody in the case TRUM input to less than 10,000 points, requiring the
of a constant opacity. To do this, we calculated spectra creation of a new radial grid. Densities, mass fractions,
for a static sphere with constant temperature, density velocities,andtemperaturesfromthefinestmeshareex-
and opacity, surrounded by a vacuum. To obtain accu- tractedfromtheRAGEprofilesandsequentiallyordered
rateblackbody spectra,the opacityanddensity must be by radius. Three different binning schemes were con-
selected so that the τ =20 surface is located a small dis- sidered: linear, logarithmic, and a hybrid scheme which
tance inside the outer surface. We initialized the sphere focuses resolutiononthe radiatingregion. Asimple bin-
in SPECTRUM with a very fine radial resolution at the ningschemewithuniform,linearlyspaced,binswillcom-
emitting surface. As we show in Figure 6, SPECTRUM pletelyexcludetheinteriorofthesimulationandonlyre-
accurately reproduces blackbody spectra for a range of solve the wind. Logarithmically spaced bins include the
temperatures. full range of radii but do not sufficiently resolve the ra-
diating regionaroundthe shock. Figure 7 shows spectra
5.2. Convergence Tests
for these two binning schemes. To obtain an accurate
As explained in Section 3, for 1D problems SPEC- binning scheme, any region where rapid changes in den-
TRUM maps data onto a 2D grid with radial and an- sity, velocity or temperature occur and from which pho-
gular bins. We test the convergenceof our SPECTRUM tons can escape, such as around the shock or any shell
calculation in both angle and radius with spectra from or torus features in the wind, must be defined and suf-
a fiducial 23 M⊙ supernova. The progenitor, which was ficiently resolved in the new grid. The criteria used to
evolved in TYCHO (Young & Arnett 2005), is a binary define this new grid will change based on the details of
starthatlosesits Henvelopeatthebaseofthe redgiant eachproblembeingstudied. Generallytheinteriorofthe
branch, evolves as a Wolf-Rayet star, and then ends its star, from which luminosity cannot escape, and the stel-
lifeasa6.4M⊙ CNOcorethatexplodesasaTypeIcSN lar wind, which is not finely structured, can be included
(Young et al. 2006). We first modeled the collapse of its at low resolution. A hybrid scheme which combines lin-
core, the launch of the shock and all nuclear burning in ear and logarithmic binning was used for all subsequent
a 1D hydrocode (Fryer 1999) and then mapped its pro- spectrainthispaperandiscomparedtosimplerschemes
file into RAGE together with a wind envelope, defined in Figure 7. This scheme includes finely spaced linear
10
Fig.7.— OutputspectraforidenticalRAGEdatafromaType Fig.8.— Radial convergence test, running SPECTRUM with
IcSNat400secondsaftercore-collapse,calculatedwith160angu- the same RAGE data as 7 and 160 angular bins. This data is
lar bins and 5200 radial bins organized by three different binning binned with varying radial resolution between the τ = 20 surface
schemes: linearbins,logbins,andthehybridschemedescribedin andtheradiationfront. Notethatthespectracalculatedwith5000
thetext. and9000radialbinsarenearlyidentical.
bins between the radiation front and the τ=20 surface
and more widely spaced log bins inside that radius and
in the wind. The τ=20 surface is calculated using the
Thomson scattering opacity, and defines the minimum
inner bound of the regionfrom which luminosity can es-
cape. The RAGE data inside each of these new radial
bins is mass-averagedto ensure that very sharp features
are included.
To find the minimum number of bins that resolve the
radiating region well, we performed a convergence test
usingbetween500and9,000radialbinsintheregionbe-
tweenthe τ=20surfaceandthe radiationfront,withthe
hybrid scheme described above. For each of these tests,
the interior region and the wind are each resolved by an
additional 100 logarithmically spaced bins. As shown in
Figure 8, this SN Ic spectrum converges at about 5000
radialbins. Thespectrainthis figurewereallcalculated
with160angularbins. InadditiontotheseSPECTRUM
tests,wecomparedeachbinnedparameterwiththecom- Fig. 9.— Angularconvergence test, runningSPECTRUMwith
thesameRAGEdataas7,5200radialbinsandavaryingnumber
plete RAGE data set to confirm that the convergence
ofangularbins.
shown here does correspond to a well resolved profile.
The numerical calculations of SN light curves per-
5.2.2. Angular Convergence formed to date have implemented varying degrees of
physics, as noted in Section 2. To investigate the im-
When the number of angular bins is increased in
portanceofradiationprocessesthathavebeenneglected
SPECTRUM,a finer gridis usedto calculateluminosity
in past studies of SN explosions, we performed five 1D
along the line of sight for a single grid point (see Figure
simulations of the Type Ic SN described in Section 5.2
1. This results in a more accurate luminosity calcula-
in which additional physics was incrementally activated
tion. A convergence study was performed with between
in RAGE. Because our goal is to examine how radia-
2 and 400 angular bins and a fixed number of radial
tionphysics affects SN profiles,spectra and lightcurves,
bins. Figure 9 shows that any number of angular bins
we make no attempt to compare hydrodynamical algo-
greater than two gives a nearly identical continuum and
rithms. The five cases we consider are as follows:
that line intensities converge at 160 angular bins, with
5200 radial bins. However, we note that higher resolu-
2T: Gray flux-limited diffusion (FLD) in which
tion might be required for detailed studies of individual •
matterandradiationtemperaturesareevolvedsep-
lines. We again emphasize that convergence in general
arately in the simulation. This physics option is
is problem-dependent because the optimum scheme for
usedinallthe simulationspresentedinsubsequent
allocating grid points depends on the structure of the
sections.
radiating flow.
1T: Gray FLD in which matter and radiation are
6. PHYSICSEFFECTSOFRADIATIONHYDRODYNAMICS • assumed to have the same temperature.