Late time kilonova light curves and implications to GW 170817

TL;DR

This paper models late-time kilonova light curves from neutron star mergers, showing a universal decay rate of approximately t^-2.8, consistent with observations of GW 170817 if ionization losses dominate after about a week.

Contribution

It introduces a universal model for late-time kilonova emission based on ionization energy loss, applicable across different compositions and radioactive decay rates.

Findings

Late-time luminosity decays approximately as t^-2.8.

IR and optical data of GW 170817 align with the model if ionization losses dominate after 7 days.

The decay rate is nearly independent of composition and decay energy release.

Abstract

We discuss the late time (tens of days) emission from the radioactive ejecta of mergers involving neutron stars, when the ionization energy loss time of beta-decay electrons and positrons exceeds the expansion time. We show that if the e$^\pm$ are confined to the plasma (by magnetic fields), then the time dependence of the plasma heating rate, $\dot{\varepsilon}_d$, and hence of the bolometric luminosity $L=\dot{\varepsilon}_d$, are given by $d\log L/d\log t\simeq-2.8$, nearly independent of the composition and of the instantaneous radioactive energy release rate, $\dot{\varepsilon}$. This universality of the late time behavior is due to the weak dependence of the ionization loss rate on composition and on e$^\pm$ energy. The late time IR and optical measurements of GW 170817 are consistent with this expected behavior provided that the ionization loss time exceeds the expansion time at $t>t_\varepsilon\approx 7$~d, as predicted based on the early (few day) electromagnetic emission.

Tables

Table 1. Table 1 Late time observations of GW 170817

Time	Frequency	$L_{\nu }$	$\Delta L_{\nu }$	Instrument	Reference
(day)	(Hz)	(erg s-1 Hz-1)	(erg s-1 Hz-1)
9.20	$7.25\times {10}^{16}$	$1.65\times {10}^{21}$	$8.42\times {10}^{20}$	Chandra	Nynka et al. 2018
9.21	$7.25\times {10}^{16}$	$1.12\times {10}^{21}$	$4.81\times {10}^{20}$	Chandra	Margutti et al. 2018
10.37	$6.20\times {10}^9$	$1.49\times {10}^{25}$	$4.98\times {10}^{24}$	VLA	Hallinan et al. 2017
15.39	$7.25\times {10}^{16}$	$1.53\times {10}^{21}$	$4.66\times {10}^{20}$	Chandra	Margutti et al. 2018
15.60	$7.25\times {10}^{16}$	$7.41\times {10}^{21}$	$2.67\times {10}^{21}$	Chandra	Nynka et al. 2018

Table 2. Table 2 Fitted parameters for the synchrotron model

mean time	time range	${\log }_{10}\left[A/\left(\mathrm{erg}{s}^{-1}{\mathrm{Hz}}^{-1}\right)\right]$	$a$	$b$	${\chi }^2$/dof	RMS	$N_X$	$N_R$
(day)	(day)
17.1	9– 30	$29.555\pm 0.782$	$1.09\pm 0.51$	$-0.56\pm 0.02$	10.9/13	0.204	5	11
94.8	20–150	$29.823\pm 0.203$	$0.75\pm 0.09$	$-0.55\pm 0.01$	16.7/31	0.106	2	32
132.7	100–180	$30.674\pm 0.576$	$0.37\pm 0.26$	$-0.55\pm 0.01$	21.1/29	0.108	3	29
171.3	150–200	$31.031\pm 1.417$	$0.30\pm 0.62$	$-0.57\pm 0.01$	4.9/14	0.074	1	16
225.8	180–275	$36.210\pm 0.579$	$-1.98\pm 0.25$	$-0.58\pm 0.01$	3.2/18	0.060	1	20
256.7	210–400	$34.491\pm 0.483$	$-1.25\pm 0.21$	$-0.58\pm 0.01$	4.8/17	0.058	2	18

Table 3. Table 3 Ratio between observed flux and interpolated synchrotron flux

time	Band	Obs/pred.	Error
(day)
43	4.5$\mu$m	39	0.05
74	4.5$\mu$m	4.1	0.2
172.12	F810W	$\sim 0.84$	0.2
172.19	F606W	$\sim 0.5$	0.2
296.74	F606W	0.7	0.2
326.61	F606W	0.6	0.2
361.7	F606W	0.4	0.4

Equations

$$\frac{dE}{dX}$$

$$\frac{d{E}^B}{dX}$$

$$\kappa_{e}\equiv E(dE/dX)^{-1}\big{|}_{E=1{\rm MeV}}\approx 1{\rm cm^{2}/g},$$

$$\frac{1}{\rho c}\frac{dE}{dt}=\frac{v_e}{c}\frac{dE}{dX}=1g(E){\kappa }_e\mathrm{MeV},$$

$$t_{\epsilon }$$

$$\frac{dE}{dX}$$

$$\frac{E_{\alpha }^{-1}\frac{d{E}_{\alpha }}{dt}}{E_e^{-1}\frac{d{E}_e}{dt}}\approx 2(E_{\alpha }/10\mathrm{MeV}{)}^{-3/2}(E_e/1\mathrm{MeV}).$$

$${\dot{\epsilon }}_d/1\mathrm{MeV}={\kappa }_e\rho c\int dE\frac{dn(E,t)}{dE}g(E).$$

$$\frac{dn(E,t)}{dE}=\dot{n}\left[t_0(E,t)\right]\frac{\partial t_0}{\partial E}.$$

$${\dot{\epsilon }}_d/1\mathrm{MeV}={\kappa }_e\rho c{\int }^{t}d{t}_0\dot{n}(t_0)g[E(t_0,t)].$$

$$\frac{d{E}_{\mathrm{MeV}}}{dt}=-x\frac{E_{\mathrm{MeV}}}{t}-{\kappa }_e\rho cg(E).$$

$$g(E)\propto E^{-{\omega }_I},$$

$$\frac{dϵ}{d\tilde{t}}=-x\frac{ϵ}{\tilde{t}}-\frac{{ϵ}^{-{\omega }_I}}{{\tilde{t}}^3}.$$

$$ϵ=({\tilde{t}}_0/\tilde{t}{)}^x.$$

$${ϵ}_m\approx {\tilde{t}}^{-x}.$$

$${\dot{\epsilon }}_d\propto t^{-3}{\int }_1^{\tilde{t}}d{\tilde{t}}_0\dot{n}(t_0)({\tilde{t}}_0/\tilde{t}{)}^{-x{\omega }_I}.$$

$$\frac{d\log {\dot{\epsilon }}_d}{d\log t}=-3+x{\omega }_I\mathrm{for}x\le 2/(1+{\omega }_I).$$

$${ϵ}_c\approx {\tilde{t}}^{-2/(1+{\omega }_I)}.$$

