Kilonovae Across the Nuclear Physics Landscape: The Impact of Nuclear Physics Uncertainties on r-process-powered Emission

Much of the uncertainty in simulations of r-process nucleosynthesis is due to the discrepant predictions of various theoretical models used to compute unmeasured properties of the nuclei involved. While a full discussion of the nuclear physics inputs we consider is provided in Zhu et al. (2021, henceforth Z20), the main points are summarized here as well.

We base our theoretical nuclear inputs on each of eight distinct nuclear models listed in Table 1. Unmeasured reaction rates, as well as α- and β-decay lifetimes and branching ratios, are adjusted to take into account Q-values corresponding to each mass model in a manner consistent with the predicted nuclear masses (Z20, and references therein, especially Mumpower et al. 2015).

Table 1. Nucleosynthesis Parameter Space for the Full Set of Simulations

Nuclear Mass Models
Mass Model	References
Finite-Range Droplet Model (FRDM2012)	Möller et al. (2016)
Duflo & Zuker (DZ33)	Duflo & Zuker (1995)
Hartree–Fock-Bogoliubov 22 (HFB22)	Goriely et al. (2013)
Hartree–Fock-Bogoliubov 27 (HFB27)	Goriely et al. (2013)
Skyrme-HFB with UNEDF1 (UNEDF1)	Kortelainen et al. (2012)
Skyrme-HFB with SLY4 (SLY4)	Chabanat et al. (1998)
Extended Thomas-Fermi plus Strutinsky Integral (ETFSI)	Aboussir et al. (1995)
Weizächer-Skyrme (WS3)	Liu et al. (2011)
Fission Prescription
Mass Model	Barrier Heights Adopted
HFB22	HFB14 (Goriely et al. 2009)
HFB27	HFB14
ETFSI	ETFSI (Mamdouh et al. 2001)
All others	Finite-Range Liquid-Droplet Model
	(FRLDM; Möller et al. 2015)
Fission Fragment	References
Distribution
Symmetric	e.g., Rauscher et al. (1994)
Gaussian	Kodama & Takahashi (1975)
Spontaneous	Reference
Fission Lifetimes
Xu & Ren (XuRen)	Xu & Ren (2005)
Karpov/Zagrebaev (KZ)	Karpov et al. (2012); Zagrebaev et al. (2011)
Astrophysical Parameters
Ye	0.02, 0.12, 0.14, 0.16, 0.18, 0.21, 0.24, 0.28
sB	40kB
${\tau }_{\exp }$	20 ms

Download table as:  ASCIITypeset image

The potential energy landscape of fission influences fission half-lives and the mass distribution of the daughter fragments, both of which can be important, particularly in the case of strong fission cycling. Our nucleosynthesis calculations include neutron-induced, β-delayed, and spontaneous fission, and so the treatment of fission adds an additional dimension of uncertainty.

Fission barrier heights have been calculated within the framework of several nuclear models, but have not been calculated for every model we consider here. When possible, we adopt barrier heights computed for a mass model that closely resembles our model of interest. In cases where there is no clearly associated mass model for which barrier heights have been determined, we use the barrier heights predicted by the Finite-Range Liquid-Droplet Model (FRLDM; Möller et al. 2015).

Our r-process simulations use phenomenological descriptions of the fission fragment mass distribution that do not take barrier heights as input. The distribution is taken to be either Gaussian, according to the formalism of Kodama & Takahashi (1975), or symmetric (A_daught = A_parent/2). The exception is ²⁵⁴Cf, for which we employ fission yields calculated with FRLDM (Mumpower et al. 2020), as in Zhu et al. (2018). We also use phenomenological prescriptions for spontaneous fission lifetimes, calculating these rates using the either the method of Xu & Ren (2005) or the barrier height-dependent prescription of Karpov et al. (2012) and Zagrebaev et al. (2011).

Because spontaneous fission plays an important role in heating, the four combinations of fission yields and spontaneous fission rates are intended to bracket the range of possible fission behavior. The rates of Xu & Ren (2005) become substantial around Z > 94, while those of Zagrebaev et al. (2011), which depend explicitly on fissibility (Z²/A) and barrier height, are typically fairly low until Z > 100. (For a more detailed discussion of these two spontaneous fission models, see Vassh et al. 2020). Meanwhile, a symmetric split produces a very narrow set of daughter products, while the yields of Kodama & Takahashi (1975) are more broadly distributed.

In addition to nuclear physics uncertainties, the r-process is also sensitive to the outflow properties of the r-processing gas. While the thermodynamic and hydrodynamic properties of the outflow (entropy and expansion rate, e.g.). can impact the r-process, in the regime of robust heavy element production (i.e., in our regime of interest), nucleosynthesis has been shown to be particularly sensitive to the initial electron fraction Y_e (Freiburghaus et al. 1999b; Goriely et al. 2011; Lippuner & Roberts 2015), defined as the ratio of protons to baryons in the gas. We therefore use Y_e to probe the dependence of the r-process on astrophysical conditions, and consider for each set of nuclear physics inputs eight values in the range 0.02 ≤ Y_e ≤ 0.28. The nuclear physics models and astrophysical quantities that define our parameter space are presented in Table 1. All possible combinations of these parameters are considered, yielding a full set of 256 models.

We use the Portable Routines for Integrated nucleoSynthesis Modeling (PRISM; Sprouse T. 2021 in preparation) to simulate the r-process for each of our 256 models for a gas trajectory with an initial entropy per baryon s_B = 40k_B (with k_B the Boltzmann constant) and an expansion timescale ${\tau }_{\exp }=20$ ms during the r-process epoch. This trajectory has also been used in other studies of post-merger nucleosynthesis (Zhu et al. 2018, Z20). The choices of s_B and ${\tau }_{\exp }$ are motivated by models of disk wind outflows from NS mergers, and are therefore consistent with our focus on the heavy r-process (red) kilonova component, which both spectral analysis (Chornock et al. 2017; Kasen et al. 2017) and theoretical simulations (Siegel & Metzger 2018) suggest is associated with a disk wind (although see Miller et al. 2019 for an alternate view). We set the initial temperature to 10 GK for all calculations, and compute the seed nuclear abundances with the SFHo equation of state (Steiner et al. 2013). The temperature evolution of the expanding gas is determined self-consistently for each simulation using the same reheating procedure as in Zhu et al. (2018) and Vassh et al. (2019).

PRISM tracks the nuclear abundances for each model as a function of time. The absolute nuclear heating is then computed from the evolving composition, using known or predicted decay modes and branching ratios for all nuclei, and measured or theoretical nuclear masses, which allow a determination of the energy Q released in each decay. The absolute heating rate, defined as the total energy emitted by radioactive decays per unit mass and time, is then a sum over decays.

The calculations performed on the full set of 256 models delineated, for the parameter space under consideration, the potential variation in total r-process heating. From the full set of parameter combinations, we selected 10 models whose absolute nuclear heating spans the range indicated by the full set.

By design, these models represent the maximum variation in absolute heating predicted by the full set, and some therefore appear as outliers. However, even the full set does not account for all possible nuclear physics uncertainties, and moreover reflects our decision to focus on a limited set of mass models and fission descriptions. The true allowed spread in heating is therefore expected to be larger even than our model set indicates.

The parameters that define our 10-model subset are recorded in Table 2. To this set of 10, we add a single model calculated using masses and fission barriers determined by the Thomas-Fermi (TF) mass model (Myers & Swiatecki 1996, 1999). TF barrier heights were calculated using the prescription applied in the GEF code (2016 version; Schmidt et al. 2016), as in Vassh et al. (2019). We include this model to explore the impact of TF barriers, which tend to be lower than the other barriers we consider in key regions in a way that limits the production of long-lived fissioning isotopes (Vassh et al. 2019). For this model, we apply the "D3C*" β-decay rates of Marketin et al. (2016), which have been used in conjunction with Thomas-Fermi barriers for some fission channels in studies of the r-process (e.g., Wu et al. 2019). While we do not incorporate TF+D3C* into the nuclear physics parameter space defined in Table 1, we adopt it in Model 11 to facilitate comparisons between our work and that of other groups.

Table 2. Properties of the Model Suite

Model Index	Nuclear Mass Model	Fission Fragment Distribution	Spontaneous Fission τnuc	Ye	Model Label	XL,A a
1	FRDM2012	Symmetric	Karpov/Zagrebaev	0.16	frdm.y16	0.22
2	FRDM2012	Symmetric	Karpov/Zagrebaev	0.28	frdm.y28	0.02
3	HFB22	Symmetric	Karpov/Zagrebaev	0.16	hfb22.y16	0.37
4	HFB27	Gaussian	Karpov/Zagrebaev	0.16	hfb27.y16	0.48
5	DZ33	Gaussian	Karpov/Zagrebaev	0.16	dz33.y16	0.31
6	UNEDF1	Gaussian	Karpov/Zagrebaev	0.16	unedf.kz.y16	0.34
7	UNEDF1	Gaussian	Xu & Ren	0.16	unedf.xr.y16	0.33
8	UNEDF1	Symmetric	Karpov/Zagrebaev	0.24	unedf.y24	0.38
9	SLY4	Symmetric	Karpov/Zagrebaev	0.18	sly4.y18	0.03
10	SLY4	Symmetric	Karpov/Zagrebaev	0.21	sly4.y21	0.12
11 b	TF+D3C*	Symmetric	Karpov/Zagrebaev	0.16	tf.y16	0.22

Notes.

^a ${X}_{{\rm{L}},{\rm{A}}}^{}$ is a result of, not an input parameter to, our nuclear physics simulations. ^b Model 11 was not selected from the full suite of simulations. It was instead constructed specifically to expedite comparison between this work and that of other groups. See Section 2.2 for further discussion.

Download table as:  ASCIITypeset image

The absolute heating rate and abundance information for each model in our 11-model subset are presented in Figure 1. Like the full suite of simulations (the results of which we plot in gray scale for comparison), the model subset exhibits significant diversity in both the energy produced and the composition synthesized. The first of these has implications for thermalization, which affects the amount of energy available to power the light curve, while the latter is important for determining opacity, which is sensitive to the abundances of lanthanides and actinides in the ejecta.

Zoom In Zoom Out Reset image size

In the interest of fully exploring the potential variability in r-process heating, we did not restrict ourselves to models whose abundance patterns conformed to the solar r-pattern, and the abundance yields for our models vary substantially. We did ensure that all our models had a combined lanthanide and actinide mass fraction, ${X}_{{\rm{L}},{\rm{A}}}^{}$, high enough to plausibly produce the high opacity required to explain the red component of the G170817 kilonova. (In our model subset, X_L,A ≥ 0.02 at t ≈ 1 day; see Table 2).

