# Three dimensional end-to-end simulation for kilonova emission from a black-
hole neutron-star merger

Kyohei Kawaguchi,1,2,3 Nanae Domoto,4 Sho Fujibayashi,1 Kota Hayashi,1,3 Hamid
Hamidani,4 Masaru Shibata,1,3 Masaomi Tanaka,4,5 and Shinya Wanajo1
1Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
Am Mühlenberg 1, Potsdam-Golm, 14476, Germany
2Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8582, Japan
3Center for Gravitational Physics and Quantum Information, Yukawa Institute
for Theoretical Physics, Kyoto University, Kyoto, 606-8502, Japan
4Astronomical Institute, Tohoku University, Aoba, Sendai 980-8578, Japan
5Division for the Establishment of Frontier Sciences, Organization for
Advanced Studies, Tohoku University, Sendai 980-8577, Japan E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We study long-term evolution of the matter ejected in a black-hole neutron-
star (BH-NS) merger employing the results of a long-term numerical-relativity
simulation and nucleosynthesis calculation, in which both dynamical and post-
merger ejecta formation are consistently followed. In particular, we employ
the results for the merger of a $1.35\,M_{\odot}$ NS and a $5.4\,M_{\odot}$ BH
with the dimensionless spin of 0.75. We confirm the finding in the previous
studies that thermal pressure induced by radioactive heating in the ejecta
significantly modifies the morphology of the ejecta. We then compute the
kilonova (KN) light curves employing the ejecta profile obtained by the long-
term evolution. We find that our present BH-NS model results in a KN light
curve that is fainter yet more enduring than that observed in AT2017gfo. This
is due to the fact that the emission is primarily powered by the lanthanide-
rich dynamical ejecta, in which a long photon diffusion time scale is realized
by the large mass and high opacity. While the peak brightness of the KN
emission in both the optical and near-infrared bands is fainter than or
comparable to those of binary NS models, the time-scale maintaining the peak
brightness is much longer in the near-infrared band for the BH-NS KN model.
Our result indicates that a BH-NS merger with massive ejecta can
observationally be identified by the bright and long lasting ($>$two weeks)
near-infrared emission.

###### keywords:

gravitational waves – stars: neutron – nucleosynthesis – radiative transfer –
hydrodynamics

††pubyear: 20XX††pagerange: Three dimensional end-to-end simulation for
kilonova emission from a black-hole neutron-star merger–B.4

##  1 Introduction

Neutron star (NS) mergers are known to be among the most promising targets of
the ground-based gravitational-wave (GW) detectors (LIGO: Aasi et al. 2015,
Virgo: Acernese et al. 2015, KAGRA: Kuroda 2010) as well as one of the most
important sources of high-energy astrophysical transients, such as gamma-ray
bursts (GRB, Paczynski, 1991; Nakar, 2007; Berger, 2014; Abbott et al.,
2017c), kilonovae (KN, Li & Paczynski, 1998; Kulkarni, 2005; Metzger et al.,
2010; Kasen et al., 2013; Tanaka & Hotokezaka, 2013), jet heated cocoons
(Nakar & Piran, 2017; Hamidani & Ioka, 2023a, b), and synchrotron flares
(Nakar & Piran, 2011; Hotokezaka & Piran, 2015; Hotokezaka et al., 2018;
Margalit & Piran, 2020). NS mergers are also considered to be important
production sites of elements heavier than iron in the universe (Lattimer &
Schramm, 1974; Eichler et al., 1989; Freiburghaus et al., 1999; Cowan et al.,
2021). The first detection of GWs from a binary neutron star (BNS) merger
(GW170817; Abbott et al. 2017a) and its multi-wavelength electromagnetic (EM)
counterparts (Abbott et al., 2017b) demonstrated that those simultaneous
observations will provide a valuable opportunity to extend our knowledge of
fundamental physics in the extreme (strongly self-gravitating, high-density,
and high-temperature) environments.

Among NS mergers, the mergers of black-hole neutron-star (BH-NS) binaries can
provide us with interesting insights that are different from BNS mergers.
While the mass ratios of the compact stars in BNS binaries are expected to
close to unity, BH-NS binaries can be more asymmetric in the mass ratio, and
hence, will provide valuable opportunity to study higher-order GW multipole
moments (Abbott et al., 2021). Also, if the NS is tidally disrupted before
reaching the innermost circular orbit of the BH, an applicable amount of NS
matter can remain outside the remnant BH and be ejected from the system. Such
ejecta formed during the NS tidal disruption as well as the matter
subsequently ejected during the evolution of the remnant BH-tours system will
be the source of various EM counterparts to the GW event. In addition, since
BH-NS mergers can potentially produce a large amount of very low ($\lesssim
0.1$) electron fraction ($Y_{e}$) ejecta, the nucleosynthetic abundances can
be different to those in the case of BNS mergers. In fact, it has been pointed
out that BH-NS mergers can provide an explanation to the observed elemental
abundances of a subclass of $r$-process-enhanced stars, so-called "actinide-
boosted" stars (Wanajo et al., 2022).

To extract the physical information from the observation of EM counterparts,
accurate modeling of the light curves and spectra consistent with the source
properties are crucial. Since the detection of GW170817, light curve modeling
of EM counterparts, particularly, for KNe has been significantly developed in
this decade. In particular, the studies by employing numerical-simulation-
based/motivated ejecta profiles and by performing radiative transfer (RT)
simulations with realistic heating rates and/or detailed opacity tables enable
us to directly connect the properties of the progenitor binary to the
observables (e.g., Kasen et al., 2013, 2015; Barnes et al., 2016; Wollaeger et
al., 2018; Tanaka et al., 2018; Wu et al., 2019; Kawaguchi et al., 2018;
Hotokezaka & Nakar, 2020; Kawaguchi et al., 2020; Korobkin et al., 2021; Bulla
et al., 2021; Zhu et al., 2021; Barnes et al., 2021; Nativi et al., 2020;
Kawaguchi et al., 2021; Wu et al., 2022; Just et al., 2022, 2023). Previous
studies showed that the complex ejecta profile in the presence of the multiple
ejecta components of different mass ejection processes induces significant
spatial dependences in radioactive heating as well as strong geometrical
effects in radiative transfer, which have great impacts on the resulting light
curves (Kasen et al., 2015; Kawaguchi et al., 2018, 2020; Bulla, 2019; Zhu et
al., 2020; Darbha & Kasen, 2020; Korobkin et al., 2021; Almualla et al., 2021;
Kedia et al., 2023). Hence, the employment of the realistic ejecta profile
consistently taking multiple ejecta components into account is essential for
the accurate prediction of KN light curves.

One of the important missing links for the accurate prediction of KNe is the
long-term hydrodynamics evolution of ejecta after the formation. While the
ejecta formation takes place on a time scale of $\lesssim 1$–10 s after the
onset of a merger (Hayashi et al., 2022, 2023), the KN emission peaks in a
much longer time scale of 0.1–10 d (Li & Paczynski, 1998; Kulkarni, 2005;
Metzger et al., 2010; Kasen et al., 2013; Tanaka & Hotokezaka, 2013), at which
the homologous expansion of ejecta has been achieved. Since ejected matter can
be accelerated by the ejecta pressure gradient and interact with different
ejecta components during these epochs, the ejecta profile at the time of KN
emission is non-trivial just from the ejecta properties at the time of
formation. In fact, Rosswog et al. (2014) and Grossman et al. (2014) performed
pseudo-Newtonian hydrodynamics simulations for BNS mergers, and studied the
long-term evolution of the dynamical ejecta component until it reached the
homologously expanding phase. They found that the thermal pressure induced by
radioactive heating in ejecta significantly changes the ejecta morphology (see
also Foucart et al. (2021)). Fernández et al. (2015) and Fernández et al.
(2017) performed long-term simulations for BH-NS mergers to investigate the
effect of the interplay between the dynamical and post-merger components and
found that the interaction of the multiple ejecta components can modify the
ejecta profile. Thus, to accurately predict KN light curves, it is also
important to follow the hydrodynamics evolution of the multiple ejecta
components until the homologously expanding phase.

Recently, the development of numerical simulation techniques and the
significant increase in the computational resources have enabled us to
consistently follow the NS mergers from the onset of the merger to the time
that ejecta formation saturates (Kiuchi et al., 2022; Fujibayashi et al.,
2023, 2020b; Shibata et al., 2021; Hayashi et al., 2022; Kiuchi et al., 2022;
Fujibayashi et al., 2023; Hayashi et al., 2023; Kiuchi et al., 2023; Just et
al., 2023; Gottlieb et al., 2023; Kiuchi et al., 2024). In this paper, we
study the KN emission associated with a BH-NS merger employing the results
obtained by the numerical-relativity (NR) simulation and nucleosynthesis
calculation consistently following the entire ejecta formation from the merger
(Hayashi et al., 2022, 2023; Wanajo et al., 2022). In particular, we focus on
the KN emission from $\approx 1\,{\rm d}$ after the onset of the merger for
the model of a large amount of dynamical ejecta with $\approx 0.04M_{\odot}$
in this paper.111During the submission process of this paper, LIGO and Virgo
have detected GWs plausibly from a BH-NS merger (The LIGO Scientific
Collaboration et al., 2024b). The alert shows that the system is likely to
contain a NS, and the mass of the other object is likely to be in between
$3\,M_{\odot}$ and $5\,M_{\odot}$ with a 50% probability. The probability for
the matter outside the remnant object to be present after the merger is also
high ($>99\%$). Hence, the system can be by chance similar to the BH-NS model
studied in this paper (a binary of a $1.35\,M_{\odot}$ NS with the radius of
$\approx 13.2\,{\rm km}$ and a $5.4\,M_{\odot}$ BH with the dimensionless spin
of 0.75). We note that the amount of the dynamical ejecta is broadly the same
($\approx 0.04\,M_{\odot}$) also for a BH-NS merger with the same NS mass, NS
radius, and dimensionless BH spin but with a larger BH mass ($8.1\,M_{\odot}$)
(Hayashi et al., 2022). The result of this paper indicates that, if it is a
BH-NS merger with significant amount of the dynamical ejecta formation, this
GW event may be associated with a kilonova of which near-infrared emission is
bright and long-lasting.

This paper is organized as follows: In Section 2, we describe the method
employed in this study. In Section 3, we describe the BH-NS model we study in
this work. In Section 4, we present the property of the ejecta obtained by the
long-term hydrodynamics evolution. In Section 5, we present the KN light curve
obtained by RT simulations. Finally, we discuss the implication of this paper
in Section 6. Throughout this paper, $c$ denotes the speed of light.

##  2 Method

###  2.1 hydrodynamics simulation

In a BH-NS merger, matter ejected by various mechanisms is expected to
experience hydrodynamics interactions between different ejecta components
before eventually reaching a homologous expansion phase at $\sim 0.1$ d
(Kawaguchi et al., 2021). In order to obtain the spatial profile of the rest-
mass density, elemental abundances, and radioactive heating rate after 0.1 d,
which are necessary for accurate prediction of KN, we perform hydrodynamics
simulations using the outflow data obtained by NR simulations as boundary
conditions, as in our previous studies. To distinguish it from the NR
simulation, the present hydrodynamics simulation is referred to as the HD
simulation in this paper.

The simulation code for the HD simulation is a 3D extension of the code
developed in our previous studies (Kawaguchi et al., 2021, 2022, 2023). This
code solves the relativistic Euler equations under a spherical coordinate
system. In order to incorporate the effect of gravity, a fixed background
metric for a non-rotating black hole expressed in isotropic coordinates is
used. See Appendix A for the formulation of the basic equations. The effect of
radioactive heating is incorporated in the same way as in the previous studies
(Kawaguchi et al., 2021, 2022, 2023). See Appendix B for the method of
particle tracing used to employ the nucleosynthesis results in the HD
simulation. We note that the equatorial symmetry is imposed for the HD
simulation following the setup of the NR simulation.

For the equation of state (EOS), we consider both contributions from gas and
radiation: the total pressure $P$ is given by $P=P_{\rm gas}+P_{\rm rad}$ with
$P_{\rm gas}=n_{\rm B}k_{\rm B}T$ and $P_{\rm rad}=a_{\rm rad}T^{4}/3$, where
$n_{\rm B}$, $T$, $k_{\rm B}$, and $a_{\rm rad}$ are the baryon number
density, temperature, Boltzmann constant, and radiation density constant,
respectively. Here, we simplified the gas pressure assuming that atoms are
fully ionized with an electron fraction of unity, and the gas pressure is
dominated by the contribution from electrons (since the average atomic mass
number is expected to be much larger than unity). We note that, although this
simplification may overestimate the gas pressure component, the contribution
of the gas pressure is found to be nevertheless subdominant. In fact, we
confirm that the resulting ejecta profiles as well as the KN light curves are
essentially unchanged even if we employ the ideal-gas EOS with the adiabatic
index of $\Gamma=4/3$, which corresponds to the case that the radiation
pressure dominates.

