Physics of Luminous Transient Light Curves: A New Relation Between Peak Time and Luminosity

TL;DR

This paper compares analytic models of supernova light curves with numerical simulations, revealing limitations of existing models and introducing new relations that better connect peak luminosity and time to physical parameters.

Contribution

It introduces new analytic relations that accurately link peak light curve features to ejecta and heating parameters, improving upon previous models.

Findings

Neglect of time-dependent diffusion limits model accuracy.

Arnett's rule often overestimates radioactive mass.

New relations incorporate effects of recombination and heating distribution.

Abstract

Simplified analytic methods are frequently used to model the light curves of supernovae and other energetic transients and to extract physical quantities, such as the ejecta mass and amount of radioactive heating. The applicability and quantitative accuracy of these models, however, have not been clearly delineated. Here we carry out a systematic study comparing certain analytic models to numerical radiation transport calculations. We show that the neglect of time-dependent diffusion limits the accuracy of common Arnett-like analytic models, and that the widely-applied Arnett's rule for inferring radioactive mass does not hold in general, with an error that increases for models with longer diffusion times or more centralized heating. We present new analytic relations that accurately relate the peak time and luminosity of an observed light curve to the physical ejecta and heating parameters. We further show that recombination and the spatial distribution of heating modify the peak of the light curve and that these effects can be accounted for by varying a single dimensionless parameter in the new relations. The results presented should be useful for estimating the physical properties of a wide variety of transient phenomena.

Figures

Tables

Table 1. Table 1: Peak Time-Luminosity Relation for Generic Heating Sources

Heating Type	Source Function $\mathscr{H}(t,t_s)$	$f(\tau ,\beta )$
Exponential	$\exp \left[-t/t_s\right]$	$1-\left(1+\beta \tau \right)\exp \left[-\beta \tau \right]$
Magnetar $n>2$	$\frac{n-1}{{\left(1+t/t_s\right)}^n}$	$\frac{1}{n-2}\left[1-\left(1+\left(n-1\right)\beta \tau \right){\left(1+\beta \tau \right)}^{1-n}\right]$
Magnetar $n=2$	${\left(1+t/t_s\right)}^{-n}$	$\ln \left(1+\beta \tau \right)-{\left(1+1/\beta \tau \right)}^{-1}$
Constant w/ Shutoff	$\Theta \left(t-t_s\right)$	$1/2$
Impulse	$E_0\delta \left(t-t_s\right)$	$1$
Power-Law $\alpha \ne 2$	${\left(\frac{t}{t_s}\right)}^{-\alpha }$	$\frac{1}{2-\alpha }\left[{\left(\beta \tau \right)}^{2-\alpha }-1\right]$
Power-Law $\alpha =2$	${\left(\frac{t}{t_s}\right)}^{-2}$	$\ln \beta \tau$

Table 2. Table 2: Values of β 𝛽 \beta for specific transients

Transient	$\beta$
Generic - Central source, constant opacity	$4/3$
Generic - Mixed source, constant opacity	$2$
Generic - Central source, Hydrogen recombination	$0.94$
Generic - Central source, Helium recombination	$0.7$
Type Ia SNe	$1.6$
Type Ibc SNe	$9/8$
Type IIb/pec SNe	$0.82$

Equations

$$\dot{\mathcal{E}}=-P\dot{\mathcal{V}}+\epsilon -\frac{\partial L}{\partial m}$$

$$R_{\mathrm{ej}}(t)=v_{\mathrm{ej}}t$$

$$\mathcal{E}=3P/\rho =a{T}^4/\rho$$

$$L(r)=-4\pi r^2\frac{c}{3\kappa \rho }\frac{\partial e}{\partial r}$$

$$e(x,t)=\frac{E_{\mathrm{int}}(t)}{V(t)}\psi (x)$$

$$L=\frac{t{E}_{\mathrm{int}}(t)}{{\tau }_d^2}$$

$${\tau }_d={\left[\frac{3}{4\pi }\frac{\kappa M_{\mathrm{ej}}}{v_{\mathrm{ej}}c}\frac{1}{\xi }\right]}^{1/2}$$

$$e(x,t)=\frac{E_{\mathrm{int}}(t)}{V(t)}\left[\frac{\pi }{3}\frac{\sin (\pi x)}{x}\right]$$

$$\frac{d{E}_{\mathrm{int}}(t)}{dt}=-\frac{E_{\mathrm{int}}(t)}{t}+L_{\mathrm{heat}}(t)-L(t)$$

$$\frac{{\tau }_d^2}{t}\frac{dL}{dt}=L_{\mathrm{heat}}(t)-L(t)$$

$$L(t)=\frac{2}{{\tau }_d^2}{e}^{-t^2/{\tau }_d^2}{\int }_0^t{t}^{\prime }{L}_{\mathrm{heat}}(t^{\prime })e^{t^{\prime 2}/{\tau }_d^2}\mathrm{d}{t}^{\prime }$$

$$t_d=\sqrt{\frac{\kappa M_{\mathrm{ej}}}{v_{\mathrm{ej}}c}}$$

$$⟨x_s⟩={\left(\frac{{\int }_0^1{x}^2{\dot{e}}_{\mathrm{heat}}(x)dx}{{\int }_0^1{\dot{e}}_{\mathrm{heat}}(x)dx}\right)}^{1/2}$$

$$L_{\mathrm{peak}}=L_{\mathrm{heat}}(t_{\mathrm{peak}})$$

$$\frac{d(tE)}{dt}=t[L_{\mathrm{heat}}(t)-L(t)]$$

$$tE(t)={\int }_0^t{t}^{\prime }{L}_{\mathrm{heat}}(t^{\prime })d{t}^{\prime }-{\int }_0^t{t}^{\prime }L(t^{\prime })d{t}^{\prime }$$

$$\frac{t^2}{2}{L}_{\mathrm{peak}}={\int }_0^t{t}^{\prime }{L}_{\mathrm{heat}}(t^{\prime })d{t}^{\prime }+ϵ(t)$$

$$ϵ(t)=\left[\frac{t^2}{2}{L}_{\mathrm{peak}}-{\int }_0^t{t}^{\prime }L(t^{\prime })d{t}^{\prime }\right]-tE(t)$$

$$L_{\mathrm{peak}}=\frac{2}{{\beta }^2{t}_{\mathrm{peak}}^2}{\int }_0^{\beta t_{\mathrm{peak}}}{t}^{\prime }{L}_{\mathrm{heat}}(t^{\prime })d{t}^{\prime }$$

