Kilonova Light Curves Research Articles

Cautionary Tales on Heating-rate Prescriptions in Kilonovae

A major ingredient for kilonova lightcurves is the radioactive heating rate and its dependence on the electron fraction and velocity of the ejecta and, in principle, on the nuclear mass formula. Heating-rate formulae commonly used as the basis for kilonova models previously employed in the literature produce substantially different outputs for high electron fractions (Y e ≳ 0.3) and at late times (t ≳ 1 day) compared to newer prescriptions. Here, we employ standard semianalytical models for kilonovae with better heating rate prescriptions valid for the full parameter space of kilonova velocities and electron fractions to explore the impact of the heating rate on kilonova lightcurves. We show the dangers of using inappropriate heating rate estimates by simulating realistic observations and inferring the kilonova parameters via a misspecified heating-rate prescription. While providing great fits to the photometry, an incorrect heating-rate prescription fails to recover the input ejecta masses with a bias significantly larger than the typical statistical uncertainty. This bias from an incorrect prescription has significant consequences for interpreting kilonovae, their use as additional components in gamma-ray burst afterglows, and understanding their role in cosmic chemical evolution or for multimessenger constraints on the nuclear equation of state. We showcase a framework and tool to better determine the impact of different modeling assumptions and uncertainties on inferences into kilonova properties.

Modeling the Light Curve of HD 222925's Progenitor Kilonova

The leading theory for heavy element nucleosynthesis is that a majority of isotopes heavier than iron are formed in neutron star mergers (NSMs), visible through the kilonovae that proceed them. In this work, we derive the potential kilonova light curve associated with the elemental abundances of the “gold standard” metal-poor star HD 222925, which is thought to be enriched by a single NSM. Previous work explores which initial kilonova conditions best match the elemental abundances in HD 222925 and finds two main solutions: either the ejecta is highly constrained, or the underlying nuclear microphysics model is incomplete. We construct the light curves for these two cases in hopes of using kilonova to break this degeneracy. We find that our merger models produce light curves that are observationally too alike to robustly distinguish between the two valid models, offering little hope for constraining kilonova ejecta or the behavior of neutron-rich nuclei.

Statistical RTA simulations of atomic data for astrophysical opacity modeling in the context of kilonova emission

When considering some complex lanthanide ions characterized by a half-filled 4f subshell, the atomic structure Hamiltonian matrix sizes are so large that their diagonalization is challenging and therefore the atomic data of these ions are only used to compute the expansion opacity of a kilonova with difficulty. To avoid this problem, we propose a statistical simulation method to compute kilonova expansion opacities based on the resolved transition array (RTA) method of Bauche et al (1991 Phys. Rev. A 44 5707). The atomic structure relativistic Hartree–Fock (HFR) method has been employed to compute the radial integrals necessary for our statistical RTA simulations where the atomic data are randomly drawn using their corresponding statistical distributions and to determine the exact expansion opacities where the atomic data are obtained by the diagonalization of the Hamiltonian matrix. The statistical RTA simulations carried out for two specific ions, i.e. Sm VIII and Eu VI, for which it is still possible to diagonalize the Hamiltonian matrix, reproduce well the expansion opacities computed using HFR atomic data. Based on this good agreements, the statistical RTA method was used to compute the expansion opacity of Dy VIII, which is determined through diagonalization with difficulty. The proposed statistical RTA simulation method allows the computation of reliable astrophysical expansion opacities which are of paramount importance for kilonova light curve modeling and spectral analysis.

Rapid Generation of Kilonova Light Curves Using Conditional Variational Autoencoder

The discovery of the optical counterpart, along with the gravitational waves (GWs) from GW170817, of the first binary neutron star merger has opened up a new era for multimessenger astrophysics. Combining the GW data with the optical counterpart, also known as AT 2017gfo and classified as a kilonova, has revealed the nature of compact binary merging systems by extracting enriched information about the total binary mass, the mass ratio, the system geometry, and the equation of state. Even though the detection of kilonovae has brought about a revolution in the domain of multimessenger astronomy, there has been only one kilonova from a GW-detected binary neutron star merger event confirmed so far, and this limits the exact understanding of the origin and propagation of the kilonova. Here, we use a conditional variational autoencoder (CVAE) trained on light-curve data from two kilonova models having different temporal lengths, and consequently, generate kilonova light curves rapidly based on physical parameters of our choice with good accuracy. Once the CVAE is trained, the timescale for light-curve generation is of the order of a few milliseconds, which is a speedup of the generation of light curves by 1000 times as compared to the simulation. The mean squared error between the generated and original light curves is typically 0.015 with a maximum of 0.08 for each set of considered physical parameters, while having a maximum of ≈0.6 error across the whole parameter space. Hence, implementing this technique provides fast and reliably accurate results.

