Exploring the Impact of Ejecta Velocity Profile on Kilonova Evolution: 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 lightcurve and spectral evolution 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 kilonova evolution due to the ejecta velocity profile.


INTRODUCTION
Mergers of two neutrons stars (NSs) or a neutron star and a stellar-mass black hole (BH) eject a small fraction of matter (10 −3 − 10 −1 M ⊙ ) into space at (subrelativistic) speeds ≲ 0.35 c (Kasen et al. 2017;Metzger 2019).The neutron-rich ejecta material rapidly decompresses from high densities, and the r-process nucleosynthesis forges heavy elements in the ejecta.Subsequently, these elements decay radioactively and become a long-term heat source for the ejecta lasting weeks, thus powering the kilonova emission.The kilonova emission is thought to be attributed to more than one ejecta component (e.g., Kasen et al. 2017;Perego et al. 2017;Siegel & Metzger 2017;Tanaka et al. 2017;Radice et al. 2018).The dynamical ejecta is expected to have an ejecta speed of v ∼ 0.2 − 0.3c with neutron-rich material (electron fraction Y e ≲ 0.1), which is ejected over a timescale of few milliseconds.On the other hand, the disk winds have v ∼ 0.1c with a higher electron fraction (Y e ∼ 0.2 − 0.3) on longer timescales of several hundreds of milliseconds.* Corresponding author: Z. Lucas Uhm In the study of supernovae, the density profile of the ejecta envelope and its composition are key ingredients to determine the evolution of a supernova at various wavelengths.For Type Ia supernovae, a quasiexponential profile known as the W7 model (Nomoto 1984;Branch et al. 1985;Thielemann et al. 1986) has found extensive application since the model effectively reproduces the distribution of elements with high accuracy (e.g., Jeffery et al. 1992;Mazzali et al. 1993;Dutta et al. 2021).While a power-law profile is widely used for Type II supernova studies (e.g., Blinnikov et al. 2000;Dessart & Hillier 2006), various other profiles have also been adopted (e.g., a combination of distinct inner and outer profiles; Nagy et al. 2014;Szalai et al. 2016).
In the case of kilonovae, a consensus regarding the ejecta density profile has not yet been reached.To explain the observational data of kilonovae, various functional forms of the density profile, heavily relying on our understanding of supernovae, have been adopted; e.g., a constant density profile with a sharp cutoff (Darbha & Kasen 2020), a broken power-law (Kasen et al. 2017), a power-law (Tanaka et al. 2017;Watson et al. 2019;Gillanders et al. 2022), and multi-dimensional anisotropic profiles (e.g., Bulla 2019;Kawaguchi et al. 2020;Breschi et al. 2021;Heinzel et al. 2021;Korobkin et al. 2021;Wollaeger et al. 2021).
Deriving the kilonova ejecta profile would be more complicated, as the mass ejection from the merger and the properties of the ejected flow depend on the binary parameters and the NS equation of state (EOS).Such details have been vigorously studied with hydrodynamic and numerical simulations (e.g., Hotokezaka et al. 2013;Tanaka & Hotokezaka 2013;Kyutoku et al. 2015;Siegel & Metzger 2017;Radice et al. 2018;Krüger & Foucart 2020;Nedora et al. 2022).For example, Kyutoku et al. (2015) simulated the merger of black hole-neutron star binaries and showed that depending on key parameters such as the EOS, mass ratio, and dimensionless spin parameter a range of complicated ejecta profiles can be obtained.On the other hand, Radice et al. (2018) performed numerical simulations on the merger of binary neutron stars, incorporating the EOS and binary parameters.From the simulations, they extracted dynamical ejecta properties and then derived the density profiles of ejecta as a function of the entropy and electron fraction of materials.
The complexity behind the merging process can lead to diversity in the ejecta velocity profile.As the ejecta propagates, the parts of ejecta with high values of velocity gradient experience significant spread-out or dilution in radial direction and become the regions of low density.In other words, the gradient or stratification in the ejecta velocity profile can give rise to various forms of the ejecta density profile.
In view of no consensus regarding the kilonova ejecta profile, in this work, we hypothesize several velocity profiles as a function of the merger's ejection time.We introduce an ejecta model that calculates the ejecta density profile for those velocity profiles with various stratification.We further investigate how the ejecta velocity profiles leave their imprints on the kilonova lightcurves by performing a radiative-transfer simulation.The ejecta model is presented in Section 2. The method of the radiative-transfer simulation is presented in Section 3. The result and the discussion are followed in Section 4 and Section 5, respectively.