$$L_{\mathrm{heat}}(t)=L_0{e}^{-t/t_s}$$

$$L_{\mathrm{peak}}=\frac{2L_0{t}_s^2}{{\beta }^2{t}_{\mathrm{peak}}^2}\left[1-(1+\beta t_{\mathrm{peak}}/t_s)e^{-\beta t_{\mathrm{peak}}/t_s}\right]$$

$$L_{\mathrm{heat}}(t)=\frac{L_0}{(1+t/t_s{)}^2}$$

$$\frac{t_{\mathrm{peak}}}{t_d}=0.11\ln \left(1+\frac{9t_s}{t_d}\right)+0.36$$

$$\beta (x_s)\approx \frac{4}{3}\left(1+x_s^4\right)$$

$$\kappa (T)={\kappa }_0+\frac{{\kappa }_0-ϵ{\kappa }_0}{2}\left[1+\tanh \left(\frac{T-T_{\mathrm{ion}}}{\Delta T_{\mathrm{ion}}}\right)\right].$$

$$L_{\mathrm{ion}}\approx 4\pi R_{\mathrm{ej}}^2{\sigma }_{\mathrm{sb}}{T}_{\mathrm{ion}}^4$$

$$L\sim \frac{E_s{t}_s}{t_d^2}$$

$$\eta$$

$$\eta \sim 0.2E_{s,51}{t}_{s,10}{\kappa }_{0.2}^{-2}{M}_{10}^{-2}{T}_4^{-4}$$

$$L_{\mathrm{peak}}=\frac{2}{{\beta }^2{t}_{\mathrm{peak}}^2}{\int }_0^{\beta t_{\mathrm{peak}}}{t}^{\prime }{L}_{\mathrm{heat}}(t^{\prime })d{t}^{\prime }$$

Full text

radiative transfer – supernovae: general

††software: Sedona (Kasen et al., 2006), Matplotlib, NumPy

1 Introduction

Wide-field surveys are gathering data on an increasing number of common supernovae (SNe) and related transients such as tidal disruption events, fast-evolving luminous transients, superluminous supernovae, and kilonovae. The general physics controlling the light curves of these events is similar – energy deposited either by a propagating shock (e.g. Type II-P SNe) or a heating source (e.g. radioactivity or a central engine) radiatively diffuses through the optically thick and expanding ejecta, undergoing adiabatic loses until the radiation reaches the surface and escapes. Analysis of the observed light curves can provide information on the ejecta properties and the nature of the powering source.

With the increasing number of observed transients, there is increased need for fast, empirical techniques to infer their physical properties and discriminate between competing theoretical explanations. Simplified analytic models are commonly used to analyze observations and make theoretical predictions (e.g. Li & Paczyński (1998); Chatzopoulos et al. (2012); Inserra et al. (2013); Villar et al. (2017); Nicholl et al. (2017); Guillochon et al. (2018)). Most common among these are the Arnett models (Arnett, 1980, 1982), in which bolometric light curves are calculated through a simple numerical integral. These models also provide several “rules of thumb” for estimating physical properties from the light curve brightness and duration. In particular, “Arnett’s rule” states that the instantaneous heating rate at peak is equal to the peak luminosity (Arnett, 1982). For Type I supernovae, this rule in principle allows one to extract the mass of radioactive 56Ni (Stritzinger et al., 2006; Valenti et al., 2007; Drout et al., 2011; Scalzo et al., 2014; Prentice et al., 2016).

Despite the frequent application of these analytic models and rules, a systematic study of their accuracy and applicability has not been carried out. Previous numerical models of Type Ia SNe have noted that Arnett’s rule is usually accurate to $\sim 20\%$ (Blondin et al., 2013; Hoeflich et al., 2017; Khokhlov et al., 1993). For Type Ib/c SNe Arnett’s rule is typically off by $\sim 50\%$ (Dessart et al., 2015, 2016), and for Type II SNe like SN1987A the error is a factor of $\sim 2$ (Woosley, 1988). The reasons for these discrepancies – and why they are more extreme for certain classes of transients – have not been fully spelled out.

In this paper we carry out a systematic investigation of certain analytic models compared to numerical light curve simulations. We find that the main limitations in the Arnett models stem from the assumption that a self-similar temperature profile is immediately established in the ejecta. This fails to account for the time-evolving propagation of a radiative “diffusion wave” from the heating source to the surface. The neglect of the diffusion wave is worse for more centrally concentrated heating sources, which is why Arnett’s rule is worse for more stratified Type Ib/c SNe and better for more thoroughly mixed Type Ia SNe.

We derive a new relation between the peak time and peak luminosity of transient light curves which accurately captures the results of numerical models. We study how the relation depends on the spatial distribution of heating as well as the effects of a non-constant opacity due to recombination. The new relation is parameterized by a dimensionless constant $\beta$, and works for a variety of assumed heating sources and ejecta characteristics.

In Section 2, we describe the assumptions and limitations of the Arnett light curve models and Arnett’s rule. In Section 3, we derive the new peak-time luminosity relation and compare it to numerical simulations. In Section 4, we investigate the relation between the peak and diffusion time. In Section 5, we examine the effects of concentration of the heating source. In Section 6, we look into the effects of a non-constant opacity due to recombination on the light curve and the new relation. Finally, in Section 7, we apply the results and relation to radioactive 56Ni-powered transients. In Appendix A, we provide a table of the new relation for a variety of luminous transients.

2 Limitations of Arnett-like Models

The analytic light curve modeling approach of Arnett (1980, 1982) (hereafter A80,A82) is widely used to analyze luminous transients. A closely related “one-zone” modeling approach (Arnett, 1979; Kasen & Bildsten, 2010; Villar et al., 2017) differs in its mathematical details but results in a similar expression for calculating the light curve.

2.1 Assumptions and Light Curve Solution

The Arnett-like models begin with the first law of thermodynamics

[TABLE]

where $\mathcal{E}$ is the specific (i.e. per unit mass) energy density, $P$ the pressure, $\mathcal{V}=1/\rho$ is the specific volume, $\varepsilon$ is the specific heating rate, and $L$ the emergent luminosity. Several simplifying assumptions are then made: (1) the ejecta is expanding homologously and so the radius evolves as

[TABLE]

where $v_{\rm ej}$ is the maximum ejecta velocity; (2) Radiation pressure dominates over the gas pressure and so we can express the specific energy density as

[TABLE]

where $T$ is the temperature; (3) The luminosity is described by the spherical diffusion equation

[TABLE]