$${\tilde{t}}_c={\tilde{t}}_0({ϵ}_c,\tilde{t})={\left({ϵ}_c{\tilde{t}}^x\right)}^{1/(x-{\omega }_E)}\approx {\tilde{t}}^{\frac{(1+{\omega }_I)x-2}{(x-{\omega }_E)(1+{\omega }_I)}}.$$

$${\dot{\epsilon }}_d\propto t^{-3}{\int }_{{\tilde{t}}_c}^{\tilde{t}}d{\tilde{t}}_0\dot{n}(t_0)({\tilde{t}}_0/\tilde{t}{)}^{-x{\omega }_I}.$$

$$\frac{d\log {\dot{\epsilon }}_d}{d\log t}=-3+\frac{2{\omega }_I}{1+{\omega }_I}-\frac{[(1+{\omega }_I)x-2](\beta -{\omega }_E-1)}{(1+{\omega }_I)(x-{\omega }_E)}.$$

$$\frac{d\log {\dot{\epsilon }}_d}{d\log t}=-3+\frac{2{\omega }_I}{1+{\omega }_I}\frac{(3-\beta )}{(2-{\omega }_E)}.$$

$$\frac{d\log {\dot{\epsilon }}_d}{d\log t}=-3+\frac{2{\omega }_I}{1+{\omega }_I}-f({\omega }_I)\frac{t_{\epsilon }}{\tau }{\left(\frac{t}{t_{\epsilon }}\right)}^{\frac{{\omega }_I}{1+{\omega }_I}}.$$

$$\frac{d\log {\dot{\epsilon }}_d}{d\log t}\simeq -3+\frac{2{\omega }_I}{1+{\omega }_I}-f({\omega }_I)\frac{t_{\epsilon }}{\tau }.$$

$$f({\omega }_I)=\frac{1}{2}{\left(\frac{2{\omega }_I}{1+{\omega }_I}\right)}^{1-\frac{1}{2(1+{\omega }_I)}}.$$

$$t_{\epsilon }$$

$${\dot{\epsilon }}_d\propto t^{-3+x{\omega }_I},$$