Note that the magnetic field effects are not taken into account in our present
HD simulations. As a consequence, and due to the coarse grid resolution in the
polar region, the relativistic jet outflow launched in the NR simulation is
not well resolved in the present HD simulations. The previous study suggests
that the presence of the jet may affect the ejecta profile and hence the KN
light curves near the jet axis (Nativi et al., 2020; Klion et al., 2021).
Since resolving the propagation of the relativistic jet in long-term three-
dimensional simulations requires high computational costs, we leave the
investigation of the effect of the jet for a future work.

We employ the same time origin for the HD simulations as in the NR
simulations. The uniform grids with $N_{\theta}$ and $N_{\phi}$ grid points
are prepared for the polar angle $\theta$ and the longitudinal angle $\phi$,
respectively. For the radial direction, the following non-uniform grid
structure is employed; for a given $j$-th radial grid point

$\displaystyle{\rm ln}\,r_{j}={\rm ln}\left(\frac{r_{\rm out}}{r_{\rm
in}}\right)\frac{j-1}{N_{r}}+{\rm ln}\,r_{\rm in},\,j=1\cdots N_{r}+1,$ (1)

where $r_{\rm in}$ and $r_{\rm out}$ denote the inner and outer radii of the
computational domain, respectively, and $N_{r}$ denotes the total number of
the radial grid points. In the present work, we employ
$(N_{r},N_{\theta},N_{\phi})=(1024,64,128)$, and $r_{\rm in}$ and $r_{\rm
out}$ are initially set to be $3,000\,{\rm km}$ and $10^{3}\,r_{\rm in}$,
respectively. We confirm that this grid resolution is sufficiently high enough
for our purpose of the study by checking the results of the ejecta profile and
KN light curves being semi-quantitatively unchanged for the HD simulation with
$(N_{r},N_{\theta},N_{\phi})=(512,32,64)$ (less than 10% and 3% difference in
the total bolometric luminosity at 1 d and 2 d, respectively).

The hydrodynamics properties of the outflow are extracted at $r=r_{\rm ext}$
in the NR simulations of Hayashi et al. (2022, 2023), and the time-sequential
data are employed as the inner boundary condition of the present HD
simulations. The outflow data obtained from the NR simulation run out at $t>1$
s, and after then, the HD simulation is continued by setting a very small
floor value to the rest-mass density of the inner boundary. To follow the
evolution of ejecta even after the high-velocity edge of the outflow reaches
the outer boundary of our HD simulation, the radial grid points are added to
the outside of the original outer boundary, while at the same time the
innermost radial grid points are removed so as to keep the total number of the
radial grid points. By this prescription, the value of $r_{\rm in}$ is
increased in the late phase of the HD simulations. The outermost radial grids
are added so that the location of the outer radial boundary, $r_{\rm out}$, is
always $10^{3}r_{\rm in}$. Note that the region of $r\gtrsim 10^{-3}ct$ is
always covered with the computational domain up to $t=0.1\,{\rm d}$ in the HD
simulations.

The so-called Courant–Friedrichs–Lewy (CFL) condition restricts the time steps
in the HD simulation to ensure the numerical stability. For our setup, the
time interval should be approximately less than the smallest value among
$\Delta r_{\rm min}/c$, $r_{\rm in}\Delta\theta_{\rm min}/c$, and $r_{\rm
in}{\rm sin}\theta_{\rm min}\Delta\phi_{\rm min}/c$ with $\theta_{\rm min}$,
$\Delta r_{\rm min}$, $\Delta\theta_{\rm min}$, and $\Delta\phi_{\rm min}$
being the minimum cell center value of the $\theta$ coordinate and the minimum
cell sizes of $r$, $\theta$, and $\phi$ directions, respectively. For the
present grid setup, the most strict constraint comes from the last condition
of $r_{\rm in}{\rm sin}\theta_{\rm min}\Delta\phi_{\rm min}/c$, and this
restricts the time interval to be so small that the computational costs
becomes practically quite high. To relax this condition, we average over the
conservative variables of hydrodynamics in the direction of $\phi$ for all the
cells located in $\theta\leq\theta_{\rm c}$ for each sub-step of the
evolution. By this prescription, the HD simulation is kept numerically stable
if the time interval is within $r_{\rm in}{\rm sin}\theta_{\rm
c}\Delta\phi_{\rm min}/c$. For the present study, we choose $\theta_{\rm c}$
to be $\pi/24$, while we confirm that the resulting LCs are essentially
unchanged even if we employ $\theta_{\rm c}=\pi/12$.

###  2.2 radiative-transfer simulation

The light curves of KNe are calculated using a wavelength-dependent RT
simulation code (Tanaka & Hotokezaka, 2013; Tanaka et al., 2017, 2018;
Kawaguchi et al., 2020, 2021). In this code, the photon transfer is simulated
by a Monte Carlo method for given ejecta profiles composed of the density,
velocity, and elemental abundance under the assumption of the homologous
expansion. The time-dependent thermalization efficiency is taken into account
following an analytic formula derived by Barnes et al. (2016). The ionization
and excitation states are determined under the assumption of the local
thermodynamic equilibrium (LTE) by using the Saha’s ionization and Boltzmann
excitation equations.

For the photon-matter interaction, bound-bound, bound-free, and free-free
transitions, and electron scattering are taken into account for the transfer
of optical and infrared photons (Tanaka & Hotokezaka, 2013; Tanaka et al.,
2017, 2018). The formalism of the expansion opacity (Friend & Castor, 1983;
Eastman & Pinto, 1993; Kasen et al., 2006) and the new line list derived in
Domoto et al. (2022) are employed for the bound-bound transitions. In this
line list, the atomic data of VALD (Piskunov et al., 1995; Kupka et al., 1999;
Ryabchikova et al., 2015) or Kurucz’s database (Kurucz & Bell, 1995) is used
for $Z=20$–$29$, while the results of atomic calculations from Tanaka et al.
(2020) are used for $Z=30$–$88$. For Sr II, Y I, Y II, Zr I, Zr II, Ba II, La
III, and Ce III, which are the ions producing strong lines, the line data are
replaced with those calibrated with the atomic data of VALD and NIST databases
(Kramida et al., 2021). Note that, since our atomic data include only up to
the triple ionization for all the ions, the early phase of the light curves
($t\leq 0.5\,{\rm d}$) may not be very reliable due to high ejecta temperature
(see Banerjee et al. 2020 for the work taking the opacity contribution from
higher ionization states into account).

The RT simulations are performed from $t=0.1\,{\rm d}$ to $30\,{\rm d}$
employing the density and internal energy profiles of the HD simulations at
$t=0.1\,{\rm d}$ and assuming the homologous expansion for $t>0.1$ d. The
spatial distributions of the heating rate and elemental abundances are
determined by the table obtained by the nucleosynthesis calculations referring
to the injected time and angle of the fluid elements. Note that, as an
approximation, the elemental abundances at $t=1\,{\rm d}$ are used during the
entire time evolution in the RT simulations to reduce the computational cost,
but this simplified prescription gives an only minor systematic error on the
resultant light curves as illustrated in Kawaguchi et al. (2021).

A three-dimensional cylindrical grid is applied for storing the local
elemental abundances and radioactive heating rate as well as for solving the
local temperature and opacity. The 50, 50, and 32 cells are set to the
cylindrical radius, vertical, and longitudinal directions, which cover the
domain with the coordinate ranges of $(0,0.6\,ct)$, $(0,0.6\,ct)$, and
$(0,2\pi)$, respectively. We confirm that the resulting light curves are
unchanged by changing each cell numbers from 50, 50, and 32 cells to $40$,
$40$, and $28$ cells or changing the maximum cylindrical radius and vertical
coordinate ranges from $0.6\,ct$ to $0.75\,ct$.

##  3 The BH-NS model

In this work, we employ the NR outflow profiles and nucleosynthetic data
obtained in Hayashi et al. (2022, 2023) and Wanajo et al. (2022) as the input
for the HD simulations. In particular, we employ the outflow data of model
Q4B5H in Hayashi et al. (2022). For this model, a BH-NS binary of which the NS
mass, BH mass and dimensionless spin are initially $1.35\,M_{\odot}$,
$5.4\,M_{\odot}$ (thus 4 times larger than the NS mass), and 0.75,
respectively, is considered with the DD2 EOS (Banik et al., 2014). The
poloidal magnetic field with the maximum strength of $5\times 10^{16}\,{\rm
G}$ is initially set in the NS, while the resulting ejecta profile is not
sensitive to the initial magnetic-field strength and configuration (Hayashi et
al., 2023). We set $6.6\,M_{\odot}$ as the BH mass of the metric employed in
the HD simulations, which approximately agrees with the summation of the
remnant BH mass and matter outside the BH measured at $t=0.1\,{\rm s}$.

For model Q4B5H, the NS experiences significant tidal disruption before it
reaches the inner-most stable circular orbit of the binary ($t\approx 10\,{\rm
ms}$). This leads to the formation of massive ejecta and torus around the
remnant BH. Ejecta formed at the time of the NS tidal disruption, which often
referred to as the dynamical ejecta, are concentrated in the vicinity of
equatorial plane and exhibit significant non-axisymmetric geometry. The
dynamical ejecta typically have low electron fraction (0.03–0.07) because
those are driven primarily by gravitational torque and do not go through
significant weak processes in the merger.

Subsequently, the magnetic field is amplified in the remnant torus, and the
effective viscosity is induced by the magnetohydrodynamical turbulence, driven
by the magnetorotational instability (Balbus & Hawley, 1998). Initially,
viscous heating in the torus is balanced with neutrino cooling. As the disk
rest-mass density and temperature drop due to the expansion driven by angular
momentum transport, neutrino cooling becomes inefficient, and viscosity-driven
mass ejection sets in ($t\approx 0.2$–$0.3\,{\rm s}$). In parallel, magneto-
centrifugal force in the central region might play a role for enhancing mass
ejection. Mass ejection in this stage, which is referred to as the post-merger
mass ejection, lasts for $\sim 1$–$10\,{\rm s}$. In contrast to the dynamical
ejecta, since thermal and weak processes play important roles during the post-
merger stages, the electron fraction of ejecta has a broad distribution in the
range of 0.1–0.4 with its peak being 0.24.

For model Q4B5H, the masses of the dynamical and post-merger ejecta are found
to be $0.045\,M_{\odot}$ and $0.028\,M_{\odot}$, respectively, at the end time
of the NR simulation. It is worth being remarked that the combination of
dynamical and post-merger ejecta approximately reproduces a solar-like
$r$-process pattern (Wanajo et al., 2022).

In this paper, we study the ejecta and KN property for one case of a BH-NS
merger among available NR results as the first step for the end-to-end
kilonova simulation. However, we should note that the disk and ejecta masses
formed in BH-NS mergers can have large variety depending on the binary
parameters, such as the BH and NS masses, BH spin, and NS radius (Rosswog,
2005; Shibata & Taniguchi, 2008; Etienne et al., 2009; Lovelace et al., 2013;
Kyutoku et al., 2015; Foucart et al., 2018), as well as the adopted EOS
(Hayashi et al., 2023). For example, the smaller amount of disk and ejecta
would be formed for the case that the NS radius is smaller due to softer EOS,
such as the SFHo EOS (Steiner et al., 2013; Hayashi et al., 2023). Hence, the
resulting property of the KN light curves can also have a large diversity.
Therefore, we emphasize that the ejecta and KN property found for model Q4B5H
with the DD2 EOS may not be universal property for every case of BH-NS
mergers, and we leave the investigation of the binary parameter and EOS
dependences for a future work.

##  4 Results: Hydrodynamics simulation

###  4.1 Ejecta mass evolution

Figure 1: Time evolution of the total rest mass in the computational domain of
the HD simulation (the purple curve). The green curve denotes the same as for
the purple curve but only for the matter which satisfies the geodesic
criterion ($u_{t}<-1$ where $u_{t}$ is the lower time component of the four
velocity). The blue dashed curve with the label “input" denotes the rest mass
obtained by integrating the mass flux of the NR outflow data which is employed
as the inner boundary condition of the HD simulation. The black dashed line
denotes the time at which the NR outflow data run out.

Figure 1 shows the total rest mass in the computational domain as a function
of time. We can consider that the ejecta has reached the homologously
expanding phase at $t=0.1\,{\rm d}$, because the total internal energy of
ejecta is smaller by 4 order of magnitudes than the total kinetic energy. In
general, two distinct ejecta components are seen in Figure 1. One found in
$t_{\rm in}\sim\,0.1$ s corresponds to the dynamical ejecta, and the other
found in $t_{\rm in}\gtrsim\,0.5$ s corresponds to the post-merger ejecta.