where $e=\rho\mathcal{E}$ is the energy density (per unit volume) and $\kappa$ the opacity; and (4) the ejecta is characterized by a constant opacity.

The Arnett models make an additional consequential, but often overlooked, assumption: (5) The energy density profile is self-similar, i.e., the spatial dependence is fixed and only the overall normalization changes with time

[TABLE]

where $E_{\rm int}$ is the total internal energy of the ejecta and $x=r/R_{\rm ej}(t)$ is the (comoving) dimensionless coordinate. The dimensionless function $\psi(x)$ describes the spatial dependence of the radiation energy density, which by assumption does not change with time. Substituting Eq. 5 into the diffusion equation (Eq. 4) gives the emergent luminosity at $r=R_{\rm ej}$

[TABLE]

where

[TABLE]

is the characteristic diffusion time through the ejecta with mass $M_{\rm ej}$. The quantity $\xi=d\psi/dx|_{x=1}$ specifies the energy density gradient at the ejecta surface and is a constant when self-similarity is assumed. The one-zone models make the ansatz $\xi=1$. Arnett (1982) uses a more sophisticated separation of variables method to derive a self-consistent solution for $e(x,t)$. This requires making a final assumption: (6) The spatial distribution of the heating is proportional to the energy density. Eq. 1 can then be solved to find

[TABLE]

which gives $\xi=\pi^{2}/3$.

To solve for the light curve, the Arnett models integrate Eq. 1 over the entire ejecta and apply the assumptions of homology and radiation energy domination to derive an equation for global energy conservation

[TABLE]

where $L_{\rm heat}(t)$ is the total input heating rate. Using Eq. 6 to replace $E_{\rm int}=L\tau_{d}^{2}/t$ and rearranging gives

[TABLE]

Arnett’s rule follows, since the condition of an extremum $dL/dt=0$ implies $L=L_{\rm heat}$.

The ordinary differential equation Eq. 10 can be solved for $E_{\rm int}(t)$ and hence the emergent luminosity

[TABLE]

Both the one-zone and the separation of variables approaches result in the same expression for the light curve; the only difference is the value of $\xi$, which reflects different assumptions about the shape of the self-similar energy density profiles. The diffusion time for the one-zone models is a factor of $\pi/\sqrt{3}\approx 2$ larger. To avoid confusion, we hereafter define a characteristic diffusion timescale without any numerical factors

[TABLE]

and so $\tau_{d}=[3/4\pi\xi]^{1/2}t_{d}$. Physically, the characteristic diffusion timescale gives the time at which the expansion timescale $t_{\rm exp}=R_{\rm ej}/v_{\rm ej}$ equals the diffusion time $t_{\rm diff}=\kappa\rho R_{\rm ej}^{2}/c$. The peak time scales with the diffusion timescale $t_{\rm peak}\propto t_{d}$, but the numerical coefficient relating them depends on the distribution of heating, nature of the opacity, and other effects.

The self-similarity assumption will be shown below to limit the accuracy of the Arnett models. However, this is not a necessary assumption as Eq.8 is only the first eigenfunction of the separated spatial equation. The full solution of the energy density can be expressed as an infinite sum of higher order eigenfunctions whose normalization will be set by the spatial distribution of heating and boundary conditions. Pinto & Eastman (2000) show how such an approach can be used to relax the assumptions (5) and (6) and produce more accurate light curves. However, due to the more complicated nature of the solution, the full solution with higher-order eigenmodes is rarely used in practice.

2.2 Comparison to Numerical Simulations

To assess the accuracy of the Arnett solutions, we compare them to numerical monte-carlo radiation transport calculations run with Sedona (Kasen et al., 2006). We adopt similar assumptions as A82: homologous expansion, uniform density, and a constant opacity. Non-constant opacity will be considered in Section 6. In this section, the ejecta has a diffusion timescale $t_{d}=100$ days and the heating source is at the center and follows $L_{\rm heat}(t)=L_{0}e^{-t/t_{s}}$, where the timescale $t_{s}=10$ days.

Fig. (1) compares the numerical light curve to the Arnett analytic solution. The numerical models have an initial “dark period” until $t\sim 0.1t_{d}$, before which the photons have not had sufficient time to diffuse from the center of the ejecta (Piro & Nakar, 2013). In contrast, the analytic solutions predict a steeper rise beginning at $t=0$, a consequence of the assumption that radiation energy is immediately distributed throughout the ejecta. The A82 solution predicts a peak time a factor of 2 shorter than the numerical result, but gives roughly the correct peak luminosity. The peak time of the one-zone model is closer to the numerical simulation, but under-predicts the peak luminosity and is overall too broad. There is no choice of $\xi$ such that the analytic solution closely matches the numerical light curve.

This inaccuracy of the analytic models is more pronounced for more centrally concentrated heating sources. Fig. (2) shows numerical models where the heating source has been uniformly mixed out to dimensionless radius $x_{s}$, with $x_{s}=0$ corresponding to a central source. The Arnett analytic solution most closely resembles a well-mixed numerical model with $x_{s}\approx 0.8$. We can define a heating-weighted radius where the bulk of heating occurs as

[TABLE]

where $\dot{e}_{\rm heat}(x)$ is the energy density heating rate at $x$. For constant heating out to radius $x_{s}$, we have the relation $\langle x_{s}\rangle=x_{s}/\sqrt{3}$. In the Arnett solution, $\dot{e}_{\rm heat}(x)\propto e(x)$ and using Eq.(8) we find that $\langle x_{s}\rangle\approx 0.4$, which indeed corresponds to $x_{s}\approx 0.7$.

The limitation of the Arnett models stems from the assumption that the spatial distribution of the radiation field is self-similar. In reality, for central sources a radiation diffusion wave initially propagates outwards, only reaching the surface and establishing a self-similar profile after a timescale $\sim t_{d}$. In Fig. (3), we show the evolution of the energy density profile defined in Eq.(8) for a central exponential heating source. At early times, a diffusion wave propagates outwards, and the self-similar assumption fails. By neglecting this diffusion wave, the Arnett-like models overestimate the luminosity at early times. For a more uniformly mixed source, self similarity is established earlier, and so the Arnett models are more applicable.

Thus, it is a common misconception that the Arnett models assume a centrally-located heating source – while the energy density increases towards the center, the heating luminosity peaks close to the surface, producing the faster rise and earlier peak compared to the central-heating numerical solution.

This likely explains why the A82 solution more closely predicts Type Ia rather than core collapse SN light curves – Type Ia SN typically have a much larger degree of mixing (Blondin et al., 2013) while core collapse SNe have more centrally concentrated 56Ni. We explore the effects of the spatial distribution of the heating source in more detail in Section 5.

