An Analytic Description of Electron Thermalization in Kilonovae Ejecta

A simple analytic description is provided of the rate of energy deposition by $\beta$-decay electrons in the homologously expanding radioactive plasma ejected in neutron star mergers, valid for a wide range of ejecta parameters -- initial entropy, electron fraction $\{s_0,Y_e\}$ and density $\rho t^3$. The formulae are derived using detailed numerical calculations following the time-dependent composition and $\beta$-decay emission spectra (including the effect of delayed deposition). The deposition efficiency depends mainly on $\rho t^3$ and only weakly on $\{s_0,Y_e\}$. The time $t_e$ at which the ratio between the rates of electron energy deposition and energy production drops to $1-e^{-1}$, is given by $t_e=t_{0e}\Big(\frac{\rho t^3}{0.5(\rho t^3)_0}\Big)^a$, where $(\rho t^3)_0=\frac{0.05M_{\odot}}{4\pi(0.2c)^3}$, $t_{0e}(s_0,Y_e)\approx17$ days and $0.4\le a(s_0,Y_e)\le0.5$. The fractional uncertainty in $t_e$ due to nuclear physics uncertainties is $\approx10\%$. The result $a\le0.5$ reflects the fact that the characteristic $\beta$-decay electron energies do not decrease with time (largely due to"inverted decay chains"in which a slowly-decaying isotope decays to a rapidly-decaying isotope with higher end-point energy). We provide an analytic approximation for the time-dependent electron energy deposition rate, reproducing the numerical results to better than $50\%$ (typically $<30\%$, well within the energy production rate uncertainty due to nuclear physics uncertainties) over a 3-4 orders-of-magnitude deposition rate decrease with time. Our results may be easily incorporated in calculations of kilonovae light curves (with general density and composition structures), eliminating the need to numerically follow the time-dependent electron spectra. Identifying $t_e$, e.g. in the bolometric light curve, will constrain the (properly averaged) ejecta $\rho t^3$.