$$L\propto\left\{\begin{array}[]{ll}(t/t_{\varepsilon})^{-\beta-2},&\hbox{free electron escape;}\\ (t/t_{\varepsilon})^{-3+x\omega_{I}\simeq-2.8},&\hbox{electron confinement.}\end{array}\right.$$

$$L_{\nu }({\nu }_1)=(2-\Gamma )\frac{{\nu }_1^{1-\Gamma }}{{\nu }_2^{2-\Gamma }-{\nu }_1^{2-\Gamma }}L.$$

$$L_{\nu }=A{t}^a{\nu }^b.$$

Full text

††slugcomment: Draft of

1. Introduction

The merger of a neutron star with its binary neutron star or black hole companion has been suggested (Lattimer & Schramm, 1974) to produce high density neutron rich ejecta, in which heavy elements beyond Iron are produced by the r-process. Heating of the expanding ejecta by radioactive decay of unstable isotopes was in turn predicted (Li & Paczyński, 1998) to produce a strong optical-UV emission, commonly referred to as a ”kilonova” (see e.g. Fernández & Metzger, 2016, for a recent review). Remarkably, the observed UV-IR emission following the neutron star merger event GW 170817 (Coulter et al., 2017; Soares-Santos et al., 2017; Valenti et al., 2017; Arcavi et al., 2017; Lipunov et al., 2018; Tanvir et al., 2017), which was detected through its gravitational wave emission (Abbott et al., 2017), is broadly consistent with these predictions (e.g. Kasen et al., 2017; Drout et al., 2017; Rosswog et al., 2018; Cowperthwaite et al., 2017; Tanaka et al., 2017; Perego et al., 2017; Kilpatrick et al., 2017; Waxman et al., 2018). However, our understanding of the UV-IR signal and of the constraints it provides on the structure and composition of the ejecta is incomplete.

The early UV/blue emission suggests the existence of a fairly massive, few $\times 10^{-2}M_{\odot}$, and fast, $v\sim 0.3c$, component of the ejecta with low opacity, $\kappa<1$ cm2 g*-1*, corresponding to a large initial electron fraction $Y_{e}$ leading to a small fraction in the ejecta of $r$-process elements beyond $A\simeq 140$, which are expected to have much larger opacity (e.g. Kasen et al., 2013). The presence of this massive component is in tension with the results of detailed numeric merger calculations, that generally predict high opacity for the massive, fast ejecta components. This has lead to a discussion of possible mechanisms for the generation of the observed low opacity ejecta (e.g. Metzger et al., 2018; Fujibayashi et al., 2018). Alternative models have also been proposed, in which the blue emission is produced by boosted relativistic material and/or shock cooling of an expanding mildly relativistic shell (e.g. Kasliwal et al., 2017; Piro & Kollmeier, 2018; Gottlieb et al., 2018).

On a few days time scale, the peak emission shifted from the UV/blue to the red/IR. This observed blue to red evolution may be explained by the existence of several ejecta components, characterized by largely differing compositions, with higher opacity components dominating at later times (e.g. Kasen et al., 2017; Drout et al., 2017; Rosswog et al., 2018; Cowperthwaite et al., 2017; Tanaka et al., 2017; Perego et al., 2017; Kilpatrick et al., 2017). Using a simple analytic model for the UV-IR emission we have shown (Waxman et al., 2018, hereafter Paper I) that an alternative explanation is possible, in which the entire ejecta is composed of low opacity, $\kappa<1$ cm2 g*-1*, material (see also Smartt et al., 2017; Rosswog et al., 2018). The inferred significant opacity at the 1-2 $\mu m$ band at a few days, $\kappa\approx 0.1$ cm2 g*-1*, provides an important constraint on the composition, the implication of which is uncertain due to the uncertainty in the opacity (see discussion in Paper I).

Our model deviates from those of earlier work in two aspects, which are key to enabling an explanation of all data (see Waxman et al., 2018; Arcavi, 2018) with a nearly uniform composition ejecta: allowing a wider velocity distribution than previously assumed, and including the suppression of radioactive heating due to adiabatic losses of beta-decay electrons and positrons. Hereafter, we use ”beta-decay electrons” to refer to both beta-decay electrons and positrons.

On a time scale of minutes to tens of days the main radioactive energy source is expected to be beta-decays (alpha-decays may provide a significant contribution on a 100 d time scale, e.g. Barnes et al., 2016). While neutrinos (and, at this stage, also $\gamma$-rays) escape the ejecta, the electrons produced in beta-decays heat it through ionization and plasma losses. As the ejecta expand, the electron ionization loss time increases in proportion to $\rho^{-1}\propto t^{3}$, exceeding the expansion time $t$ at $t>t_{\varepsilon}$ (the expansion and loss time are equal at $t=t_{\epsilon}$). We have pointed out in Paper I that while the beta-decay electrons are likely confined (by magnetic fields111Magnetic confinement appears likely since it would be facilitated by a relatively weak magnetic field. Assuming $B\propto r^{-2}$ in the expanding ejecta, as may be appropriate for a tangled field, $B\gtrsim 0.1\mu$G is expected at $t\sim 10$ d, $r\sim 10^{16}$ cm (the equipartition field at 10 d is $\sim 1$ mG), implying an electron Larmor radius of $\sim 10^{10}$ cm$\ll r\sim 10^{16}$.) to the plasma, the fraction of their energy that is converted to heat at $t>t_{\varepsilon}$ is only $(t/t_{\varepsilon})^{-2}$, since they lose energy to adiabatic expansion (applying pressure against the expanding plasma) faster than to ionization. This leads to a steepening of the decline of the bolometric luminosity at $t\sim t_{\varepsilon}$. Note that adiabatic losses become important for supernovae much later than in the kilonova case due to the larger mass and slower velocity of the ejecta of supernovae, $t_{\varepsilon}\propto M^{1/2}v^{-3/2}$, see Eq. (5).

Analyzing the UV-IR data of GW 170817, we have shown in Paper I that the analytic model for the emission by a uniform composition ejecta is over-constrained. All its parameters are determined by the observations at $t<6$ d, predicting a transition to inefficient energy deposition and transparency at roughly the same time, $t=t_{\varepsilon}\approx 7$ d, and $L\propto t^{-3}$ at $t>t_{\varepsilon}$. This is consistent with the observed steepening of the light curve at $\approx 7$ d, and with the strong deviation from thermal spectra at t\mathrel{\raise 1.29167pt\hbox{>}\mkern-14.0mu\lower 2.58334pt\hbox{\sim}}5 d.

We return here to the issue of the late time kilonova light curve for two reasons. First, new IR and optical observations extend the range of time over which data are available to $>100$ d (Kasliwal et al., 2018; Lamb et al., 2019), compared to $\sim 15$ d at the time Paper I was written. This enables further tests of the models. Second, in Paper I we have neglected plasma heating by low energy electrons, that may accumulate in the ejecta if they are confined by magnetic fields. Kasen & Barnes (2018) have pointed out that such accumulation leads to a shallower decline of the bolometric luminosity at $t>t_{\varepsilon}$, in particular if the ionization energy loss rate increases with decreasing electron energy. We improve our analytic model to include these effects. Our results differ from those of Kasen & Barnes (2018) (who find $L=\dot{\varepsilon}_{d}\propto\dot{\varepsilon}t^{-1}$), for reasons explained in § 2 (see end of § 2.3).

This paper is organized as follows. The processes affecting the late time, $t\gtrsim t_{\varepsilon}$, radioactive heating of the ejecta are discussed in § 2.1. A simple analytic derivation of the asymptotic, $t\gg t_{\varepsilon}$, light curve behavior for the case where electrons are confined to the plasma is given in § 2.2. The range of validity of the asymptotic behavior is discussed based on a comparison to complete solutions of the equations (given in the appendix). The main results of the analysis of § 2.1 and § 2.2 are given in § 2.3. Simple analytic formulae describing the late-time kilonova light curves, updating Eq. (12) of Paper I to include the effects of the accumulation of low-energy electrons and of the energy dependence of the ionization losses, are given. The new, late time, observations of GW 170817 are analyzed in § 3. Our conclusions are discussed in § 4.

2. Late time kilonova light curves

2.1. Late radioactive heating

We consider a fluid element within the expanding ejecta, with velocity $v$ and density $\rho\propto t^{-3}$. As noted in the introduction, on a time scale of minutes to tens of days the main radioactive energy source is expected to be beta-decays, while alpha-decays may provide a significant contribution on a 100 d time scale. We consider first heating by beta-decay, and return to alpha-decay at the end of this sub-section. Since neutrinos and, at the late times in which we are interested here, also gamma-rays escape the plasma, we define $\dot{\varepsilon}$ as the rate at which radioactive energy is released in the form of electron kinetic energy.

The electrons and positrons lose energy to plasma heating mainly by ionization (e.g. Longair, 1992),

[TABLE]

where $dE/dX$ is the energy loss per unit column density (grammage) traversed by the electron/positron, $dX=\rho dx$, $Z$ and $A$ are the atomic and mass numbers of the plasma nuclei, and $\bar{I}$ is the effective average ionization energy of the atoms, which is empirically determined and approximately given by $\bar{I}=10Z$ eV for heavy nuclei.

Electrons lose energy also by scattering free plasma electrons (”plasma losses”). The plasma loss rate is given by an equation similar to Eq. (1), the main differences being replacing $Z$ with the number of free electrons per atom and $\bar{I}$ with $\hbar\omega_{p}$, where $\omega_{p}$ is the plasma frequency. At the times under discussion, the ejecta temperature is a fraction of 1 eV and its number density is $n\sim 10^{6}(t/10{\rm d})^{-3}{\rm cm^{-3}}$ of $A\sim 100$ atoms. Under these conditions, the losses are dominated by ionization for electron kinetic energy larger than $\bar{I}\sim 1$ keV (for $\hbar\omega_{p}\sim 10^{-7}$ eV and ionization of a few, plasma losses at $E\sim 1$ MeV amount to $\lesssim 20\%$ of the ionization losses for high $Z$ nuclei). Finally, at highly relativistic energy Bremsstrahlung losses dominate over ionization (e.g. Longair, 1992),

[TABLE]

Note that since gamma-rays escape the plasma, only a fraction of the Bremsstrahlung energy is deposited in the plasma.

Figure 1 shows the energy loss rate of electrons as a function of their kinetic energy, $E=(\gamma_{e}-1)m_{e}c^{2}$. The figure demonstrates that the ionization loss rate is nearly independent of $E$ and of the plasma composition. Following Paper I, we define an effective ionization ”opacity”,

[TABLE]

such that $\kappa_{e}^{-1}$ is the grammage over which an electron of $E=1$ MeV loses its energy by ionization and plasma losses. Using $\kappa_{e}$, we may write

[TABLE]

where $g(E)\approx 1$ is a weak function of energy. For the relevant energy range (see below) of $\sim 0.01$ to $\sim 1$ MeV, we may approximate $g(E)=(E/1{\rm MeV})^{-0.15}$, reflecting a factor 2 change in the loss rate as the energy decreases from 1 to 0.01 MeV.

The time at which the ionization loss time is equal to the expansion time is given by $t_{\varepsilon}dE/dt=E$, or $\kappa_{e}\rho(t_{\varepsilon})ct_{\varepsilon}g(E)=(E/1{\rm MeV})$. Since $tdE/dt\propto\rho t\propto t^{-2}$, the fraction of the energy of electrons produced at time $t\gg t_{\varepsilon}$, that is deposited in and heat the plasma, is $\approx(t/t_{\varepsilon})^{-2}$. As pointed out in paper I, this is valid both for electrons that escape the plasma freely, in which case they spend a time $\sim r/c\propto t$ in the plasma, and for electrons that are confined to the plasma by magnetic fields. In the latter case, the electrons lose energy adiabatically, i.e contribute to the kinetic energy of the plasma by applying pressure against the expanding flow. Assuming that the confined energetic electrons behave as an ideal gas, we have $p\propto\rho^{4/3}$ and $E\propto\rho^{1/3}\propto t^{-1}$ for highly relativistic electrons, and $p\propto\rho^{5/3}$, $E\propto\rho^{2/3}\propto t^{-2}$ in the highly non-relativistic limit. Both relativistic and non-relativistic electrons lose most of their energy to adiabatic expansion over time $t$, hence the fraction of their energy that is converted to thermal energy is $\approx(t/t_{\varepsilon})^{-2}$.

For an ejecta mass $M$ and characteristic velocity $v_{M}$, neglecting the very weak dependence of $g$ on $E$ (in particular at $E\sim 1$ MeV), we have

[TABLE]

Here, $E=1E_{\rm MeV}$ MeV, $v_{M}=10^{10}v_{10}{\rm cm/s}$ . $f_{\rho}$ is a dimensionless (order unity) coefficient that depends on the density distribution of the ejecta, defined such that the mass averaged value of the fraction of electron energy lost to ionization is given at late time by $(t/t_{\varepsilon})^{-2}$.

Let us consider next the energy $E$ with which electrons are produced by beta-decay. In general, the decay-time is longer for lower $Q$-values, i.e. for lower energy release of the beta-decay. This behavior is due mainly to the larger phase space available for the emitted electron and neutrino. A rough lower limit to the decay time as a function of $Q$-value for $\beta^{-}$ decays is given by $\tau=10^{4}(Q/1\,{\rm MeV})^{-5}$ s (e.g. Arnould et al., 2007). The scatter in this relation is however very large- approximately two orders of magnitude at $Q=10$ MeV, and six orders of magnitude at few MeV. This spread is due to selection rules of the transitions, and it decreases with increasing $Q$-value as the number of final states increases. Thus, while we expect in general that the $Q$-value of the decays at the time scale in which we are interested would be $\sim 1$ MeV, one cannot adopt a simple relation between $Q$ (hence $E$) and $t$. A simple relation between $Q$ ($E$) and $t$ is unlikely to be valid also due to the fact that the decay of an isotope with a low $Q$-value and a long life time may produce an unstable isotope with a short life time and a high $Q$-value222For example, the 3 d, $0.5$ MeV decay of 132Te produces 132I with 2 hr, 4 MeV decay..

The very large spread in decay times for $Q\sim$ few MeV, and the possible production of short lived isotopes by the decay of longer lived ones, imply that a time independent electron release energy $E$ may be more appropriate for our discussion. A more accurate treatment would require detailed calculations of the beta-decay chains. However, the accuracy of such calculations would at present be questionable, as experimental data for relevant heavy isotopes is partial at best, and theoretical estimates may deviate from experimental data by orders of magnitude (see e.g. Arnould et al., 2007). For completeness, we nevertheless discuss in the derivation of the late time deposition rate (and bolometric luminosity), given in § 2.2, also in the case where $E$ is uniquely related to $t$. We find that the late time behavior is not sensitive to the time dependence of $E$.

Finally, let us consider the possible contribution of alpha-decays to the heating of the plasma. For low values of the initial electron fraction, $Y_{e}\lesssim 0.1$, some nuclear network models, used in calculations of the isotopic evolution of merger ejecta, predict a significant, few tens of percent, contribution of alpha-decays to the radioactive energy release at $\sim 100$ d (e.g. Barnes et al., 2016). The ionization loss rate of energetic massive particles (mass $\gg m_{e}$) is given by (e.g. Longair, 1992),

[TABLE]

where $ze$, $v$ and $\gamma$ are the particle’s charge, velocity and Lorentz factor, respectively. alpha particles are expected to be produced with $\sim 10$ MeV kinetic energy. At this energy range, the ratio between the loss rate of alpha particles and electrons is approximately given by

[TABLE]

Thus, the release of a fraction of the energy in $\alpha$-particles instead of in electrons will lead to some increase in the energy deposition rate by particles produced at $\sim 100$ d. This will, however, have only a minor effect on the energy deposition rate since, as we show in § 2.2, assuming that electrons are confined to the plasma, all electrons produced at $t\gtrsim t_{\varepsilon}$ (rather than only those produced at $t\sim 100$ d) contribute to the ionization heating at $t\sim 100$ d$\gg t_{\varepsilon}$.

2.2. Asymptotic energy deposition for confined electrons

We consider the ionization energy deposition by electrons, which are produced in beta-decays at a rate $\dot{n}(t)$ with initial energy $E_{i}(t)$. Let us denote by $E(t_{0},t)$ the time dependent energy of an electron produced at time $t_{0}$. The ionization heating rate at time $t$ is given by

[TABLE]

Here, $dn/dE$ is the number of electrons per unit energy, and we have used Eq. (4) for the ionization loss rate. The differential electron number is given by

[TABLE]

Note that $t_{0}(E,t)$ is the time at which an electron should be produced in order to have an energy $E$ at time $t$. We may therefore write

[TABLE]

In order to determine $\dot{\varepsilon}_{d}$, we need to determine $E(t_{0},t)$. The evolution of the energy of the electrons is determined by

[TABLE]

The first term on the r.h.s. accounts for adiabatic energy losses, and the second term accounts for ionization losses, following Eq. (4). As explained in § 2.1, the term describing the adiabatic losses is obtained assuming that the high energy electrons behave as an ideal gas, in which case we have $x\approx 1$ for highly relativistic electrons, and $x\approx 2$ in the highly non-relativistic limit. This is of course only an approximate description. Moreover, the electrons are mildly relativistic, with $\gamma_{e}v_{e}/c\sim 1$, and as they lose energy the value of $x$ that best describes the evolution is changing. We will solve below the evolution with a fixed value of $x$, and will show that due to the weak dependence of $g$ on $E$, the dependence of the late time evolution of the energy deposition on $x$ is small.

For simplicity, we will first discuss the case of time independent $E_{i}$, and will later provide the (straightforward) generalization to a time dependent $E_{i}$. For our analytic solutions, we will use a power-law approximation of $g(E)$,

[TABLE]

which holds with $\omega_{I}=0.15$ for $0.01<E_{\rm MeV}<1$ (see § 2.1 and Fig. 1). We show below that this is a good approximation by comparing the results of the analytic solutions to those obtained using the exact dependence of $g$ on $E$, as given by Eq. (1).

Measuring $t$ in units of $t_{\varepsilon}$, $\tilde{t}=t/t_{\varepsilon}$, and $E$ in units of $E_{i}$, $\epsilon=E/E_{i}$, Eq. (11) with the approximation of Eq. (12) becomes

[TABLE]

The evolution of the energy of an electron produced at time $t_{0}$, $E(t_{0},t)$, is determined by this equation with the initial condition $\epsilon(\tilde{t}_{0},\tilde{t}_{0})=1$.

Let us consider an electron produced at $\tilde{t_{0}}\gg 1$ (i.e. $t_{0}\gg t_{\varepsilon}$). At production, the adiabatic losses are larger than the ionization losses by a factor $\tilde{t}_{0}^{2}$, and the energy evolution of the electron may be approximated using only the adiabatic loss term in Eq. (11), yielding

[TABLE]

As the electron loses energy, ionization losses may become significant. For the energy evolution of Eq. (14), the ratio of the ionization to adiabatic losses (second to first term in Eq. (13)) is $\epsilon^{-1-\omega_{I}}\tilde{t}^{-2}\propto\tilde{t}^{-2+x(1+\omega_{I})}$.

Let us first consider the case $x\leq 2/(1+\omega_{I})$. In this case, ionization losses never become important (for electrons produced at $\tilde{t}_{0}>1$). Eq. (14) holds at all times for all electrons produced at $\tilde{t}>1$, i.e. down to energy given by Eq. (14) with $\tilde{t}_{0}\approx 1$,

[TABLE]

At time $\tilde{t}$, the electron distribution is strongly suppressed below $\epsilon_{m}(\tilde{t})$ since all electrons produced at $\tilde{t}<1$ lost most of their energy by ionization by $\tilde{t}<1$. We may therefore approximate

[TABLE]

Since the integral over $\dot{n}(t_{0})$ must converge at large $t$, we have for $\tilde{t}\gg 1$

[TABLE]

The asymptotic behavior is independent of the evolution of the radioactive energy release. This is due to the fact that, as can be seen from Eq. (16), all electrons produced at $t_{\varepsilon}<t_{0}<t$ contribute to the ionization heating at time $t$, weakly weighted by $t_{0}^{-x\omega_{I}}$ (recall that $\omega_{I}\ll 1$). Since $x\omega_{I}\ll 1$, the convergence to the asymptotic behavior is determined by the convergence of $\int dt\dot{n}$.

Fig. 2 shows the evolution of $\dot{\varepsilon}_{d}$ as a function of time, obtained by numerically solving the complete implicit analytic solution of Eq. (13) for $t_{0}(E,t)$ given in the appendix (and using it to numerically integrate Eq. (10)). The solutions converge to the asymptotic behavior given by Eq. (17) for a wide range of functional forms of $\dot{n}(t)$.

It is straightforward to generalize the above discussion for the case of a decreasing $E_{i}(t)\propto t^{-\omega_{E}}$. Normalizing the energy $E$ to $E_{i}(t_{\varepsilon})$, Eq. (13) remains unchanged, with initial conditions modified to $\epsilon(\tilde{t}_{0},\tilde{t}_{0})=\tilde{t}_{0}^{-\omega_{E}}$. Eq. (14) is replaced with $\epsilon=\tilde{t}_{0}^{-\omega_{E}}(\tilde{t}_{0}/\tilde{t})^{x}$, while Equation (15) remains unchanged (as long as $\omega_{E}<2/(1+\omega_{I})$ so that at $t\gg t_{\varepsilon}$ the ionization-loss time of electrons produced at time $t$ is longer than $t$). The result of Eq. (17) therefore also remains unchanged (as long as $\int dt_{0}\dot{n}(t_{0})t_{0}^{-(x-\omega_{E})\omega_{I}}$ converges).

Noting that $\omega_{I}\ll 1$, the regime $x>2/(1+\omega_{I})\approx 2$ is not relevant for our discussion, since $x=2$ is obtained only in the highly non-relativistic limit, while in our case the beta-decay electrons remain mildly relativistic over the relevant time. The energy of electrons produced at $\tilde{t}\sim 1$ ($t\sim t_{\varepsilon}$) with $E\sim 1$ MeV, drops to $\sim 0.1$ MeV with $v_{e}/c=0.6$ at $t=10t_{\varepsilon}$, and to $\sim 0.01$ MeV with $v_{e}/c=0.2$ at $t=100t_{\varepsilon}$. Nevertheless, we give below the results for this case as well, both for completeness and for demonstrating that the results are not strongly modified also for $x=2$.

For $x>2/(1+\omega_{I})$, as the electron cools adiabatically following Eq. (14), it eventually reaches an energy at which ionization losses become important. Comparing the two loss terms of Equation (13), we find that at time $\tilde{t}$ ionization losses dominate below

[TABLE]

Using Eq. (14), we find that electrons reaching this energy at time $\tilde{t}$ were produced at,

[TABLE]

At time $\tilde{t}$, the electron distribution is strongly suppressed below $\epsilon_{c}(\tilde{t})$ since all electrons produced at $\tilde{t}<\tilde{t}_{c}$ lost most of their energy by ionization. We may therefore approximate

[TABLE]

Since $\tilde{t}_{c}$ diverges as $\tilde{t}$ diverges, and since $\dot{n}(t)$ drops faster than $1/t$, the integral is dominated by the contribution from $t_{0}\sim t_{c}$. In this case, the asymptotic behavior does depend on the temporal evolution of $\dot{n}$ and of $E_{i}$, as the ionization heating is determined at time $t$ by electrons produced at time $t_{0}\gtrsim t_{c}(t)$.

For a power-law behavior, $\dot{n}\propto t^{-(\beta-\omega_{E})}$ and $E_{i}\propto t^{-\omega_{E}}$, and hence $\dot{\varepsilon}\propto t^{-\beta}$, we find

[TABLE]

This is the result obtained by Kasen & Barnes (2018) (see their Eq. (26)). For $x=2$,

[TABLE]

For $x=2$ and an exponential $\dot{\varepsilon}\propto\dot{n}\propto e^{-t/\tau}$, with time independent $E_{i}$, we have

[TABLE]

Since $\omega_{I}\ll 1$, we may approximate for the times of interest

[TABLE]

In order to determine $f(\omega_{I})$ an exact relation between $\tilde{t}_{c}$ and $\tilde{t}$ is required. Using Eq. (4) of the Appendix we find

[TABLE]

Fig. 3 shows the evolution of $\dot{\varepsilon}_{d}$ as a function of time, for $x=2$ and a wide range of functional forms of $\dot{n}(t)$. The solutions converge to the asymptotic behavior given by Eqs. (22) and  (23).

The convergence of the solutions to the asymptotic behavior appears slow, both for $x=1$ (Fig. 2) and for $x=2$. In Fig. 4 we show, however, that simple interpolations between the $t\gg t_{\varepsilon}$ asymptotic behavior and the $t\ll t_{\varepsilon}$ asymptotic behavior, where $\dot{\varepsilon}_{d}=\dot{\varepsilon}$, provide excellent approximations to the complete solutions. This figure also compares $\dot{\varepsilon}_{d}$ obtained using the approximation $g(E)\propto E^{-0.15}$ and that obtained using the exact $g(E)$ given by Eq. (1). The results are nearly identical.

2.3. Summary and analytic light curves

The ionization energy-loss rate of electrons is nearly independent of their energy $E$ and of the plasma composition (see Fig. 1). It is given by Eq. (4) with an effective ionization ”opacity”, $\kappa_{e}\approx 1$ cm2g*-1*, and a weak energy dependence given by $g(E)$. For the relevant energy range of $\sim 0.01$ to $\sim 1$ MeV, we may approximate $g(E)\propto E^{-\omega_{I}}$ with $\omega_{I}=0.15$, reflecting a factor 2 change in the loss rate as the energy decreases from 1 to 0.01 MeV.

For an ejecta mass $M$ and characteristic velocity $v_{M}$, neglecting the very weak dependence of $g$ on $E$ (in particular at $E\sim 1$ MeV), we have (Paper I)

[TABLE]

Here, $E=1E_{\rm MeV}$ MeV, $v_{M}=10^{10}v_{10}{\rm cm/s}$ . $f_{\rho}$ is a dimensionless (order unity) coefficient that depends on the density distribution of the ejecta, defined such that the mass averaged value of the fraction of electron energy lost to ionization is given at late time by $(t/t_{\varepsilon})^{-2}$. In Paper 1 we have considered an ejecta with a power-law dependence of mass on velocity, $dm/dv\propto v^{-(\alpha+1)/\alpha}$ for $v\geq v_{M}$, for which $f_{\rho}=[\alpha(2+3\alpha)]^{-1/2}$ (for GW 170817 observations imply $\alpha\approx 0.6$ and hence $f_{\rho}\approx 0.7$).

Eq. (26) generalizes Eq. (10) of Paper I, where $E=1$ MeV was assumed. As discussed in § 2.1, while in general the decay-time is longer for lower energy ($Q$-value) decays, with a rough lower limit to the decay time as a function of $Q$-value given by $\tau=10^{4}(Q/1\,{\rm MeV})^{-5}$ s (e.g. Arnould et al., 2007), one cannot adopt a simple relation between $Q$ (hence $E$) and $t$. The scatter in this relation is very large- approximately two orders of magnitude at $Q=10$ MeV, and six orders of magnitude at few MeV, and the decay of an isotope with a low $Q$-value and a long life time may produce an unstable isotope with a short life time and a high $Q$-value. This suggests that a time independent electron release energy $E\sim 1$ MeV may be more appropriate for the times under consideration.

If the electrons freely escape the plasma, then at $t>t_{\varepsilon}$ we simply have $\dot{\varepsilon}_{d}\approx\dot{\varepsilon}(t/t_{\varepsilon})^{-2}$. For confined electrons, we have shown that

[TABLE]

independent of the evolution of the radioactive energy release, as well as of the possible evolution of the energy $E$ with which electrons are produced. Here $x$ is the energy decay exponent describing adiabatic expansion losses, $dE/dt=-xE/t$. This description of adiabatic losses holds if the high energy electrons behave as an ideal gas, in which case we have $x\approx 1$ for highly relativistic electrons and $x\approx 2$ for highly non-relativistic electrons. This is obviously a rough approximation. However, since $\omega_{I}\ll 1$, the result is not sensitive to the exact value of $x$, and for the mildly relativistic electrons, $\gamma_{e}v_{e}/c\sim 1$, in which we are interested we may use $x=1.5$, which implies $-3+x\omega_{I}\simeq-2.8$.

The convergence of the solutions to the asymptotic behavior derived in § 2.2 appears slow, see Fig. 2 and Fig. 3. In Fig. 4 we show, however, that simple interpolations between the $t\gg t_{\varepsilon}$ asymptotic behavior and the $t\ll t_{\varepsilon}$ asymptotic behavior, where $\dot{\varepsilon}_{d}=\dot{\varepsilon}$, provide excellent approximations to the complete solutions. This figure also shows that $\dot{\varepsilon}_{d}$ obtained using the approximation $g(E)\propto E^{-0.15}$ and that obtained using the exact $g(E)$ given by Eq. (1) are nearly identical.

At times later than $t_{M}$, defined as the time at which the photon diffusion time out of the ejecta equals the expansion time $t$, the bolometric luminosity is given by $L=\dot{\varepsilon}_{d}$. For $t_{M}<t_{\varepsilon}$, as is the case for GW 170817 (see Paper I), we thus have for $t>t_{\varepsilon}$ and $\dot{\varepsilon}\propto t^{-\beta}$

[TABLE]

This corrects equations (12, A23) of Paper I to include the effects of the accumulation of electrons and of the energy dependence of the ionization loss rate.

Our results differ from those of Kasen & Barnes (2018) for two reasons. First, they approximate the energy dependence of the ionization loss rate as $g(E)\propto E^{-\omega_{I}}$ with $\omega_{I}=0.5$ (they do mention that the value may be slightly smaller, but adopt $\omega_{I}=0.5$ as the ”default” value used to determine their predictions for a power-law energy deposition, and as the only value for which results are derived for an exponential deposition; this value was later adopted by Kasliwal et al. (2018) for the analysis of the late time observations of GW 170817). Choosing $\omega_{I}=0.5$ was motivated by noting that $dE/dt\propto v_{e}dE/dX\propto v_{e}^{-1}\log(E)$ (see Eq.(1)), neglecting the $\log(E)$ dependence and using the non-relativistic approximation $v_{e}\propto E^{1/2}$. These approximations are not accurate for the energy range of interest, producing a much steeper energy dependence than that given by Eq. (1), as can be clearly seen also from examining Fig. (1). This large value of $\omega_{I}$ leads to the prediction of a much shallower decline of $L=\dot{\varepsilon}_{d}$. Second, they describe adiabatic losses using $x=2$, which is appropriate only in the highly non relativistic limit, and therefore not applicable for the current discussion. As explained in § 2.2, the result of Eq. (27) holds for $x<2/(1+\omega_{I})$, while for $x>2/(1+\omega_{I})$ the asymptotic behavior is given by Eq. (21), which reduces to Eq. (22) for $x=2$ (which is the result obtained by Kasen & Barnes, 2018). It should be noted that for $\omega_{I}\ll 1$ the asymptotic decay index is $\approx-3$, independent of $x$ and $\beta$. Hence, using Eq. (22) with $\omega_{I}=0.15$ would yield a late time decay index which is close to that of Eq. (28), in the range of $-2.8$ to $-3$ for $1.1<\beta<3$ and $\omega_{E}=0$.

3. The case of GW 170817

In this section, we discuss the implications of the late-time observations of GW 170817. The observations are discussed in § 3.1, and the analysis and implications is discussed in § 3.2.

3.1. Observations

For the early time observations ($t<20$ days) we adopt the bolometric light curve of Paper I. This bolometric light curve was estimated by polynomialy interpolating all the available observations in all the bands, and integrating over wavelength (see details in Paper I). We note that using bolometric light curves from other sources does not change our conclusions. This is demonstrated by using also the bolometric light curve presented in Kasliwal et al. (2017) (see Figure 6, § 3.2).

The main sources of the late time visible-light and IR observations ($t>20$ days) are Kasliwal et al. (2018) and Lamb et al. (2019). At late time, of order tens of days, synchrotron emission from the collisionless shock driven by the ejecta into the surrounding medium, which is responsible for the observed X-ray and radio emission, may contribute a significant fraction of the flux observed at the UV-IR bands. The removal of this component from the optical observations is discussed in § 3.1.1.

On tens of days time scale, the synchrtron luminosity increases with time while the radioactive kilonova emission decreases with time. We show that while the kilonova emission dominates at the time of the first Spitzer observations, $t=43$ day, by $t=74$ days the synchrotron contribution is estimated as $\approx 25$% with large uncertainty (due to uncertainty in modelling the synchrotron emission, which may deviate from a simple interpolation between radio and X-ray measurements). At still later times (i.e., HST observations), the observed flux is consistent with being dominated by the synchrotron component, rendering credible estimates of the kilonova emission at these late times impossible.

3.1.1 The synchrotron component

Table 1 lists all available observations taken at $t>20$ days, along with their references. Since the X-ray observations are usually obtained over a very wide band, we calculated the specific luminosity using the best fit photon power-law index $\Gamma$, defining the specific luminosity at the lower-limit of the band ($\nu_{1}$) as

[TABLE]

Here $nu_{2}$ is the band frequency upper limit.

In order to be able to remove the synchrotron component from the optical observations, we fitted the X-ray and radio specific luminosity, $L_{\nu}$ given in Table 1, with a physically motivated simple model of the form

[TABLE]

Here $\nu$ is the frequency, and $A$, $a$ and $b$ are free parameters, chosen independently for different time intervals. The best fitted $A$, $a$, and $b$ in some selected time bins are listed in Table 2. The time bins in which we performed the fit were selected by trial and error to represent periods in which the observations are consistent with this simple model.

As can be seen from Table 2, the ratio of X-ray and radio specific luminosities is independent of time, with $b\simeq-0.57$. This is consistent with synchrotron emission from a power-law distribution of electrons, $dN/dE\propto E^{-p}$, which do not radiate a significant fraction of their energy over the observed time. In this case $p=1-2b\cong 2.14$. A power-law distribution with $p\gtrsim 2$ is expected for collisionless shock acceleration (e.g. Marcowith et al., 2016), supporting the slow-cooling synchrotron emission interpretation.

The single power-law fit of Equation (30) is valid assuming that the ”cooling frequency”, $\nu_{c}$, above which radiation is produced by electrons that lose most of their energy over the observed time, is above the X-ray frequency, $\nu_{X}$. For $\nu_{c}<\nu_{X}$, a broken power-law should be fitted, with $L_{\nu}\propto\nu^{-p/2}$ for $\nu>\nu_{c}$. Such a fit would yield a smaller value of $p$, and a larger estimate of the flux in the optical/IR bands. We note that $\nu_{c}<\nu_{X}$ is unlikely since we generally expect $\nu_{c}$ to decrease with time, which will lead to variability in $b$ over time for $\nu_{c}<\nu_{X}$. This is not observed (as can be seen in Table 2). Moreover, for a smaller value of $p$ the predicted optical flux will exceed the flux measured by HST (see Fig. 5). We also note that the X-ray spectral power-law index ($\Gamma\cong 1.6$) is consistent with the inferred value of $p$.

Table 3 presents the ratio between the observed IR and visible-light flux and the synchrotron flux interpolated into the optical wavebands. The error term includes only the reported relative photometric errors without any uncertainty related to the interpolation of the synchrotron flux. At the first Spitzer observations ($t=43$ day), the kilonova emission dominates the observed flux. However, by $t=74$ days the synchrotron is contributing about 25% of the observed flux, and the exact fraction is highly uncertain due to the uncertainty in the interpolation of the synchrotron emission. At later epochs (i.e., HST observations), the observed flux is consistent with being dominated by the synchrotron emission.

3.2. Implications of the late time kilonova emission

Figure 5 presents the kilonova bolometric light curve as well as the late time Spitzer observations and some fitted models. The black circles represent the bolometric light curve derived in Paper I. The Kasliwal et al. (2018) IR observations are shown as magenta triangles, where the upper-pointed tip triangles show the integrated luminosity in the Spitzer 4.5$\mu$m band (i.e., lower limit on the bolometric luminosity), while the lower-tip triangles show $\nu L_{\nu}$. The black line shows the simple analytic model of Paper I, for which $t_{\varepsilon}=7$ d, with $L\propto t^{-2.8}$ at $t>t_{\varepsilon}$ instead of $L\propto t^{-3}$ in paper I, to correct for the accumulation of low energy electrons. The gray line shows the best fitted exponential model for $t>3$ days with exponential time scale of $3.6\pm 0.2$ days. Figure 6 shows the analytic model of Paper I overlaid on the bolometric light curve, including the late time Spitzer points, as derived in Kasliwal et al. (2018).

Prior to discussing the implication of the late time observations, the following point should be noted. On a time scale of tens of days the plasma is optically thin, far from Local Thermodynamic Equilibrium (LTE), and the emission is far from thermal. The luminosity will be produced at wavebands where atomic transitions with sufficiently strong oscillator strengths and sufficiently large population of excited states exit. Without a detailed understanding, which is currently lacking, of the oscillator strengths and excitation cross sections of the relevant atoms, the luminosity measured at a given band provides only a lower limit to the bolometric luminosity, rather than a robust estimate of it. Thus, while models that predict a bolometric luminosity lower than that observed within a single band can be ruled out, models predicting larger bolometric luminosity would be consistent with the data.

The late time IR observations are consistent with the analytic model of Paper I, for which $t_{\varepsilon}\approx 7$ d. In models with $t_{\varepsilon}\gg 7$ d, the steepening of the bolometric light curve must be explained as due to a suppression of the radioactive energy release rate, $\dot{\varepsilon}$. As pointed out in Paper I, and shown in Fig. 5, the 3-17 d data are consistent with an exponential suppression, $\dot{\varepsilon}\propto e^{-t/\tau}$ with $\tau\approx 3.6$ d. Such a suppression could occur if the radioactive energy release is dominated by the decay of isotopes with the appropriate life time. In this case, $t_{\varepsilon}\lesssim 20$ d would be required, with $L\propto t^{-2.8}$ at $t>20$ d, in order for the bolometric luminosity not to fall below the observed IR luminosity (see Fig. 5). We note that if the bolometric luminosity is larger than the observed luminosity in the 4.5$\,\mu$m band, then the upper limit on $t_{\varepsilon}$ is lower.

4. Discussion

The inefficient ionization energy loss of beta-decay electrons at $t>t_{\varepsilon}$ leads to a steepening of the bolometric light curve. $t_{\varepsilon}$ is given by Eq. (26), where $E$ is the characteristic energy with which beta decay electrons are released, and $\kappa_{e}$ is defined in Eq. (3). Its value is nearly independent of the composition of the ejecta, $\kappa_{e}\approx 1{\rm cm^{2}/g}$.

For the case where beta-decay electrons are confined to the expanding ejecta, the asymptotic behavior is given by Eq. (17), $d\log\dot{\varepsilon}_{d}/d\log t=-3+x\omega_{I}$, where $\omega_{I}\approx 0.15$ describes the dependence of the ionization loss rate on the electron energy, $dE/dt\propto E^{-\omega_{I}}$ (see § 2.1) and $x$ characterizes the adiabatic energy loss rate, $dE/dt=-xE/t$ (with $x=1$ for highly relativistic electrons and $x=2$ in the highly non-relativistic limit). For the mildly relativistic electrons of interest, we have $d\log\dot{\varepsilon}_{d}/d\log t\simeq-2.8$. The asymptotic behavior is independent of the time dependence of the radioactive energy release rate, and of the possible temporal dependence of the kinetic energy with which beta decay electrons are produced. This implies that the late time bolometric light curve does not provide strong constraints on the composition of the ejecta.

The late time IR observations of GW 170817 are consistent with the analytic model of Paper I, for which $t_{\varepsilon}\approx 7$ d, see Figs. (5) and (6). Models with larger $t_{\varepsilon}$ values may also be constructed to explain the data. In such models the steepening of the bolometric light curve at $\approx 7$ d may be explained as due to a suppression of the radioactive energy release rate, $\dot{\varepsilon}$, which may be obtained if the radioactive energy release is dominated (at $\sim 7$ d) by the decay of isotopes with $\approx 4$ d life time (see Fig. 5). The late time IR measurements may be consistent with such models if $t_{\varepsilon}\lesssim 20$ d, or, for $t_{\varepsilon}>20$ d, if isotopes with $\approx 17$ d life time and appropriate abundance dominate the energy release at late time (see Fig. 5). For $t_{\varepsilon}\gg 7$ d, the agreement of the light curve with the predicted steepening due to adiabatic losses at 7 d would be a coincidence.

At late time, the plasma is optically thin and far from LTE, and the luminosity is produced at wavebands where atomic transitions with sufficiently strong oscillator strengths and sufficiently large population of atoms in excited states exist. Inferring the total bolometirc luminosity from single band observations is therefore difficult. The fact that a considerable part of the emission at $t\sim 50$ d is at the 4.5$\mu$m band (and not detected in the 3.6$\mu$m band), provides an important handle on the composition of the plasma. Its interpretation requires a detailed understanding, which is currently lacking, of the oscillator strengths, the excitation cross sections of the relevant atoms, and the energy distribution of the beta-decay electrons.

Let us finally comment on the reasons due to which our results differ from those of earlier work. As explained in some detail at the end of § 2.1, adopting a strong energy dependence of the energy loss rate $\omega_{I}=0.5$ (and $x=2$ appropriate only in the highly non relativistic limit), Kasen & Barnes (2018) find that the late time bolometric light curve decreases slower than obtained here, and does depend on the radioactive energy release rate (as given by Eq. (22) for $\dot{\varepsilon}\propto t^{-\beta}$). Adopting these results, Kasliwal et al. (2018) use the fact that the late-time IR flux drops faster than would be predicted in this case for $\dot{\varepsilon}\propto t^{-\beta}$ with $\beta\simeq 1.3$, to constrain $\dot{\varepsilon}$ and hence the composition. As we have shown here, for weak energy dependence of the energy loss rate ($\omega_{I}\ll 1$), the late time behavior of the bolometric light curve is independent of $\dot{\varepsilon}$ and consistent with the IR observations. Moreover, as noted in the preceding paragraph, estimating the bolometric luminosity based on the IR band luminosity is highly uncertain, and cannot be used therefore to robustly constrain models.

E.W. acknowledges support by grants form the Israel Science Foundation, Minerva, and the Israeli Ministry of Technology and Science. E.O.O. is grateful for the support by grants from the Israel Science Foundation, Minerva, Israeli Ministry of Technology and Science, the US-Israel Binational Science Foundation, Weizmann-UK, and the I-CORE Program of the Planning and Budgeting Committee and the Israel Science Foundation. The solution of Eq. (13) for $\epsilon(\tilde{t}_{0},\tilde{t})$, with initial conditions $\epsilon(\tilde{t}_{0},\tilde{t}_{0})=\tilde{t}_{0}^{-\omega_{E}}$, is straightforwardly obtained by writing a differential equation for the temporal evolution of $(t^{x}\epsilon)$. The solution is

[TABLE]

For this solution,

[TABLE]

Let us consider now the asymptotic form of $\tilde{t}_{0}(\epsilon,\tilde{t})$ at $\tilde{t}\gg 1$. For $\epsilon\gg t^{-\min[x,2/(1+\omega_{I})]}$ we have

[TABLE]

with $\partial\tilde{t}_{0}/\partial\epsilon\propto\epsilon^{-1+1/(x-\omega_{E})}$. For $\epsilon\ll t^{-\min[x,2/(1+\omega_{I})]}$ we have $\tilde{t}_{0}=Const.\approx 1$ for $x<2/(1+\omega_{I})$, and

[TABLE]

for $x>2/(1+\omega_{I})$. Using Eq. (2), we find $\partial\tilde{t}_{0}/\partial\epsilon\propto\epsilon^{\omega_{I}}$ in the limit $\epsilon\ll t^{-\min[x,2/(1+\omega_{I})]}$.