2.3 Arnett’s Rule

A specific prediction of the Arnett models is that the peak luminosity is equal to the heating rate at peak, i.e.

[TABLE]

This is commonly referred to as Arnett’s rule (or Arnett’s law) and is widely used to infer e.g. the nickel mass in radioactive SNe (Stritzinger et al., 2006; Valenti et al., 2007).

We see from Fig.(2) that Arnett’s rule does not hold in general and that its accuracy depends on the heating source concentration. For centrally concentrated sources ($x_{s}\lesssim 0.8$) Arnett’s law gives an under-estimate of the true peak luminosity, with the error being systematically worse for more centralized heating. For nearly fully mixed heating sources ($x_{s}\gtrsim 0.8$) Arnett’s law is an over-estimate.

The failure of Arnett’s rule again stems from the assumption of self-similarity, which implies a proportionality between the luminosity and the total ejecta internal energy, $L\propto Et$. Under this assumption, the light curve has to peak simultaneously with $E(t)$ – i.e., at the time when the rate of energy loss, $L(t)$ equals the rate of energy gain, $L_{\rm heat}(t)$. In reality, the time-dependent propagation of a diffusion wave means that the luminosity does not strictly track the internal energy. For central sources, $L(t)$ generally lags $E(t)$ and the light curve peaks at a time when $L>L_{\rm heat}$. For fully mixed cases, $L(t)$ leads $E(t)$ and the light curve peaks when $L<L_{\rm heat}$. Arnett’s law holds only for the case of a specific concentration ($x_{s}\approx 0.8$) for which the light curve coincidentally peaks at the same time as does the internal energy.

3 A New Relation Between Peak Time and Luminosity

Given the limits of Arnett’s rule, we look for a more robust relationship between the peak time and peak luminosity of a transient light curve. We proceed by considering the evolution of the global internal energy, $E$, and rewrite Eq. 9 as

[TABLE]

which integrates to

[TABLE]

Eq. (16) is similar to the analysis presented in (Katz et al., 2013), which considers times $t\gg t_{\rm peak}$ when $E(t)=0$. Here we instead consider times around peak $t\sim t_{\rm peak}$. Furthermore, we assume the initial energy content in the ejecta is zero and ignore the initial stellar radius in our assumption of homology. Thus, the analysis presented here does not necessarily apply to Type IIP/L SNe, whose light curves are dominated by the initial shock-deposited energy.

We rewrite Eq. (16) as

[TABLE]

where

[TABLE]

The first term in brackets can be shown to be positive (since $L(t)\leq L_{\rm peak}$) and monotonically increasing (see Appendix B). The second term $tE(t)$ is also positive and is a decreasing function when $L>L_{\rm heat}$, which is typically obtained for $t\gtrsim t_{\rm peak}$. We therefore anticipate there may be a time when the two functions cross and cancel to give $\epsilon(t)=0$.

We express this time as $t=\beta t_{\rm peak}$ and rearrange Eq.(17) to get

[TABLE]

which is our desired expression for $L_{\rm peak}$. In Appendix B we show that for common heating functions there is indeed a time when $\epsilon(t)=0$ for a value of $\beta\sim 1$ that can be calibrated from the numerical simulations and is essentially independent of the heating timescale or functional form.

Eq.(19) can be analytically evaluated for certain heating functions. For example, for an exponential heating function

[TABLE]

the peak time-luminosity relation can be evaluated to get

[TABLE]

This can be contrasted with Arnett’s rule, which predicts $L_{\rm peak}=L_{0}e^{-t_{\rm peak}/t_{s}}$. The two expressions make similar predictions when $t_{\rm peak}\ll t_{s}$ but increasingly diverges for $t_{\rm peak}\gg t_{s}$.

In Fig. (4) we compare our expression for $L_{\rm peak}$ to those of numerical light curve calculations for a central exponential heating source with $t_{s}=10$ days. The numerical models span a wide range of ejecta masses, velocities, and opacities, and hence result in a range of peak times. Eq. 21 with $\beta=4/3$ gives a near-perfect match to the numerical simulations, independent of the ejecta properties. In comparison, Arnett’s rule predicts systematically too low values of $L_{\rm peak}$ and becomes progressively worse for larger values of $t_{\rm peak}/t_{s}$.

The peak time-luminosity relation Eq. 19 with $\beta=4/3$ applies for most central heating functions, as long as the opacity is constant and the density uniform. In Fig. 5, we show the peak time-luminosity relations for both an exponential source and a power-law source appropriate for magnetar energy injection

[TABLE]

which can also be analytically evaluated (see Appendix A). Fig. 5 shows that the Eq. 19 with $\beta=4/3$ accurately reproduces the numerical calculations. In later sections, we show that the value of $\beta$ does change if the source heating is spatially mixed or the opacity is non-constant due to recombination, but that $\beta$ remains largely independent of the heating function, source timescale $t_{s}$, or ejecta diffusion timescale $t_{d}$.

4 Relation Between Peak Time and Diffusion Time

Analyses of observed light curves often attempt to constrain the ejecta mass and velocity by setting the observed time of peak, $t_{\rm peak}$, equal to the diffusion timescale $\tau_{d}$ (e.g., Drout et al., 2011; Prentice et al., 2018). Here we study that relation for constant opacity models, and show that $t_{\rm peak}$ depends not only on $t_{d}$, but also on the heating timescale $t_{s}$.

In Fig.(6), we show the dependence of $t_{\rm peak}$ on the ratio $t_{s}/t_{d}$, for a large number of numerical models with uniform density ejecta and two different central heating sources. The models have a range of masses, velocities, and constant opacities, although only the combination $t_{d}$ is relevant for the light curve behavior. For $t_{s}/t_{d}\ll 1$, the peak time asymptotes to $t_{\rm peak}\approx 0.4t_{d}$ independent of $t_{s}$. In this limit, the source can thus be approximated as an instantaneous “pulse” of energy deposited at $t_{s}$. The energy from such a pulse diffuses out and peaks at around $\sim 0.4t_{d}$. In comparison, the Arnett models predict $t_{\rm peak}\approx 0.2t_{d}$ (see Fig. (1)).

For $t_{s}/t_{d}\gtrsim 0.1$, the continuing source deposition begins to lengthen the peak time. The dependence is fairly weak – $t_{\rm peak}$ only increases by a factor of $\sim 2$ as $t_{s}$ changes over three orders of magnitude, implying that for the sources considered the light curve peak is mostly powered by heating deposited at early times.