After the NR outflow data run out at $t\approx 1\,{\rm s}$, we impose a floor
rest-mass density value to the inner boundary. It is clearly seen in Figure 1
that the total mass in the computational domain decreases after that time,
indicating that the matter is artificially falling back and escaping through
the inner boundary. This happens because the pressure support from the inner
boundary vanishes after the outflow data run out. However, the total mass of
the matter with gravitationally unbound orbits remains increasing even after
the time when the NR outflow data run out, as the consequence of the
acceleration of the matter in the presence of the thermal pressure gradient.

After $t\approx 100\,{\rm s}$, approximately all the ejecta matter remaining
in the computational domain becomes gravitationally unbound, and the value of
the total mass in the computational domain converges to $0.063\,M_{\odot}$.
This value is smaller than the ejecta mass estimated in Hayashi et al. (2023)
by $\approx 0.01\,M_{\odot}$. We interpret this discrepancy as a consequence
of the mismatch in the employed EOS between the NR and HD simulations and the
inconsistency of the matter flux at the inner boundary. In fact, Fernández et
al. (2015) show similar results: they performed BH-disk simulations to follow
the formation of the post-merger wind ejecta and used the extracted ejecta
property as the inner boundary condition of the subsequent simulation for
long-term ejecta evolution in the same manner as our present work. They found
that the mass of the post-merger wind ejecta which becomes gravitationally
unbound in the subsequent simulation decreases by a factor of $\approx 2$ from
the values estimated in the BH-disk simulations. They interpreted this
difference as a consequence of the discrepancy between the stresses at the
inner boundary and those that would be obtained in a self-consistent
simulation.

Nevertheless, by performing the HD simulation with artificially modified inner
boundary conditions, we confirmed that our main results are essentially the
same and the modification to the resulting KN light curves is only minor: we
perform a HD simulation in which the ejecta injection is sustained with the
final value of the mass flux at $1\,{\rm s}$ after the NR outflow data run
out. By this prescription, the total ejecta mass in the HD simulation at the
homologously expanding phase increases by $0.01\,M_{\odot}$, but the
bolometric luminosity increases only at most $\approx 10\%$ since the unbound
matter increased by this prescription has the velocity only less than
$\lesssim 0.05\,c$ and hence has a long diffusion time scale, which gives a
minor contribution to the brightness of the emission.

###  4.2 Ejecta profiles at homologously expanding phase

Figure 2: Rest-mass density and electron fraction ($Y_{e}$) profiles on the
equatorial plane at $t=0.1\,{\rm d}$. The yellow dotted lines denote the
angles for which the meridional ejecta profiles are shown in Figures 3 and 4.
The white dotted curves denote the longitudinal angle ranges in which the KN
light curves shown in Figures 8 and 9 are obtained. The value employed as the
initial condition of the nucleosynthesis calculation is shown in the $Y_{e}$
profile (see Appendix B and Wanajo et al. (2022) for the detail.)

Figure 3: Rest-mass density profiles on the meridional planes at $t=0.1\,{\rm
d}$. The top left, top right, bottom left, and bottom right panels denote the
profiles on the $\phi\approx 66^{\circ}$, $156^{\circ}$, $246^{\circ}$,
$336^{\circ}$ planes, respectively (see also the left panel of Figure 2 for
the location of each plane). $R$ denotes the cylindrical radius.

Figure 4: Same as Figure 3 but for electron fraction $Y_{e}$ (see also the
left panel of Figure 2 for the location of each plane).

Figures 2, 3, and 4 show the rest-mass density and electron fraction ($Y_{e}$)
profiles of ejecta with two-dimensional various cross sections at $t=0.1\,{\rm
d}$ obtained by the HD simulation. Here, the value employed as the initial
condition of the nucleosynthesis calculation is shown in the $Y_{e}$ profile
(see Appendix B and Wanajo et al. (2022) for the detail). The center of mass
for the matter with $Y_{e}<0.1$ is located in the direction of $\phi\approx
141^{\circ}$ with $\phi$ being the longitudinal angle measured from the $+x$
axis. The longitudinal angles of the meridional planes shown in Figures 3 and
4 are selected to show the profiles in which the dynamical ejecta are
approximately mostly (‘b)’: $\phi\approx 156^{\circ}$), moderately (‘a)’:
$\phi\approx 66^{\circ}$ and ’c)’: $\phi\approx 246^{\circ}$), and least
(‘d)’: $\phi\approx 336^{\circ}$) present. As we mentioned above, the entire
ejecta have reached the homologously expanding phase at this epoch. Broadly
speaking, the dynamical and post-merger ejecta are present around the regions
where the cylindrical radius is larger and smaller than $\approx
0.05$–$0.1ct$, respectively. Those two components are clearly distinguishable
with the value of $Y_{e}$. The value of $Y_{e}$ for the dynamical ejecta is
typically below $0.1$, which primarily reflects the original $Y_{e}$ values of
the disrupted NS. On the other hand, the post-merger ejecta have wider range
of $Y_{e}$ values from $0.1$ to $0.4$.

The rest-mass density profile of the dynamical ejecta exhibits clear non-
axisymmetric geometry, with its mass mostly distributed in the fan-like shape
in $70^{\circ}\lesssim\phi\lesssim 250^{\circ}$. The dynamical ejecta are
extended up to $\approx 0.5ct$ in the cylindrical radius direction, while
their vertical extent is $\approx 0.2ct$. The aspect ratio of the cylindrical
and vertical extents for the dynamical ejecta is close to unity. This is in
contrast to the fact that the dynamical ejecta are launched initially confined
around the equatorial plane within the latitudinal opening angle of $\sim
10^{\circ}$ (Kyutoku et al., 2013; Foucart et al., 2014). As we show below,
this ejecta expansion is due to thermal pressure enhanced by radioactive
heating.

On the other hand, the post-merger ejecta exhibit approximately an
axisymmetric shape. It has two distinct components with one having
approximately a spherical shape and the other having the cone-like shape. The
former is concentrated in the region within $\approx 0.05ct$ while the latter
is more extended in the vertical direction with the polar opening angle of
$\approx 10^{\circ}$ and the vertical extent reaches $\approx 0.25ct$. As we
show below, this complex geometry of the post-merger ejecta is realized by the
interaction with the dynamical ejecta which significantly expand due to
thermal pressure enhanced by radioactive heating.

Figure 5: Rest-mass density and electron fraction ($Y_{e}$) profiles on the
equatorial plane at $t=0.1\,{\rm d}$ for the HD simulation in which
radioactive heating is turned off. Figure 6: Mass weighted average of the
total specific radioactive heating rate in our HD simulation. The specific
heating rate of model H4 in Darbha et al. (2021) is also shown.

Figure 5 shows the rest-mass density and electron fraction profiles of ejecta
on the equatorial and meridional planes at $t=0.1\,{\rm d}$ obtained by the HD
simulation but switching off radioactive heating. Under the presence of
radioactive heating, the dynamical ejecta expand significantly due to the
increase in thermal pressure and the inhomogeneities in the rest-mass density
are also smoothed out, as clearly seen in Figures 2 and 3. These results are
consistent with the finding of Rosswog et al. (2014); Grossman et al. (2014)
in the context of BNSs, and of Fernández et al. (2015); Darbha et al. (2021)
in the context of BH-NSs. In fact, the resulting aspect ratio of the dynamical
ejecta is found to be close to unity as the model H4 in Darbha et al. (2021).
The radioactive heating rate of the dynamical ejecta in our model also agrees
with that in Darbha et al. (2021) (see Figure 6).

The comparison between Figure 4 and Figure 5 shows that the profile of the
post-merger ejecta is affected by the modification of the dynamical eject
profile. Figure 5 shows that, in the absence of radioactive heating, the post-
merger ejecta exhibit a prolate shape with the extension of $0.1ct$ and
$0.25ct$ in the equatorial and vertical directions, respectively. On the other
hand, the radioactive heating significantly expands the dynamical ejecta,
which compress the post-merger ejecta in $0.05ct\leq z\leq 0.15ct$ and confine
the ejecta in the region of $\lesssim 0.05ct$ as found in Figure 4. This
happens because of the higher typical electron fraction of the post-merger
ejecta. The higher electron fraction leads to the relatively small radioactive
heating rate and hence small enhancement of the pressure of the ejecta
compared to the dynamical component.

Significant expansion of the dynamical ejecta and enforced confinement of the
post-merger ejecta in the presence of radioactive heating are not found in the
BNS models in our previous studies (Kawaguchi et al., 2021, 2022, 2023). This
is because the dynamical ejecta of the present BH-NS model and the BNS models
studied in Rosswog et al. (2014); Grossman et al. (2014) are massive compared
to the post-merger ejecta and also much more confined around the equatorial
plane compared to the BNS models in our previous studies. As a result, higher
internal energy density and high thermal pressure are realized. Hence, the
importance of the radioactive heating will depend on the density and isotopic-
abundance profiles of ejecta, which can have a variety even among BH-NS
mergers depending on the binary parameters or the adopted EOS.

##  5 Results: KN light curves

###  5.1 bolometric light curves

Figure 7: Total bolometric luminosity and total energy deposition rate for
model Q4B5H. The isotropically equivalent bolometric luminosity observed in
AT2017gfo with the distance of 40 Mpc is shown by the filled circles adopting
the data in Waxman et al. (2018), which assume a black-body fit to the photo-
metric observations. Note that the bolometric light curve before 1 d is hidden
since it is not reliable due to the lack of opacity data in the high
temperature regime ($\gtrsim 20,000\,{\rm K}$).

Figure 7 shows the bolometric luminosity calculated by the RT simulation
employing the ejecta rest-mass density, elemental abundance, and radioactive
heating rate profiles obtained by the combination of the results of the HD
simulation and nucleosynthesis calculation (Wanajo et al., 2022). The total
energy deposition rate taking the thermalization efficiency into account is
also plotted in Figure 7. As we mentioned in Section 2, our atomic data
include only up to the triple ionization for all the ions, and the opacity of
ejecta in the early phase ($t\leq 0.5\,{\rm d}$) may be underestimated due to
high temperature ($\gtrsim 20,000\,{\rm K}$). Hence, hereafter, we only focus
on the light curves after 1 d.

For 1–10 d, the bolometric luminosity is approximately constant with the value
of 1–2$\times 10^{\rm 41}\,{\rm erg/s}$. It decreases only slowly and the
change is only by a factor of 2 during this epoch. However, after 10 d, the
bolometric luminosity starts decreasing more rapidly, and it decreases by a
factor of 5 during 10–30 d. This faint and long-lasting emission is caused by
the fact that it is primarily powered by the lanthanide-rich dynamical ejecta,
in which a long photon diffusion time scale is realized by the large mass and
high opacity. This behaviour of the bolometric luminosity is qualitatively the
same as that found in the models with massive dynamical ejecta studied in the
previous study (MS1Q3a75 and H4Q3a75 in Tanaka et al. (2014)). The bolometric
luminosity converges to the total deposition rate after 20 d, which suggests
that the entire thermal photons created in the ejecta immediately diffuse out
from the ejecta after this epoch. As we show below, however, the viewing-angle
dependence of the emission due to the aspherical profile of the ejecta opacity
is still playing a role up to 30 d.

Our present BH-NS KN model shows significantly distinct light curves from the
observation of AT2017gfo. Specifically, our BH-NS model is fainter by a factor
of 2 around $\sim 1\,{\rm d}$ than AT2017gfo. However, due to the slow
decrease in the bolometric luminosity, compared to AT2017gfo, our BH-NS KN
model becomes comparably bright at 4 d, and brighter by a factor of 5 at 10 d.
This result clearly shows that a BH-NS binary which we study particularly in
this work is not likely to be the progenitor of AT2017gfo.

Figure 8: Isotropically equivalent bolometric luminosities observed from
various viewing angles for model Q4B5H. The top panels denote the comparison
among the results for different longitudinal directions, while the middle and
bottom panels denote the comparison among the results for different
latitudinal directions. The isotropically equivalent bolometric luminosity
observed in AT2017gfo with the distance of 40 Mpc is also shown by the filled
circles adopting the data in Waxman et al. (2018), which assume a black-body
fit to the photo-metric observations. Note that the light curves before 1 d
are hidden since they are not reliable due to the lack of opacity data in the
high temperature regime ($\gtrsim 20,000\,{\rm K}$).

Figure 8 shows the results of the isotropically equivalent bolometric light
curves observed from various viewing angles for the present model. Focusing on
the observer in the polar direction with $\theta\leq 28^{\circ}$, the KN
emission is brightest in b): $135^{\circ}\leq\phi<180^{\circ}$. This direction
approximately matches to the longitudinal direction in which the dynamical
ejecta have the most of their mass (see Figures 2 and 3). On the other hand,
the faintest emission is observed from the direction in which the dynamical
ejecta least present (‘d)’: $315^{\circ}\leq\phi<360^{\circ}$). Nevertheless,
the variation in the bolometric luminosity is not so large and it is always
within 40%. This is reasonable because the observers with different
longitudinal angles are in similar directions for the polar view $\theta\leq
28^{\circ}$.