In fact, some models have ${X}_{{\rm{L}},{\rm{A}}}^{}$ greater than typically considered in studies of kilonovae (${X}_{{\rm{L}},{\rm{A}}}^{\mathrm{typ}}\lesssim 0.1$). The effects of higher ${X}_{{\rm{L}},{\rm{A}}}^{}$ are explored in Section 5. Here, we simply affirm that consideration of high-${X}_{{\rm{L}},{\rm{A}}}^{}$ models is warranted given (a) their potential relevance for NSBH kilonovae, (b) indications from r-process-enriched stars that a sizeable fraction of r-process events should have ${X}_{{\rm{L}},{\rm{A}}}^{}$ greater than is usually attributed to the kilonova of GW170817 (Ji et al. 2019), and (c) observations of an actinide "boost" in ∼30% metal-poor stars (Mashonkina et al. 2014), seeming to require copious actinide production by the r-process (Holmbeck et al. 2019). We also note that the imperfect (though still impressive) agreement between models and observations (e.g., Chornock et al. 2017; Kasliwal et al. 2017; Smartt et al. 2017) of the GW170817 kilonova—the lone event of its class—make it premature to conclude that high ${X}_{{\rm{L}},{\rm{A}}}^{}$ lie outside the realm of possibility even for NS mergers.

Finally, we note that X_L,A in our models does not uniformly increase as Y_e decreases (e.g., sly4.y18 and sly4.y21). For relatively high Y_e, a complete r-process is not obtained, and an abundance drop-off occurs just after the point at which the system runs out of neutrons. If this happens before the material reaches the third peak, there will be only small amounts of lanthanides and no actinides. If, on the other hand, the system has just enough neutrons to make a full third peak but no heavier material, then the lanthanide mass fraction is usually at a maximum, although actinide production is minimal. (This occurs around Y_e = 0.21 for SLY4). At only slightly lower Y_e, when the material pushes past the third peak and just begins to substantially fission, the lanthanide distribution reaches a local minimum as a function of Y_e. The production of actinides may be enhanced, but this is unlikely to compensate for the loss of lanthanides, leading to an overall decrease in X_L,A. (This occurs at around Y_e = 0.18 for SLY4). Continuing to still lower Y_e, the lanthanide distribution starts to build up again, and actinide production becomes increasingly robust. For a more detailed discussion, see Z20 (their Section 3 and Figure 5).

We carry out detailed numerical calculations of thermalization efficiencies and kilonova light curves for these 11 models.

Each set of nuclear physics models we consider predicts the decay modes and branching ratios of any nuclear decay for which these quantities have not been measured. Together with experimental data, these predictions allow us to determine the fraction of the absolute heating rate contributed by the three decay channels important on kilonova timescales: α-decay, β-decay, and spontaneous fission. (PRISM also tracks β-delayed and neutron-induced fission, but these do not operate after very early times). The decay products from each channel have distinct energy scales and cross sections for energy-loss processes, and so are treated individually in thermalization calculations.

In α-decay (fission), the decay energy mainly accrues to the α-particle (fission fragments), and so the kinetic energy of the emitted α-particle (fission fragments) is taken to equal the energy Q of the decay. In contrast, the energy from a β-decay is shared by a β-particle, a neutrino, and an arbitrary number of γ-rays. The γ-rays thermalize only weakly, and neutrinos do not thermalize at all, so additional analysis is needed to determine the division of β-decay energy among these species.

We determined at each time step the 50 nuclei producing the most energy through β-decay, and calculated from this subset the fraction of the total β-decay energy emerging as β-particles, γ-rays, and neutrinos using β-decay endpoint energies, average kinetic energies, and intensities from the Nuclear Science References database (Pritychenko et al. 2011), which we accessed through the website of the International Atomic Energy Agency. We assumed parent nuclei decay from the nuclear ground state (though see Fujimoto & Hashimoto (2020) and Misch et al. (2021) for a discussion of isomeric decay in the r-process). When possible, we estimated unknown quantities (e.g., unrecorded average β-particle energies were assumed to be 1/3 of the corresponding β endpoint energy). If there were insufficient data to build a reasonable model of a given decay, the nucleus in question was excluded from our analysis. However, by the time the kilonova begins to rise (t ∼ few hours), the composition has already stabilized relative to earlier times, and exotic unmeasured isotopes are the exception rather than the rule. Thus, few nuclei were omitted.

The distribution of the absolute heating among decay products is shown for each model as a function of time in the left-hand panels of Figure 2. The models exhibit considerable diversity, particularly in regard to the importance of α-decay and fission, which have been shown to enhance thermalization (B16). Roughly half the models show the contributions from these channels growing slowly, accounting for ≳30% of the absolute heating by t = 30 days. This increase is due to the production of certain actinide species (in the case of α-decay) and the spontaneous fissioning isotope ²⁵⁴Cf (in the case of fission), whose long half-lives allow them to dominate the energy release at late times after the background β-decay radioactivity has subsided. Other models (frdm.y28, unedf.y24, sly4.y21, and tf.y16) have only negligible contributions from α-decay or fission, and still others (hfb22.y16 and hfb27.y16) have heating rates that are dominated by fission over the entire kilonova timescale.

Zoom In Zoom Out Reset image size

We used the particle-specific absolute heating rates extracted from the PRISM calculations in our thermalization simulations. This improves on the method of B16, which assumed the absolute heating from all decay modes and particle types obeyed a power law. Recent analytic work (Kasen & Barnes 2019; Waxman et al. 2019) showed that the interplay between the absolute heating associated with a given decay product and that product's time-dependent emission spectrum influences thermalization, motivating the more rigorous treatment of ${\dot{\epsilon }}_{\mathrm{abs}}$ in this work. The range of behavior seen in our models—in terms of both of the shape and magnitude of the absolute heating curve, and the relative importance of the different decay products to the total heating—suggests that thermalization efficiencies will be variable.

Because the cross sections of thermalizing processes are energy dependent, the emission spectra of radioactive decay products are required for thermalization calculations. We construct the time-dependent emission spectra for α-particles, β-particles, and γ-rays with a procedure similar to that of B16. We assign to emitted α-particles the entire Q-value of the decay that produced them (α-decay typically proceeds directly to the daughter ground state, so the incidence of emitted γ-rays is vanishingly low). The α-decay spectrum then simply reflects the fraction of the total α-decay heating contributed by every α-decaying nucleus at a given time.

The spectra of β-decay products are less straightforward due to the three-body nature of the decay and the fact that β-particle and neutrino emission generally land the daughter nucleus in one of many excited nuclear states, with de-excitation to the ground state occurring γ-ray emission. We compute the β-particle spectrum for each nucleus from available endpoint energies and intensities (Pritychenko et al. 2011), using the simplified fit to the Fermi formula for allowed transitions proposed by Schenter & Vogel (1983). In cases where forbidden transitions are important, the true spectrum may deviate from our calculations. The use of forbidden-transition shape factors in the determination of β-decay spectra could improve thermalization calculations in the future.

The γ-ray spectrum was determined from recorded γ-ray energies and intensities. The β-particle (γ-ray) spectrum for the entire composition was then calculated by weighting the spectrum from each isotope by the fraction of the total energy in β-particles (γ-rays) contributed by that isotope as a function of time.

Snapshots of the α-particle, β-particle, and γ-ray emission spectra for each model at t = 1 day are provided in the right-hand panels of Figure 2. The α-spectrum is dominated by about a dozen nuclei, each of which emits an α-particle with a single characteristic energy. The β-particle and γ-ray spectra are smoother because of the continuous nature of the emission from individual nuclei (β-particles) or the much larger number of discrete emission energies (γ-rays).

Because of the inverse relationship between Q and lifetime for β-decaying nuclei (Colgate & White 1966; Metzger et al. 2010; Hotokezaka et al. 2017), β-decay spectra tend to shift toward lower energies as time progresses. Even by t = 1 day, many of the less stable (higher Q) isotopes have decayed away, concentrating the β- and γ-ray spectra at E ≲ 1 MeV. However, late-time fission could produce a population of unstable β-decaying nuclei with high Q, in defiance of this trend (Z20).

Unlike for other decay channels, we do not attempt a full calculation of the fission fragment spectrum, which would need to be described in terms of both daughter mass A_ff and initial kinetic energy E_ff,0. However, test calculations of fission fragment thermalization that adopted δ-function distributions of A_ff and E_ff,0 over the ranges 100 ≤ A_ff ≤ 200 and 100 ≤ E_ff,0/MeV ≤ 200 yielded very similar results, suggesting that fission fragment thermalization is not highly sensitive to these quantities. Therefore, although fission daughter distributions are tracked internally by PRISM to model nucleosynthesis, we do not extract from PRISM simulations a detailed description of fission fragment production for our thermalization calculations. Rather, we adopt a fiducial fragment mass A_ff = 150 and a flat emission spectrum that ranges from 100–150 MeV.

The rate at which radioactive decay products transfer their energy to the ejecta depends on the ejecta composition. Because our models produce different abundance patterns, these rates vary slightly from model to model, and we calculate for each model a set of energy-loss rates consistent with its composition. The variation is generally small, however, as can be seen in Figure 3, which shows the range of energy-loss rates predicted by our models. The variance is usually within a factor of 2, and is often much less. We summarize the relevant interactions below, but refer to B16 (their Section 3) for a detailed description.

Zoom In Zoom Out Reset image size

γ-rays: The γ-rays from radioactive decay lose their energy in photoionization and Compton-scattering events. The associated opacities are plotted in the top panel of Figure 3. We constructed photoionization cross sections for each model's composition from the element-specific, energy-dependent cross sections available through the Photon Cross Section Database (XCOM; Berger et al. 2010), published by the National Institute of Standards and Technology (NIST). Compton-scattering cross sections were calculated analytically using the Klein-Nishina formula. Photoionization opacity is much higher than Compton opacity at E_γ ≲ 0.5 MeV. Both processes produce an electron, which often has an energy greater than the thermal energy of the background gas. These secondary electrons are propagated through the ejecta until they have been effectively thermalized.

Although not included in our thermalization model, pair production in the electric field of the nucleus is an additional energy-loss channel for γ-rays with energies above a few MeV. Our models produce relatively few such γ-rays (e.g., Figure 2), so the cost of this omission is not as high as it otherwise would be. However, this channel deserves to be investigated in the future. This is especially true because pair production could take on additional importance for thermalization in models featuring particularly energetic β-decays.

β -particles: Energetic β-particles transfer their energy by exciting or ionizing bound electrons, or through Coulomb interactions with free electrons. The highest energy β-particles may produce Bremsstrahlung radiation that can itself thermalize or escape from the ejecta. The importance of each of these energy-loss channels is presented in the second panel of Figure 3.

We model ionization and excitation losses—which dominate β-particle thermalization—with the well-established Bethe-Bloch formula (Heitler 1954; Berger & Seltzer 1964; Gould 1975; Blumenthal & Gould 1970), including relativistic corrections. We use mass-weighted averages for composition-dependent quantities, such as the bound-electron number density and the average excitation energy. Losses due to Coulomb interactions with unbound thermal electrons are described by the formulae of Huba (2013). We adopt the formula appropriate in the limit that the energy of the thermalizing particle greatly exceeds the energy of thermal electrons (as in Equation (4) of B16). Finally, the Bremsstrahlung stopping is computed following Seltzer & Berger (1986), which requires the use of Z-dependent empirical fitting constants. We use constants corresponding to the mass-weighted average atomic number of each model's composition.