Interpreting kilonovae observations requires understanding the ★ E-mail: ben.shenhar@weizmann.ac.il processes by which secondary particles produced in radioactive decay deposit energy, i.e., "thermalize," in the ejected plasma.Radioactivity produces energetic particles (photons, electrons, alphas, and fission fragments).Initially, the particles' energy is wholly absorbed in the plasma.However, as the ejecta expand and dilute, the energy is only partially absorbed, leading to an energy deposition rate that deviates from the total energy release rate.-rays initially dominate the radioactive heating, up to  ∼ 1 day at which time they escape the plasma (Barnes et al. 2016;Guttman et al. 2024).On a timescale of  ∼ 1 − 10 days, the main heating source is expected to be -decay electrons, with possible  and fission heating contributions depending on initial   and nuclear mass model: For initial   ≳ 0.2, -decay electron heating is expected to dominate the heating for all  ≳ 1 day, with a possible comparable contribution from fission heating at  ∼ 100 days for some nuclear mass models (e.g Zhu et al. 2021); For initial   ≲ 0.2, mass model uncertainties yield larger variance in the heating contributions--particle heating was found to be of the same order as -decay electron heating for  ≲ 10 days for some mass models (e.g.Zhu et al. 2021; Barnes et al. 2021) and fission fragments may be equally important on timescales of 10 ≲  ≲ 100 for most mass models investigated (However, the fission contribution is highly uncertain due to uncertainties in theoretical mass models, fission barriers, and fission rates, e.g.Zhu et al. 2021; Barnes et al. 2021).Therefore, notwithstanding the effects of fission, late-time kilonovae light-curve interpretation requires an understanding of the inefficient thermalization timescale on which electrons efficiently deposit their energy, particularly for compositions with initial   ≳ 0.2.In this work, we calculate the time-dependent energy deposition rate per unit mass,  dep , by electrons and positrons (in what follows, we refer to both as "electrons") emitted by -decays in homologously expanding plasma ejected in binary neutron star (BNS) mergers for a wide range of ejecta parameters, { 3 ,  0 ,   }.Using detailed numeric nucleosynthesis calculations, we obtain the time-dependent composition, stopping power, and -decay energy spectra, and follow the evolution of the electron energy distribution, assuming that the electrons are confined by magnetic fields to the fluid element within which they were produced 1 (i.e., including "delayed energy deposition", the deposition at time  by electrons that were emitted at  ′ <  and did not lose all their energy by the time ).The secondary electrons that are ionized by -decay electrons have very low kinetic energies and thermalize rapidly, we thus assume that they fully thermalize upon emission.Finally, we test the sensitivity of our results to uncertainties in nuclear physics.
We carry out calculations over a broad range of values of ejecta parameters { 0 ,   ,  3 }, that generously covers the expected ejecta parameter range from binary neutron star merger simulations (e.g.Radice et al. 2018;Nedora et al. 2021) : 10 −3 ≤  3 /( 3 ) 0 ≤ 102 , where ( 3 ) 0 = 0.05 ⊙ 4  (0.2) 3 is our benchmark value motivated by the inferred mass and characteristic velocity of the ejecta of GW170817 (Drout et al. 2017;Hotokezaka et al. 2018;Kasen et al. 2017;Kasliwal et al. 2017;Waxman et al. 2018) (for homologous expansion,  3 is a constant of the ejecta determined by its mass and velocity spread); 1  /baryon ≤  0 ≤ 100  /baryon; 0.05 ≤   ≤ 0.45.Our results may be readily applied to ejecta with general density and composition structures by properly integrating over the ejecta components.Under the assumption that electrons are confined to the plasma by magnetic fields, such integration would be valid also at late times, >   , at which delayed deposition may become important.
We provide a simple analytic description of the thermalization efficiency that is valid for a wide range of ejecta parameters, { 3 ,  0 ,   }, extending the results of earlier works, that compile tables of parameterized fits to the thermalization efficiency for specific combinations of ejecta mass and velocity (e.g.Barnes et al. 2016).Although our calculated   values are broadly consistent with those of earlier works (Metzger et al. 2010;Barnes et al. 2016;Hotokezaka & Nakar 2020;Kasen & Barnes 2019), our analytic approximation differs from earlier analytic results.This is partly due to our calibration against a systematic numeric analysis, and partly due to assumptions adopted in earlier work, which we find invalid (see § 4).In particular, it is commonly assumed that the characteristic energy of electrons released in -decays drops with time (e.g.Kasen & Barnes 2019;Hotokezaka & Nakar 2020), leading to a dependence of   on  3 which is stronger than the   ∝ ( 3 ) 1/2 scaling expected for timeindependent characteristic energy (see discussion in § 4.2).While the -value of -decays decreases with decay half-time , this does not imply (as pointed out in Waxman et al. 2019) a decrease of the characteristic -decay electron energy with time , since there is a very large spread in the relation between  and  and since the decay of an isotope with a low -value and a long lifetime may produce an unstable isotope with a short lifetime and a high -value.We find here 1 Magnetic confinement appears likely since it would be facilitated by a relatively weak magnetic field.Assuming  ∝  −2 in the expanding ejecta, as may be appropriate for a tangled field,  ∼ 0.1G is expected at  ∼ 10 days,  ∼ 10 16 cm (the equipartition field at 10 days is ∼ 1 mG), implying an electron Larmor radius of ∼ 10 10 cm. that such "inverted decay chains" lead to a nearly time-independent, or slightly rising with time, characteristic -decay electron energy, and   ∝ ( 3 )  with  ≲ 0.5.
The plasma stopping power depends on the ejecta composition and on the energy spectrum of the -decay electrons, which in turn depend on   (and  0 ).  is therefore often expected to depend strongly on   and nuclear physics inputs (Barnes et al. 2016(Barnes et al. , 2021;;Zhu et al. 2021).We find that the deposition efficiency, particularly   , depends primarily on  3 , with only a weak dependence on {  ,  0 } and on nuclear physics uncertainties.This is due to the fact that the ionization stopping power.which dominates the energy loss, depends on the composition mainly through the ratio of atomic number and mass, /, which is nearly composition independent, and since the characteristic -decay electron energy depends weakly on composition.
This paper is organized as follows.In § 2, we describe our nucleosynthesis and electron spectra calculations.In § 3, we give the equations we used for (numerically) calculating the electron energy loss and deposition, and the definition of   .§ 4 presents our results.In § 4.1 we present our analytical formulae for   .In § 4.2 we discuss the characteristic energy of -decay electrons.In § 4.3 we present our analytical formulae for the energy deposition rate     ().§ 4.4 presents an analysis of the sensitivity of our results to nuclear physics uncertainties.We summarize and discuss our conclusions in § 5.
We compute time-dependent yield abundances using the publicly available SkyNet nuclear network, including 7843 isotopes up to 337 Cn (Lippuner & Roberts 2017).SkyNet requires time-dependent trajectories of Lagrangian fluid elements to predict the temporal evolution of the abundances.Matter density evolves first through an exponential phase which is largely constant, before transitioning to a homologous expansion (Lippuner & Roberts 2015): where  is the expansion timescale of the ejecta.NSE is assumed to hold until  = 7GK, at which time SkyNet switches to a full network evolution.Most works that examine SkyNet nucleosynthesis yields initialize the calculation with given ,  0 ,   ,  0 (Perego et al. 2022;Lippuner & Roberts 2015), where  0 is the initial temperature of the ejecta.For these given values, the EoS determines the initial density  0 (see Lippuner & Roberts (2017) regarding EoS and NSE solver used by SkyNet) and then uses the density history provided by Eq.
(1) to evolve the network.
In contrast, we initialize our nucleosynthesis calculations as follows: All runs are initialized with  0 = 10 GK, as NSE is expected for this temperature.We supply initial values of  0 ,   , and then find  0 using SkyNet's EoS solver.We determine  from  0 using Eq. ( 1), such that the ejecta has the desired  3 value.We stress that in this approach  is not a meaningful physical parameter; we assume that the ejecta has a particular  3 ,  0 and was initially in NSE, in which case  specifies the time before expansion.Initializing with different  0 would result in a different , however as NSE conditions hold for both, the same composition is achieved.This alternative approach allows us to directly probe the dependence on  3 for fixed  0 within the range obtained in BNS simulations.
We employ the latest JINA REACLIB database (Cyburt et al. 2010), with specific corrections as described in Appendix B of Guttman et al. (2024), using the same setup as in Lippuner & Roberts (2015) for the other input nuclear physics.More specifically, strong inverse rates are computed assuming detailed balance.Spontaneous and neutron-induced fission rates are taken from Frankel & Metropolis (1947); Panov & Janka (2009), adopting fission barriers and fission fragment distributions from Mamdouh et al. (2001);Wahl (2002).Nuclear masses are taken from the REACLIB database, which includes experimental masses and theoretical masses calculated from the FRDM model (Möller et al. 2016).
We use experimental data from the latest ENDF database (Brown et al. 2018) to calculate the time-dependent energy released in electrons.We ensure that the total energy release we compute is equivalent to  2 as calculated by SkyNet.We use BetaShape (Mougeot 2017) to calculate spectra of individual -decays.We then compute the full time-dependent emission spectra using the calculated activities of -decaying isotopes.

RADIOACTIVE HEATING
In this section, we first present the energy-loss mechanisms for electrons, and then describe the energy deposition calculation method and the definition of   .

Electron energy loss
-decay electrons emitted in the 10 −1 − 1 MeV range lose energy primarily through ionization losses.However, for initial entropies  0 ≈ 10 2   /baryon, the ejecta composition is dominated by H and/or He (Perego et al. 2022).In this case, plasma losses dominate over other energy loss mechanisms.
We calculate mass-weighted, time-dependent stopping powers based on the instantaneous composition of the ejecta, Where   iso is the energy loss per unit column density for a particular isotope,  iso is the isotope's mass number and  iso is its number fraction, defined such that iso  iso = 1.For electron ionization losses, we use (Longair 2011) where   / is the energy loss per unit column density traversed by the electron,  = ,   is the electron's velocity,   is the Lorentz factor of the electron,  and  are the atomic and mass numbers of the plasma nuclei, and Ī is the effective average ionization potential which can be approximated for an element of atomic number Z as (Segre 1977): Ī = 9.1 1 + 1.9  2/3 eV.For plasma losses, we use (Solodov & Betti 2008) where   is the number of free electrons per atom,   = (4 2   /  ) 1/2 is the plasma frequency, and   =   /(   ) is the number density of free electrons.We use   = 1 and ℏ  = 10 −7 eV (we keep   fixed at this value since it affects the stopping power only logarithmically).Finally, at highly relativistic energies, Bremsstrahlung losses significantly contribute to electron losses (Longair 2011) Figure 1 shows the electron stopping power contributions for Xe composition.

Deposition Equations
The energy deposited by electrons at time  is given by where is the number of electrons per unit energy, determined by the continuity equation in energy space Here,    is the number of electrons released per time per unit energy (as determined from the time-dependent abundances), and  is the full energy-loss rate for electrons, given by The first term on the right-hand side accounts for adiabatic energy losses and the second term accounts for energy transferred to the plasma as described in § 3.1.The adiabatic losses are obtained assuming that the high-energy electrons behave as an ideal gas, in which case we have  ≈ 1 for highly relativistic electrons and  ≈ 2 in the highly non-relativistic limit.-decay electrons are mildly relativistic, with     / ∼ 1, and as they lose energy, the value of  that best describes the evolution is changing.We take a fixed  = 1 throughout all our calculations.We numerically solve Eq. ( 7) for   (, )  and use the result in order to integrate Eq. (6).

Instantaneous Deposition and 𝑡 𝑒
We define the fraction of energy instantaneously deposited at time  as where   is the energy production rate in -decay electrons, and the fraction of energy instantaneously deposited by an electron produced at time  with energy   is approximated as Here is the deposition energy loss timescale.Under this approximation, deposition deviates from full efficiency when the deposition energy loss timescale becomes larger than the dynamical time .Note that at late times   (  , ) ∝  −2 because   ∝  −1 ∝  3 .  is defined as In Figure 2 we plot  inst (/  ) for different ejecta parameters.
At late times  inst () approaches a 1/ 2 behavior, as expected from Eq. ( 10), although deviations from this behavior may occur due to variations in the emitted electron spectra.The rather common behavior of instantaneous deposition for different ejecta parameters suggests that an approximate analytic description may be obtained.
In Figure 3, we compare the rate of electron energy release and deposition.At early times, when  <   ,  dep <  inst .This is due to energy that is lost to adiabatic expansion (first term on righthand side of Eq. ( 8)).At  =   ,  inst ≃ 1.5  dep , with the exact factor depending on ejecta parameters.We are therefore confident in using our definition of   as the inefficient deposition timescale, even though it is defined with respect to instantaneous deposition.At  ∼ 2  ,  dep >  inst as delayed deposition starts to dominate.At  ≫   ,  dep ∼  −2.85 regardless of   (), as was obtained analytically in Waxman et al. (2019).
Table 1 presents the best-fit parameters for Eq. ( 12) for the 3 regions considered.In Figure 6 we plot   ( 3 ) for all three regions, together with their respective broken power-law fits from Table 1.Error bars show the standard error of the estimate, defined as where   is the value of the fit, calculated according to Eq. ( 12), and  , is calculated according to Eq. ( 11) for the different   and  0 values considered.The fit for Region I exhibits the best uniformity, with a fractional standard error (defined as   ≡     ) of   < 0.1 for all  3 values excluding  3 ( 3 ) 0 = 10 −3 , for which   ≈ 0.2.For Regions II and III,   ≈ 0.1 − 0.2, with the exact value depending on the value of  3 .
We find that  1 = 1/2 for our entire parameter space, while  2 = 0.37, 0.42, 0.5 for Regions I, II and III, respectively.As noted in the introduction, the fact that the power-law indices  1,2 are close to 1/2 implies that the characteristic energy of the -decay electrons does not vary significantly with time, and the lower than 1/2 values suggest a minor increase in this energy with time.We discuss this further in § 4.2.

The characteristic energy of 𝛽-decay electrons
The relation between the dependence of   on  3 and the dependence of the characteristic -decay electron energy on time may be qualitatively understood as follows.  is approximately the time at which the expansion time, , equals the deposition (mainly by ionization) energy loss time,   =   (/) −1 where   is the initial electron energy.Noting that   / is only weakly dependent on energy at the relevant energy range (see Figure 1), and For characteristic -decay electron energy that varies with time as   ∝  − , we have   ∝ ( 3 ) 1/2/(1−/2) .Figure 7 shows the characteristic -decay energies vs. time for several  3 values (with  0 = 20   /baryon).Motivated by Eq. ( 14), we define the characteristic energy as  −1/2 −2 , rather than ⟨⟩, where the brackets denote an average over the electron energy spectra (weighted by the electron energy).Several conclusions may be drawn from this figure.Firstly, only for   = 0.4 does the characteristic energy monotonously decrease with time up until  ∼ 100 days.This is consistent with Figure 4, where only the   = 0.4 curve exhibits a logarithmic slope (of the dependence of   on  3 ) that is > 1/2 at late times.For   = 0.05, 0.1, 0.15, we see a sharp rise in characteristic energy around  ∼ 20 days.For these   values, we find that the radioactive heating at  ≳ 30 days becomes dominated by the following decay-chain 194 Os This is an example of an "inverted decay-chain" in which a slowly decaying isotope decays to a rapidly decaying isotope with higher energy release.We find that at these times, other inverted decay-chains are also active and contribute significantly to the heating.These include  = 140, 132, 106 and others, depending on   .We note that these occur only for even values of , and may be understood using the semi-empirical mass formula (Weizsäcker 1935).The pairing term in the formula creates two mass parabolas for a given : one for even-even and another for odd-odd nuclei.As the nucleus decays, it transitions between these parabolas, resulting in uneven energy  releases until reaching stability.Overall, there are 40 inverted decaychains for which a slowly decaying isotopes decays into a rapidly decaying isotopes with half-life > 10 2  1/2 of its parent isotope, 26 of these are for decay chains with  ≥ 90.
In Region I (  < 0.22 & ∀ 0 , see Table 1), third peak −process elements are robustly produced.Therefore, the inverted decay chains for  = 132, 140, 194 are active, leading to an increase in late-time characteristic energy release and, thus, a lower value of the exponent,  2 = 0.37.In Region II, the situation is more subtle.We note that in NSE, the entropy  0 scales monotonically with the photon-to-baryon ratio .The abundance of an isotope (, ) in NSE is proportional to ∼  1−  and thus NSE with high entropy favors more free neutrons and less heavy-nuclei (Meyer 1993).As the ejecta expand and exit Figure 6.  ( 3 ) for the 3 regions in parameter space as defined in Figure 5, together with broken power-law fits as defined in Eq. ( 12) and Table 1.Different colors correspond to different  0 values and different markers correspond to different   values.Error-bars show the standard-error of the estimate.
NSE, the many free neutrons are captured on the few seed nuclei, leading to heavy-element production.Thus for   = 0.25, 0.3, third peak elements are produced, albeit less than in Region I by a factor ∼ 10 −3 , despite the ejecta being relatively neutron-poor.The inverted decay chains for  = 132, 140, 194 are similarly active.For   = 0.35, 0.4, -process elements are produced up to the second peak (up to  ≈ 125).Despite the high entropy, the ejecta is too neutronpoor to create second and third-peak -process elements.For these values of   in Region II, the inverted decay chains for  = 106, 125 dominate the heating from  ∼ 10 days, leading to an increase in characteristic initial electron energy.For these same high   values but for low values of  0 (Region III), isotopes up to  ≈ 90 are synthesized, and thus, none of the significant inverted decay chains are active.In Region III, the characteristic initial electron energy, for the most part, either stays constant or slightly declines over time, yielding  2 = 1/2 as is expected from Eq. ( 14).
We conclude that inverted -decay chains are not insignificant anomalies that can be disregarded when considering characteristic energies of -decay electrons.Their presence in nearly all the nucleosynthesis runs considered, especially those with significant heavy element production, invalidates the assumption that the characteristic electron energy steadily declines with time.Notable rises in characteristic energies occur for   = 0.05, 0.10, 0.15 in which 3rd-peak elements are produced.This is due to the activity of "inverted decay-chains" in which slowly decaying isotopes decay to quickly decaying isotopes with higher  −1/2 −2 .Their presence impacts the characteristic electron energy for nearly all nucleosynthesis calculations considered.

Analytic approximation for 𝑄 dep
In this section we provide an interpolating function  int dep that describes the electron energy deposition rate both at early ( <   ) and late ( >   ) times, for different values of { 3 ,   ,  0 }.Since  dep ∼  −2.85 at late times regardless of   (see Waxman et al. (2019)), we provide an interpolation for  dep rather than for the fractional deposition,  dep /   .We suggest the following interpolation: Here, ,  1 ,  2 ,   are free parameters to be fitted for and   is given by Eq. ( 12). governs the transition from the early-time to late-time interpolation. 1 governs the shift from full to partial deposition until   . 2 governs the late-time asymptotic deposition, which should approach a  −2.85 behavior.  is the time at which the delayed deposition overtakes the instantaneous deposition, and should be roughly ∼ 1 − 2 ×   .The first term in  late ensures that for  ≪    early <  late such that  int dep interpolates correctly between the early/late-time regimes.
We fit {,  1 ,  2 ,   } values for the three different regions in {  ,  0 } as described in § 4.1.Our fitting employed a gradient descent over these four parameters while minimizing the cost function defined as cost = 1   =1 log  int dep (  )/  dep (  ) , where  = 40,   are logarithmically spaced from   = 0 up to   = 6  .The best-fit parameters are compiled in Table 1.The values obtained for  2 and   are  2 ≈ 2.85,   = 1.3  , as is expected for the correction transition to the late-time asymptotic heating.
In Figure 8, we compare  int dep and  dep for the same sample of parameters as shown in Figure 3.In Figure 9 we check the accuracy of the analytic description of  dep for all nucleosynthesis parameters considered.Our interpolation is most accurate in Region I, within which the nucleosynthesis yield is mostly uniform.In Regions II and III, different parameters result in different nucleosynthesis yields, and therefore different deposition curves, which reduces the accuracy of the fit.The quality of the fit and its spread may be gauged by considering the standard deviation of   dep /  dep at each time, for each region.Averaging the standard deviation over logarithmic time for each region, we find that σ ≈ 0.1, 0.29, 0.18 for Regions I, II, and III, respectively.Overall, our interpolation is nearly accurate to within 50% over 3 − 4 orders-of-magnitude change in  dep for our entire parameter space.Note that Eq. (12, 16) provide an analytic description that requires   () as input in order to determine the full electron deposition up to  ∼ 10   .

Nuclear Physics Robustness Checks
Modeling -process nucleosynthesis from BNS mergers requires calculating nuclear masses and nuclear reaction rates for which no experimental data are available.For   < 0.25, different nuclear mass models may result in orders-of-magnitude differences in final composition of post 1st-peak elements ( ≳ 100) and up to an orderof-magnitude difference in radioactive energy release (Mumpower et al. 2016;Zhu et al. 2021;Barnes et al. 2021).In this section we study the sensitivity of our results to the nuclear physics models' uncertainties.
To gauge the impact of these uncertainties, we reran our calculations for two different nuclear mass models together with randomized nuclear rates.We edited our local REACLIB file which contains all the rates and cross-sections that go into SkyNet.We modify every strong or weak theoretical interaction rate  to where  ∈ [10 −2 , 10 2 ] was randomly chosen from a uniform distribution in log().This amounts to changing ∼70,000 rates (approximately 90% of all rates included).With our new reaction rates file, we reran   for a subset of initial parameters:  3 /( 3 ) 0 = 0.1, 1, 10,  0 = 20, 70   /baryon,   = 0.1, 0.25, 0.4 for both FRDM and UNEDF1 (Kortelainen et al. 2012, http: //massexplorer.frib.msu.edu,)mass models.We repeat the above process 100 times and calculate   for each run and compare the values obtained to our nominal values,  , presented in Figures 4. We stress that these changes lead to orders-of-magnitude differences in final composition between different runs.
In Figure 10, we show histograms of   / , for different values of  3 and   , where   is calculated according to Eq. ( 11) for each SkyNet run with randomized reaction rates using both FRDM and UNEDF1 mass models.We see that the mean falls at ≈ 1 for all parameter values considered.The spread in   values is small, with an average fractional standard deviation of  , ≈ 0.12 for   < 0.4.For   = 0.4, the spread is even smaller, with  , ≈ 0.03.
These tests demonstrate that our results are only weakly sensitive to nuclear physics uncertainties.Different nuclear mass models and different reaction rates may lead to large variations in composition  6)) respectively.The dotted line shows our analytic approximation for  dep ,  int dep (Eq.( 16) and Table 1). 0 = 20, (60)   /baryon for top (bottom) row.Left, middle, and right panels correspond to   = 0.05, 0.2, 0.45, respectively.Blue, green, and red lines correspond to  3 ( 3 ) 0 = 10 −1 , 1, 10 respectively.Left panel: Ejecta parameters for each region in parameter space of {  ,  0 }.Middle panel: Fitted broken power-law parameters for   ( 3 ) as defined in Eq. ( 12).Right panel: Fitted interpolation parameters as defined in Eq. ( 16).The regions are those outlined in Figure 5.The accuracy of the interpolation is examined in Figure 9. and heating rates.However, these changes do not significantly affect the inefficient thermalization timescales.In retrospect, this was to be expected -Figure 4 shows that   is relatively insensitive to the   value and thus to the detailed electron energy spectra and ejecta composition.