EJECTA MODEL
We adopt here the ejecta model given in Uhm (2011), which describes a radially-stratified relativistic ejecta of gamma-ray bursts (GRBs).This model makes use of a Lagrangian description and views the ejecta as a sequence of infinitesimal shells that are ejected from the merger.Each shell is prescribed an ejection time τ as its Lagrangian coordinate.After collisions between adjacent shells, every shell can stationarily coast with its velocity v(τ ) such that any two adjacent shells no longer collide with each other.1 Therefore, it is legitimate to consider a non-increasing function of v(τ ): i.e., dv/dτ ≤ 0. The relativistic continuity equation for a mass flow in spherical polar coordinates can be written as where ρ is the mass density measured in the rest frame of ejecta, c is the speed of light, t is the time in laboratory frame, r is the radius measured from the merger, u α = Γ(c, v, 0, 0), Γ = (1 − β 2 ) −1/2 , and β = v/c.Note that the radius r of a τ -shell at time t is given by where t ≥ τ , 0 ≤ τ ≤ τ d , and τ d is the duration of the ejection time.By solving the continuity equation, Uhm (2011) provides an exact solution of the ejecta density ρ for each Lagrangian shell τ at radius r, where v ′ (τ ) ≡ dv/dτ and Ṁ (τ ) represents the mass ejection rate from the merger.In the three-dimensional expansion, the ejecta cross-sectional area is proportional to r 2 , which gives the density dilution with r −2 .The term with v ′ shows the contribution arising from radial spread-out due to the velocity stratification.
Once the merger shapes its mass ejection rate Ṁ (τ ) and velocity profile v(τ ) as a function of the ejection time τ , Equation 2 gives the radius of every τ -shell and Equation 3 determines the density of all these τ -shells.Therefore, Equation 3 offers a model to fully describe the density profile of the kilonova ejecta and its evolution over time.This model is very general such that Ṁ (τ ) and v(τ ) can be arbitrary functions of τ as long as v ′ ≤ 0. When t ≫ τ d and v ′ ̸ = 0, Equation 2 and Equation 3 display the nature of homologous expansion, r = vt and ρ ∝ t −3 , respectively.It is worth mentioning that the density profile, ρ ∝ v −3 t −3 , employed in the modeling of AT2017gfo (e.g., Watson et al. 2019;Gillanders et al. 2022), can be achieved by setting β = sech(τ /τ 0 + θ 0 ), where τ 0 and θ 0 are arbitrary constants.We introduce six velocity-profile models, where the outermost shell velocity is fixed to β outer = 0.35, and inner shells have a velocity of β(τ ) = β outer + β ′ (τ )dτ , where β ′ = dβ/dτ ≤ 0. We assume τ d = 100 ms since the mass ejection from the merger can persist for several tenths of milliseconds (Kasen et al. 2017).The minimum velocity is set to 0.08 c when τ = τ d .Figure 1 shows the velocity (left) and derivative (right) profiles of our six models.Models 1 to 3 represent relatively simple cases, where the derivative is either constant or changing monotonically with time.Models 4 to 6 encompass more complex scenarios, which could reflect irregular ejecting activities of the merger.We remark that these models are not based on any prior information or numerical studies regarding the ejecta velocity profile.Instead, we conceive these models to explore various scenarios, including some complicated cases, as several hydrodynamic simulations have shown irregularities in the density profile (e.g., Kyutoku et al. 2015;Radice et al. 2018).
The ejected shells outside the optically thick photosphere constitute the so-called line-forming (lf ) region, where photons interact with the expanding ejecta material, resulting in the line-transition signature in the kilonova spectrum (for a schematic illustration, see e.g., Figure 2 in Gillanders et al. 2022).The innermost boundary of the line-forming region is fixed to a radius of r min = 10 15 cm from the merger.The velocity of the innermost boundary at a given time is obtained by v min = r min /t since t ≫ τ d for the timescale of days.This implies that, at a given time, only a portion of the velocity profile between the innermost and outermost boundaries constitutes the density profile of the line-forming region.Note that the minimum velocity reduces to about 0.08 c at 4.8 days after the merger time.For each model, the mass ejection rate, Ṁ , in Equation 3 is assumed to remain constant over the ejection time τ and is determined by setting the mass within the line-forming region to be M lf = 2 × 10 −4 M ⊙ at 2.4 days after the explosion.In other words, each model has a distinct value of Ṁ while maintaining the same mass in M lf at 2.4 days.Note that M ej = M ph + M lf , where M ej is the total ejecta mass, and M ph is the mass beneath the photosphere.
Figure 2 shows the density profile as a function of its velocity and radius at 4.8 days after the explosion, as a reference.Due to the radial stratification, various density profiles can be obtained, especially in Models 4 to 6 (see right panel of Figure 2).To explain in more detail, the relative velocity difference between the inner and outer boundaries of a shell (equivalent to β ′ ) affects the density dilution within the shell; the higher |β ′ |, the higher the dilution.For example, in Model 6, we have an extremely high |β ′ | at β ∼ 0.15 (or τ ∼ 67 ms), leading to an exceedingly low density at that velocity.