α -particles: Like β-particles, α-particles thermalize via interactions with background thermal electrons. In the case of free electrons, these interactions can be modeled using the formula of Huba (2013), describing the energy lost by a suprathermal ion to thermal electrons in a plasma. The energy lost to bound electrons is calculated from α-particle stopping powers retrieved from NIST's Stopping Power and Range Tables for Helium Ions (ASTAR; Berger et al. 2005). As in B16, we map the full composition of each model to a simplified composition comprised of elements for which α-stopping data are available. Alpha-particle energy-loss rates are plotted in the third panel of Figure 3. Plasma losses assume the ejecta is singly ionized, a condition that may not be met at very late times.

Fission fragments: In addition to free and bound electrons, fission fragments thermalize through Coulomb interactions with atomic nuclei. Fission fragment thermalization is complicated by the question of ion state; fragments are not fully ionized, and their charge state Z_ff,eff depends on their kinetic energy. We determine Z_ff,eff using the formula of Schiwietz & Grande (2001).

We can then calculate losses due to free electrons using the same ion-electron formula as for α-particles, again approximating the ejecta as singly ionized, and with terms involving ion mass and charge updated appropriately. As in B16, in the case of free electrons only, we set a floor Z_ff,eff ≥ 7, motivated by the fact that low ionic charge states permit thermal electrons to impact the fragment at distances smaller than the fission fragment radius. To model interactions of fission fragments with bound electrons, we follow the procedure of Ziegler (1980) and scale proton stopping powers, which we take from NIST's PSTAR database (Berger et al. 2005), by ${Z}_{\mathrm{ff},\mathrm{eff}}^{2}$. The energy loss from interactions between fission fragments and atomic nuclei is given by the nuclear stopping formula of Ziegler (1980); it is subdominant for the energies of interest, as is apparent in the bottom panel of Figure 3.

Radioactive particles are emitted in the ejecta at a rate proportional to the local mass density. We assume that the ejecta has a uniform composition, so radioactive decay does not depend on position within the ejecta. (In reality, thermodynamic conditions are likely to vary within an outflow, resulting in spatially dependent patterns of nucleosynthesis and radioactive decay. We defer an exploration of these effects to later work). As decay products traverse the ejecta, either following magnetic field lines (in the case of charged fission fragments, α-, and β-particles) or undergoing discrete scattering and absorption events (in the case of γ-rays), the energy they transfer to the ejecta is calculated and tallied up.

Unlike in B16, we track the deposition of radioactive energy in both time and space, enabling a more detailed description of thermalization. It has long been recognized that because thermalization depends on density, its efficiency should vary within the ejecta (e.g., Wollaeger et al. 2018). In particular, because the densest regions will emit more radioactive energy, thermalize that energy more efficiently, and remain optically thick out to later epochs than less dense regions, spatially dependent thermalization can impact the temperature of the ejecta's photosphere, and thus the overall brightness and effective temperature of the kilonova (see Korobkin et al. 2020 for a nice discussion).

To better highlight the variation in light-curve predictions arising from the choice of nuclear physics parameters, we consider in this study only one ejecta mass, velocity, and density profile. Our model is spherically symmetric with a broken power-law density profile, ρ(v) ∝ v^−η , with η = 1 (10) in the inner (outer) ejecta layers, and the transition point set by the requirement that the mass distribution yields the specified ejecta mass and kinetic energy. Motivated by the inferred mass and velocity of the red component of the GW170817 kilonova (Kasen et al. 2017; Drout et al. 2017; Villar et al. 2017; Tanaka et al. 2017; Perego et al. 2017), we set M_ej = 0.04 M_⊙, and ${v}_{\mathrm{ej}}=2{({E}_{\mathrm{kin}}/{M}_{\mathrm{ej}})}^{1/2}=0.1c$, with c the speed of light.

Magnetic field configuration was shown in B16 to affect thermalization efficiency. Here, we assume the magnetic field lines are "tangled" as a result of turbulent processes early in the ejecta's expansion history. From a thermalization standpoint, this represents a middle ground between radial fields (which escort charged particles efficiently out of the ejecta, reducing the time during which the particles experience thermalizing interactions) and toroidal fields (which hold the charged particles in the ejecta interior, maximizing thermalization).

Much of the variability in f_α (t) and f_β (t) can be explained by differences in the radioactive energy generation rates, ${\dot{\epsilon }}_{{\rm{i}}}$, which are shown in the top middle panels of Figure 4, normalized to ${\dot{\epsilon }}_{{\rm{i}}}=1$ at t = 0.1 days to aid comparison. We find that a more steeply declining ${\dot{\epsilon }}_{{\rm{i}}}$ corresponds to a shallower decrease in f_i(t), and therefore to higher levels of thermalization. (Of course, because less energy is emitted in the steeply declining cases, the total amount of thermalized energy is still generally lower).

The effect of ${\dot{\epsilon }}_{{\rm{i}}}$ is strongest at later times, after thermalization has started to become inefficient. Prior to this, particles thermalize easily and f_i(t) hovers near unity, regardless of the rate of decline. In realistic models of the r-process, where the energy generation rate changes with time, the slope at later times will have a stronger effect on thermalization than the early-time decline.

This trend is illustrated by the β-particle thermalization of sly4.y18 (model 9). At early times, the β-heating of this model, ${\dot{\epsilon }}_{\beta ,9}$, falls rapidly, and f_β,9(t) is high compared to the full set of models. But by t ∼ 10 days, ${\dot{\epsilon }}_{\beta ,9}$ has flattened, becoming less steep than the β-heating rate of frdm.y28 (model 2). Around the same time, f_β,9(t) diverges from f_β,2(t), falling to lower values while the latter retains its shallow slope.

This behavior is expected from analytical models (e.g., Kasen & Barnes 2019, Equation (26)), and can be understood by considering the energy evolution of a population of particles generated by radioactive decays. The energy thermalized at any time t comes from a collection of suprathermal decay products emitted at times ≤t. The oldest unthermalized particles will have the lowest energies and the lowest thermalization times on average because to zeroth order, thermalization time scales with particle energy. When ${\dot{\epsilon }}_{{\rm{i}}}$ declines steeply, the fraction of the total unthermalized energy carried by older, rapidly thermalizing particles increases, which flattens the slope of f_i(t).

The decay rate for single-isotope heating, ${\dot{\epsilon }}_{{\rm{i}}}\propto {e}^{-t/{\tau }_{{\rm{r}}}}$, with τ_r the isotope lifetime, will generically be steeper than the rate from an ensemble of nuclei (Li & Paczyński 1998). Analytic thermalization models suggest that for exponential decay, f_i(t) will pass through a local minimum and begin to increase with time. While the sheer number of nuclei produced make truly exponential decay unlikely for the r-process, in some instances, the heating may be dominated by one or a handful of nuclei, and become approximately exponential in nature.

Nucleosynthesis in high-Y_e conditions burns a narrower distribution of nuclei, and so is more likely to become dominated by a single decay chain and transition into the regime of quasi-exponential decay. Indeed, the three models with the highest initial Y_e (models frdm.y28, unedf.y24, and sly4.y21, with Y_e = 0.28, 0.24, and 0.21, respectively) derive their energy from a much smaller set of nuclear decays than other models. However, the decrease and late-time increase in f_i(t) characteristic of exponential heating is found only in the β-particle thermalization of frdm.y28 (model 2). As shown in Figure 4, f_β,2(t) reaches a minimum of 0.84 at t = 36 days, then rises to 0.87 by t = 50 days, a value ≳70% greater than the next-highest model. This behavior is the result of the nearly exponential form of ${\dot{\epsilon }}_{\beta ,2}$.

This is shown in Figure 5, which plots the total absolute β-particle heating (dashed lines) for the three high-Y_e models, as well as the heating from individual isotopes (solid lines), all normalized to the total β-heating rate at t = 1.0 day. By t = 15 days, the heating for frdm.y28 (model 2) tracks the decay of ¹³¹I, which is responsible for ∼70%–80% of ${\dot{\epsilon }}_{\beta ,2}$ thereafter. The decay timescale of ¹³¹I is 8.03 days, and its production, in our simulations, occurs at t ≤ 1 day. Thus, by the time ¹³¹I dominates the total heating, it is well into its decay phase, and can impart to ${\dot{\epsilon }}_{\beta ,2}$ its steep exponential decline (see also the top third panel of Figure 4).

Zoom In Zoom Out Reset image size

The relationship between ${\dot{\epsilon }}_{{\rm{i}}}$ and f_i(t) slightly attenuates the variability in effective heating. Models with more steeply declining ${\dot{\epsilon }}_{{\rm{i}}}$ will produce less energy through radioactivity, but will thermalize that energy more effectively. However, as we show in Section 4.4, the effect is not strong enough to overcome the underlying differences in absolute heating.

While the form of ${\dot{\epsilon }}_{{\rm{i}}}$ explains some of the diversity in f_i(t), the magnitude and shape of the thermalization curves also depend on the energy spectra of the emitted particles, which varies among models and changes with time. The importance of the spectrum is particularly apparent for β-particle and γ-ray thermalization.

β -particles: The role of the emission spectrum in β-particle thermalization is most obvious for frdm.y28 (model 2). While the quasi-exponential shape of ${\dot{\epsilon }}_{\beta ,2}$ explains the unique falling and rising behavior of f_β,2(t), the sustained efficient thermalization is additionally due to an unusually low-energy β-spectrum. Figure 6 compares the later-time β-emission spectra of the high-Y_e models (frdm.y28, unedf.y24, and sly4.y21) and presents for all models the fraction of the total β-particle heating supplied by particles with E_β,0 ≤ 0.5 MeV as a function of time. The concentration of frdm.y28's β-spectrum at lower energies makes it an outlier even relative to the other high-Y_e models, which generally produce lighter nuclei with lower decay energies. This again is due to the predominance of ¹³¹I. Almost 90% of ¹³¹I decays proceed through the emission of β-particles whose maximum energy is 606 keV, which places the average energy of a β-particle from this decay at 〈E_β 〉 = 191 keV. These low-energy particles thermalize more quickly, enhancing f_β (t).

Zoom In Zoom Out Reset image size

Compared to other decay products, β-particles exhibit a high degree of variability in their thermalization. Our findings indicate that even in the regime of high-Y_e r-process nucleosynthesis, radioactivity and thermalization can vary substantially, and f_β (t) is sensitive to the particular constellation of isotopes synthesized.

γ -rays: Gamma-ray thermalization occurs via a series of distinct scattering and absorption events, rather than the continual slowing-down processes that mediate the thermalization of massive particles. The probability of scattering or being absorbed is a function of the optical depth, τ_γ = ∫ρ κ_γ ds, with ρ the mass density, κ_γ the γ-ray opacity, and s the path length.