CONCLUSIONS
We calculated the time-dependent energy deposition rate per unit mass,  dep (see Eq. ( 6)), by electrons (and positrons) emitted by -decays in homologously expanding ejecta produced in BNS mergers for a wide range of ejecta parameters, { 3 ,  0 ,   }.Using detailed numeric nucleosynthesis calculations, we obtained the timedependent composition, stopping power, and -decay energy spectra, and followed the evolution of the electron energy distribution, assuming that the electrons are confined by magnetic fields to the fluid element within which they were produced (see Eq. ( 8)).
The deposition efficiency, the ratio  inst between the rates of instantaneous electron energy deposition and energy production rates, Eq. ( 10), depends mainly on  3 and only weakly on { 0 ,   }.   , the time at which  inst drops to 1 −  −1 , is well described by a broken   is calculated according to the broken power-law fits as defined in Eq. ( 12) and Table 1. int dep is calculated according to Eq. ( 16) with the best-fit parameters for each region as compiled in Table 1.
The result  1,2 ≤ 0.5 reflects the fact that contrary to naive expectations, the characteristic energy of the -decay electron spectrum does not decrease with time (see Eq. ( 14)).This is largely due to "inverted decay chains" in which a slowly-decaying isotope decays to a rapidly-decaying isotope with higher end-point energy.A detailed discussion of the relevant decay chains at different {  ,  0 } values is given in § 4.2.We find that inverted decay chains play a significant role at  ≳ 10 days, affecting the thermalization timescale at all relevant {  ,  0 } values.
Using our analytic results for   , we provided an analytic description of the energy deposition of -decay electrons (and positrons),  int dep () given by Eq. ( 16) and Table 1.Our formulae provide an analytic description of the thermalization efficiency, which we find to be only weakly sensitive to { 0 ,   } (see § 4. 1,4.2,and Figures 4 ,5 ,6) and nuclear physics uncertainties (see § 4.4, Figure 10).This enables one to determine  int dep () for a given total radioactive energy release rate,   ().Our analytic description reproduces the numerical results to better than ≈ 50% (typically < 30%) accuracy  In contrast to previous works that compile tables of parameterized fits for the thermalization of specific combinations of ejecta mass and velocity in a specific velocity structure (e.g.Barnes et al. (2016)), we provide a simple analytic description that can be utilized for a wide range of ejecta parameters and density profiles.Our results may thus be easily incorporated in calculations of kilonovae light curves (with general composition, density and velocity structures), eliminating the need to follow the time-dependent electron spectra numerically.
The weak dependence of the deposition efficiency on {  ,  0 } and on nuclear physics uncertainties implies that identifying   in the bolometric kilonova light curve will constrain the (properly averaged) ejecta  3 .