RADIATIVE-TRANSFER SPECTRAL ANALYSIS
We perform the radiative-transfer spectral analysis for the ejecta models, using the well-established code, tardis (Kerzendorf & Sim 2014;Kerzendorf et al. 2023).tardis is a widely-used Monte Carlo radiative-transfer spectral synthesis software for supernova (e.g., Dutta et al. 2021;Williamson et al. 2021) and kilonova simulations (e.g., Smartt et al. 2017;Gillanders et al. 2022).Assuming a spherically-symmetric explosion with an opthick region beneath the innermost boundary of the line-forming region, tardis computes what the resultant spectrum will be after the primary photospheric spectrum propagates through the homologouslyexpanding outer regions of ejecta.
We follow the tardis setup implemented by Gillanders et al. ( 2022 2022).In total, we include 1,029,577 transition lines from 49,117 energy levels in elements of 1≤Z≤92 (for details, see Gillanders et al. 2022).
Since we aim to explore the effect of the stratified velocity model, we refrain from varying the abundance and instead invoke the solar-system r-process composition here, adopted from Prantzos et al. (2020).
For the photospheric spectrum, we assume a steadystate photosphere: its temperature T decreases linearly in time, where t d is the time in days since the merger.This formula is approximately derived using the results of Gillanders et al. (2022), where their model temperatures of AT2017gfo decrease by ∼ 200 K per day, from 2.4 days to 4.4 days after the explosion. 3This simple assumption provides the luminosity of the photosphere as a function of time, L(t) = 4πr2 min σ R T 4 , where σ R is the Stefan-Boltzmann constant.The ejecta temperature is also fixed to the photospheric temperature to avoid complex and inaccurate effects, stemming from the uncertain temperature profile (see Gillanders et al. 2022, for details).We acknowledge that these assumptions are not entirely realistic, but they are sufficient to explore the impact of the ejecta velocity profile on kilonova evolution.
Since the kilonova ejecta is expanding with a relativistic speed during early epochs (v ≳ 0.1 − 0.2 c), we use the full treatment of special relativity in the radiative transfer, including angle aberration, as implemented in tardis by Vogl et al. (2019).Furthermore, we employ the local thermal equilibrium (LTE) approximation to handle ionization and the dilute-LTE approach for excitation.To account for line interactions, we adopt the macro-atom scheme (Lucy 2002(Lucy , 2003)), which accounts for the effects of fluorescence and multiple internal line transitions.