At energies E_γ ≳ 0.5 MeV, κ_γ is set by the Compton-scattering opacity, which has a fairly low value (κ_Comp ≲ 0.1 cm² g⁻¹) and is relatively flat with increasing E_γ . For E_γ ≲ 0.5 MeV, photoionization becomes important, providing an opacity κ_PI orders of magnitude greater than κ_Comp, and rising sharply as E_γ decreases (Figure 3). As a result, γ-ray thermalization depends strongly on the emission spectrum of γ-rays; if most γ rays are emitted at energies for which κ_γ is high, f_γ (t) will be augmented.

To illustrate this, we calculate as a function of time a spectrum-weighted average γ-ray opacity,

$$\begin{eqnarray*}&&\langle \kappa \rangle =\sum {\kappa }_{\gamma }(E)\,\phi (E),\end{eqnarray*}$$

where ϕ(E) is the fraction of γ-ray energy emitted in the range E ≤ E_γ ≤ E + δ E. We show 〈κ〉 for all models in the top left panel of Figure 4.

While some of the model-to-model differences are due to composition (e.g., different abundance patterns resulting in slightly different net opacities), most of the variation among models, and all of the variation with time for any particular model, is due to differences in the γ-emission spectrum. Because κ_γ increases so steeply as E_γ falls below ∼0.5 MeV, 〈κ〉 is much more sensitive to the spectrum at low E_γ ; low-energy γ-rays can dominate 〈κ〉 and therefore thermalization, even when they represent only a small fraction of the total energy from γ-ray emission.

As anticipated, 〈κ〉 is strongly correlated with γ-ray thermalization (lower-left panel of Figure 4). To the extent that 〈κ〉 is not perfectly predictive of f_γ (t), the discrepancy is most likely due to a confluence of second-order effects, such as the rate of decline of ${\dot{\epsilon }}_{\gamma }$, or the strength of the γ-emission spectrum at energies just above the point where κ_PI begins to dominate the opacity. Photons emitted at these energies could down-scatter into the photoionization regime, and thereafter thermalize more efficiently than their initial energies would suggest. Because the ejecta is largely transparent to γ-rays, as indicated by the overall low efficiencies, even a small number of such down-scattering events could represent a relatively large fraction of all thermalizing interactions, and therefore impact f_γ (t).

Unlike γ-rays, β-particles, and α-particles, fission fragment thermalization is uniform across models, seeming to counter the conclusion that f_i(t) is sensitive to the details of radioactivity. However, fission fragments are in fact the exception that proves the rule: the results are consistent because the predictions of our models of fission radioactivity—despite the range of fission treatments considered—are quite uniform, and because we assume a single fission fragment mass and flat energy distribution for all models.

At early times, thermalization is generally efficient, particularly for rapidly thermalizing fission fragments. But at times late enough for f_ff(t) to be sensitive to the slope of ${\dot{\epsilon }}_{\mathrm{ff}}$, only two isotopes contribute to fission heating. The first, ²⁵⁴Cf, has a well-established half-life of 60.5 days (Phillips et al. 1963). It is present in all models save frdm.y28, which has no measurable fission. Its production by the r-process and its potential impact on kilonova light curves via thermalization were explored by Zhu et al. (2018). The second nucleus, ²⁶⁹Rf, appears for hfb22.y16 and hfb27.y16. There have been no experimental measurements of ²⁶⁹Rf, and its predicted fission half-life is sensitive to the theoretical models of nuclear masses and spontaneous fission lifetimes adopted. The HFB models for which ²⁶⁹Rf is found to contribute substantively to fission heating have fission barriers computed for HFB14 and half-lives calculated according to Karpov et al. (2012) and Zagrebaev et al. (2011). With these inputs, the spontaneous fission half-life of ²⁶⁹Rf is predicted to be very close to that of ²⁵⁴Cf: 58.95 days. (The longevity of these isotopes, along with the overall efficient thermalization of fission fragments, keep f_ff(t) from reaching a local minimum within the time limit of our simulation, despite ${\dot{\epsilon }}_{\mathrm{ff}}$ being effectively exponential).

Because nearly all fission after early times is due to these two isotopes with similar half-lives (Figure 7) and because we lack a detailed description of the fission fragment distribution and so use the same fission spectrum for each case (Section 3.2), late-time fission in every model is, modulo normalization, effectively identical. As a result, our results for f_ff(t) are similarly consistent.

Zoom In Zoom Out Reset image size

However, this consistency depends on the theoretical calculation of the fission half-life and branching ratios of allowed decay modes of ²⁶⁹Rf. For example, some models populate ²⁶⁹Rf, but find α-decay, rather than fission, to be its dominant decay mode. A thorough experimental investigation of potentially long-lived fissioning isotopes could change this picture and allow more accurate models of fission fragment thermalization and its effect on kilonova light curves.

The total thermalization efficiency from r-process decay depends on the efficiencies associated with individual particle types, as well as on the partition of radioactive energy into various decay channels. In this section we present the net thermalization efficiency of each model, as well as the corresponding effective heating rates for the associated kilonova model. Because thermalization depends on ejecta parameters, our results are valid only for the ejecta model described in Section 3.5.

Figure 8 shows the total thermalization efficiency for all models, highlighting the contribution of each particle type to the overall thermalization. Across all models, thermalization efficiencies range from 0.53 to 0.9 at early times, and decrease (though not necessarily monotonically) to between 0.1 and 0.5 at late times.

Zoom In Zoom Out Reset image size

The most conspicuous difference is between models that feature significant heating by α-decay and fission—which have higher overall levels of thermalization—and models that do not. Even frdm.y28's extremely high β-particle thermalization efficiency cannot compensate for its lack of rapidly thermalizing α-particles and fission fragments. Because the radioactivity in frdm.y28 (and unedf.y24 and sly4.y21) is dominated by β-decay, most of the energy is lost to neutrinos and γ-rays, and even when β-particles themselves thermalize efficiently, this sets a restrictive upper limit on total thermalization. In many cases where fission and α-decay grow in importance, the net thermalization efficiency flattens dramatically, or even begins to rise, around t = 10–20 days, as a result of the more efficient thermalization associated with these decay modes.

In the right-hand panel of Figure 9, we plot the rate at which radioactive energy is deposited in the ejecta (i.e., ${\dot{\epsilon }}_{\mathrm{eff}}(t)\,={M}_{\mathrm{ej}}\,{\dot{\epsilon }}_{\mathrm{abs}}(t)\times f(t)$) for all models. These curves exhibit a degree of variation comparable to that of the absolute heating rates (Figure 1), although the effects of model-specific thermalization have adjusted the shapes and magnitudes of the curves. Critically, for 0.5 days ≤ t ≤ 15 days, the time over which most kilonovae are expected to peak and fall, a difference of more than an order of magnitude separates the models with the lowest and highest heating. This indicates that the kilonovae corresponding to the models in our subset will have a similarly wide range of luminosities.

Zoom In Zoom Out Reset image size

We compute light curves for the ejecta model described in Section 3.5. We have modified Sedona so that the heating rate can be specified in time and space by a discrete data set $\{{\dot{\epsilon }}_{\mathrm{eff}}(v,t)\}$, where v is the ejecta velocity coordinate. This allows us to import our effective heating rates directly into Sedona, instead of relying on best-fit functional forms to approximate ${\dot{\epsilon }}_{\mathrm{eff}}$. It also enables us to capture the variation of ${\dot{\epsilon }}_{\mathrm{eff}}$ in space. This is important because, as discussed in Korobkin et al. (2020), regions of higher mass density will emit more radioactive energy and better thermalize that energy. This results in further stratification of the ejecta temperature: hot areas stay hotter, and relatively cool areas cool faster.

The composition of the outflow is a critical input for kilonova modeling (Even et al. 2020) because it sets the bound-bound opacity of the ejecta. Because our atomic data set does not include data for all species in the periodic table, we derive from our PRISM results a modified composition for use in our Sedona calculations. The major determinant of an element's opacity is its position in the periodic table: f-block elements have the highest opacity, followed by d-block, p-block, and s-block elements (Kasen et al. 2013; Tanaka & Hotokezaka 2013). However, opacity also varies within blocks. Elements near the center of a periodic table row contribute a higher opacity due to their greater atomic complexity (Tanaka et al. 2020).

The simplified compositions are constructed to preserve the distribution of the highest-opacity species within their periodic table row as closely as possible. We use the same atomic data set as in Kasen et al. (2017). It includes synthetic atomic data for all lanthanides except Lu (Z = 71), but lacks actinide data. In our model compositions, the mass fraction of any lanthanide is the mass fraction X^P of the element predicted by PRISM, plus the predicted mass fraction of the corresponding actinide. For example, ${X}_{\mathrm{Nd}}^{\mathrm{mod}}={X}_{\mathrm{Nd}}^{{\rm{P}}}+{X}_{{\rm{U}}}^{{\rm{P}}}$. Because we have no data for Lu, we send both Lu and its atomic structure homolog Lr (Z = 103) to Yb (Z = 70). The mass fractions of the d-block elements are split evenly among Z = 21–28, while Z = 20 stands in for the s- and p-block remainder. For these lighter elements, we use atomic data from the CMFGEN database (Hillier & Lanz 2001), as in Kasen et al. (2017).

Simple analytic considerations reveal how the peak time, t_pk, of a light curve depends on key properties of the outflow,

$$\begin{eqnarray}&&{t}_{\mathrm{pk}}\propto {\left({M}_{\mathrm{ej}}\kappa /{v}_{\mathrm{ej}}c\right)}^{1/2},\end{eqnarray} \tag{ 1 }$$

where κ is an effective (gray) opacity. In r-process ejecta, opacity is dominated by lanthanide and actinide elements (Kasen et al. 2013; Tanaka & Hotokezaka 2013), and κ increases with ${X}_{{\rm{L}},{\rm{A}}}^{}$. The scalings of Equation (1), along with the rule of thumb ${L}_{\mathrm{pk}}\approx {\dot{\epsilon }}_{\mathrm{eff}}({t}_{\mathrm{pk}})$ ("Arnett's Law"; Arnett 1980, 1982), define a straightforward relation between ejecta parameters, peak time, and peak luminosity.

Were ${\dot{\epsilon }}_{\mathrm{eff}}$ consistent across our model suite, this relationship would account entirely for variability in our models' light curves: with M_ej and v_ej held constant, light-curve evolution would depend only on composition, and models with higher ${X}_{{\rm{L}},{\rm{A}}}^{}$ would have later-peaking, dimmer light curves (e.g., Barnes & Kasen 2013).

As can be seen in the right-hand panel of Figure 9, our models do exhibit a range of peak times and luminosities. The diversity in predicted ${\dot{\epsilon }}_{\mathrm{eff}}$ (left-hand panel) propagates through to the light curves; luminosities for the single ejecta model we consider differ by ∼1 order of magnitude at peak and by even more on the light-curve tail. However, because our models vary both ${\dot{\epsilon }}_{\mathrm{eff}}$ and ${X}_{{\rm{L}},{\rm{A}}}^{}$ (re-plotted in the inset axes of Figure 9), they do not reproduce the expected correlation between t_pk, L_pk, and κ.