Figure 2 .
Figure 2. Instantaneous deposition fraction  inst (/  ) as defined in Eq. (9) for different ejecta parameters.Different colours correspond to different values of  3 and different line-styles correspond to different   values.Dotted, dashed, dashed-dotted, and solid lines correspond to   = 0.05, 0.15, 0.25, 0.4, respectively. 0 = 20   /baryon for all calculations shown.The instantaneous deposition fractions follow similar trajectories, especially before   .

Figure 7 .
Figure 7. Top, middle, and lower panel show  −1/2 −2 plotted vs. time for  3 /( 3 ) 0 = 2, 10, 100, respectively.Different colors correspond to different   values, and  0 = 20   /baryon for all runs considered.We see that only for   = 0.4 does  −1/2 −2 monotonously decrease with time.Notable rises in characteristic energies occur for   = 0.05, 0.10, 0.15 in which 3rd-peak elements are produced.This is due to the activity of "inverted decay-chains" in which slowly decaying isotopes decay to quickly decaying isotopes with higher  −1/2 −2 .Their presence impacts the characteristic electron energy for nearly all nucleosynthesis calculations considered.

Figure 8 .
Figure 8. Electron energy release and energy deposition rates as a function of time.Solid and dashed lines correspond to the numerically calculated   and  dep (Eq.(6)) respectively.The dotted line shows our analytic approximation for  dep ,  int dep (Eq.(16) and Table1). 0 = 20, (60)   /baryon for top (bottom) row.Left, middle, and right panels correspond to   = 0.05, 0.2, 0.45, respectively.Blue, green, and red lines correspond to

Figure 9 .
Figure 9.  int dep /  dep vs /  for entire parameter space of { ,  3 ,  0 }, for the three regions in parameter space as defined in Figure 5. Colors indicate different   values. is calculated according to the broken power-law fits as defined in Eq. (12) and Table1. int dep is calculated according to Eq. (16) with the best-fit parameters for each region as compiled in Table1.

Table 1 .
Fitted parameters for an analytic description of   and  int dep .