An equation that captures the peak time of numerical models with constant opacity and central heating is

[TABLE]

In the limit that $t_{s}\ll t_{d}$, Eq.(23) goes to $t_{\rm peak}\approx 0.4t_{d}$, while for $t_{s}\gg t_{d}$ it grows logarithmically with $t_{s}$. The relation is relatively insensitive to the functional form of the heating source (e.g. exponential vs. power-law) as long as the function is smoothly and gradually declining.

5 Spatial Distribution of Heating

Another important effect in shaping the light curve is the spatial distribution of heating within the ejecta, e.g. different amounts of “mixing” of 56Ni in Type I SNe or assumed distribution of magnetar heating (see Dessart (2018)). Indeed, In Fig. (2), we showed how the spatial distribution impacts both the peak time and luminosity of the light curve, and in particular found that the Arnett solution and Arnett’s rule are most appropriate for less concentrated/more uniform heating.

To account for the spatial distribution of heating, we take the heating rate to be uniform out to a (scaled) radius $x_{s}$. Fig.(7) shows how the concentration affects the time of peak. The overall effect is to systematically drop the relation, i.e. for a given $t_{s}$ and $t_{d}$, concentration causes the light curve to peak earlier. This was shown for the case of $t_{s}$ = 10 days in Fig.(2). Interestingly, the peak time does not differ much unless the concentration radius is greater than $x_{s}>1/3$.

In all cases, there is similar behavior of a “flattening” in the relation for $t_{s}\ll t_{d}$. For the most mixed case $x_{s}=0.9$, the relation flattens as $t_{\rm peak}\approx 0.1t_{d}$. This lends further caution to using $t_{d}$ as a proxy for $t_{\rm peak}$; in addition to depending on $t_{s}$, there is also another dependence on $x_{s}$.

In Fig.(8), we show the peak time-luminosity relation for the different spatial distributions of heating, for an exponential source with $t_{s}=10$ days. Interestingly, for different concentrations, the relation Eq.(19) still holds. The only difference is in the value of $\beta$. For $x_{s}=1/3$, the numerical simulations lie on the $\beta=4/3$ relation, which was found to be appropriate for a central source. This is in agreement with the results shown in Figs. (2) and (7), where the $x_{s}=1/3$ does not differ significantly from simply assuming a central source.

More centrally concentrated heating acts to increase the value of $\beta$. For the most uniform heating, $x_{s}=0.9$, $\beta$ increases by about a factor of $2$ compared to a central source. For the central exponential source used in Fig.(7) and a constant opacity, we find that $\beta$ depends on $x_{s}$ approximately as

[TABLE]

Note that we assume local deposition of the heating source. In reality, for the case of e.g. 56Ni decay, there is the additional effect of gamma-ray deposition, which introduces a non-locality to the heating. In particular, gamma rays emitted closer to the center may deposit their energy farther out (or may escape entirely). Exploring this effect is outside the scope of this work (although see e.g. Dessart et al. (2016)).

6 Non-Constant Opacity and Recombination

While the previous results assumed a constant opacity, for certain compositions the opacity drops sharply when the ejecta cools and ions recombine. As the ejecta is typically hotter at the center, a cooling “recombination front” propagates from the surface inward (Grassberg et al., 1971; Grassberg & Nadyozhin, 1976; Popov, 1993). The photosphere is nearly coincident with the recombination front, with a temperature set by the ionization/recombination temperature $T_{\rm ion}$.

To account for recombination effects in our numerical calculations, we prescribe a temperature dependence that mimics the behavior of the opacity in hydrogen and helium-rich compositions, for which electron scattering dominates for $T>T_{\rm ion}$

[TABLE]

The opacity $\kappa=\kappa_{0}$ for temperatures $T>T_{\rm ion}$ but drops to $\kappa=\epsilon\kappa_{0}$ for $T<T_{\rm ion}$. The tanh function ensures a smooth transition over a temperature range $\Delta T_{\rm ion}$. We use $\epsilon=10^{-3}$ and $\Delta T_{\rm ion}=0.1T_{\rm ion}$, although our results are not sensitive to the exact values. We take the temperature $T$ to be equal to the radiation temperature $T_{\rm rad}=(E/a)^{1/4}$, where $E$ is the energy density and $a$ the radiation constant. This is appropriate for a radiation-dominated ejecta and wavelength-independent opacity.

Fig.(9) shows the effect of changing the recombination temperature $T_{\rm ion}$, while keeping the heating source and ejecta properties fixed. These runs use a central exponential source with timescale $t_{s}=10$ days and energy $E_{s}=10^{50}$ ergs, and uniform ejecta with mass $M_{\rm ej}=5M_{\odot}$, velocity $v_{\rm ej}=10^{9}$ cm s*-1*, and opacity $\kappa_{0}=0.1$ cm2 g*-1*. For low $T_{\rm ion}$, most of the ejecta remains ionized at and after peak and the light curve resembles the constant opacity case, with the exception of a late-time “bump” that occurs when recombination sets in and allows radiation to escape more easily. For higher $T_{\rm ion}$ recombination occurs earlier; for $T_{\rm ion}\gtrsim 6000$ K recombination results in a brighter and earlier light curve peak.

When recombination is included, the total heating energy scale, $E_{s}$, impacts the light curve morphology. This is in contrast to constant opacity models, where $E_{s}$ simply sets the normalization of the light curve but leaves the shape the same. Fig. (10) shows a set of numerical light curves where only $E_{s}$ is varied. At sufficiently large values of $E_{s}$, the heating source keeps the ejecta ionized until very late times, and the light curve shape resembles a constant opacity light curve. As $E_{s}$ decreases, the ejecta recombines earlier, resulting in a “bump” at late times. For sufficiently low $E_{s}$, the heating source is unable to keep the ejecta temperature above $T_{\rm ion}$, and so recombination impacts the light curve peak.

To determine whether recombination effects are important or not, we can compare the heating rate to the luminosity necessary to keep the ejecta ionized. The ionizing luminosity is set by the ejecta radius and the recombination temperature as

[TABLE]

From Section 3, the luminosity will roughly scale as

[TABLE]

We define a ratio of the luminosity to the critical ionizing luminosity

[TABLE]

We calibrate the proportionality based on numerical simulations of an exponential heating source to find

[TABLE]