Most obviously, diversity in ${\dot{\epsilon }}_{\mathrm{eff}}$ disrupts the relation between t_pk and L_pk. The effect of Arnett's Law can still be seen in certain cases, although its influence is diminished relative to a scenario with uniform ${\dot{\epsilon }}_{\mathrm{eff}}$. For example, the comparatively high luminosity of sly4.y18 (model 9) at t ≲ 3 days is a product less of its effective heating than of its low opacity, which allows a larger fraction of ${\dot{\epsilon }}_{\mathrm{eff},9}$ to escape the ejecta on short timescales. However, for the most part, the light curves with the most luminous peaks are those with the most energetic ${\dot{\epsilon }}_{\mathrm{eff}}$, rather than the shortest rise times.

A more subtle effect is that of the shape of ${\dot{\epsilon }}_{\mathrm{eff}}$ on the kilonova rise time. While low-${X}_{{\rm{L}},{\rm{A}}}^{}$ cases (frdm.y28, sly4.y18, and sly4.y21) generally rise faster and sustain shorter photospheric phases than their high-${X}_{{\rm{L}},{\rm{A}}}^{}$ counterparts (e.g., unedf.kz.y16 and unedf.xr.y16), not all models adhere to the expectation that peak time increases with ${X}_{{\rm{L}},{\rm{A}}}^{}$. Model 5 (dz33.y16), for instance, peaks later and transitions from peak to tail more slowly than unedf.y24 (model 8), despite having a lower ${X}_{{\rm{L}},{\rm{A}}}^{}$ (${X}_{{\rm{L}},{\rm{A}}}^{5}=0.31\lt {X}_{{\rm{L}},{\rm{A}}}^{8}=0.38$).

This is due to the shallower slope of ${\dot{\epsilon }}_{\mathrm{eff},5}$. When ${\dot{\epsilon }}_{\mathrm{eff}}$ declines more slowly, the energy supplied by radioactivity more effectively offsets the loss, due to adiabatic expansion and diffusion of radiation, of energy deposited earlier, causing the light curve to evolve more slowly. (Heat from a flatter ${\dot{\epsilon }}_{\mathrm{eff}}$ also allows the ejecta to stay singly ionized longer, which enhances the opacity; see Section 5.3). If ${\dot{\epsilon }}_{\mathrm{eff}}$ remains shallow as the system becomes optically thin, it will also reduce the luminosity difference between the light-curve peak and tail. (The largest anomaly in the ${X}_{{\rm{L}},{\rm{A}}}^{}$–t_pk relation appears to be hfb22.y16, with ${X}_{{\rm{L}},{\rm{A}}}^{}=0.37$, and t_pk ≈ 3 days. However, the quick rise of hfb22.y16 results from a unique set of circumstances, as explained in Section 5.3).

To more clearly disentangle the effects of ${\dot{\epsilon }}_{\mathrm{eff}}$ and ${X}_{{\rm{L}},{\rm{A}}}^{}$ on light curves, we ran two sets of additional calculations. In the first, we adopted a standard r-process heating rate, ${\dot{\epsilon }}_{\mathrm{eff},\mathrm{std}}$, similar to that of Kasen et al. (2017), but retained the compositions predicted by PRISM. In the second, we fixed the composition and set ${X}_{{\rm{L}},{\rm{A}}}^{}=0.15$, but used ${\dot{\epsilon }}_{\mathrm{eff}}$ from our thermalization calculations. The light curves of the uniform heating (composition) models are shown in the left-hand (right-hand) panel of Figure 10, while the middle panel displays the fully consistent light curves of Figure 9. For Figure 10, we employ an alternate color scheme in which the colors of the curves reflect each model's original ${X}_{{\rm{L}},{\rm{A}}}^{}$.

Zoom In Zoom Out Reset image size

As the left-hand panel of Figure 10 demonstrates, when ${\dot{\epsilon }}_{\mathrm{eff}}$ is fixed, ${X}_{{\rm{L}},{\rm{A}}}^{}$ is tightly correlated with t_pk and L_pk. However, adjusting ${\dot{\epsilon }}_{\mathrm{eff}}$ (middle panel) confuses this pattern, and produces a more diverse set of light curves than would be expected from variation in ${X}_{{\rm{L}},{\rm{A}}}^{}$ alone. A comparison of the middle and right-hand panels shows that the bulk of the diversity in our model light curves—in terms of both brightness and light-curve morphology—is due to differences in ${\dot{\epsilon }}_{\mathrm{eff}}$, rather than ${X}_{{\rm{L}},{\rm{A}}}^{}$. The light curves of the middle panel, which reflect the full diversity in nucleosynthesis and radioactivity predicted by our PRISM and thermalization simulations, clearly exhibit the greatest variety and most complex relation between ${X}_{{\rm{L}},{\rm{A}}}^{}$, ${\dot{\epsilon }}_{\mathrm{eff}}$, and the light-curve evolution.

Effective heating ${\dot{\epsilon }}_{\mathrm{eff}}$ thus joins M_ej, v_ej, and κ as a critical determinant of kilonova light curves, which introduces new complexities into the modeling and interpretation of kilonova emission. For example, a plateau-shaped light curve could reflect a shallow-sloped ${\dot{\epsilon }}_{\mathrm{eff}}$ as much as the high M_ej, low v_ej, and/or high opacity that Equation (1) would suggest. Likewise, if ${\dot{\epsilon }}_{\mathrm{eff}}$ is higher than is typically assumed, a given luminosity may be achieved by a lower M_ej than would otherwise be required. In our models, higher ${\dot{\epsilon }}_{\mathrm{eff}}$ usually signals significant heating from fission, which is correlated with both greater ${X}_{{\rm{L}},{\rm{A}}}^{}$ (because fission requires robust heavy element nucleosynthesis) and flatter ${\dot{\epsilon }}_{\mathrm{eff}}$ due to the efficient thermalization of fission fragments and the long half-lives of the dominant fission reactions. One or both of these mechanisms could extend the light curve, compensating for the decrease in optical depth and t_pk caused by a lower M_ej. (However, the initial high heating rate of frdm.y28, which has a low ${X}_{{\rm{L}},{\rm{A}}}^{}$ and no appreciable fission or α-decay, bucks this trend, reflecting the dependence of radioactivity on a complex interplay of factors. We caution that no single heuristic can account for the full variety of nucleosynthesis outcomes).

The effective heating for some of our models is higher than in many earlier studies of kilonovae (e.g., Barnes & Kasen 2013; Barnes et al. 2016; Kasen et al. 2017), and the resulting hotter temperatures impact the ionization structure of the ejecta. The evolution of the ionization structure—and the increase in opacity that occurs when doubly ionized lanthanide species recombine—in some cases leaves an imprint on the kilonova emission.

Where bound-bound absorption dominates opacity, opacity is sensitive to temperature, peaking at ionization threshold temperatures and falling to local minima above these transition points. This is because bound-bound opacity depends strongly on the number of distinct bound-bound transitions ("lines") a photon encounters. For a system in local thermodynamic equilibrium (LTE; a reasonable assumption for kilonovae during the photospheric phase), atoms in a given ionization state populate a broader distribution of energy levels as temperature increases from the ionization threshold. At temperatures just above the threshold, atoms are confined to a handful of low-lying energy states, and there are comparatively few distinct lines contributing to the opacity (Kasen et al. 2013; Tanaka & Hotokezaka 2013; Tanaka et al. 2020).

Kilonova opacity is almost entirely from lanthanide and actinide ions, and therefore varies more strongly with their ionization state than with that of lower-opacity species. To illustrate this relation, we calculate as a function of temperature the average lanthanide ionization state, χ_ion,lanth., ⁷ and the mean opacity for conditions typical of our ejecta model (ρ = 10⁻¹⁴ g cm⁻³ and t = 5 days). We use the Rosseland mean, ${\bar{\kappa }}_{{\rm{R}}}$, and average over a wavelength-dependent opacity that includes free–free and bound-free absorption, electron-scattering, and bound-bound transitions. The latter we treat using the expansion opacity formalism (Karp et al. 1977; Eastman & Pinto 1993).

The results for frdm.y28 and hfb27.y16 (models 2 and 4, respectively) are presented in Figure 11. Although the two models have the highest and lowest ${X}_{{\rm{L}},{\rm{A}}}^{}$ of our subset (${X}_{{\rm{L}},{\rm{A}}}^{4}=0.48$, while ${X}_{{\rm{L}},{\rm{A}}}^{2}=0.02$), the evolution of ${\bar{\kappa }}_{{\rm{R}}}$ tracks χ_ion,lanth. in both cases. (However, ${X}_{{\rm{L}},{\rm{A}}}^{}$ does affect the magnitude of ${\bar{\kappa }}_{{\rm{R}}}$—note the different y-axis scales for each subplot). The opacity attains a maximum just as the lanthanide elements shift from singly to doubly ionized, but falls significantly once the transition is complete.

Zoom In Zoom Out Reset image size

Where lanthanides are nearly exactly doubly ionized, a low-opacity window opens up in the ejecta. Whether this window affects the observed luminosity depends on its location in the ejecta. In many models of kilonovae, particularly those with lower ${\dot{\epsilon }}_{\mathrm{eff},}$ the region where χ_ion,lanth. = 2 lies deep within the ejecta after very early times. The effective opacity, and therefore the diffusion timescale that controls the evolution of the light curve, is then set by the higher opacity of the singly ionized material outside the low-opacity region. On the other hand, if a higher ${\dot{\epsilon }}_{\mathrm{eff}}$ sustains until later times a spatially extended region where χ_ion,lanth. = 2, the increase in opacity that occurs as the gas cools and lanthanide ions recombine to χ_ion,lanth. = 1 can influence the effective opacity of the system as a whole.

While it is a simplification to describe optically thick, multiwavelength radiation transport in terms of any single parameter, we find the opacity at the surface where the Rosseland mean optical depth, ${\tau }_{{\rm{R}}}(r)\equiv -{\int }_{\infty }^{r}\rho {\bar{\kappa }}_{{\rm{R}}}{{dr}}^{{\prime} }$, is unity to be a useful proxy for the effective opacity of the ejecta. Specifically, the opacity at τ_R = 1 is correlated with the timescale on which the light curve evolves; an abrupt reduction in this opacity accelerates the light-curve evolution, while a sudden increase will slow it.

This can be understood if the τ_R = 1 surface is interpreted as a (rough) boundary separating an outer region where radiation escapes easily from an inner region where it diffuses out more slowly. A rising opacity necessarily moves this boundary outward, forcing a larger fraction on the thermalized energy to diffuse out more slowly, while falling opacity has the opposite effect.