On the other hand, the longitudinal variation is larger for the equatorial
view ($82^{\circ}\leq\theta\leq 90^{\circ}$). For this case, the KN emission
is also brightest in b): $135^{\circ}\leq\phi<180^{\circ}$ in which direction
the dynamical ejecta is most present. By contrast, the bolometric luminosity
is the faintest in a): $45^{\circ}\leq\phi<90^{\circ}$. This is because, in
this direction, the relatively thin part of the dynamical ejecta present in
$R\gtrsim 0.3ct$ around the equatorial plane (see Figure 3) suppresses
radiation from the ejecta center (note that such suppression is not
significant from the direction of c): $225^{\circ}\leq\phi<270^{\circ}$ due to
the absence of the low density ejecta above $R\approx 0.3ct$). Meanwhile, the
emission from the post-merger ejecta enhances the luminosity in this view
($82^{\circ}\leq\theta\leq 90^{\circ}$) in which the dynamical ejecta are
least present (‘d)’: $315^{\circ}\leq\phi<360^{\circ}$). The variation in the
bolometric luminosity is larger for the equatorial view than that for the
polar view, and it is larger than a factor of 2 for 1–7 d.

The latitudinal viewing-angle dependence of the bolometric luminosity is not
significant and the variation in the bolometric luminosity is always less than
a factor of 2 in our present model except for the equatorial view in a):
$45^{\circ}\leq\phi<90^{\circ}$. The dependence of the KN brightness on the
latitudinal direction is weak in particular from the viewing angle of b):
$135^{\circ}\leq\phi<180^{\circ}$, in which direction the KN emission is
brightest. As we show below, the latitudinal viewing-angle dependence is much
weaker than that for BNS KNe. This is due to the fact that, for the present
BH-NS KN model, the emission is dominated by the dynamical ejecta of which the
aspect ratio is close to unity. In fact, compared with a previous study
(Tanaka et al., 2014), our present BH-NS KN model shows a less significant
viewing-angle dependence on the latitudinal direction. This is because the
ejecta in this study have larger aspect ratio than those in the models of
their previous study. This comparison indicates that the modification of the
ejecta morphology due to radioactive heating has a great impact on the
viewing-angle dependence of the KN emission. Hence, this work demonstrates the
importance of modeling KN light curves taking the ejecta long-term evolution
into account.

Interestingly, the latitudinal viewing-angle dependence is still present in
d): $315^{\circ}\leq\phi<360^{\circ}$ even after 20 d at which the total
bolometric luminosity converges to the total deposition rate (see Figure 7).
This can be understood by the fact that the post-merger ejecta are present in
between the high-density part of the dynamical ejecta and the observer in the
direction of d): $315^{\circ}\leq\phi<360^{\circ}$ (see Figures 2 and 3).
While the post-merger ejecta have a minor contribution to the luminosity in
the late epoch due to their relatively low heating rate, they still contribute
as an opacity source (because of the relatively high density) to prevent
photons emitted in the high-density part of the dynamical ejecta diffusing out
to the direction of d): $315^{\circ}\leq\phi<360^{\circ}$. Although it is not
quantitatively significant, this long-lasting viewing-angle dependence due to
the non-axisymmetric geometry of the ejecta might have an impact on estimating
the total deposition rate in the ejecta from the late-time observation.

###  5.2 broad-band magnitudes

Figure 9: gzK-band light curves for model Q4B5H observed from various viewing
angles with the distance of 40 Mpc. The top panels denote the comparison among
the results for different longitudinal directions, while the middle and bottom
panels denote the comparison among the results for different latitudinal
directions. The data points denote the AB magnitudes of AT2017gfo taken from
Villar et al. (2017).

Figure 9 shows the optical (the g and z bands) and near-infrared (the K-band)
light curves observed from various viewing angles. As is the case for the
bolometric luminosity, for the polar view ($\theta\leq 28^{\circ}$), the gzK-
bands are the brightest and the faintest in b):
$135^{\circ}\leq\phi<180^{\circ}$ and in d): $315^{\circ}\leq\phi\leq
360^{\circ}$, in which the dynamical ejecta are mostly and least present,
respectively. The viewing-angle dependence of the emission is weak from the
polar view, and the variation is always within 0.5 mag around the peak
magnitudes.

For the equatorial view ($\theta\geq 28^{\circ}$), the gzK-band emission is
the brightest in b): $135^{\circ}\leq\phi<180^{\circ}$, while the emission in
a): $45^{\circ}\leq\phi\leq 90^{\circ}$ becomes the faintest. The longitudinal
viewing-angle dependence of the emission is significant in the gz-bands, and
it is always larger than 1 mag among different longitudinal directions. The
variation in the K-band magnitude among different longitudinal directions is
relatively small compared to that in the gz-bands, and it is always
approximately within 1 mag among the all viewing angles.

The g-band emission is fainter and declines faster than AT2017gfo even if it
is observed from the brightest direction (‘b)’:
$135^{\circ}\leq\phi<180^{\circ}$). The peak brightness in the z-band is
comparable to the observation of AT2017gfo, but it declines much faster. In
contrast, the K-band emission is comparable to that of AT2017gfo in a few
days, and then, it becomes brighter after 4 d. The K-band magnitude finally
reaches its peak at $\approx 10\,{\rm d}$ after the onset of the merger with
its emission brighter than AT2017gfo by more than 1 mag.

Interestingly, the gz-band emission observed from b):
$135^{\circ}\leq\phi<180^{\circ}$ becomes slightly brighter in the equatorial
direction than in the polar direction after 1 d for the present BH-NS model.
This brightness dependence on the latitudinal direction is opposite compared
to the BNS KN models, for which the emission becomes brighter in the polar
direction222However, it should be noted that here we do not consider the
impact that the short GRB jet might have on the polar ejecta and on the KN
emission (see Hamidani et al. (2024)).. The same latitudinal-angle dependence
is also found for the emission observed in the direction of the ejecta bulk
motion for the model in Darbha et al. (2021) in which the radioactive heating
rate agrees with our model (model H4). The brighter emission in the equatorial
plane is explained by the enhancement of the radiation energy flux due to the
Doppler effect induced by the bulk motion of the ejecta. It is also important
to have the aspect ratio of the dynamical ejecta close to unity to realize the
present latitudinal-angle dependence of the emission brightness: otherwise the
Doppler effect can be obscured by the suppression of the emission due to the
decrease in the projected area toward the observer for the case that the
dynamical ejecta have a more oblate shape (see Darbha et al. (2021) for the
discussion).

###  5.3 radiative-transfer effect of non-axisymmetric ejecta geometry

Figure 10: Comparison of the isotropically equivalent bolometric luminosities
(left) and gzK-band light curves (right) between the fiducial (the same as in
Figures 8 and 9) and axisymmetrized (labeled as “2D model”) models. The top
and bottom panels denote the light curves observed from
$0^{\circ}\leq\theta\leq 28^{\circ}$ and $82^{\circ}\leq\theta\leq
90^{\circ}$, respectively. The solid and dashed curves denote the light curves
of the axisymmetrized model and those for the fiducial model observed from b):
$135^{\circ}\leq\phi<180^{\circ}$, respectively. The dotted curves in the
upper and bottom panels denote the light curves of the fiducial model observed
from d): $315^{\circ}\leq\phi<360^{\circ}$ and a):
$45^{\circ}\leq\phi<90^{\circ}$, respectively.

To clarify the RT effect of the non-axisymmetric ejecta geometry, we perform a
RT simulation for an axisymmetrized ejecta profile. The axisymmetrized ejecta
profile is generated by averaging over the rest-mass density, specific
internal energy, elemental abundances, and radioactive heating rate profiles
obtained by the HD simulation at $t=0.1\,{\rm d}$ with respect to the
longitudinal direction. Note that the volume and mass in each grid cell are
used as the weights of the average for the rest-mass density and the latter
three quantities, respectively. Figure 10 compares the isotropically
equivalent bolometric luminosities and gzK-band light curves observed from the
polar ($0^{\circ}\leq\theta\leq 28^{\circ}$) and equatorial
($82^{\circ}\leq\theta\leq 90^{\circ}$) directions between the axisymmetrized
and fiducial models (the same as in Figures 8 and 9). For the fiducial model,
the light curves observed from b): $135^{\circ}\leq\phi<180^{\circ}$, d):
$315^{\circ}\leq\phi<360^{\circ}$ (for the polar view), and a):
$45^{\circ}\leq\phi<90^{\circ}$ (for the equatorial view) are shown.

The upper panels of Figure 10 show that the KN emission observed from the
polar direction becomes slightly brighter for the axisymmetrized model than
for the fiducial model except for the g-band emission. This reflects the fact
that the area of the dynamical ejecta projected toward the observer increases
for the axisymmetrized model due to the longitudinal average. The bolometric
light curves declines slightly earlier than the original fiducial model
because the optical depth decreases due to the decrease in the rest-mass
density of the dynamical ejecta for the axisymmetrized model. Nevertheless,
the effect of the longitudinal average is found to be minor for the polar
view, particularly, in the gzK-band magnitudes, for which the differences
between the axisymmetrized model and fiducial model are always smaller than
$0.5$ mag for $t\geq 0.5\,{\rm d}$.

The difference in the brightness between the axisymmetrized model and fiducial
model is more pronounced for the emission observed from the equatorial
direction. For the equatorial view, the bolometric luminosity observed from
the longitudinal direction of the brightest region (‘b)’:
$135^{\circ}\leq\phi<180^{\circ}$) is brighter approximately by a factor of
1.5–2 for the fiducial model than that for the axisymmetrized model. The gzK-
band magnitudes observed from the same direction are also brighter than those
for the axisymmetrized model by $\sim 1\,{\rm mag}$. On the other hand, the KN
brightness observed from the longitudinal direction of the faintest region
(‘a)’: $45^{\circ}\leq\phi<90^{\circ}$) is comparable to or slightly fainter
for the fiducial model than that for the axisymmetrized model. This
discrepancy in the equatorial brightness between the fiducial and
axisymmetrized models is due to the difference in the ejecta aspect ratio: as
a consequence of the longitudinal average, the polar projected area of the
ejecta is larger than that in the equatorial direction for the axisymmetrized
model, which makes photons to preferentially diffuse in the polar direction
and thus the equatorial brightness to be fainter.

###  5.4 Comparison with various BNS KN models

Figure 11: Comparison of the bolometric and gzK-band light curves among the
present BH-NS KN model and various BNS KN models. For BNS KN models, three
cases are shown: a case in which the remnant MNS survives for a short time
(the dashed curves; SFHo-125145, Kiuchi et al. (2022); Fujibayashi et al.
(2023); Kawaguchi et al. (2023)), a case in which the remnant MNS survives for
a long time (the dash-dotted curves; DD2-135135, Fujibayashi et al. (2020b);
Kawaguchi et al. (2022)), and a case in which large-scale magnetic field
significantly plays a role in the long-surviving remnant MNS (the dotted
curves; MNS75a, Shibata et al. (2021); Kawaguchi et al. (2022)). Note that the
light curve for model SFHo-125145 in the top-right panel is below the plot
range.

Figure 11 compares the gzK-band light curves among the present BH-NS KN model
and various BNS KN models obtained in our previous studies (Kawaguchi et al.,
2021, 2022, 2023). For BNS KN models, three cases are shown as representative:
a case in which the remnant massive NS (MNS) survives for a short time (the
dashed curves; SFHo-125145, Kiuchi et al. (2022); Fujibayashi et al. (2023);
Kawaguchi et al. (2023)), a case in which the remnant MNS survives for a long
time (the dash-dot curves; DD2-135135, Fujibayashi et al. (2020b); Kawaguchi
et al. (2022)), and a case in which large-scale magnetic field significantly
plays a role in the long-surviving remnant MNS (the dotted curves; MNS75a,
Shibata et al. (2021); Kawaguchi et al. (2022)). We note that the BNS KN
models are obtained by imposing axisymmetry in all the post-merger NR
simulations, subsequent HD simulations, and RT simulations. For the BH-NS
model, we show the light curves observed from b):
$135^{\circ}\leq\phi<180^{\circ}$, d): $315^{\circ}\leq\phi<360^{\circ}$ (for
the polar view), and a): $45^{\circ}\leq\phi<90^{\circ}$ (for the equatorial
view), which represent the longitudinal directions for the brightest and
faintest emission, respectively.