Exploring the Impact of the Ejecta Velocity Profile on the Evolution of Kilonova: Diversity of the Kilonova Lightcurves

A kilonova is a short-lived explosive event in the Universe, resulting from the merger of two compact objects. Despite its importance as a primary source of heavy elements through r-process nucleosynthesis, its nature is not well understood due to its rarity. In this work, we introduce a model that determines the density of a radially stratified relativistic ejecta. We apply the model to kilonova ejecta and explore several hypothesized velocity profiles as a function of the merger’s ejection time. These velocity profiles result in diverse density profiles of the ejecta, for which we conduct radiative transfer simulations using tardis with the solar r-process composition. Consequently, we investigate the impact of the ejecta velocity profile on the resulting evolution of the lightcurve and spectra through the line transitions of heavy elements. The change in the rate at which these elements accumulate in the line-forming region leaves its imprint on the kilonova lightcurve at specific wavelengths, causing the lightcurves to decay at different rates. Furthermore, in several profiles, plateau-like behaviors (slow and/or flat decline) are also observed. In conclusion, this work proposes potential scenarios of the evolution of kilonova due to the ejecta velocity profile.

Radioactive decay of specific heavy elements as an energy source for late-time kilonovae and potential James Webb Space Telescope observations

ABSTRACT Revealing the temporal evolution of individual heavy elements synthesized in the merger ejecta from binary neutron star mergers not only improves our understanding of the origin of heavy elements beyond iron but also clarifies the energy sources of kilonovae. In this work, we present a comprehensive analysis of the temporal evolution of the energy fraction of each nuclide based on r-process nucleosynthesis simulations. The heavy elements dominating the kilonova emission within ∼100 days are identified, including 127Sb, 128Sb, 129Sb, 130Sb, 129Te, 132I, 222Rn, 223Ra, 224Ra, and 225Ac. It is found that the late-time kilonova light curve (t ≳ 20 days) is highly sensitive to the presence of the heavy element 225Ac (with a half-life of 10.0 days). Our analysis shows that the James Webb Space Telescope (JWST), with its high sensitivity in the near-infrared band, is a powerful instrument for the identification of these specific heavy elements.

A multimessenger model for neutron star–black hole mergers

ABSTRACT We present a semi-analytic model for predicting kilonova light curves from the mergers of neutron stars with black holes (NSBH). The model is integrated into the mosfit platform, and can generate light curves from input binary properties and nuclear equation-of-state considerations, or incorporate measurements from gravitational wave (GW) detectors to perform multimessenger parameter estimation. The rapid framework enables the generation of NSBH kilonova distributions from binary populations, light curve predictions from GW data, and statistically meaningful comparisons with an equivalent binary neutron star (BNS) model in mosfit. We investigate a sample of kilonova candidates associated with cosmological short gamma-ray bursts, and demonstrate that they are broadly consistent with being driven by NSBH systems, though most have limited data. We also perform fits to the very well sampled GW170817, and show that the inability of an NSBH merger to produce lanthanide-poor ejecta results in a significant underestimate of the early (≲2 d) optical emission. Our model indicates that NSBH-driven kilonovae may peak up to a week after merger at optical wavelengths for some observer angles. This demonstrates the need for early coverage of emergent kilonovae in cases where the GW signal is either ambiguous or absent; they likely cannot be distinguished from BNS mergers by the light curves alone from ∼2 d after the merger. We also discuss the detectability of our model kilonovae with the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST).

Atomic data and expansion opacity calculations in two representative 4d transition elements, niobium and silver, of interest for kilonovae studies