where $E_{s,51}=E_{s}/10^{51}$ erg, $t_{s,10}=t_{s}/10$ days, $\kappa_{0.2}=\kappa/0.2$ cm2 g*-1*, $M_{10}=M_{\rm ej}/10M_{\odot}$, and $T_{4}=T_{\rm ion}/10^{4}$ K. For $\eta\lesssim 1$, the heating luminosity is too low to keep the ejecta sufficiently ionized and so recombination effects become important.

In Fig.(11), we show the results of a set of numerical simulations with $T_{\rm ion}=6000$K and various ejecta properties and $E_{s}$. Interestingly, the numerical simulations with recombination still fall on the relation Eq. (19). The only difference is that the value of $\beta$ changes. Specifically, recombination tends to decrease the value of $\beta$. Also shown in Fig. (11) are the respective values of $\eta$ for the numerical simulations. Points with $\eta>1$ are not affected by recombination at peak, and so fall on the relation with $\beta=4/3$, appropriate for a constant opacity. For $\eta<1$, recombination is important and the points fall on a smaller $\beta=0.94$ curve.

7 Discussion and Conclusions

We have shown how the light curve peak time and luminosity are related for luminous transients, and derived analytic relations that can be used to infer the physical properties of the heating mechanism. In particular, Eq.(19)

[TABLE]

captures the relationship between $t_{\rm peak}$ and $L_{\rm peak}$ where the light curve physics (i.e. recombination and concentration) is contained in the $\beta$ parameter. Furthermore, Eq.(24)

[TABLE]

gives the approximate dependence of $\beta$ on the spatial distribution of heating. Another useful result is given in Eq.(23) for central sources,

[TABLE]

which shows how $t_{\rm peak}$ depends not only on the diffusion time $t_{d}$, but also the heating timescale $t_{s}$. In addition, recombination will change the value of $\beta$ compared to a constant opacity. In Table 1, we give approximate values of $\beta$ for a variety of transients. In Appendix A, we evaluate the peak time-luminosity relation for specific heating sources.

For example, one of the common sources of heating in luminous transients is the radioactive decay chain of 56Ni (Stritzinger et al., 2006; Valenti et al., 2007). In particular, the decay chain of 56Ni $\rightarrow^{56}$Co $\rightarrow^{56}$Fe dominates the heating at timescales of interest. Fig.(13) shows the peak time-luminosity parameter space of 56Ni-powered transients, where the relation is given in Appendix A. For a given peak time and luminosity, one can thus infer an approximate value for the 56Ni mass for an appropriate choice in $\beta$.

We now consider the implications of our findings to analyzing observed SNe. In particular, Arnett’s rule is often used to infer the nickel mass of Type Ia SNe. In Fig. (14) we show the effects of varying the nickel mass $M_{\rm Ni}$ in a set of toy Ia models with a total mass $M_{\rm ej}=1.4M_{\odot}$ and constant density and opacity. The inner layers of ejecta in these models are composed of pure 56Ni, and so higher $M_{\rm Ni}$ corresponds to a larger nickel core and less centrally concentrated heating. As expected, Arnett’s rule works better for larger $M_{\rm Ni}$ and becomes progressively worse for the more centrally concentrated low $M_{\rm Ni}$ models. This suggests that analyses of SNe Ia using Arnett’s rule may be systematically biased, with the nickel mass of sub-luminous Ia’s being overestimated.