The gz-band emission observed from the polar direction
($0^{\circ}\leq\theta\leq 20^{\circ}$) for the present BH-NS KN model is by
0.5–1 mag brighter than that for the BNS models in which the remnant MNS
survives only for a short time ($<10\,{\rm ms}$, SFHo-125145), but is by
$\approx$ 1 mag fainter than that for the BNS models in which the remnant MNSs
survive for a long time ($>1\,{\rm s}$, DD2-135135 and MNS75a). On the other
hand, the gz-band emission observed from the equatorial direction
($86^{\circ}\leq\theta\leq 90^{\circ}$) is comparably bright or brighter than
those for the BNS models in which remnant MNSs survive for a long time except
for the z-band emission of model MNS75a. This is due to the fact that the BNS
KN models show stronger latitudinal viewing-angle dependence than the BH-NS KN
model and become significantly faint in the equatorial view. The difference in
the latitudinal viewing-angle dependence reflects the fact that the dynamical
ejecta are the primary source of the emission in the optical wavelength for
the BH-NS model, while for the BNS models, the post-merger ejecta are main
source of the emission and the dynamical ejecta are mostly acting as the
opacity source rather than the emission source (lanthanide curtain effect;
Kasen et al. 2015; Kawaguchi et al. 2018, 2020; Bulla 2019; Zhu et al. 2020;
Darbha & Kasen 2020; Korobkin et al. 2021). We note that, in the BNS cases,
enhancement of the brightness due to the Doppler effect is obscured by the
latitudinal-angle dependence of the emission induced by the angle-dependent
opacity of the dynamical ejecta.

The K-band emission for the present BH-NS KN model has comparable peak
brightness to that for the BNS models without significant large-scale magnetic
field effect in the remnant NS (SFHo-125145 and DD2-135135). However, it is
only the BH-NS model that maintains the K-band brightness within 1 mag of its
peak for a two-week period. The BNS models in which the remnant MNSs survive
for short and long periods of time become fainter than the BH-NS model after
$1$–$2\,{\rm d}$ and $5$–$7\,{\rm d}$, respectively. The BNS model in which
large-scale magnetic field significantly plays a role in the remnant NS shows
bright K-band emission for a week, but the brightness declines much faster
than that for the BH-NS model. Hence, the observation of a KN with long-
lasting near-infrared emission which is bright for more than two weeks will
indicate that the progenitor of a KN is a BH-NS merger with massive ejecta (in
particular dynamical ejecta) formation.

##  6 Summary and Discussions

In this paper, we studied the long-term evolution of the matter ejected in a
BH-NS merger by employing the results of the NR simulation and nucleosynthesis
calculation, in which both dynamical and post-merger ejecta are followed
consistently. In particular, we employed the results for the merger of a
$1.35\,M_{\odot}$ NS and $5.4\,M_{\odot}$ BH with the dimensionless spin of
0.75. We confirmed the finding in the previous studies that, thermal pressure
induced by radioactive heating in the ejecta could significantly modify the
morphology of the ejecta. In our studied case of a BH-NS binary, the dynamical
ejecta expand significantly and the aspect ratio becomes close to unity with
the fine structure being smeared out in the presence of radioactive heating.
On the other hand, the post-merger ejecta were compressed and confined in the
region with the radial velocity $\lesssim 0.05\,c$ due to the significant
expansion of the dynamical component.

We then computed the KN light curves employing the ejecta profile obtained by
the HD simulation of the ejecta matter. We found that our present BH-NS model
results in KN light curves that are fainter but longer lasting than those
observed in AT2017gfo, reflecting the fact that the emission is primarily
powered by the lanthanide-rich massive dynamical ejecta. The optical-band
emission is comparable to or fainter than those for the various BNS models
obtained in our previous studies. While the peak brightness of the near-
infrared emission is also comparable to the BNS models, the time-scale
maintaining the brightness is much longer, and the emission comparable to the
peak brightness within 1 mag is sustained for more than two weeks for the BH-
NS model. The wide-field infrared observations with the ground-based
telescopes, such as VISTA (Ackley et al., 2020), WINTER (Frostig et al., 2022)
and PRIME (Kondo et al., 2023), can detect such bright infrared KN emission up
to $\approx 14$ d if the distance to the event is within 150 Mpc since the
K-band emission will be apparently brighter than 21 mag for all the viewing
angles. However, the field of views of infrared telescopes are typically not
as large as those for the optical telescopes for a given sensitivity (Nissanke
et al., 2013). Therefore, a tight constraint of the localization area by the
GW data analysis or the follow-up observation within $\approx 1$ d in the
optical bands is crucial to detect the KN emission unless the event occurs as
close as in the case of AT2017gfo. Once a KN with long-lasting near-infrared
emission is found, follow-up observations in the radio band may also be useful
to support the presence of massive dynamical ejecta, by finding the
synchrotron radio flares with the relatively delayed peak time of $\sim
10\,{\rm yr}$ (Kyutoku et al., 2013).

We found that the non-axisymmetric geometry of the ejecta induces various
interesting radiative-transfer effects in the viewing-angle dependence of the
KN emission. In particular, we found the Doppler effect induced by the bulk
velocity of the ejecta to the emission, which is pointed out by Fernández et
al. (2017) and Darbha et al. (2021), is in fact present. Due to this effect,
the optical light curves observed from the direction of the bulk ejecta motion
show a slightly inverted latitudinal angle dependence to those found in the
BNS models: The optical-band emission observed from b):
$135^{\circ}\leq\phi<180^{\circ}$ becomes slightly brighter in the equatorial
direction than in the polar direction for the present BH-NS model. Since the
KNe emission becomes fainter in the equatorial direction than in the polar
direction for BNS mergers, our results suggest that, for the edge-on view, the
KN emission for BH-NS mergers can be brighter in the optical-band than that of
BNS mergers.

Our results indicate that the long-lasting near-infrared emission is the key
to distinguish the types of progenitors by the KN observation. If the K-band
emission of which brightness comparable to its peak is maintained for more
than two weeks, it may indicate that the progenitor is a BH-NS merger with
massive ejecta formation. This is consistent with our finding in the previous
study (Kawaguchi et al., 2020). On the other hand, only from the optical
emission, the BH-NS KN light curves can be similar to those associated with
BNS mergers, and hence, it may be difficult to infer the information of the
progenitor. We should note that the ejecta mass and hence the brightness of
the KN of BH-NS mergers can have large variety depending on the binary
parameters, such as the BH and NS masses, BH spin, and NS radius (Rosswog,
2005; Shibata & Taniguchi, 2008; Etienne et al., 2009; Lovelace et al., 2013;
Kyutoku et al., 2015; Foucart et al., 2018), as well as on the adopted EOS
(Hayashi et al., 2023). We also note that the assumption of LTE employed in
our RT simulation will not be valid in the region where the rest-mass density
has been significantly dropped. It is implied in Hotokezaka et al. (2021);
Pognan et al. (2022) that the matter temperature can be higher than that
estimated under the assumption of LTE if non-LTE effects take place. In such a
case, the infrared emission can be dimmer with the combination of the
suppression of the neutral and poorly ionized atoms (Kawaguchi et al., 2021,
2022, 2023). Hence, we should note that the absence of the long-lasting bright
near-infrared emission does not necessarily rule out the possibility that the
progenitor of the observed KN is a BH-NS merger.

So far, four candidates have been reported for BH-NS GW events: GW190814
(Abbott et al., 2020), GW200105/GW200115 (Abbott et al., 2021), and GW230529
(The LIGO Scientific Collaboration et al., 2024a). Among them, according to
the inferred masses and spins of the binary, the latest GW event, GW230529,
was most likely to be accompanied by EM counterparts. Unfortunately, the EM
counterpart was not found in GW230529 due to the poorly constrained sky
localization, although the luminosity distance to the event was relatively
close ($201^{+102}_{-96}\,{\rm Mpc}$ with the error bar being the 90% credible
intervals). Nevertheless, the discovery of this system increases the expected
rate of the GW detection of a BH-NS merger with EM counterparts in the future.
For this event, (under the assumption that this event was a BH-NS merger) the
ratio of the BH mass to the NS mass and the dimensionless BH spin were less
than $\approx 4$ and larger than $\approx 0$, respectively. For such a case of
a BH-NS merger, the ratio of the post-merger ejecta mass to the dynamical one
can be larger compared to the BH-NS model which we studied in this paper
(Hayashi et al., 2021). This indicates that the resulting KN may become bluer
than the present result, while it is not always trivial since the long-term
hydrodynamics evolution of ejecta may also differ. Hence, the systematic study
on the KNe for various configurations of BH-NS binaries would be crucial to
quantitatively interpret the EM observational data in the future.

Figure 12: Comparison between the present BH-NS KN model and GRB KN
candidates. The solid and dashed curves denote the polar light curves in the
observer frame ($0^{\circ}\leq\theta\leq 28^{\circ}$ ) for the present BH-NS
KN model observed from b): $135^{\circ}\leq\phi<180^{\circ}$ and d):
$315^{\circ}\leq\phi<360^{\circ}$, respectively. The square and triangle
symbols denote, respectively, the observed magnitudes and upper-limits of the
optical and near-infrared counterparts of GRBs taken from Jin et al. (2016,
2015); Yang et al. (2015); Berger et al. (2013); Tanvir et al. (2013); Lamb et
al. (2019); Troja et al. (2019); Rastinejad et al. (2022); Troja et al.
(2022); Gompertz et al. (2023); Levan et al. (2024). The afterglow models
which broadly reproduce the models in the literature are also plotted in the
dotted curves.

There are a number of KN candidates reported which are associated with the
observation of GRBs: GRB050709 (Jin et al., 2016), GRB060614 (Jin et al.,
2015; Yang et al., 2015), GRB130603B (Berger et al., 2013; Tanvir et al.,
2013), GRB160821B (Lamb et al., 2019; Troja et al., 2019), GRB211211A
(Rastinejad et al., 2022; Troja et al., 2022; Gompertz et al., 2023), and
GRB230307A (Levan et al., 2024). In Figure 12, we compare our present BH-NS KN
model with these observational data. The optical and near-infrared brightness
of KN candidates found in GRB211211A and GRB230311 are comparable to that of
AT2017gfo. We find that our present BH-NS KN model is too faint to explain the
optical brightness of these KN candidates at a few days, while the K-band
emission of our present BH-NS KN model after 4 d is too bright to be
consistent with the later time upper limits. Our present BH-NS KN model is
also too faint in the optical bands to explain the KN candidates found in
GRB050709, GRB060614, and GRB160821B after $\approx 2$ d. The K-band
brightness of GRB160821B at 4.3 d is comparable to that of our present BH-NS
KN model, while the BNS model in which the remnant MNS survives for a long
time ($>1\,{\rm s}$, DD2-135135 in Figure 11) also has comparable K-band
brightness at that epoch. Interestingly, despite the bright and long-lasting
K-band emission, our present BH-NS model has fainter H-band emission at $\sim
10\,{\rm d}$ than that observed in GRB130603B. This is due to the fact that
the KN of our present BH-NS model is very red and the peak of the red-shifted
spectrum is located in the wavelength longer than the H-band at that epoch. In
summary, currently we do not find a KN candidate that can be only explained by
our present BH-NS model with significant dynamical ejecta formation. However,
as we mentioned above, we cannot rule out the possibility that some of those
KN candidates are KNe associated with BH-NS mergers, since BH-NS KNe can have
a large diversity reflecting the variety of binary parameters as well as the
adopted EOS.

We found more than a factor of 2 variation in the KN brightness depending on
the viewing angle in our present BH-NS model. Such a variation in the
brightness can induce the same degree of the systematic error in conducting
the ejecta parameter estimation for the ejecta properties, such as the mass,
velocity, and effective ejecta opacity, from the observational data. We should
also note that we only focus on one single case of a BH-NS merger with the DD2
EOS, and it is not clear whether the property of KN is always the same for
other setups of BH-NS mergers. For example, if the longitudinal opening angle
of the ejecta is close to $2\pi$, the BH-NS KN can have the viewing-angle
dependence in the brightness comparably strong to those of BNS mergers as we
indeed see in the results of the axisymmetrized model (see Fig. 10). We
further should note that uncomprehended systematic errors in the opacity and
heating rate can induce large systematic errors in the ejecta parameter
inference. In particular, the latter can be severe for KNe from BH-NS mergers
since the uncertainty is more significant for the ejecta with low values of
$Y_{e}$ ($<0.24$, see Barnes et al. (2021); Zhu et al. (2021)). Hence, it is
essential to consider that these systematic errors can significantly alter the
results of the ejecta parameter inference, and those estimated values should
be used with a great caution.

## Acknowledgements

Numerical computation was performed on Yukawa21 at Yukawa Institute for
Theoretical Physics, Kyoto University and the Yamazaki, Sakura, Cobra, and
Raven clusters at Max Planck Computing and Data Facility, and the Cray XC50 at
CfCA of the National Astronomical Observatory of Japan. ND acknowledges
support from Graduate Program on Physics for the Universe (GP-PU) at Tohoku
University. This work was supported by Grant-in-Aid for Scientific Research of
JSPS/MEXT (20H00158, 21H04997, 22KJ0317, 23H00127, 23H04894, 23H04900,
23H05432, and 23H01172) and JST FOREST Program (JPMJFR212Y).

## Data availability

Data and results underlying this article will be shared on reasonable request
to the corresponding author.