Aims. Neutron star (NS) mergers are thought to be a source of heavy trans-iron element production. The latter can be detected in the spectra of the ejected materials, from which bright electromagnetic radiation is emitted. This latter is due to the radioactive decay of the produced heavy r-process nuclei and is known as kilonova. Because of their complex atomic structures – characterized by configurations involving unfilled nd or nf subshells – the heavy elements of the kilonova ejecta often give rise to numerous absorption lines generating significant opacities. The determination of the latter, which are of paramount importance for the analysis of kilonova light curves, requires knowledge of the radiative parameters of the spectral lines belonging to the ions expected to be present in the kilonova ejecta. The aim of the present work is to provide new atomic opacity data for two representative 4d elements, niobium (Nb) and silver (Ag), in their first four charge states, namely for Nb I–IV and Ag I–IV. Methods. Large-scale calculations based on the pseudo-relativistic Hartree-Fock (HFR) method were performed to obtain the atomic structure and radiative parameters while the expansion formalism was used to estimate the opacities. Results. Wavelengths and oscillator strengths were computed for several million spectral lines in Nb I–IV and Ag I–IV ions. The reliability of these parameters was estimated by comparison with the few previously published experimental and theoretical results. The newly obtained atomic data were then used to calculate expansion opacities for typical kilonova conditions expected one day after NS merger, a density of ρ = 10−13 g cm−3, and temperatures ranging from T = 5000 K to T =15 000 K.

Constraining Inputs to Realistic Kilonova Simulations through Comparison to Observed r-process Abundances

Kilonovae, one source of electromagnetic emission associated with neutron star mergers, are powered by the decay of radioactive isotopes in the neutron-rich merger ejecta. Models for kilonova emission consistent with the electromagnetic counterpart to GW170817 predict characteristic abundance patterns, determined by the relative balance of different types of material in the outflow. Assuming that the observed source is prototypical, this inferred abundance pattern in turn must match r-process abundances deduced by other means, such as what is observed in the solar system. We report on analysis comparing the input mass-weighted elemental compositions adopted in our radiative transfer simulations to the mass fractions of elements in the Sun, as a practical prototype for the potentially universal abundance signature from neutron star mergers. We characterize the extent to which our parameter inference results depend on our assumed composition for the dynamical and wind ejecta and examine how the new results compare to previous work. We find that a dynamical ejecta composition calculated using the FRDM2012 nuclear mass and FRLDM fission models with extremely neutron-rich ejecta (Y e = 0.035) along with moderately neutron-rich (Y e = 0.27) wind ejecta composition yields a wind-to-dynamical mass ratio of M w /M d = 0.47, which best matches the observed AT2017gfo kilonova light curves while also producing the best-matching abundance of neutron capture elements in the solar system, though, allowing for systematics, the ratio may be as high as of order unity.

Magnetic field effects on nucleosynthesis and kilonovae from neutron star merger remnants

ABSTRACT We investigate the influence of parametric magnetic field configurations of a hypermassive neutron star (HMNS) on the outflow properties, nucleosynthesis yields, and kilonova light curves. We perform three-dimensional dynamical space–time general-relativistic magnetohydrodynamic simulations, including a neutrino leakage scheme, microphysical finite-temperature equation of state, and an initial poloidal magnetic field. We find that varying the magnetic field strength and falloff impacts the formation of magnetized winds or mildly relativistic jetted outflows, which in turn has profound effects on the outflow properties. All of the evolved configurations collapse to a black hole ∼38–40 ms after coalescence, where the ones forming jetted outflows seem more effective at redistributing angular momentum, which result in earlier collapse times. Larger mass ejecta rates and radial velocities of unbound material characterize the systems that form jetted outflows. The bolometric light curves of the kilonovae and r-process yields that are produced by the post-merger remnant system change considerably with different magnetic field parameters. We conclude that the magnetic field strength and falloff have robust effects on the outflow properties and electromagnetic observables. This can be particularly important as the total ejecta mass from our simulations (≃10−3 M⊙) makes the ejecta from HMNS a compelling source to power kilonova through radioactive decay of r-process elements.

Tellurium emission line in kilonova AT 2017gfo