Within our model subset, we identify four distinct ways that ionization state impacts light curves. These are summarized in Figure 12, which maps out χ_ion,lanth. and the position of the τ_R = 1 surface as functions of the normalized interior mass coordinate m_enc, defined as the fraction of M_ej that lies interior to a given ejecta layer. The bolometric light curves for the models presented have been re-plotted for convenience. For an analytic discussion of how changing opacity influences light curves, see Appendix B.

Zoom In Zoom Out Reset image size

Case A: Standard evolution. In the standard case, which is exemplified by unedf.xr.y16 in the first panel of Figure 12, sufficient cooling and a fairly high ${X}_{{\rm{L}},{\rm{A}}}^{}$ allow material with χ_ion,lanth. = 1 to set the effective opacity over essentially the entire light curve. The τ_R = 1 surface may skirt the χ_ion,lanth. = 2 region at very early times, when the entire ejecta is dense and multiply ionized, but at all times ≳a few hours, it tracks the singly ionized ejecta layer. In addition to unedf.xr.y16, unedf.kz.y16, unedf.y24, sly4.y21, and tf.y16 belong to this category. All (see Figure 9) have low ${\dot{\epsilon }}_{\mathrm{eff}}$ and ${X}_{{\rm{L}},{\rm{A}}}^{}\geqslant 0.12.$

Case B: Minor early excess. In models with slightly greater ${\dot{\epsilon }}_{\mathrm{eff}}$, represented by dz33.y16 in Figure 12 (panel 2), higher temperatures produce a doubly ionized layer near the ejecta edge, reducing the effective opacity of the system and allowing the light curve to rise sharply. However, ${\dot{\epsilon }}_{\mathrm{eff}}$ is not large enough to maintain this state for long, and the opacity increases as lanthanides recombine, helped along by a high ${X}_{{\rm{L}},{\rm{A}}}^{}$ that exacerbates the opacity difference between singly- and doubly ionized gas. The τ_R = 1 surface quickly jumps to the singly ionized region, signaling an increase in effective opacity that flattens the light-curve rise and creates the appearance of early excess luminosity. In addition to dz33.y16, early excesses are seen for frdm.y16 and especially hfb27.y16. All have higher ${\dot{\epsilon }}_{\mathrm{eff}}$ and greater ${X}_{{\rm{L}},{\rm{A}}}^{}$ ($0.22\leqslant {X}_{{\rm{L}},{\rm{A}}}^{}\leqslant 0.48$) than the models of Case A.

Case C: Extended decline phase. In the third case, a low ${X}_{{\rm{L}},{\rm{A}}}^{}$, perhaps assisted by a large ${\dot{\epsilon }}_{\mathrm{eff}}$, keeps the effective opacity low through and beyond the light-curve peak, even as the outer layers of the ejecta start to recombine. This is illustrated by frdm.y28 in the third panel of Figure 12. The lower opacity of the χ_ion,lanth. = 2 region causes the luminosity to peak sharply and fall off swiftly. However, the boundary where χ_ion,lanth. transitions from 2 to 1 recedes inward with time, and eventually there is enough singly ionized ejecta to drive the opacity up. The increasing opacity, which manifests in the abrupt displacement of the τ_R = 1 surface, slows the light-curve decline, producing a "shoulder" feature in the light-curve tail. (An analogous effect is responsible for the infrared rebrightening of Type Ia supernovae; Kasen 2006). This is seen for sly4.y18 as well as for frdm.y28.

Case D: Dual-phase light curve. Under certain conditions, recombination allows a light curve to reach two distinct bolometric peaks. This is the case for hfb22.y16, shown in the fourth panel of Figure 12 as the sole exemplar of Case D.

Early on, strong heating produces an extended highly ionized region, which sets a low initial effective opacity that supports a rapid light-curve rise. As indicated by the τ_R = 1 surface, this lower opacity is preserved through the first bolometric peak at t ≈ 3 days.

While the opacity of the doubly ionized region is low relative to singly ionized material, it is still fairly high in an absolute sense because of the high ${X}_{{\rm{L}},{\rm{A}}}^{}$ (=0.37) of hfb22.y16 (model 3). However, early peaks are possible even for high-opacity systems if ${\dot{\epsilon }}_{\mathrm{eff}}$ drops sharply (Section 5.2), and the decline of ${\dot{\epsilon }}_{\mathrm{eff},3}$ around the time of the first peak is particularly steep ($\mathrm{dlog}{\dot{\epsilon }}_{\mathrm{eff},3}/\mathrm{dlog}t\approx -2$ for 1 ≤ t/day ≤ 5; see Figure 9).

The interplay of ${X}_{{\rm{L}},{\rm{A}}}^{}$ and ${\dot{\epsilon }}_{\mathrm{eff}}$ can be understood by contrasting the light curves of the two HFB models. Both have high ${\dot{\epsilon }}_{\mathrm{eff}}$, but compared to hfb22.y16, hfb27.y16 has a slightly shallower early decline and a higher ${X}_{{\rm{L}},{\rm{A}}}^{}$ = 0.48. Thus, while hfb22.y16 forms a peak, hfb27.y16 merely exhibits a pronounced early excess. The situation changes if ${X}_{{\rm{L}},{\rm{A}}}^{}$ is reduced. This can be seen in the right-hand panel of Figure 10, which shows light curves calculated using our numerically determined ${\dot{\epsilon }}_{\mathrm{eff}}$ and a standard composition with ${X}_{{\rm{L}},{\rm{A}}}^{}=0.15$. The brightest light curve in that panel corresponds to hfb22.y16, and the second-brightest to hfb27.y16. The reduced opacity for this variant of hfb27.y16 allows an early peak to form, despite the fact that ${\dot{\epsilon }}_{\mathrm{eff}}$ declines less steeply for hfb27.y16 than for hfb22.y16. As these examples show, an early peak depends on ${\dot{\epsilon }}_{\mathrm{eff}}$ and opacity; the higher ${X}_{{\rm{L}},{\rm{A}}}^{}$, the steeper the ${\dot{\epsilon }}_{\mathrm{eff}}$ an early peak requires.

The second peak of hfb22.y16 is due to an increase in its opacity that coincides with a dramatic flattening of its ${\dot{\epsilon }}_{\mathrm{eff}}$. As can be inferred from Figure 12, for hfb22.y16, the opacity change from recombination is significant; it pushes the τ_R = 1 surface out almost to the edge of the ejecta, effectively resetting the ejecta to a condition similar to that of the pre-peak phases of Case A (standard) kilonovae.

At the same time, the heating rate becomes much shallower ($\mathrm{dln}{\dot{\epsilon }}_{\mathrm{eff},3}/\mathrm{dln}t\to -0.5$). This is the result of the growing importance of the long-lived fissioning isotopes ²⁵⁴Cf and ²⁶⁹Rf (Figure 7), which increases thermalization (Figure 8) and softens the decline of ${\dot{\epsilon }}_{\mathrm{eff}}$. The additional energy injection from the flatter ${\dot{\epsilon }}_{\mathrm{eff}}$ replenishes the internal energy, providing power for a second peak at t_pk,2nd ≈ 16 days.

As these four cases demonstrate, the effect of high heating rates and highly ionized regions in the ejecta exist on a continuum. While more than half our models are impacted in at least a minor way, the strongest effects are observed only for hfb22.y16 and hfb27.y16, which stand out for their high ${X}_{{\rm{L}},{\rm{A}}}^{}$ and extremal ${\dot{\epsilon }}_{\mathrm{eff}}$.

To illustrate the range of outcomes, we show in Figure 13 the position of the τ_R = 1 surface and the bolometric luminosity for all models. Increases in effective opacity, which appear in the plots as abrupt increases in the mass coordinate where τ_R = 1, are indeed correlated with apparent discontinuities in the light curves.

Zoom In Zoom Out Reset image size

Select optical and NIR broad-band light curves for each of our models are presented in Figure 14. There is some variation in the magnitude and morphology of the light curves, which parallels variability in bolometric luminosity. However, the evolution of the SEDs are overall consistent from model to model. Emission at optical wavelengths is suppressed after very early times, with much of the radiation instead emerging in the NIR. This behavior conforms to expectations established by earlier theoretical studies of lanthanide-rich kilonovae (Barnes & Kasen 2013; Tanaka & Hotokezaka 2013) and to observations of the red kilonova of GW170817 (Arcavi et al. 2017; Drout et al. 2017; Villar et al. 2017, and references therein).

Zoom In Zoom Out Reset image size

While the distinct spectra of kilonovae that are rich (${X}_{{\rm{L}},{\rm{A}}}^{}\gtrsim {10}^{-2}$) versus poor (${X}_{{\rm{L}},{\rm{A}}}^{}\lesssim {10}^{-3}$) in lanthanides and actinides has been well documented (e.g., Metzger & Fernández 2014; Tanaka et al. 2020), the effect on emission colors of very high ${X}_{{\rm{L}},{\rm{A}}}^{}$, like those considered here, has received less attention. Consistent with earlier studies, we find that generally speaking, lower ${X}_{{\rm{L}},{\rm{A}}}^{}$ is correlated with bluer colors. However, this correlation weakens with increasing ${X}_{{\rm{L}},{\rm{A}}}^{}$ and is further complicated by the effects—for some models—of an early low-opacity phase associated with a ∼doubly ionized ejecta. As a result, there is not necessarily an obvious observational fingerprint that can distinguish between moderately and extremely lanthanide-rich compositions, at least if uncertainties in ${\dot{\epsilon }}_{\mathrm{eff}}$ are taken into account.

The assumption that low ${X}_{{\rm{L}},{\rm{A}}}^{}$ leads to bluer colors is least reliable at early times, when even lanthanide-rich outflows may shine blue if high levels of ionization lower the effective opacity and allow the photosphere to form deep in the ejecta where temperatures are hotter. This is shown in Figure 15, which presents R–I and I–H colors for all models at t ≤ 5 days. We also show ${X}_{{\rm{L}},{\rm{A}}}^{}$ for each model in the top panel, while the bottom panel estimates how long the photosphere of each model, defined as the surface where τ_R = 1, remains in the χ_ion,lanth. = 2 region. The bluest models at this epoch include those with the lowest ${X}_{{\rm{L}},{\rm{A}}}^{}$ (e.g., frdm.y28, with ${X}_{{\rm{L}},{\rm{A}}}^{}$ = 0.02), and those with some of the highest (hfb22.y16, with ${X}_{{\rm{L}},{\rm{A}}}^{}$ = 0.37, and hfb27.y16, with ${X}_{{\rm{L}},{\rm{A}}}^{}$ = 0.48). In the latter case, ionization by large ${\dot{\epsilon }}_{\mathrm{eff}}$ compensates for the high ${X}_{{\rm{L}},{\rm{A}}}^{}$, lowering the opacity enough to place the photosphere at temperatures hotter than the first ionization threshold of lanthanides.

Zoom In Zoom Out Reset image size