## References

  * Aasi et al. (2015) Aasi J., et al., 2015, Class. Quant. Grav., 32, 074001
  * Abbott et al. (2017a) Abbott B., et al., 2017a, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2017b) Abbott B. P., et al., 2017b, Astrophys. J., 848, L12
  * Abbott et al. (2017c) Abbott B. P., et al., 2017c, Astrophys. J. Lett., 848, L13
  * Abbott et al. (2020) Abbott R., et al., 2020, Astrophys. J. Lett., 896, L44
  * Abbott et al. (2021) Abbott R., et al., 2021, Astrophys. J. Lett., 915, L5
  * Acernese et al. (2015) Acernese F., et al., 2015, Class. Quant. Grav., 32, 024001
  * Ackley et al. (2020) Ackley K., et al., 2020, Astron. Astrophys., 643, A113
  * Almualla et al. (2021) Almualla M., Ning Y., Salehi P., Bulla M., Dietrich T., Coughlin M. W., Guessoum N., 2021, arXiv e-prints, p. arXiv:2112.15470
  * Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70, 1
  * Banerjee et al. (2020) Banerjee S., Tanaka M., Kawaguchi K., Kato D., Gaigalas G., 2020, Astrophys. J., 901, 29
  * Banik et al. (2014) Banik S., Hempel M., Bandyopadhyay D., 2014, Astrophys. J. Suppl., 214, 22
  * Barnes et al. (2016) Barnes J., Kasen D., Wu M.-R., Martínez-Pinedo G., 2016, Astrophys. J., 829, 110
  * Barnes et al. (2021) Barnes J., Zhu Y. L., Lund K. A., Sprouse T. M., Vassh N., McLaughlin G. C., Mumpower M. R., Surman R., 2021, Astrophys. J., 918, 44
  * Berger (2014) Berger E., 2014, Ann. Rev. Astron. Astrophys., 52, 43
  * Berger et al. (2013) Berger E., Fong W., Chornock R., 2013, Astrophys. J. Lett., 774, L23
  * Bulla (2019) Bulla M., 2019, Mon. Not. Roy. Astron. Soc., 489, 5037
  * Bulla et al. (2021) Bulla M., et al., 2021, Mon. Not. Roy. Astron. Soc., 501, 1891
  * Cowan et al. (2021) Cowan J. J., Sneden C., Lawler J. E., Aprahamian A., Wiescher M., Langanke K., Martínez-Pinedo G., Thielemann F.-K., 2021, Rev. Mod. Phys., 93, 15002
  * Darbha & Kasen (2020) Darbha S., Kasen D., 2020, ApJ, 897, 150
  * Darbha et al. (2021) Darbha S., Kasen D., Foucart F., Price D. J., 2021, Astrophys. J.
  * Domoto et al. (2022) Domoto N., Tanaka M., Kato D., Kawaguchi K., Hotokezaka K., Wanajo S., 2022, Astrophys. J., 939, 8
  * Eastman & Pinto (1993) Eastman R. G., Pinto P. A., 1993, ApJ, 412, 731
  * Eichler et al. (1989) Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126
  * Etienne et al. (2009) Etienne Z. B., Liu Y. T., Shapiro S. L., Baumgarte T. W., 2009, Phys. Rev., D79, 044024
  * Fernández et al. (2015) Fernández R., Quataert E., Schwab J., Kasen D., Rosswog S., 2015, Mon. Not. Roy. Astron. Soc., 449, 390
  * Fernández et al. (2017) Fernández R., Foucart F., Kasen D., Lippuner J., Desai D., Roberts L. F., 2017, Class. Quant. Grav., 34, 154001
  * Foucart et al. (2014) Foucart F., et al., 2014, Phys. Rev., D90, 024026
  * Foucart et al. (2018) Foucart F., Hinderer T., Nissanke S., 2018, Phys. Rev., D98, 081501
  * Foucart et al. (2021) Foucart F., Moesta P., Ramirez T., Wright A. J., Darbha S., Kasen D., 2021, Phys. Rev. D, 104, 123010
  * Freiburghaus et al. (1999) Freiburghaus C., Rosswog S., Thielemann F. K., 1999, ApJ, 525, L121
  * Friend & Castor (1983) Friend D. B., Castor J. I., 1983, ApJ, 272, 259
  * Frostig et al. (2022) Frostig D., et al., 2022, Astrophys. J., 926, 152
  * Fujibayashi et al. (2020a) Fujibayashi S., Shibata M., Wanajo S., Kiuchi K., Kyutoku K., Sekiguchi Y., 2020a, Phys. Rev. D, 101, 083029
  * Fujibayashi et al. (2020b) Fujibayashi S., Wanajo S., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2020b, Astrophys. J., 901, 122
  * Fujibayashi et al. (2023) Fujibayashi S., Kiuchi K., Wanajo S., Kyutoku K., Sekiguchi Y., Shibata M., 2023, Astrophys. J., 942, 39
  * Gompertz et al. (2023) Gompertz B. P., et al., 2023, Nature Astron., 7, 67
  * Gottlieb et al. (2023) Gottlieb O., et al., 2023, Astrophys. J. Lett., 953, L11
  * Grossman et al. (2014) Grossman D., Korobkin O., Rosswog S., Piran T., 2014, Mon. Not. Roy. Astron. Soc., 439, 757
  * Hamidani & Ioka (2023a) Hamidani H., Ioka K., 2023a, Mon. Not. Roy. Astron. Soc., 520, 1111
  * Hamidani & Ioka (2023b) Hamidani H., Ioka K., 2023b, Mon. Not. Roy. Astron. Soc., 524, 4841
  * Hamidani et al. (2024) Hamidani H., Kimura S. S., Tanaka M., Ioka K., 2024, Astrophys. J., 963, 137
  * Hayashi et al. (2021) Hayashi K., Kawaguchi K., Kiuchi K., Kyutoku K., Shibata M., 2021, Phys. Rev. D, 103, 043007
  * Hayashi et al. (2022) Hayashi K., Fujibayashi S., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2022, Phys. Rev. D, 106, 023008
  * Hayashi et al. (2023) Hayashi K., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2023, Phys. Rev. D, 107, 123001
  * Hotokezaka & Nakar (2020) Hotokezaka K., Nakar E., 2020, ApJ, 891, 152
  * Hotokezaka & Piran (2015) Hotokezaka K., Piran T., 2015, Mon. Not. Roy. Astron. Soc., 450, 1430
  * Hotokezaka et al. (2018) Hotokezaka K., Kiuchi K., Shibata M., Nakar E., Piran T., 2018, Astrophys. J., 867, 95
  * Hotokezaka et al. (2021) Hotokezaka K., Tanaka M., Kato D., Gaigalas G., 2021, MNRAS, 506, 5863
  * Jin et al. (2015) Jin Z.-P., Li X., Cano Z., Covino S., Fan Y.-Z., Wei D.-M., 2015, Astrophys. J. Lett., 811, L22
  * Jin et al. (2016) Jin Z.-P., et al., 2016, Nature Commun., 7, 12898
  * Just et al. (2022) Just O., Kullmann I., Goriely S., Bauswein A., Janka H. T., Collins C. E., 2022, MNRAS, 510, 2820
  * Just et al. (2023) Just O., et al., 2023, Astrophys. J. Lett., 951, L12
  * Kasen et al. (2006) Kasen D., Thomas R. C., Nugent P., 2006, Astrophys. J., 651, 366
  * Kasen et al. (2013) Kasen D., Badnell N. R., Barnes J., 2013, Astrophys. J., 774, 25
  * Kasen et al. (2015) Kasen D., Fernandez R., Metzger B., 2015, Mon. Not. Roy. Astron. Soc., 450, 1777
  * Kawaguchi et al. (2018) Kawaguchi K., Shibata M., Tanaka M., 2018, Astrophys. J., 865, L21
  * Kawaguchi et al. (2020) Kawaguchi K., Shibata M., Tanaka M., 2020, ApJ, 889, 171
  * Kawaguchi et al. (2021) Kawaguchi K., Fujibayashi S., Shibata M., Tanaka M., Wanajo S., 2021, Astrophys. J., 913, 100
  * Kawaguchi et al. (2022) Kawaguchi K., Fujibayashi S., Hotokezaka K., Shibata M., Wanajo S., 2022, Astrophys. J., 933, 22
  * Kawaguchi et al. (2023) Kawaguchi K., Fujibayashi S., Domoto N., Kiuchi K., Shibata M., Wanajo S., 2023, Mon. Not. Roy. Astron. Soc.
  * Kedia et al. (2023) Kedia A., et al., 2023, Phys. Rev. Res., 5, 013168
  * Kiuchi et al. (2022) Kiuchi K., Held L. E., Sekiguchi Y., Shibata M., 2022, Phys. Rev. D, 106, 124041
  * Kiuchi et al. (2023) Kiuchi K., Fujibayashi S., Hayashi K., Kyutoku K., Sekiguchi Y., Shibata M., 2023, Phys. Rev. Lett., 131, 011401
  * Kiuchi et al. (2024) Kiuchi K., Reboul-Salze A., Shibata M., Sekiguchi Y., 2024, Nature Astron., 8, 298
  * Klion et al. (2021) Klion H., Duffell P. C., Kasen D., Quataert E., 2021, Mon. Not. Roy. Astron. Soc., 502, 865
  * Kondo et al. (2023) Kondo I., Sumi T., Koshimoto N., Rattenbury N. J., Suzuki D., Bennett D. P., 2023, AJ, 165, 254
  * Korobkin et al. (2021) Korobkin O., et al., 2021, ApJ, 910, 116
  * Kramida et al. (2021) Kramida A., Yu. Ralchenko Reader J., and NIST ASD Team 2021, NIST Atomic Spectra Database (ver. 5.9). Available: https://physics.nist.gov/asd. National Institute of Standards and Technology, Gaithersburg, MD.
  * Kulkarni (2005) Kulkarni S. R., 2005, arXiv e-prints, pp astro–ph/0510256
  * Kupka et al. (1999) Kupka F., Piskunov N., Ryabchikova T. A., Stempels H. C., Weiss W. W., 1999, A&AS, 138, 119
  * Kurganov & Tadmor (2000) Kurganov A., Tadmor E., 2000, Journal of Computational Physics, 160, 241
  * Kuroda (2010) Kuroda K., 2010, Class. Quant. Grav., 27, 084004
  * Kurucz & Bell (1995) Kurucz R. L., Bell B., 1995, Atomic line list
  * Kyutoku et al. (2013) Kyutoku K., Ioka K., Shibata M., 2013, Phys. Rev. D, 88, 041503
  * Kyutoku et al. (2015) Kyutoku K., Ioka K., Okawa H., Shibata M., Taniguchi K., 2015, Phys. Rev., D92, 044028
  * Lamb et al. (2019) Lamb G. P., et al., 2019, ApJ, 883, 48
  * Lattimer & Schramm (1974) Lattimer J. M., Schramm D. N., 1974, Astrophys. J., 192, L145
  * Levan et al. (2024) Levan A. J., et al., 2024, Nature, 626, 737
  * Li & Paczynski (1998) Li L.-X., Paczynski B., 1998, Astrophys. J., 507, L59
  * Lovelace et al. (2013) Lovelace G., Duez M. D., Foucart F., Kidder L. E., Pfeiffer H. P., Scheel M. A., Szilegyi B., 2013, Class. Quant. Grav., 30, 135004
  * Margalit & Piran (2020) Margalit B., Piran T., 2020, MNRAS, 495, 4981
  * Metzger et al. (2010) Metzger B. D., et al., 2010, Mon. Not. Roy. Astron. Soc., 406, 2650
  * Nakar (2007) Nakar E., 2007, Phys. Rept., 442, 166
  * Nakar & Piran (2011) Nakar E., Piran T., 2011, Nature, 478, 82
  * Nakar & Piran (2017) Nakar E., Piran T., 2017, Astrophys. J., 834, 28
  * Nativi et al. (2020) Nativi L., Bulla M., Rosswog S., Lundman C., Kowal G., Gizzi D., Lamb G. P., Perego A., 2020, Mon. Not. Roy. Astron. Soc., 500, 1772
  * Nishimura et al. (2015) Nishimura N., Takiwaki T., Thielemann F.-K., 2015, ApJ, 810, 109
  * Nissanke et al. (2013) Nissanke S., Kasliwal M., Georgieva A., 2013, Astrophys. J., 767, 124
  * Paczynski (1991) Paczynski B., 1991, Acta Astron., 41, 257
  * Piskunov et al. (1995) Piskunov N. E., Kupka F., Ryabchikova T. A., Weiss W. W., Jeffery C. S., 1995, A&AS, 112, 525
  * Pognan et al. (2022) Pognan Q., Jerkstrand A., Grumer J., 2022, MNRAS, 510, 3806
  * Rastinejad et al. (2022) Rastinejad J. C., et al., 2022, Nature, 612, 223
  * Rosswog (2005) Rosswog S., 2005, Astrophys. J., 634, 1202
  * Rosswog et al. (2014) Rosswog S., Korobkin O., Arcones A., Thielemann F., Piran T., 2014, Mon. Not. Roy. Astron. Soc., 439, 744
  * Ryabchikova et al. (2015) Ryabchikova T., Piskunov N., Kurucz R. L., Stempels H. C., Heiter U., Pakhomov Y., Barklem P. S., 2015, Phys. Scr., 90, 054005
  * Shibata (2015) Shibata M., 2015, Numerical Relativity. WORLD SCIENTIFIC, doi:10.1142/9692, https://doi.org/10.1142/9692
  * Shibata & Taniguchi (2008) Shibata M., Taniguchi K., 2008, Phys. Rev., D77, 084015
  * Shibata et al. (2021) Shibata M., Fujibayashi S., Sekiguchi Y., 2021, Phys. Rev. D, 104, 063026
  * Steiner et al. (2013) Steiner A. W., Hempel M., Fischer T., 2013, Astrophys. J., 774, 17
  * Tanaka & Hotokezaka (2013) Tanaka M., Hotokezaka K., 2013, Astrophys. J., 775, 113
  * Tanaka et al. (2014) Tanaka M., Hotokezaka K., Kyutoku K., Wanajo S., Kiuchi K., Sekiguchi Y., Shibata M., 2014, Astrophys. J., 780, 31
  * Tanaka et al. (2017) Tanaka M., et al., 2017, Publ. Astron. Soc. Jap., 69, Publications of the Astronomical Society of Japan, Volume 69, Issue 6, 1 December 2017, 102, https://doi.org/10.1093/pasj/psx121
  * Tanaka et al. (2018) Tanaka M., et al., 2018, Astrophys. J., 852, 109
  * Tanaka et al. (2020) Tanaka M., Kato D., Gaigalas G., Kawaguchi K., 2020, Mon. Not. Roy. Astron. Soc., 496, 1369
  * Tanvir et al. (2013) Tanvir N. R., Levan A. J., Fruchter A. S., Hjorth J., Wiersema K., Tunnicliffe R., de Ugarte Postigo A., 2013, Nature, 500, 547
  * The LIGO Scientific Collaboration et al. (2024a) The LIGO Scientific Collaboration the Virgo Collaboration the KAGRA Collaboration 2024a, arXiv e-prints, p. arXiv:2404.04248
  * The LIGO Scientific Collaboration et al. (2024b) The LIGO Scientific Collaboration the Virgo Collaboration the KAGRA Collaboration 2024b, GCN Circular, 36240
  * Troja et al. (2019) Troja E., et al., 2019, Mon. Not. Roy. Astron. Soc., 489, 2104
  * Troja et al. (2022) Troja E., et al., 2022, Nature, 612, 228
  * Villar et al. (2017) Villar V. A., et al., 2017, Astrophys. J., 851, L21
  * Wanajo et al. (2022) Wanajo S., Fujibayashi S., Hayashi K., Kiuchi K., Sekiguchi Y., Shibata M., 2022, arXiv e-prints, p. arXiv:2212.04507
  * Waxman et al. (2018) Waxman E., Ofek E. O., Kushnir D., Gal-Yam A., 2018, Mon. Not. Roy. Astron. Soc., 481, 3423
  * Wollaeger et al. (2018) Wollaeger R. T., et al., 2018, Mon. Not. Roy. Astron. Soc., 478, 3298
  * Wu et al. (2019) Wu M.-R., Barnes J., Martinez-Pinedo G., Metzger B. D., 2019, Phys. Rev. Lett., 122, 062701
  * Wu et al. (2022) Wu Z., Ricigliano G., Kashyap R., Perego A., Radice D., 2022, MNRAS, 512, 328
  * Yang et al. (2015) Yang B., et al., 2015, Nature Commun., 6, 7323
  * Zhu et al. (2020) Zhu J.-P., Yang Y.-P., Liu L.-D., Huang Y., Zhang B., Li Z., Yu Y.-W., Gao H., 2020, Astrophys. J., 897, 20
  * Zhu et al. (2021) Zhu Y. L., Lund K., Barnes J., Sprouse T. M., Vassh N., McLaughlin G. C., Mumpower M. R., Surman R., 2021, Astrophys. J., 906, 94