ABSTRACT The late-time spectra of the kilonova AT 2017gfo associated with GW170817 exhibit a strong emission line feature at $2.1\, {\rm \mu m}$. The line structure develops with time and there is no blueshifted absorption feature in the spectra, suggesting that this emission line feature is produced by electron collision excitation. We attribute the emission line to a fine structure line of Tellurium (Te) III, which is one of the most abundant elements in the second r-process peak. By using a synthetic spectral modelling including fine structure emission lines with the solar r-process abundance pattern beyond the first r-process peak, i.e. atomic mass numbers A ≳ 88, we demonstrate that [Te iii] $2.10\, \rm \mu m$ is indeed expected to be the strongest emission line in the near infrared region. We estimate that the required mass of Te iii is ∼10−3 M⊙, corresponding to the merger ejecta of 0.05 M⊙, which is in agreement with the mass estimated from the kilonova light curve.

Nucleosynthesis in outflows of compact objects and detection prospects of associated kilonovae

ABSTRACT We perform a comparative analysis of nucleosynthesis yields from binary neutron star (BNS) mergers, black hole-neutron star (BHNS) mergers, and core-collapse supernovae (CCSNe) with the goal of determining which are the most dominant sources of r-process enrichment observed in stars. We find that BNS and BHNS binaries may eject similar mass distributions of robust r-process nuclei post-merger (up to third peak and actinides, A ∼ 200−240), after accounting for the volumetric event rates. Magnetorotational (MR) CCSNe likely undergo a weak r-process (up to A ∼ 140) and contribute to the production of light element primary process (LEPP) nuclei, whereas typical thermal, neutrino-driven CCSNe only synthesize up to first r-process peak nuclei (A ∼ 80−90). We also find that the upper limit to the rate of MR CCSNe is $\lesssim 1~{{\ \rm per\ cent}}$ the rate of typical thermal CCSNe; if the rate was higher, then weak r-process nuclei would be overproduced. Although the largest uncertainty is from the volumetric event rate, the prospects are encouraging for confirming these rates in the next few years with upcoming surveys. Using a simple model to estimate the resulting kilonova light curve from mergers and our set of fiducial merger parameters, we predict that ∼7 BNS and ∼2 BHNS events will be detectable per year by the Vera C. Rubin Observatory (LSST), with prior gravitational wave (GW) triggers.

Superheavy Elements in Kilonovae

As LIGO-Virgo-KAGRA enters its fourth observing run, a new opportunity to search for electromagnetic counterparts of compact object mergers will also begin. The light curves and spectra from the first “kilonova” associated with a binary neutron star merger (NSM) suggests that these sites are hosts of the rapid neutron capture (“r”) process. However, it is unknown just how robust elemental production can be in mergers. Identifying signposts of the production of particular nuclei is critical for fully understanding merger-driven heavy-element synthesis. In this study, we investigate the properties of very neutron-rich nuclei for which superheavy elements (Z ≥ 104) can be produced in NSMs and whether they can similarly imprint a unique signature on kilonova light-curve evolution. A superheavy-element signature in kilonovae represents a route to establishing a lower limit on heavy-element production in NSMs as well as possibly being the first evidence of superheavy-element synthesis in nature. Favorable NSM conditions yield a mass fraction of superheavy elements X Z≥104 ≈ 3 × 10−2 at 7.5 hr post-merger. With this mass fraction of superheavy elements, we find that the component of kilonova light curves possibly containing superheavy elements may appear similar to those arising from lanthanide-poor ejecta. Therefore, photometric characterizations of superheavy-element rich kilonova may possibly misidentify them as lanthanide-poor events.

The luminosity functions of kilonovae from binary neutron star mergers under different equation of states

ABSTRACTKilonovae produced by mergers of binary neutron stars (BNSs) are important transient events to be detected by time domain surveys with the alerts from the ground-based gravitational wave detectors. The observational properties of these kilonovae depend on the physical processes involved in the merging processes and the equation of state (EOS) of neutron stars (NSs). In this paper, we investigate the dependence of kilonova luminosities on the parameters of BNS mergers, and estimate the distribution functions of kilonova peak luminosities (KLFs) at the u, g, r, i, y, and z bands as well as its dependence on the NS EOS, by adopting a comprehensive semi-analytical model for kilonovae (calibrated by the observations of GW170817), a population synthesis model for the cosmic BNSs, and the ejecta properties of BNS mergers predicted by numerical simulations. We find that the kilonova light curves depend on both the BNS properties and the NS EOS, and the KLFs at the considered bands are bimodal with the bright components mostly contributed by BNS mergers with total mass $\lesssim 3.2\, \mathrm{M}_\odot$/$2.8\, \mathrm{M}_\odot$ and fainter components mostly contributed by BNS mergers with total mass $\gtrsim 3.2\, \mathrm{M}_\odot$/$2.8\, \mathrm{M}_\odot$ by assuming a stiff/soft (DD2/SLy) EOS. The emission of the kilonovae in the KLF bright components is mostly due to the radiation from the wind ejecta by the remnant discs of BNS mergers, while the emission of the kilonovae in the KLF faint components is mostly due to the radiation from the dynamical ejecta by the BNS mergers.