As another case study, we show in Fig. (12) the observed bolometric light curve of SN1987A, a Type II supernova whose primary peak is powered by radioactive 56Ni (Woosley, 1988; Suntzeff & Bouchet, 1990). The late-time light curve behavior gives a constraint on the 56Ni mass to be $M_{\rm Ni}\approx 0.07M_{\odot}$. Arnett’s rule predicts a $M_{\rm Ni}$ a factor of $2$ too large, whereas using the new relation Eq.(19) with $\beta=0.82$ (appropriate for hydrogen recombination $T_{\rm ion}\approx 6000$K and a largely centrally-located Ni-Co heating source, as inferred from numerical simulations (see Fig.(11)), gives $M_{\rm Ni}\approx 0.07M_{\odot}$, in agreement with the late-time determination.

As another example, we show in Fig.(15) the peak time-luminosity relation for the Type Ib/c SNe models presented in Dessart et al. (2016). As noted in their work, Arnett’s rule seems to overestimate the 56Ni mass of their models. Using Eq.(19), we find that the models lie on a $\beta=9/8$ relation. Given that the models do not have much mixing, we can assume centrally located heating and attribute any deviation from $\beta=4/3$ to recombination effects. Interestingly, a $\beta=9/8$ is in agreement with a recombination temperature of $T_{\rm ion}=4000$K, which is roughly that for a C- and O-rich composition. On the other hand, helium has a much higher recombination temperature and would imply a much smaller $\beta$; this indicates that the 56Ni in the Dessart et al. (2016) models are primarily diffusing out from the much denser carbon/oxygen inner ejecta rather than the outer helium ejecta. This is in agreement with the results in Piro & Morozova (2014), who similarly showed that light curve modeling is a better constraint on the C/O core rather than the helium.

The above examples demonstrate that, in principle, the peak time-luminosity relation may allow one to infer the composition of the ejecta solely from photometric observations. Suppose we know from observations of the radioactive tail of SN1987A that it is powered by $0.07M_{\odot}$ of 56Ni, and we assume the nickel to be largely centrally concentrated. From the peak time and luminosity we can solve for $\beta\approx 0.94$. Since each recombination temperature has its own unique value of $\beta$, we can then infer $T_{\rm ion}\sim 6000$K, suggesting a hydrogen-rich composition.

If the composition, and hence recombination temperature, of an observed supernova is constrained (e.g. from spectroscopic observations), then the derived value of $\beta$ may indicate the spatial distribution of the heating source. For example, one can assume to good approximation a constant opacity for Type Ia SNe. Assuming central heating, this would point to a $\beta=4/3$, yet SNe-Ia seem to obey Arnett’s rule fairly well, which corresponds to a larger value of $\beta$ if using Eq.(19). This is in agreement with the results presented in Fig.(8), since we expect Ia SNe to have a more uniform distribution of heating and hence fall on a larger $\beta$.

The main result, Eq.(19), is general enough to be applied to an arbitrary heating source, e.g. central-engine accretion, a magnetar, kilonovae, etc. Thus, for an observed transient peak time/luminosity that might be powered by other means than 56Ni, one must simply choose a different $L_{\rm in}(t)$ (which need not be analytic). Next, by choosing an appropriate $\beta$ (e.g. $4/3$ for constant opacity and a central source), one can constrain the heating source parameters. Note that there still remains a degeneracy in the heating source timescale $t_{s}$ and energy $E_{s}$. One is not able to break this degeneracy from the peak time and luminosity alone. Such constraints require additional information/observations or by putting physical limits on allowed values.

Several physical effects were neglected in our analysis here so as to isolate the basic behavior of supernova light curves. The models presented assume spherical symmetry and adopt a grey opacity. Asymmetries in the ejecta/heating as well as non-grey effects likely play a role in the overall shape of the light curve, and on the inferred $\beta$ in the new relation. Dessart et al. (2018) show that clumping affects the recombination rate, which would impact the inferred $\beta$. We also used a simple parameterization for the spatial distribution of heating, which was taken to be uniform out to some radius; more complicated distributions (e.g. from accounting for gamma-ray deposition throughout the ejecta) warrant further investigation. We further made the assumption of homologous expansion, which may be violated for supernovae interacting with a circumstellar material. Because our relations apply to the bolometric peak, errors can also arise from observational effects, such as uncertainties in the distance, reddening, or bolometric correction.

Conclusions drawn from the relation presented here are of course conditional on the specific form of the heating source assumed (e.g. radioactive decay vs. magnetar spindown). It is thus important to perform consistency checks on the nature of the heating source using additional information aside from the properties at peak, such as examining the slope of the late-time light curve tail.

While our comparisons here have demonstrated the limitations of the Arnett-like models, there exist other analytic models of transient light curves that attempt to account for the time-dependent diffusion effects that are important for setting the luminosity before and around peak (e.g. Piro & Nakar (2013, 2014); Waxman et al. (2018)). In future work, we will investigate a broader range of analytic models, and look for improved methods for calculating analytic light curves. Additionally, there exist other analytic techniques in estimating properties of the light curve, in particular the integral relation of Katz et al. (2013). This requires knowing the full shape of the light curve out to times well after peak, rather than the properties at peak. The Katz integral approach and the new relation presented here are thus complementary in inferring the physical properties of the light curve.

There is still much work to be done to understand how/why the peak time-luminosity relation as well as it does, and in particular to better calibrate its value for specific heating mechanisms and ejecta properties. Nonetheless, the framework presented here will be useful for more detailed modeling, as well as providing a fast way to characterize the large number of transients to be discovered in current and upcoming surveys.

D.K.K was supported in part by the National Science Foundation Graduate Research Fellowship Program. D.K.K and D.N.K are supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract number DE-AC02-05CH11231 and DE-SC0017616; SciDAC award DE-SC0018297; the Gordon and Betty Moore Foundation through Grant GBMF5076; and the Exascale Computing project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy and the National Nuclear Security Administration. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC0205CH11231.

Appendix A Expressions for the Peak Time-Luminosity Relation

From Section 3, we found an expression for the peak time-luminosity relation as

[TABLE]

where $L_{\rm peak}$ is the observed peak luminosity at the peak time $t_{\rm peak}$, $L_{\rm heat}(t)$ is the time-dependent heating rate, and $\beta$ is a constant that depends on opacity/concentration effects. Eq.(A.1) can be evaluated analytically for several functional forms of $L_{\rm heat}(t)$.

Suppose the heating source can be written most generally as

[TABLE]

where $\mathcal{H}(t,t_{s})$ is the time-dependent component of $L_{\rm heat}(t)$, $t_{s}$ is the heating source timescale, and $E_{s}$ the characteristic heating energy. Equivalently, the heating source can be expressed in terms of a characteristic luminosity by setting $L_{0}=E_{s}/t_{s}$.

Let $\tau=t_{\rm peak}/t_{s}$ be the ratio between the peak time and source timescale. Then the peak luminosity can be evaluated to get

[TABLE]

where $f(\tau,\beta)$ depends on the exact functional form of $\mathcal{H}(t,t_{s})$. For an exponential source,

[TABLE]

the integral can be evaluated to get

[TABLE]

In Table 1, we give the analytic expressions of $f(\tau,\beta)$ for a variety of heating functions $\mathcal{H}(t,t_{s})$. The choice of $\beta$ again depends on opacity/recombination effects, as well as the spatial distribution of heating. In Table 2, we give approximate values of $\beta$ based on numerical results.

A.1 Radioactive Two-Decay Chain

Numerous transients are powered by the radioactive decay of synthesized elements, e.g. 56Ni in Type I and IIb/pec supernovae.

Consider a decay chain consisting of $0\rightarrow 1\rightarrow 2$ with decay timescales $t_{0}$ and $t_{1}$, respectively (ignoring the contribution to heating of species 2). The total number of the species at time $t$ is expressed as

[TABLE]

where $N=N_{0}(0)$, and it is assumed that $N_{1}(0)=0$. Let $Q_{0}$, $m_{0}$, $Q_{1}$, and $m_{1}$ be the decay energies and species mass. Define heating rates per unit mass as

[TABLE]

Then the heating luminosity can be expressed as

[TABLE]

where $M=Nm_{0}$ is the total initial mass of species [math].

With this expression we can derive the peak time-luminosity relation as

[TABLE]

A.1.1 Radioactive Nickel Decay

Supernovae of Type I and IIb/pec are powered primarily by the radioactive decay of 56Ni followed by 56Co (Stritzinger et al., 2006; Valenti et al., 2007). The heating function can be written in terms of the nickel mass $M_{\rm Ni}$ as

[TABLE]

where $\varepsilon_{\rm Ni}=3.9\cdot 10^{10}$ erg g*-1* s*-1* and $\varepsilon_{\rm Co}=6.8\cdot 10^{9}$ erg g*-1* s*-1* are the specific heating rates of Ni- and Co-decay, and $t_{\rm Ni}=8.8$ days and $t_{\rm Co}=111.3$ days are the decay timescales.

We can evaluate the peak time-luminosity relation in terms of $\beta$ as

[TABLE]

Let $t_{p}=t_{\rm peak}/$day be the peak time in days. Using the numerical values of $\varepsilon_{\rm Ni}$, $\varepsilon_{\rm Co}$, $t_{\rm Ni}$, and $t_{\rm Co}$ we get a more compact numerical expression (accurate to within $\sim 1\%$) as

[TABLE]

A.1.2 Magnetar-powered Supernovae

The spindown luminosity of a magnetar is generally described by (Bodenheimer & Ostriker, 1974; Gaffet, 1977; Woosley, 2010; Kasen & Bildsten, 2010)

[TABLE]

where $l=2$ for magnetic dipole spin-down,

[TABLE]

is the magnetar energy with $P_{10}=P/10$ms is the spindown period, and

[TABLE]

is the spindown timescale with $B_{14}=B/10^{14}$G the magnetic field strength. For $l=2$, we get

[TABLE]

were $\tau=t_{\rm peak}/t_{\rm mag}$.

A.2 Accretion-Powered Transients

Another interesting heating source is that of an accreting compact object (Dexter & Kasen, 2013). Let $\dot{M}\sim M_{\rm acc}/t_{\rm acc}$ be the accretion rate of mass $M_{\rm acc}$ and timescale $t_{\rm acc}$. We here consider two functional forms of accretion luminosity. The first is of constant heating that “shuts off” after a time $t_{\rm acc}$,

[TABLE]

where $\Theta(t-t_{\rm acc})$ is the Heaviside step function, and $\epsilon$ is the radiative efficiency. Substituting this into Eq.(19) we get

[TABLE]

Another functional form of interest is that of an $n=-5/3$ power law, appropriate for fallback accretion

[TABLE]

Again evaluating this source in Eq.(19) we get

[TABLE]

A.3 Kilonovae

The kilonova heating rate from the radioactive decay of r-process elements can be parameterized as (Li & Paczyński, 1998; Metzger et al., 2010; Roberts et al., 2011; Lippuner & Roberts, 2015).

[TABLE]

where $\varepsilon_{0}\approx 10^{11}$ erg g*-1* s*-1* is the specific heating rate, $M_{\rm ej}$ is the ejecta mass, $t_{0}\approx 1$ day, and $\eta\approx 1.3$. The function $f(t)$ gives the thermalization efficiency with which the radioactive decay energy is able to deposit as heat in the ejecta.

(Kasen & Barnes, 2018) suggest an analytic approximation of the thermalization efficiency of electrons

[TABLE]

where

[TABLE]

is the electron thermalization timescale, $M_{0.01}=M_{\rm ej}/0.01M_{\odot}$, and $v_{0.2}=v_{\rm ej}/0.2c$ the ejecta velocity.

Evaluating the integral for $\eta=1.3$ and the above approximation for the thermalization efficiency, we get

[TABLE]

where ${}_{2}F_{1}(a,b,c,x)$ is the hypergeometric function which we approximate by a functional fit of

[TABLE]

Thus, the peak time-luminosity relation for kilonova is approximately

[TABLE]

Assuming the r-process heating is uniformly mixed throughout the ejecta, we choose an approximate $\beta\approx 2$ based on numerical simulations.

Appendix B Derivation of the Peak Time-Luminosity Relation

In Section 3, we showed that a simple relation holds between the peak time and luminosity of a light curve

[TABLE]

assuming there exists some time $t=\beta t_{\rm peak}$ such that $\epsilon(t)=0$, where we defined

[TABLE]

Our goal, then, is to motivate that such a time when $\epsilon(t)=0$ exists. We define two quantities

[TABLE]

and

[TABLE]

Thus, the time when $\epsilon(t)=0$ also implies

[TABLE]

Initially, $\mathcal{E}(t)$ will rise as heat is deposited and trapped in the optically thick ejecta. Eventually, the expanding ejecta becomes optically thin, allowing radiation to freely escape. As a result, any heating goes directly into the light curve, $L(t)=L_{\rm heat}(t)$. This implies that $\mathcal{E}(t)=0$ at late times $t\gg t_{\rm peak}$.

Since $\mathcal{E}(t)$ is continuous, if $\mathcal{F}(t)$ is a monotonically increasing function of time and $\mathcal{F}(t)<\mathcal{E}(t)$ for some $t$, then it follows that $\mathcal{F}(t)$ and $\mathcal{E}(t)$ must intersect.

Taking the derivative in time of $\mathcal{F}(t)$ we get

[TABLE]

Since $L(t)\leq L_{\rm peak}$ by definition, it follows that

[TABLE]

and so $\mathcal{F}(t)$ is indeed a monotonically increasing function of time.

Next, we need to show that there exists a time such that $\mathcal{F}(t)<\mathcal{E}(t)$. At $t=0$ we have $\mathcal{F}(t)=\mathcal{E}(t)=0$. For small $t$, we can expand the derivative of $\mathcal{E}(t)$ to get

[TABLE]

where we make use of the fact that $L(0)=0$. Similarly, for $\mathcal{F}(t)$, we get

[TABLE]

If the condition $L_{\rm heat}(0)>L_{\rm peak}$ is satisfied, then

[TABLE]

for small $t$. Since $\mathcal{F}(0)=\mathcal{E}(0)=0$, we have that, at early times, $\mathcal{E}(t)>\mathcal{F}(t)$. Combined with the monotonicity of $\mathcal{F}(t)$ and the fact that $\mathcal{E}(t)\rightarrow 0$ at late times, it follows that $\mathcal{E}(t)$ and $\mathcal{F}(t)$ must intersect. In other words, there exists a time such that $\epsilon(t)=0$ and the peak time-luminosity relation holds.

There is a-priori mathematical justification for $L_{\rm heat}(0)>L_{\rm peak}$, however in physical cases of interest (e.g. radioactive decay, magnetar spindown, etc.), this seems to be a valid assumption, whereby $L_{\rm heat}(t)$ is monotonically decreasing in time. Furthermore, both diffusion and adiabatic degradation act to spread out and decrease the heating luminosity in time. This is confirmed in our numerical simulations for a wide variety of heating functions, wherein all the light curves seem to indicate $L_{\rm heat}(0)>L_{\rm peak}$ (see e.g. Figs. (1) and 2).

Another mathematical possibility is that $\mathcal{F}(t)$ and $\mathcal{E}(t)$ intersect more than once. Assuming $L_{\rm heat}(0)>L_{\rm peak}$, since $\mathcal{F}(t)$ is monotonically increasing, such behavior requires $\mathcal{E}(t)$ to “oscillate”, i.e. there exists more than one time that $d\mathcal{E}/dt=0$. It is unclear whether such behavior is physical.

In summary, it is difficult to prove definitively that the peak time-luminosity relation holds for an arbitrary heating function. For a variety of monotonically decreasing heating functions, we have confirmed numerically that $L_{\rm heat}(0)>L_{\rm peak}$. Furthermore, this intersection appears to occur only once. In Fig.(16), we show the behavior of $\mathcal{E}(t)$ and $\mathcal{F}(t)$ for a subset of numerical simulations, assuming a central source and constant opacity. We pick different functional forms (exponential and magnetar-like) as well as different source timescales. In all cases, the time when $\epsilon(t)=0$ is nearly identical, with $\beta\approx 1.3$.