##  Appendix A Formulation

In this appendix, we describe the formulation of axisymmetric hydrodynamics
equations in the spherical coordinates employed for the long-term evolution of
ejecta. Throughout this appendix, the units of $c=1=G$ are employed where $G$
is the gravitational constant, unless otherwise mentioned.

###  A.1 Basic equations

The basic equations for the numerical hydrodynamics employed in this work are
formulated in the framework of the 3+1 decomposition of the spacetime (see,
e.g., Shibata, 2015). In the 3+1 form, the metric tensor $g_{\mu\nu}$ is
decomposed as

$\displaystyle
ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=-\alpha^{2}dt^{2}+\gamma_{ij}\left(dx^{i}+\beta^{i}dt\right)\left(dx^{j}+\beta^{j}dt\right),$
(2)

where $\mu$ and $\nu$ denote the spacetime indices, $i$ and $j$ denote the
spatial indices, $\alpha$, $\beta^{i}$, and $\gamma_{ij}$ denote the lapse,
shift, and spatial metric, respectively. We treat the matter as a perfect
fluid and the energy-momentum tensor is given by

$\displaystyle T_{\mu\nu}=\rho hu_{\mu}u_{\mu}+Pg_{\mu\nu},$ (3)

where $\rho$, $h$, $u^{\mu}$, and $P$ denote the rest-mass density, specific
enthalpy, four velocity, and pressure, respectively. The equations of the
Euler equation, energy equation, and continuity equation are given,
respectively, by

$\displaystyle\gamma_{\nu i}\nabla_{\mu}T^{\mu\nu}$
$\displaystyle=-\rho{\dot{\epsilon}}_{\rm esc}u_{i}$ (4) $\displaystyle
n_{\nu}\nabla_{\mu}T^{\mu\nu}$ $\displaystyle=-\rho{\dot{\epsilon}}_{\rm
esc}n_{\nu}u^{\nu}$ (5) $\displaystyle\nabla_{\mu}\left(\rho u^{\mu}\right)$
$\displaystyle=0,$ (6)

with the covariant derivative, $\nabla_{\mu}$. Here,
$n_{\nu}=-\alpha\nabla_{\nu}t$, $\gamma_{\mu\nu}=g_{\mu\nu}+n_{\mu}n_{\nu}$,
and ${\dot{\epsilon}}_{\rm esc}$ is the specific radioactive heating rate
deposited in the form of neutrinos which entirely escape from the system
without being thermalized. Note that the contribution of other radioactive
decay channels that instantaneously thermalize (beta, alpha, and spontaneous
fission) is taken into account by the same treatment as described in Appendix
A of Kawaguchi et al. (2022).

Equations (4), (5) and (6) are rewritten in the forms

$\displaystyle\partial_{t}S_{i}+\partial_{k}\left(S_{i}v^{k}+P\alpha\sqrt{\gamma}\delta^{k}_{i}\right)=$
$\displaystyle-S_{0}\partial_{i}\alpha+S_{k}\partial_{i}\beta^{k}$
$\displaystyle-\frac{1}{2}\alpha\sqrt{\gamma}S_{jk}\partial_{i}\gamma^{jk}-\frac{\alpha{\dot{\epsilon}_{\rm
esc}}}{hw}S_{i},$ (7) $\displaystyle\partial_{t}S_{0}$
$\displaystyle+\partial_{k}\left[S_{0}v^{k}+P\sqrt{\gamma}\left(v^{k}+\beta^{k}\right)\right]=$
$\displaystyle-\gamma^{ij}S_{i}\partial_{j}\alpha+\alpha\sqrt{\gamma}S_{ij}K^{ij}-\alpha\rho_{*}{\dot{\epsilon}_{\rm
esc}},$ (8)
$\displaystyle\partial_{t}\rho_{*}+\partial_{k}\left(\rho_{*}v^{k}\right)=0,$
(9)

respectively. Here, $K_{ij}$ denotes the extrinsic curvature, and the other
variables which newly appear in the above equations are defined as follows:

$\displaystyle\sqrt{\gamma}$ $\displaystyle={\rm
det}\left(\gamma_{ij}\right),$ $\displaystyle\rho_{*}$ $\displaystyle=\rho
w\sqrt{\gamma},$ $\displaystyle w$ $\displaystyle=\alpha u^{t},$
$\displaystyle S_{i}$ $\displaystyle=\rho_{*}{\hat{u}}_{i}=\rho_{*}hu_{i},$
$\displaystyle S_{0}$
$\displaystyle=\rho_{*}{\hat{e}}=\rho_{*}\left(hw-\frac{P}{\rho w}\right),$
$\displaystyle S_{ij}$ $\displaystyle=\rho hu_{i}u_{j}+P\gamma_{ij},$
$\displaystyle v^{i}$ $\displaystyle=\frac{u^{i}}{u^{t}}.$ (10)

In the non-rotating black-hole spacetime with the isotropic coordinates, where
${\dot{\alpha}}=0$, $\beta^{i}=0$, and $K_{ij}=0$, the basic equations are
simplified as

$\displaystyle\partial_{t}\left(r^{2}{\rm sin}\theta{\tilde{\rho}}_{*}\right)$
$\displaystyle+\partial_{r}\left(r^{2}{\rm
sin}\theta{\tilde{\rho}}_{*}v^{(r)}\right)$
$\displaystyle+\partial_{\theta}\left(r{\rm
sin}\theta{\tilde{\rho}}_{*}v^{(\theta)}\right)$
$\displaystyle+\partial_{\phi}\left(r{\tilde{\rho}}_{*}v^{(\phi)}\right)=0,$
(11) $\displaystyle\partial_{t}\left(r^{2}{\rm
sin}\theta{\tilde{S}}_{(r)}\right)$ $\displaystyle+\partial_{r}\left[r^{2}{\rm
sin}\theta\left({\tilde{S}}_{(r)}v^{(r)}+P\alpha\sqrt{{\tilde{\gamma}}}\right)\right]$
$\displaystyle+\partial_{\theta}\left(r{\rm
sin}\theta{\tilde{S}}_{(r)}v^{(\theta)}\right)$
$\displaystyle+\partial_{\phi}\left(r{\tilde{S}}_{(r)}v^{(\phi)}\right)$
$\displaystyle=r^{2}{\rm
sin}\theta\left[-{\tilde{S}}_{0}\partial_{r}\alpha-\frac{1}{2}\alpha\sqrt{{\tilde{\gamma}}}{\tilde{S}}_{(i)(j)}\partial_{r}{\tilde{\gamma}}^{(i)(j)}\right.$
$\displaystyle+\frac{2}{r}\alpha\sqrt{{\tilde{\gamma}}}P+\frac{1}{r}\left({\tilde{S}}_{(\theta)}v^{(\theta)}+{\tilde{S}}_{(\phi)}v^{(\phi)}\right)$
$\displaystyle\left.-\frac{\alpha}{hw}{\tilde{S}}_{(r)}{\dot{\epsilon}_{\rm
esc}}\right],$ (12)

$\displaystyle\partial_{t}\left(r^{3}{\rm
sin}\theta{\tilde{S}}_{(\theta)}\right)$
$\displaystyle+\partial_{r}\left(r^{3}{\rm
sin}\theta{\tilde{S}}_{(\theta)}v^{(r)}\right)$
$\displaystyle+\partial_{\theta}\left[r^{2}{\rm
sin}\theta\left({\tilde{S}}_{(\theta)}v^{(\theta)}+P\alpha\sqrt{{\tilde{\gamma}}}\right)\right]$
$\displaystyle+\partial_{\phi}\left(r^{2}{\tilde{S}}_{(\theta)}v^{(\phi)}\right)$
$\displaystyle=r^{3}{\rm
sin}\theta\left[\frac{1}{r}\alpha\sqrt{{\tilde{\gamma}}}P{\rm
cot}\theta+\frac{1}{r}{\tilde{S}}_{(\phi)}v^{(\phi)}{\rm cot}\theta\right.$
$\displaystyle\left.-\frac{\alpha}{hw}{\tilde{S}}_{(\theta)}{\dot{\epsilon}_{\rm
esc}}\right],$ (13) $\displaystyle\partial_{t}\left(r^{3}{\rm
sin}^{2}\theta{\tilde{S}}_{(\phi)}\right)$
$\displaystyle+\partial_{r}\left(r^{3}{\rm
sin}^{2}\theta{\tilde{S}}_{(\phi)}v^{(r)}\right)$
$\displaystyle+\partial_{\theta}\left(r^{2}{\rm
sin}^{2}\theta{\tilde{S}}_{(\phi)}v^{(\theta)}\right)$
$\displaystyle+\partial_{\phi}\left[r^{2}{\rm
sin}\theta\left({\tilde{S}}_{(\phi)}v^{(\phi)}+P\alpha\sqrt{{\tilde{\gamma}}}\right)\right]$
$\displaystyle=-r^{3}{\rm
sin}^{2}\theta\frac{\alpha}{hw}{\tilde{S}}_{(\phi)}{\dot{\epsilon}_{\rm
esc}},$ (14)