RESULT
We compute the kilonova spectra from 1.4 to 4.8 days after the explosion with a time step of 0.1 days for Models 1 to 6.
Figure 3 shows the evolution of the spectral energy distributions (SEDs).An SED with a yellowish color cor- responds to the earlier time step, transitioning to a reddish color for later times.Although they originate from identical conditions except for the velocity profile, they undergo different evolutionary paths.Consequently, the SEDs at 4.8 days for the six models are highly distinct.For example, in Model 5 and 6, one can find that the spectra fall rapidly at ∼ 8500 Å and form a peak at ∼ 10,000 Å.This feature is the P-Cygni line profile arising from the Sr II near-infrared triplet transitions, as seen in the early epochs of AT2017gfo (Watson et al. 2019;Gillanders et al. 2022).Depending on the models and times, this P-Cygni feature can be prominent (or hidden) because its strength relies on the density of Sr II at a specific velocity and temperature, as well as the relative line-transition rates of other elements.Further details will be studied in a future work.
Figure 4 shows the lightcurves at 16 different wavelengths.We generate them between 4000 Å and 12000 Å with the bandwidth of 500 Å, resulting in 16 curves.We find that Models 4 to 6 show unusual features in some wavelengths, while the lightcurves from Models 1 to 3 decay smoothly in all wavelengths.
To understand the origin of the distinct evolutions of Models 1-6, we firstly delve into the simple scenarios (Models 1-3) to examine how the concave or convex shape in the velocity profile affects the overall decay rate of the bolometric luminosity.In the wavelengths from 500 Å to 25000 Å, the spectral energy distributions of the three models exhibit a gradual decay over time (Figure 3).The left panel of Figure 5 shows the scaled bolometric luminosity lightcurves for Models 1 to 3. We normalize the bolometric luminosity to its maximum value as we are primarily interested in the relative changes rather than the absolute magnitude of the luminosity.Model 2 shows a relatively steep decay, compared to the others.To investigate the reason behind the difference in decaying rates, we compute the total mass of lanthanide elements (57 ≤ Z ≤ 71) in the line-forming region, where they contribute significantly to the line-transition opacity at wavelengths below ∼ 7000 Å (the right panel of Figure 5).Note that these three models have similar lanthanide mass at 2.4 days after the explosion because each model is configured to have the same mass in the line-forming region at this epoch (see Section 2).The accumulated amount of lanthanide materials of Model 2 increases sharply in time, compared to the others.This leads to an increase in the total absorption rate, resulting in a steeper decay in the lightcurve.
Compared to the simple scenarios, Models 4 to 6 exhibit diverse features in their lightcurves, attributed to the various density profiles (Figure 2). Figure 6 shows the lightcurves of Models 4 to 6 in UBgVIz filters.In the three models, a significant change in the decay rate is evident in all filters, except for the "I" filter.Specifically, Model 4 tends to show a steepening at around 2.5 days, which is particularly prominent in the filters ranging from 4300 to 5500 Å.Conversely, Model 5 demonstrates an early sharp decay followed by a more gradual decline.This break between the two phases occurs at around 2 days.Lastly, Model 6 shows a dynamic evolution in its lightcurve: an initial steep decay, followed by a period of slower decline, and then a subsequent resteepening.Especially in the "B" and "g" filter, the aforementioned evolutions in Models 4 to 6 are conspicuous.We note that the lightcurves from "u" filter (3596 Å) are similar to that of "U" filter (3633 Å).In case of "r" (6122 Å), "R" (6407 Å), and "i" (7439 Å) filters, the luminosity follows a gradual decline over time, resembling the behavior observed in the "I" filter (lowermiddle panel of Figure 6).