Surrogate light curve models for kilonovae with comprehensive wind ejecta outflows and parameter estimation for AT2017gfo

The electromagnetic emission resulting from neutron star mergers have been shown to encode properties of the ejected material in their light curves. The ejecta properties inferred from the kilonova emission has been in tension with those calculated based on the gravitational wave signal and numerical relativity models. Motivated by this tension, we construct a broad set of surrogate light curve models derived for kilonova ejecta. The four-parameter family of two-dimensional anisotropic simulations and its associated surrogate explore different assumptions about the wind outflow morphology and outflow composition, keeping the dynamical ejecta component consistent. We present the capabilities of these surrogate models in interpolating kilonova light curves across various ejecta parameters and perform parameter estimation for AT2017gfo both without any assumptions on the outflow and under the assumption that the outflow must be representative of solar r-process abundance patterns. Our parameter estimation for AT2017gfo shows these surrogate models help alleviate the ejecta property discrepancy while also illustrating the impact of systematic modeling uncertainties on these properties, urging further investigation.

The critical role of nuclear heating rates, thermalization efficiencies, and opacities for kilonova modelling and parameter inference

ABSTRACT We present an improved version of the 3D Monte Carlo radiative transfer code possis to model kilonovae from neutron star mergers, wherein nuclear heating rates, thermalization efficiencies, and wavelength-dependent opacities depend on local properties of the ejecta and time. Using an axially symmetric two-component ejecta model, we explore how simplistic assumptions on heating rates, thermalization efficiencies, and opacities often found in the literature affect kilonova spectra and light curves. Specifically, we compute five models: one (FIDUCIAL) with an appropriate treatment of these three quantities, one (SIMPLE-HEAT) with uniform heating rates throughout the ejecta, one (SIMPLE-THERM) with a constant and uniform thermalization efficiency, one (SIMPLE-OPAC) with grey opacities, and one (SIMPLE-ALL) with all these three simplistic assumptions combined. We find that deviations from the FIDUCIAL model are of several (∼1–10) magnitudes and are generally larger for the SIMPLE-OPAC and SIMPLE-ALL compared to the SIMPLE-THERM and SIMPLE-HEAT models. The discrepancies generally increase from a face-on to an edge-on view of the system, from early to late epochs and from infrared to ultraviolet/optical wavelengths. This work indicates that kilonova studies using either of these simplistic assumptions ought to be treated with caution and that appropriate systematic uncertainties ought to be added to kilonova light curves when performing inference on ejecta parameters.

Binary neutron star merger simulations with neutrino transport and turbulent viscosity: impact of different schemes and grid resolution

ABSTRACT We present a systematic numerical relativity study of the impact of different physics input and grid resolution in binary neutron star mergers. We compare simulations employing a neutrino leakage scheme, leakage plus M0 scheme, the M1 transport scheme, and pure hydrodynamics. Additionally, we examine the effect of a sub-grid scheme for turbulent viscosity. We find that the overall dynamics and thermodynamics of the remnant core are robust, implying that the maximum remnant density could be inferred from gravitational wave observations. Black hole collapse instead depends significantly on viscosity and grid resolution. Differently from recent work, we identify possible signatures of neutrino effects in the gravitational waves only at the highest resolutions considered; new high-resolution simulations will be thus required to build accurate gravitational wave templates to observe these effects. Different neutrino transport schemes impact significantly mass, geometry, and composition of the remnant’s disc and ejecta; M1 simulations show systematically larger proton fractions, reaching maximum values larger than 0.4. r-process nucleosynthesis yields reflect the different ejecta compositions; they are in agreement and reproduce residual solar abundances only if M0 or M1 neutrino transport schemes are adopted. We compute kilonova light curves using spherically-symmetric radiation-hydrodynamics evolutions up to 15 d post-merger, finding that they are mostly sensitive to the ejecta mass and electron fraction; accounting for multiple ejecta components appears necessary for reliable light curve predictions. We conclude that advanced neutrino schemes and resolutions higher than current standards are essential for robust long-term evolutions and detailed astrophysical predictions.