This early phase is transient, however, and for nearly all models, the photosphere forms at χ_ion,lanth. = 1 for most of the kilonova duration. Because of this, the kilonova has fairly consistent colors at a given phase in its evolution. This is illustrated in Figure 16, which shows V–R, R–I, and I–H colors for all models as a function of time normalized to bolometric peak time. For hfb22.y16, we have normalized to the second bolometric peak, t_pk,2nd = 15.75 days.

Zoom In Zoom Out Reset image size

The models frdm.y28 (${X}_{{\rm{L}},{\rm{A}}}^{}$ = 0.02) and sly4.y18 (${X}_{{\rm{L}},{\rm{A}}}^{}$ = 0.03) appear as outliers in the left-hand column of Figure 16. The low opacity of these models produces an early light-curve peak and keeps the photosphere in the χ_ion,lanth. = 2 region through maximum light. As a result, their emission stays bluer out to later t/t_pk, although it eventually reddens.

The signature of this early highly ionized phase—blue emission that reddens as the ejecta recombines—may recall the photometric evolution of the kilonova associated with GW170817. The shift from bluer to redder emission observed in that event (e.g Arcavi et al. 2017; Chornock et al. 2017; Nicholl et al. 2017) was ascribed by many (Kasen et al. 2017; Kasliwal et al. 2017; Tanaka et al. 2017) not to ionization state, but to the presence of two distinct outflows with different levels of lanthanide enrichment. A full discussion of the implications of high heating for the interpretation of GW170817-like kilonovae is planned for future work. Here, we note just that none of the single-component models considered in this work reproduce all aspects of the electromagnetic emission observed in conjunction with GW170817.

Models frdm.y28 and sly4.y18 notwithstanding, the remaining models exhibit only subtle differences in color at a given phase. Model sly4.y21, shown in the middle column of Figure 16, has an intermediate ${X}_{{\rm{L}},{\rm{A}}}^{}$ (=0.12) and slightly bluer colors at t ≤ t_pk than the higher-${X}_{{\rm{L}},{\rm{A}}}^{}$ models of the right-hand column. However, its color evolution is consistent with that of tf.y16 (middle column), although tf.y16 has ${X}_{{\rm{L}},{\rm{A}}}^{}=0.22$, nearly twice that of sly4.y21. Conversely, frdm.y16 has ${X}_{{\rm{L}},{\rm{A}}}^{}$ equal to that of tf.y16, but shares the marginally redder colors of the models in the right-hand column, which span $0.22\,\leqslant {X}_{{\rm{L}},{\rm{A}}}^{}\leqslant 0.48$.

Despite this variance in ${X}_{{\rm{L}},{\rm{A}}}^{}$, the models of the right-hand column have R–I (I − H) colors at bolometric peak concentrated in the range 2.0–2.2 (3.8–4.1). The effects of ${X}_{{\rm{L}},{\rm{A}}}^{}$ on color are stronger for low ${X}_{{\rm{L}},{\rm{A}}}^{}$; the increase in ${X}_{{\rm{L}},{\rm{A}}}^{}$ from sly4.y18 (${X}_{{\rm{L}},{\rm{A}}}^{}=0.03$) to sly4.y21 (${X}_{{\rm{L}},{\rm{A}}}^{}=0.12$) raises R–I (I − H) at peak by 0.9 (≳1.6). This mainly reflects the formation of the photosphere in the χ_ion,lanth. = 1 region for sly4.y21, versus the χ_ion,lanth. = 2 region for sly4.y18. Broadly similar behavior over a wide range of high ${X}_{{\rm{L}},{\rm{A}}}^{}$ underscores the challenge of using color as a diagnostic of composition in the lanthanide-rich regime.

The root of this similarity is the strong dependence of opacity on temperature and the diminishing marginal impact of increasing ${X}_{{\rm{L}},{\rm{A}}}^{}$. As has been pointed out elsewhere (Barnes & Kasen 2013; Tanaka & Hotokezaka 2013, see also Figure 11), the sharp decrease in opacity as lanthanides and actinides transition from singly ionized to neutral ensures that under most conditions, the photosphere forms over a narrow range of temperatures around the first lanthanide ionization threshold. Changes to the opacity can nudge the photospheric temperature one way or the other, and the emerging pseudo-blackbody spectrum is modified by line-blanketing and absorption features. However, the effect of ${X}_{{\rm{L}},{\rm{A}}}^{}$ on the opacity, and thus on kilonova SEDs, is strongest for low ${X}_{{\rm{L}},{\rm{A}}}^{}$ (Kasen et al. 2013). Once the lanthanide concentration reaches some threshold, the impact of increasing ${X}_{{\rm{L}},{\rm{A}}}^{}$ is minimal because the opacity has been saturated. Strong lines are already absorbing so effectively that additional strength does not meaningfully increase their optical depth.

To illustrate this, we plot in Figure 17 the spectrum at bolometric peak for each of our models. Models frdm.y28 (2) and sly4.y18 (9) have the bluest peak spectra, which is not surprising given their low opacity (${X}_{{\rm{L}},{\rm{A}}}^{2}$ = 0.02 and ${X}_{{\rm{L}},{\rm{A}}}^{9}$ = 0.03) and the fact that their photospheric emission at peak originates from a hotter, more highly ionized ejecta layer. The other models, despite the large range of lanthanide concentrations represented ($0.12\leqslant {X}_{{\rm{L}},{\rm{A}}}^{}\leqslant 0.48$), show only minor spectral variation. (And even frdm.y28 and sly4.y18 evolve redward with time).

Zoom In Zoom Out Reset image size

Of course, the spectrum depends not only on ${X}_{{\rm{L}},{\rm{A}}}^{}$, but also on how ${X}_{{\rm{L}},{\rm{A}}}^{}$ is distributed among individual lanthanides and actinides. While a detailed accounting of the lines most important for spectral formation in kilonovae is a long-term project, theoretical work (Kasen et al. 2017; Even et al. 2020) suggests that certain elements may have more influence than others. Similarities in spectra (observed or simulated) may then point to comparable abundances of these dominant species, rather than identical total ${X}_{{\rm{L}},{\rm{A}}}^{}$. Both opacity saturation and the potential spectral dominance of few species suggest that a large range in ${X}_{{\rm{L}},{\rm{A}}}^{}$ can nonetheless produce spectra with similar shapes.

We start with the equation for conservation of energy in a simple one-zone, homologously expanding ejecta with a time-varying gray opacity, κ(t),

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{dE}}_{\mathrm{int}}}{{dt}}=\dot{Q}-\displaystyle \frac{{E}_{\mathrm{int}}}{t}-L\\ \,=\,\dot{Q}-\displaystyle \frac{{E}_{\mathrm{int}}}{t}-\displaystyle \frac{4\pi {vc}}{3M\kappa (t)}{E}_{\mathrm{int}}t,\end{array}\end{eqnarray} \tag{ B1 }$$

where E_int is the internal energy of the ejecta, $\dot{Q}$ is the heating from radioactivity (analogous to ${\dot{\epsilon }}_{\mathrm{eff}}$ in previous sections), and L is the luminosity. In the second line, which makes use of the diffusion equation, v is the ejecta velocity, M is the ejecta mass, and κ(t) is a time-dependent opacity.

We assume that opacity starts at some initial value κ_i and asymptotically approaches a value κ_f > κ_i. For simplicity, we express the opacity as κ(t) = κ_f k(t) and define a light-curve timescale t_lc in terms of κ_f: ${t}_{\mathrm{lc}}^{2}\equiv 3M{\kappa }_{{\rm{f}}}/4\pi {vc}$ (e.g., Arnett 1982). This allows us to rewrite L as

$$\begin{eqnarray}&&L=\displaystyle \frac{{E}_{\mathrm{int}}t}{{t}_{\mathrm{lc}}^{2}k(t)}.\end{eqnarray} \tag{ B2 }$$

Taking the derivative of L with respect to time provides a new expression for the time derivative of E_int. Accounting for the time-dependence of k, L evolves as

$$\begin{eqnarray}&&\displaystyle \frac{{dL}}{{dt}}=\displaystyle \frac{1}{{t}_{\mathrm{lc}}^{2}}\left[\displaystyle \frac{{{dE}}_{\mathrm{int}}}{{dt}}\displaystyle \frac{t}{k}+\displaystyle \frac{{E}_{\mathrm{int}}}{k}-\displaystyle \frac{{E}_{\mathrm{int}}t}{{k}^{2}}\displaystyle \frac{{dk}}{{dt}}\right]\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}&&\to \displaystyle \frac{{{dE}}_{\mathrm{int}}}{{dt}}=\displaystyle \frac{{t}_{\mathrm{lc}}^{2}k}{t}\displaystyle \frac{{dL}}{{dt}}-\displaystyle \frac{{E}_{\mathrm{int}}}{t}\left(1-\displaystyle \frac{t}{k}\displaystyle \frac{{dk}}{{dt}}\right).\end{eqnarray} \tag{ B4 }$$

When substituted into Equation (B1), this yields, after some rearranging,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{t}_{\mathrm{lc}}^{2}k}{t}\displaystyle \frac{{dL}}{{dt}}+\displaystyle \frac{{E}_{\mathrm{int}}}{k}\displaystyle \frac{{dk}}{{dt}}+L=\dot{Q}\\ \quad \displaystyle \frac{{dL}}{{dt}}+L\left(\displaystyle \frac{1}{k}\displaystyle \frac{{dk}}{{dt}}+\displaystyle \frac{1}{{t}_{\mathrm{lc}}^{2}}\displaystyle \frac{t}{k}\right)=\displaystyle \frac{\dot{Q}t}{{t}_{\mathrm{lc}}^{2}k}.\end{array}\end{eqnarray} \tag{ B5 }$$

Equation (B5) is an ordinary linear differential equation and can be solved with an integrating factor u(t) once the form of k(t) is known. (This is the motivation for making k a function of t directly, rather than E_int). We choose

$$\begin{eqnarray}&&k(t)=\displaystyle \frac{{t}^{2}+{{bt}}_{\kappa }^{2}}{{t}^{2}+(1+b){t}_{\kappa }^{2}},\end{eqnarray} \tag{ B6 }$$

which has the desired limiting behavior. In this formulation, the constant b is determined by the requirement that k(t = 0) = κ_i/κ_f, while t_κ sets the timescale of the transition from κ_i to κ_f. We show k(t) for a few sets of parameters b and t_k in Figure 19.

Zoom In Zoom Out Reset image size

With this choice, the second term in parentheses in the left-hand side of Equation (B5) becomes

$$\begin{eqnarray*}&&\displaystyle \frac{1}{{t}_{\mathrm{lc}}^{2}}\displaystyle \frac{t}{k}=\displaystyle \frac{1}{{t}_{\mathrm{lc}}^{2}}\left[t+\displaystyle \frac{t}{{\left(t/{t}_{\kappa }\right)}^{2}+b}\right].\end{eqnarray*}$$

This allows the integrating factor to be expressed analytically,