DISCUSSION AND CONCLUSION
In this work, we explored the impact of ejecta velocity profile on kilonova evolution.For various velocity profiles as a function of the merger's ejection time, we calculated the ejecta density profiles using the ejecta model by Uhm (2011) and then computed the resultant lightcurves with the radiative-transfer spectral synthesis software, tardis.The results indicate that the lightcurve of kilonovae can exhibit not only a change in its decaying rate but also a plateau-like behavior (slow and/or flat decline).
The plateau-like lightcurve is seemingly similar to that of Type II-P supernovae but they have different origins.The plateau in Type II-P supernovae is attributed to a change in the electron-scattering opacity as the ionization and recombination of hydrogen in the supernova's envelope progresses (Grassberg et al. 1971).In the case of kilonovae, the electron scattering contributes much less to the total opacity compared to line transitions by heavy elements (Kasen et al. 2013).In other words, the plateau-like behavior in kilonovae would be observed at specific wavelengths where the heavy elements play an important role in absorbing and re-emitting photons through line transitions.As shown in Figure 6, the filters and their corresponding wavelengths that exhibit a notable change in the decay rate are indeed connected to the energy levels at which the dominant line transitions of the heavy elements (e.g., Sr, Ce, Nd, and Sm) mainly take place.As reported in many studies on AT2017gfo (e.g., Kasen et al. 2017;Watson et al. 2019;Domoto et al. 2021;Gillanders et al. 2022), these elements play a key role in generating absorption and re-emission features in the kilonova spectrum at early epochs.
Generally, as the total mass of r-process elements increases, the lightcurve steepens due to a greater amount of absorption (see Figure 5).This simple picture can help to explain the evolution of Models 4 to 6 as well.In the case of Model 4, between 1.4 days and 2.5 days, the mass of heavy elements within the line-forming region does not show a notable increase due to a well-like dip in the density profile (see the red lines in Figure 2).From 2.5 days, the mass accumulation rate increases sharply, resulting in a break in the lightcurve.Conversely, in Model 5, the mass accumulation rate in the line-forming region increases during the early phase up to 2 days and then gradually decreases at later times (see the green lines in Figure 2), hence the lightcurve is steep up to 2 days and following which it flattens.The lightcurve of Model 6 can be understood by combining the slope changes observed in Model 4 and Model 5: i.e., transitioning gradually from the behavior of Model 5 to that of Model 4. In Model 6, the light curve up to approximately 2 days, resulting from the density profile within 0.2 ≲ β ≤ 0.35, resembles those of Model 5, but squeezed.However, as time progresses, the density profile of Model 6 exhibits a well-like dip, similar to that of Model 4, resulting in a change in the decay rate at around 3 days.
Since our work is to explore the change in kilonova emission patterns by varying the ejecta velocity profile, we do not narrow our focus to one or two ejecta components that have distinct properties such as ejecta velocity, timescale of ejection, and electron fraction Y e (e.g., Kasen et al. 2017;Siegel & Metzger 2017;Radice et al. 2018).In this work, we utilized ejecta speeds ranging from 0.08 to 0.35 c, which covers the speeds of dynamical and disk-wind components.The changes in minimum and maximum ejecta velocities can stretch or compress the lightcurve features.The duration of ejection, τ d , does not significantly impact the results because the timescale relevant to the observation, t, is sufficiently large (i.e., t ≫ τ d ).
However, we acknowledge that the location and duration of these variations in the kilonova lightcurve heavily depend on the element composition (or Y e of the ejecta).By varying the composition, we would expect to observe similar features in the lightcurve, but to have slight changes in the characteristic wavelengths and durations.
In addition, temperature can play an important role in shaping the lightcurve behavior.In this work, we adopt a continuously decreasing temperature profile (see Equation 4) to investigate the causality between the velocity profile and the resulting lightcurve features.However, the overall rate of kilonova luminosity decay is heavily contingent upon how the temperature of a heating source cools over time.This implies that a rapid cool or re-heating of the source could lead to obscure or amplify some lightcurve features caused by the irregular density profile.
Lastly, numerical-relativity simulations for the merger of binary compact objects (e.g., Kyutoku et al. 2015;Radice et al. 2018) present the mass ejection rate as a function of time, Ṁ (τ ).We may obtain an intricate density profile in this case as well (like in Models 4 to 6), leading to similar behaviors in the kilonova lightcurve.
To date, only a few kilonovae have been observed, mainly because of their short lifetime and relatively low luminosity.However, with the enhancement of followup strategies using current telescopes and also with the development of advanced future telescopes like the 7-Dimensional Telescope (7DT; Im 2021), the chances of detecting kilonovae are expected to increase significantly (e.g., Cowperthwaite et al. 2019;Chase et al. 2022;Ekanger et al. 2023).In this circumstance, this work proposes a model to explain peculiar kilonova lightcurves that might be observed in the near future.

Figure 1 .Figure 2 .
Figure1.Velocity and its derivative profile of the six ejection scenarios (Models 1 to 6) as a function of the ejection time τ .The left panel shows the velocity profile, and the right panel shows its derivative profile.We assume that the ejection continues up to 100 milliseconds with the outermost shell velocity of 0.35 c and the minimum velocity of 0.08 c.
), which successfully reproduced the AT2017gfo spectra during the photospheric epochs.Here we also use the same atomic data set, which, in addition to the default tardis atomic data set, 2 contains the extended Kurucz atomic data (Sr I-III , Y I-II , and Zr I-III ; Kurucz 2018), the Database on Rare Earths At Mons university (DREAM; 57≤Z≤71; Quinet & Palmeri 2020), and the Pt I-III and Au I-III atomic data presented by Gillanders et al. (2021) and McCann et al. (

Figure 3 .Figure 4 .Figure 5 .Figure 6 .
Figure 3. Evolution of the spectral energy distributions of Models 1 to 6.Each panel shows the spectra of each model at times ranging from 1.4 days to 4.8 days.Deeper colors on the plot represent later times in the evolution.