$\displaystyle\partial_{t}\left(\alpha r^{2}{\rm
sin}\theta{\tilde{S}}_{0}\right)$ $\displaystyle+\partial_{r}\left[\alpha
r^{2}{\rm
sin}\theta\left({\tilde{S}}_{0}v^{(r)}+P\sqrt{{\tilde{\gamma}}}v^{(r)}\right)\right]$
$\displaystyle+\partial_{\theta}\left[\alpha r{\rm
sin}\theta\left({\tilde{S}}_{0}v^{(\theta)}+P\sqrt{{\tilde{\gamma}}}v^{(\theta)}\right)\right]$
$\displaystyle+\partial_{\phi}\left[\alpha
r\left({\tilde{S}}_{0}v^{(\phi)}+P\sqrt{{\tilde{\gamma}}}v^{(\phi)}\right)\right]$
$\displaystyle=-\alpha^{2}r^{2}{\rm
sin}\theta{\tilde{\rho}}_{*}{\dot{\epsilon}_{\rm esc}}.$ (15)

Here, the definitions for the new variables are given as follows:

$\displaystyle\Lambda_{(r)}^{r}=1,~{}\Lambda_{(\theta)}^{\theta}=\frac{1}{r},$
$\displaystyle~{}\Lambda_{(\phi)}^{\phi}=\frac{1}{r{\rm
sin}\theta},~{}\Lambda_{(i)}^{j}=0\,(i\neq j),$
$\displaystyle{\tilde{\gamma}}_{(i)(j)}$
$\displaystyle=\Lambda_{(i)}^{k}\Lambda_{(j)}^{l}\gamma_{kl}=\psi^{4}\delta_{(i)(j)},$
$\displaystyle\sqrt{{\tilde{\gamma}}}$ $\displaystyle=\frac{1}{r^{2}{\rm
sin}\theta}\sqrt{\gamma}=\psi^{6},$ $\displaystyle{\tilde{\rho}}_{*}$
$\displaystyle=\frac{1}{r^{2}{\rm sin}\theta}\rho_{*},$ $\displaystyle
v^{(r)}=v^{r},~{}v^{(\theta)}$ $\displaystyle=rv^{\theta},~{}v^{(\phi)}=r{\rm
sin}\theta\,v^{\phi},$ $\displaystyle{\tilde{S}}_{(r)}=\frac{1}{r^{2}{\rm
sin}\theta}S_{r},~{}{\tilde{S}}_{(\theta)}$ $\displaystyle=\frac{1}{r^{3}{\rm
sin}\theta}S_{\theta},~{}{\tilde{S}}_{(\phi)}=\frac{1}{r^{3}{\rm
sin}^{2}\theta}S_{\phi},$ $\displaystyle{\tilde{S}}_{(i)(j)}$
$\displaystyle=\Lambda_{(i)}^{k}\Lambda_{(j)}^{l}S_{kl}.$ (16)

$\psi$ denotes the conformal factor. The indices without the parenthesis
denote the tensor components with respect to the coordinate basis, and the
variables of indices with the parenthesis correspond to the orthonormal basis
components in the spherical coordinates of the flat spacetime.

In this work, we numerically solve the set of these equations by employing a
Kurganov-Tadmor scheme (Kurganov & Tadmor, 2000) with a piecewise parabolic
reconstruction for the quantities of cell interfaces and the minmod-like
filter introduced in Kurganov & Tadmor (2000) for the flux-limitter.

##  Appendix B Particle tracing

To perform nucleosynthesis calculations, the Eulerian data obtained with a
simulation have to be translated into Lagrangian evolution of the
thermodynamical variables. For this purpose, the tracer particle method is
employed. In this appendix, we summarize our method of particle tracing and
how we use the nucleosynthesis results for the HD simulation.

###  B.1 Time evolution of particles

Suppose we have a time series of three-velocity field $v^{i(n)}$, where $i$ is
the index of the spatial coordinates, at a time slice $t=t^{(n)}$. Here, the
velocity field is defined at the spatial points
$(x^{1}_{j},x^{2}_{k},x^{3}_{l})$ discretely, where $j$, $k$, $l$ denote the
grid points in the three-dimensional space. For a given spatial position of a
particle $x^{i(n)}$ at the $n$th time slice, the coordinates of the particle
in the $(n+1)$th time slice is solved with the so-called semi-implicit
trapezoidal method (Fujibayashi et al., 2020a; Nishimura et al., 2015). In
this method, we solve

$\displaystyle\frac{x^{i(n+1)}-x^{i(n)}}{t^{(n+1)}-t^{(n)}}=\frac{1}{2}\big{(}\dot{x}^{i(n)}+\dot{x}^{i(n+1)}\big{)}$
(17)

for $x^{i(n+1)}$. Here, $\dot{x}^{i(n)}$ is the velocity of the particle at
$n$th time slice defined with the tri-linear interpolation from the velocity
field to the particle position. For the back-tracing of the ejecta particles,
the time step $\Delta t=t^{(n+1)}-t^{(n)}$ is limited by the frequency of the
outputs. Numerical accuracy of the particle tracing may be diagnosed with a
timescale defined as

$\displaystyle\Delta
t^{\prime}=\min\bigg{(}\frac{x^{1}_{j+1}-x^{1}_{j}}{\dot{x}^{1(n)}},\frac{x^{2}_{k+1}-x^{2}_{k}}{\dot{x}^{2(n)}},\frac{x^{3}_{l+1}-x^{3}_{l}}{\dot{x}^{3(n)}}\bigg{)},$
(18)

where the spatial coordinates of the particle are in the intervals
$x^{1}_{j+1}$–$x^{1}_{j}$, $x^{2}_{k+1}$–$x^{2}_{k}$, and
$x^{3}_{l+1}$–$x^{3}_{l}$ at $t=t^{(n)}$. $\Delta t^{\prime}<\Delta t$
indicates that the particle position changes significantly within the time
step $\Delta t$. In such a case, we perform several sub-steps between the
$(n+1)$th and $n$th time slices. The velocity fields at a time $t$ between
$t^{(n+1)}$ and $t^{(n)}$ are linearly interpolated from the two time slices.

###  B.2 Initial position of the particle

The particles are distributed on a sphere with the radius
$r_{\mathrm{ext}}=3\times 10^{8}$ cm. At a given time, the ejecta criterion is
checked at $r=r_{\mathrm{ext}}$ with a set of the polar and azimuthal angles
$(\theta_{i},\phi_{i})$, where $i$ is the index of the direction on which
particles are placed. If the ejecta criterion is satisfied at a direction, a
particle is placed at the angle. The particles are then traced backward in
time.

The polar angles in the set $\\{(\theta_{i},\phi_{i})\\}$ are on a uniform
grid in 0–$\pi/2$ with a spacing $\Delta\theta$. The azimuthal angles in the
set are also on a uniform grid in 0–$2\pi$ but with a $\theta$-dependent
spacing $\Delta\phi$, which decreases with the polar angle. The function of
$\Delta\phi$ is determined to distribute the particles in an approximately
uniform manner in the sphere.

The particles are placed repeatedly with a time interval $\Delta
t_{\textrm{p-set}}$, which is determined by the average velocity of the ejecta
at $r=r_{\mathrm{ext}}$ as

$\displaystyle\Delta
t_{\textrm{p-set}}=\frac{r_{\mathrm{ext}}\Delta\theta}{\langle v^{r}\rangle},$
(19)

where $\Delta\theta$ is the polar angle spacing of the particle-set locations
and $\langle v^{r}\rangle$ is the average radial velocity of the particles
located at a given time. In this way, the particles are located uniformly in
space, i.e., in similar radial and lateral distances.

The mass of each particle is assigned as

$\displaystyle m={r_{\mathrm{ext}}}^{2}\Delta\Omega\rho u^{r}\sqrt{-g}\Delta
t_{\textrm{p-set}}.$ (20)

This definition of the particle mass is consistent with the conserved mass
flux at $r=r_{\mathrm{ext}}$. Therefore, the total mass of the particles
converges to the ejecta mass of the NR simulation for increasing the number of
the particles set on the sphere (i.e., $\Delta\theta\rightarrow 0$ and thus
$\Delta t_{\textrm{p-set}}\rightarrow 0$).

###  B.3 Nucleosynthesis

In total, about 9600 tracer particles are generated with this method. Along
the thermodynamical histories of the tracer particles, nucleosynthesis
calculations are performed in the same manner as in Wanajo et al. (2022). The
initial $Y_{e}$ for the nucleosynthesis calculation is taken from the final
value of each tracer particle, which is that shown in Figs. 2, 4, and 5. We
define the particles with $Y_{e}\leq 0.08$ as those of dynamical ejecta. They
are found to have very similar nucleosynthesis results with only a small
difference in actinide abundances. Therefore, for the particles of the
dynamical ejecta, we assign the same elemental abundance and radioactive
heating rate, which are the mass-weighted average of those in Wanajo et al.
(2022). On the other hand, for the post-merger component ($Y_{e}>0.08$), we
assign each particle with the abundance and heating rate obtained by the
nucleosynthesis calculation along its thermodynamical history.

###  B.4 Mapping particle data for HD simulation

To take the radioactive heating into account in the HD simulation, and to
obtain the final spatial distribution of elements, the nucleosynthesis results
along with the tracer particles are used. For this purpose, the time and
angular positions at which fluid elements are injected through
$r=r_{\mathrm{ext}}$ are traced by advecting the three passive scalar
variables. The injection time and angles are then converted to the tracer
particles which reach the extraction radius at the similar times and angles.

In practice, we first construct a table of heating rates as functions of the
time ($t$), and the injection time ($t_{\mathrm{inj}}$) and angles
($\theta_{\mathrm{inj}}$, $\phi_{\mathrm{inj}}$). The tracer particles are
placed sparsely in terms of time and angles. Therefore, the properties of the
tracer particles (abundance and radioactive heating rates) have to somehow be
mapped onto the injection time and angle space.

For a given injection time and angles ($t_{\mathrm{inj}}$,
$\theta_{\mathrm{inj}}$, $\phi_{\mathrm{inj}}$), we select a few tracer
particles spatially closest to the (Cartesian) point
$\vec{x}_{\mathrm{inj}}=(r_{\mathrm{ext}}\sin\theta_{\mathrm{inj}}\cos\phi_{\mathrm{inj}},r_{\mathrm{ext}}\sin\theta_{\mathrm{inj}}\sin\phi_{\mathrm{inj}},r_{\mathrm{ext}}\cos\theta_{\mathrm{inj}})$
at the time $t_{\mathrm{inj}}$. The location of a tracer particle at the time
$t_{\mathrm{inj}}$, $\vec{x}_{i}(t_{\mathrm{inj}})$, is estimated as

$\displaystyle x_{i}(t_{\mathrm{inj}})$
$\displaystyle=r_{i}(t_{\mathrm{inj}})\sin\theta_{i}\cos\phi_{i},$ (21)
$\displaystyle y_{i}(t_{\mathrm{inj}})$
$\displaystyle=r_{i}(t_{\mathrm{inj}})\sin\theta_{i}\sin\phi_{i},$ (22)
$\displaystyle z_{i}(t_{\mathrm{inj}})$
$\displaystyle=r_{i}(t_{\mathrm{inj}})\cos\theta_{i},$ (23)

where $\theta_{i}$ and $\phi_{i}$ are the angles at which the particle is
located on $r=r_{\mathrm{ext}}$, and $r_{i}(t_{\mathrm{inj}})$ is the radius
of the $i$th particle estimated at $t=t_{\mathrm{inj}}$. The radius is
estimated as

$\displaystyle
r_{i}(t_{\mathrm{inj}})=r_{\mathrm{ext}}+v^{r}_{i}(t_{\mathrm{inj}}-t_{i}),$
(24)

where $v^{r}_{i}$ is the radial three velocity of the particle when it is
located on $r=r_{\mathrm{ext}}$ at the time $t_{i}$. In short, we assume that
tracer particles move only radially with the same velocity as what they have
at $r=r_{\mathrm{ext}}$. The radioactive heating rate of the fluid elements
injected at ($t_{\mathrm{inj}}$, $\theta_{\mathrm{inj}}$,
$\phi_{\mathrm{inj}}$) is then estimated as

$\displaystyle\dot{q}(t)=\sum_{i\in\mathcal{P}}w_{i}\dot{q}_{i}(t),$ (25)

where $\mathcal{P}$ is the set of the selected tracer particles, and $w_{i}$
is the weight that satisfies $\sum_{i}w_{i}=1$. We define the weight as
$w_{i}\propto 1/|\vec{x}_{\mathrm{inj}}-\vec{x}_{i}(t_{\mathrm{inj}})|$, i.e.,
the weight is inversely proportional to the distance between the injection
point and the particle at the same time $t_{\mathrm{inj}}$.