$$\begin{eqnarray}\begin{array}{rcl}u(t) & = & \exp [\int P(t){dt}]\\ P & = & \displaystyle \frac{1}{k}\displaystyle \frac{{dk}}{{dt}}+\displaystyle \frac{1}{{t}_{\mathrm{lc}}^{2}}\left[t+\displaystyle \frac{t}{{\left(t/{t}_{\kappa }\right)}^{2}+b}\right]\end{array}\end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray}&&\to u(t)=k(t){\left({t}^{2}+{{bt}}_{\kappa }^{2}\right)}^{({t}_{\kappa }^{2}/2{t}_{\mathrm{lc}}^{2})}\exp [{t}^{2}/2{t}_{\mathrm{lc}}^{2}].\end{eqnarray} \tag{ B8 }$$

We can now write an expression for the time-dependent luminosity,

$$\begin{eqnarray}&&\begin{array}{l}L(t)=k{\left(t\right)}^{-1}\,{\left({t}^{2}+{{bt}}_{\kappa }^{2}\right)}^{(-{t}_{\kappa }^{2}/2{t}_{\mathrm{lc}}^{2})}\,\exp [-{t}^{2}/2{t}_{\mathrm{lc}}^{2}]\\ \,\times \,\left[{L}_{0}({\kappa }_{{\rm{i}}}/{\kappa }_{{\rm{f}}}){\left({{bt}}_{\kappa }^{2}\right)}^{({t}_{\kappa }^{2}/2{t}_{\mathrm{lc}}^{2})}\right.\\ \,+\ \left.{\displaystyle \int }_{0}^{t}\displaystyle \frac{\dot{Q}({t}^{{\prime} }){t}^{{\prime} }}{{t}_{\mathrm{lc}}^{2}}\,{\left({{t}^{{\prime} }}^{2}+{t}_{\kappa }^{2}b\right)}^{({t}_{\kappa }^{2}/2{t}_{\mathrm{lc}}^{2})}\,\exp [{{t}^{{\prime} }}^{2}/2{t}_{\mathrm{lc}}^{2}]\,{{dt}}^{{\prime} }\right].\end{array}\end{eqnarray} \tag{ B9 }$$

Equation (B9) becomes slightly more palatable if expressed in terms of dimensionless units. Defining τ = t/t_lc and τ_κ = t_κ /t_lc, we arrive at a final expression for luminosity,

$$\begin{eqnarray}&&\begin{array}{l}L(\tau )=k{\left(\tau \right)}^{-1}\,{\left({\tau }^{2}+b{\tau }_{\kappa }^{2}\right)}^{(-{\tau }_{\kappa }^{2}/2)}\,\exp [-{\tau }^{2}/2]\\ \,\times \ \left[{L}_{0}({\kappa }_{{\rm{i}}}/{\kappa }_{{\rm{f}}}){\left(b{\tau }_{\kappa }\right)}^{({\tau }_{\kappa }^{2}/2)}\right.\\ \,+\,\left.{\displaystyle \int }_{0}^{\tau }\dot{Q}({\tau }^{{\prime} }){\tau }^{{\prime} }{\left({\tau }^{{\prime} }+{\tau }_{\kappa }^{2}b\right)}^{({\tau }_{\kappa }^{2}/2)}\ \exp [{{\tau }^{{\prime} }}^{2}/2]d{\tau }^{{\prime} }\right].\end{array}\end{eqnarray} \tag{ B10 }$$

In practice, the initial luminosity L₀ will be negligible compared to the luminosity at t > 0 once radioactive heating has begun, and the approximation L₀ = 0 can simplify Equation (B10) further.

The light curves described by Equation (B10) will depend on the form of $\dot{Q}$, as well as t_κ and the ratio κ_i/κ_f. We start by considering simple power-law heating, $\dot{Q}\propto {\tau }^{-1}$, to illustrate the effects of the opacity parameters. Because our interest is in the shape of the light curves, we do not normalize $\dot{Q}$, and the luminosity scales in our models are therefore not physically meaningful. In all calculations in this section, we have assumed that the contributions of radioactive heating dominate the initial luminosity of the system, and set L₀ = 0.

In Figure 20 we show the effect on light curves of κ_i/κ_f (left column) and the timescale τ_κ (right column). Light curves are shown in the top row, while the associated k(τ) are plotted below. For comparison, we show as a gray line the light curve produced for a constant k(t) = 1 (i.e., κ(τ) = κ_f). As can be seen in the left-hand panels, the impact on light curves of an early low opacity is enhanced luminosity during the light-curve rise. The smaller κ_i/κ_f, the more pronounced the enhancement, and the more obvious the effect of the increasing opacity at τ ∼ τ_κ = 0.3. From the model with κ_i/κ_f = 0.1, we see that double-peaked light curves are possible for this choice of $\dot{Q}$, but require a fairly extreme opacity increase. All models in Figure 20 peak globally at τ ≈ 1.

Zoom In Zoom Out Reset image size

The models in the right-hand column tell a similar story. The larger t_κ , the longer the system maintains its early low opacity, and as a result, the brighter the light-curve rise relative to the case with constant κ = κ_f. When the opacity transition occurs prior to peak, the new, longer diffusion time induced by the higher κ_f reduces the slope of the light-curve rise. We note that for τ_κ ≳ 1, κ_i is more characteristic of the light-curve opacity than κ_f, and scaling to t_lc (defined in terms of κ_f) via the τ parameter becomes less logical. (This is the reason the τ_κ = 1 model in the top right panel of Figure 20 appears to have an earlier, more luminous peak than the other light curves).

As discussed above (Section 2 and Section 4.4), r-process heating across large swaths of the nuclear physics parameter space is not well represented by a single power-law. We have also argued (Section 5.3) that changes in the heating rate may support the formation of double-peaked light curves. This motivates us to consider broken power-law forms of $\dot{Q}$,

$$\begin{eqnarray*}\begin{array}{rcl}\dot{Q}(\tau ) & \propto & {\tau }^{-\zeta }\\ \quad \zeta & = & \left\{\begin{array}{ll}{\zeta }_{1}, & \tau \leqslant {\tau }_{{\rm{Q}}}\\ {\zeta }_{2}, & \tau \gt {\tau }_{{\rm{Q}}}\end{array}\right..\end{array}\end{eqnarray*}$$

Figure 21 shows model light curves for a handful of broken power laws $\dot{Q}$, all with τ_Q = 0.3 and ζ₁ = 1, for constant and time-evolving opacity. All models with a time-dependent κ(τ) have κ_i/κ_f = 0.1 and τ_κ = 0.5. As can be seen in the middle panel, when κ is constant, changes to $\dot{Q}$ do not lead to any early peaks or other surprising features, although the flattening of $\dot{Q}$ does lead to brighter, later peaks. The former is expected based on the additional energy supplied by a flatter $\dot{Q}$, while the latter is due to the improved ability of a shallow $\dot{Q}$ to compensate for energy lost to diffusion and adiabatic expansion (Section 5.2).

Zoom In Zoom Out Reset image size

When κ increases in time, changes to $\dot{Q}$ strongly affect the light curve, as shown in the right-hand panel of Figure 21. All models in this panel have an early peak, enabled by a low κ_i/κ_f. However, whether that peak is the only peak produced, a global peak in a dual peaked light-curve, or a local peak followed by a global peak depends on how dramatically $\dot{Q}$ flattens. For ζ₂ = 1/2 (3/4), enough energy is injected into the system to support a second rise to a global (local) maximum. For ζ₂ = ζ₁ = 1, no secondary peak is present, and the rising opacity merely produces a shoulder feature, similar to that observed in some of our Sedona models (Section 5.3).

Allowing $\dot{Q}$ to vary further expands the parameter space under consideration and introduces a new timescale, τ_Q, into our analysis. Our model light curves now depend on the change in opacity, κ_i/κ_f; τ_κ , the timescale of the opacity transition; the slope(s) of $\dot{Q}$; and the timescale τ_Q on which any changes to that slope occur.

While we do not attempt to map out the full parameter space described above, we explore a limited slice of it in Figure 22 in order to illustrate the most important trends and highlight how certain model parameters interact. All of the models in Figure 22 have a broken power law $\dot{Q}$. For most $\dot{Q}$, we have set the first power-law index ζ₁ = 2. This is steeper than in Figure 21, but is a good approximation to the effective heating rate of the model (hfb22.y16) that produces a double-peaked light curve in our numerical calculations. However, this steep heating rate is not expected to characterize very early-time r-process decay (see Figure 1 of this work, or e.g., Korobkin et al. 2012; Lippuner & Roberts 2015; or Wu et al. 2019). To avoid overheating at early times as $\dot{Q}$ is extrapolated back toward τ = 0, we force $\dot{Q}$ at τ < 0.05 to decay no more steeply than τ⁻¹.

Zoom In Zoom Out Reset image size

The $\dot{Q}$s for each model are plotted in the top panels of Figure 22, while the corresponding light curves are presented below. For all models, we have adopted a moderate opacity ratio κ_i/κ_f = 0.2, although (Figure 20) a lower value would allow double-peaked light curves in a larger region of the parameter space. The left and middle columns show the effect of varying the timescales τ_κ and τ_Q. In the left-hand column, τ_Q = 0.3, while τ_κ is a free parameter, while in the middle column, τ_κ is set to 0.3, and we vary τ_Q.

Both panels suggest that the appearance of two well-defined peaks depends on having 0.1 ≲ τ_κ ≈ τ_Q ≲ 0.6. If either timescale is too low, the system does not spend enough time in the regime of steeply declining $\dot{Q}$ and/or low opacity for these early conditions to leave a strong imprint on the light curve, which simply evolves in accordance with κ_f and ζ₂. This is the case for the violet curve in the lower-left panel and the orange curve in the lower-middle panel. On the other hand, if κ or $\dot{Q}$ transitions too slowly, the light-curve evolution, determined by κ_i and ζ₁, is largely complete by the time the shifts take place. In these instances, represented by the brown and gray curves in the lower-left and lower-middle panels, respectively, the increasing opacity and/or the flattening of $\dot{Q}$ impact the tail of the light curve, but do not lead to the formation of a new peak.

In the right-hand panels, we explore the impact of ζ₁ and ζ₂ on the light curves while holding τ_κ = τ_Q = 0.3. While all models in the lower right panel show some sort of early feature, the data in this panel suggest that the combination of large ζ₁ and small ζ₂ favor the formation of two peaks. This is expected. A steeper decline early on accelerates the timescale for light-curve evolution (Section 5.2), allowing the first peak to fully form, while the much flatter decline at τ > τ_Q provides the additional energy needed to power a second peak.

While one-zone, gray-opacity models cannot capture the complexities of realistic numerical simulations of transients, such as those we present in Section 5, they are useful for diagnosing trends and understanding important determinants of emission. We have shown that double-peaked light curves are expected analytically when opacity increases in time and the heating rate deviates from a straightforward power law. As we argued above, the r-process—at least within the limits of current nuclear physics uncertainties—can in some cases meet both of these criteria. This raises the possibility that kilonova light curves may be much more diverse than previously expected.