Hubble constant and nuclear equation of state from kilonova spectro-photometric light curves

The merger of two compact objects of which at least one is a neutron star is signalled by transient electromagnetic emission in a kilonova (KN). This event is accompanied by gravitational waves and possibly other radiation messengers such as neutrinos or cosmic rays. The electromagnetic emission arises from the radioactive decay of heavy r-process elements synthesized in the material ejected during and after the merger. In this paper we show that the analysis of KNe light curves can provide cosmological distance measurements and constrain the properties of the ejecta. In this respect, MAAT, the new Integral Field Unit in the OSIRIS spectrograph on the 10.4 m Gran Telescopio CANARIAS (GTC), is well suited for the study of KNe by performing absolute spectro-photometry over the entire 3600 − 10 000 Å spectral range. Here, we study the most representative cases regarding the scientific interest of KNe from binary neutron stars, and we evaluate the observational prospects and performance of MAAT on the GTC to do the following: (a) study the impact of the equation of state on the KN light curve, and determine to what extent bounds on neutron star (NS) radii or compactness deriving from KN peak magnitudes can be identified and (b) measure the Hubble constant, H0, with precision improved by up to 40%, when both gravitational wave data and photometric-light curves are used. In this context we discuss how the equation of state, the viewing angle, and the distance affect the precision and estimated value of H0.

Actinide opacities for modelling the spectra and light curves of kilonovae

ABSTRACT We extend previous ab initio calculations of lanthanide opacities to include a complete set of actinide opacities for use in the modelling of kilonova (KN) light curves and spectra. Detailed, fine-structure line features are generated using the configuration-interaction approach. These actinide opacities display similar trends to those observed for lanthanide opacities, such as the lighter actinides producing higher opacity than the heavier ones for relevant conditions in the dynamical ejecta. A line-binned treatment is employed to pre-compute opacity tables for 14 actinide elements (89 ≤ Z ≤ 102) over a grid of relevant temperatures and densities. These tabular opacities will be made publicly available for general usage in KN modelling. We demonstrate the usefulness of these opacities in KN simulations by exploring the sensitivity of light curves and spectra to different actinide abundance distributions that are predicted by different nuclear theories, as well as to different choices of ejecta mass and velocity. We find very little sensitivity to the two considered distributions, indicating that opacities for actinides with Z ≥ 99 do not contribute strongly. On the other hand, a single actinide element, protactinium, is found to produce faint spectral features in the far-infrared at late times (5–7 d post merger). More generally, we find that the choice of ejecta mass and velocity have the most significant effect on KN emission for this study.

Engine-fed kilonovae (mergernovae) – I. Dynamical evolution and energy injection/heating efficiencies

ABSTRACT A binary neutron star merger is expected to be associated by a kilonova, transient optical emission powered by radioactive decay of the neutron-rich ejecta. If the post-merger remnant is a long-lived neutron star, additional energy injection to the ejecta is possible. In this first paper of a series, we study the dynamical evolution of the engine-fed kilonova (mergernova) ejecta in detail. We perform a semi-analytical study of the problem by adopting a modified mechanical blastwave model that invokes interaction between a Poynting-flux-dominated flow and a non-magnetized massive ejecta. Shortly after the engine is turned on, a pair of shocks would be excited. The reverse shock quickly reaches the wind-acceleration region and disappears (in a few seconds), whereas the forward shock soon breaks out from the ejecta (in 102–103 s) and continues to propagate in the surrounding interstellar medium. Most of the energy injected into the blastwave from the engine is stored as magnetic energy and kinetic energy. The internal energy fraction is fint < 0.3 for an ejecta mass equal to 10−3 M⊙. Overall, the energy injecting efficiency ξ is at most ∼0.6 and can be as small as ∼0.04 at later times. Contrary to the previous assumption, efficient heating only happens before the forward shock breaks out of the ejecta with a heating efficiency ξt ∼ (0.006 − 0.3), which rapidly drops to ∼0 afterwards. The engine-fed kilonova light curves will be carefully studied in Paper II.