# Kilonovae of binary neutron star mergers leading to short-lived remnant
neutron star formation

Kyohei Kawaguchi,1,2,3 Sho Fujibayashi,1 Nanae Domoto,4 Kenta Kiuchi,1,3
Masaru Shibata,1,3 Shinya Wanajo,1
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
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We study kilonova emission from binary neutron star (BNS) mergers for the case
that a remnant massive neutron star (MNS) forms and collapses to a black hole
within $20$ ms after the onset of the merger (which we refer to as “a short-
lived case") by consistently employing numerical-relativity and
nucleosynthesis results. We find that such kilonovae are fainter and last
shorter than those for BNSs resulting in the formation of long-lived ($\gg
1\,{\rm s}$) MNSs, in particular in the optical band. The resulting light
curves are too faint and last for a too short duration to explain the kilonova
observation for the BNS associated with GW170817, indicating that the merger
remnant formed in GW170817 is unlikely to have collapsed to a black hole
within a short period of time ($\sim 20$ ms) after the onset of the merger.
Our present result implies that early observation is necessary to detect
kilonovae associated with BNSs leading to short-lived MNS formation in
particular for the optical blue band as well as that kilonovae could be hidden
by the gamma-ray burst afterglow for nearly face-on observation. We provide a
possible approximate scaling law for near-infrared light curves with the given
reference time and magnitude when the decline power of the z-band magnitude,
$dM_{\it z}/d{\rm log}_{10}t$, reaches $2.5$. This scaling law suggests that
the HK-band follow-up observation should be at least $1$ mag deeper than that
for the z-band reference magnitude and earlier than 4 times the reference
time.

###### keywords:

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

††pubyear: 20XX††pagerange: Kilonovae of binary neutron star mergers leading
to short-lived remnant neutron star formation–13

##  1 Introduction

Binary neutron star (BNS) mergers are among the most efficient gravitational-
wave emitters in the universe and the most important sources of multi-
messenger high-energy astrophysical phenomena, such as gamma-ray bursts (GRB,
Paczynski, 1991; Nakar, 2007; Berger, 2014; Abbott et al., 2017c), kilonovae
(Li & Paczynski, 1998; Kulkarni, 2005; Metzger et al., 2010; Kasen et al.,
2013; Tanaka & Hotokezaka, 2013), and synchrotron flares (Nakar & Piran, 2011;
Hotokezaka & Piran, 2015; Hotokezaka et al., 2018; Margalit & Piran, 2020).
Furthermore, BNS mergers are 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). All these facts
imply that BNS mergers are unmissable research subjects from an astronomical
point of view. They are also among the unique systems in the universe in which
the most extreme (strongly self-gravitating, high-density, and high-
temperature) environments in the universe are realized. Hence, the multi-
messenger observation of BNS mergers is also an indispensable tool to extend
our knowledge of fundamental physics.

Quantitative prediction of the merger dynamics and outcomes is crucial to
correctly interpret the observed signals. Since the first simultaneous
detection of gravitational waves and electromagnetic (EM) signals from a BNS
(GW170817/AT2017gfo; Abbott et al. 2017a), remarkable progress has been
achieved in the theoretical understanding, particularly, in the studies based
on numerical simulations. For example, recent numerical studies revealed the
quantitative nature of mass ejection from BNS mergers, for which the processes
can be broadly divided into two phases: At the onset of the merger, a fraction
of neutron-rich matter is ejected by tidal force and collisional shock heating
(e.g., Rosswog et al., 1999; Ruffert et al., 2001; Hotokezaka et al., 2013).
After the merger, a massive neutron star (MNS) or a black hole (BH) surrounded
by a strongly magnetized hot and dense accretion torus is formed (e.g., Price
& Rosswog, 2006; Kiuchi et al., 2018; Kiuchi et al., 2022a). The magnetized
central objects and accretion tori are considered to launch relativistic jets
and outflows by magnetic pressure and tension, viscous heating due to magneto-
hydrodynamical turbulence, and neutrino irradiation. Quantitative properties
of the ejecta and the nucleosynthetic element abundances for each phase are
studied by various groups together with their dependence on binary parameters,
such as NS masses and NS equations of state (EoS, Hotokezaka et al. 2013;
Bauswein et al. 2013; Wanajo et al. 2014; Sekiguchi et al. 2015; Foucart et
al. 2016; Sekiguchi et al. 2016; Radice et al. 2016; Dietrich et al. 2017;
Bovard et al. 2017; Kiuchi et al. 2018; Dessart et al. 2009; Metzger &
Fernández 2014; Perego et al. 2014; Just et al. 2015; Wu et al. 2016; Siegel &
Metzger 2017; Shibata et al. 2017; Lippuner et al. 2017; Fujibayashi et al.
2018; Siegel & Metzger 2018; Ruiz et al. 2018; Fernández et al. 2019; Christie
et al. 2019; Perego et al. 2019; Miller et al. 2019; Fujibayashi et al. 2020a,
b, c; Bernuzzi et al. 2020; Ciolfi & Kalinani 2020; Nedora et al. 2021;
Foucart et al. 2020; Fernández et al. 2020; Mösta et al. 2020; Shibata et al.
2021a, b; Curtis et al. 2022; Fujibayashi et al. 2023; Kiuchi et al. 2022a;
Foucart et al. 2022; Just et al. 2023; Curtis et al. 2023; see Shibata &
Hotokezaka 2019 for a review). The light curve modeling of EM counterparts,
particularly for kilonovae, are also developed in this decade by employing
numerical-simulation-based/motivated ejecta profiles and by performing
radiative transfer simulations with realistic heating rates and/or detailed
opacity tables (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).

However, there are still various open questions remaining. For example,
whether the remnant NS has gravitationally collapsed into a BH or not is still
being an open question for GW170817 due to the lack of the detection of post-
merger gravitational waves in GW170817 (Abbott et al., 2017d). Such
information is important, because it is connected to the underlying physics of
the uncomprehended NS EoS (e.g., Margalit & Metzger, 2017; Rezzolla et al.,
2018; Shibata et al., 2019). While we expect that the observation of the EM
counterparts can provide a great hint to address this issue, it is still
unclear from what observational features we can know about the fate of the
remnant. Focusing particularly on the kilonova emission, a general consensus
has not been yet reached for the property and origin of the ejecta in GW170817
(e.g., Kasliwal et al., 2017; Cowperthwaite et al., 2017; Kasen et al., 2017;
Villar et al., 2017; Waxman et al., 2018; Kawaguchi et al., 2018; Kawaguchi et
al., 2020; Bulla, 2019; Almualla et al., 2021; Kedia et al., 2023; Bulla,
2023). Determination of the ejecta property is crucial for understanding the
post-merger evolution of the system and whether BNS mergers could be the major
production site of $r$-process elements in the universe.

To address these questions, quantitative understanding of the relation between
the initial condition and/or underlying physics, and EM signals is important.
For this purpose, conducting a study based on numerical simulations
consistently starting from the merger to the phase of EM emission is a useful
approach to link the observables that should be related to each other. In
particular, for the kilonova modeling, it is important to accurately determine
the ejecta profile for the rest-mass density and compositions at the time of
kilonova emission ($>0.1$ d). Previous studies showed that the ejecta profile
induces significant spacial dependence in radioactive heating as well as
strong geometrical effects in radiative transfer, which have great impact on
the resultant light curves (Kasen et al., 2015; Wollaeger et al., 2018;
Kawaguchi et al., 2020; Bulla, 2019; Zhu et al., 2020; Darbha & Kasen, 2020;
Korobkin et al., 2021; Almualla et al., 2021; Kedia et al., 2023). However,
there are still limited number of studies which provide the end-to-end
modeling from the merger to observational outputs following the hydrodynamics
evolution of all the ejecta components up to the time of kilonova emission
(Kawaguchi et al. (2021, 2022); Just et al. (2023); see, however, Rosswog et
al. (2014); Grossman et al. (2014); Collins et al. (2023); Neuweiler et al.
(2023) for the studies focusing on the dynamical ejecta components, and
Fernández et al. (2015, 2017); Foucart et al. (2021) in the context of BH-NS
mergers). Given the situation that a number of BNS mergers will be observed in
the next decades, the EM counterpart prediction based on the consistent
simulations by taking the BNS diversity into account is an urgent task for
correctly interpreting the observed data.

In this paper, we study the kilonova light curves of BNS mergers for the case
that a remnant MNS forms and subsequently collapses to a BH within $20$ ms
after the onset of the merger (which we refer to as “a short-lived case")
consistently employing numerical-relativity (NR) results of Kiuchi et al.
(2022b); Fujibayashi et al. (2023). This paper is organized as follows: In
Section 2, we describe the method employed in this study. In Section 3, we
describe the BNS models we study in this work. In Section 4, we present the
property of the ejecta obtained by the long-term hydrodynamics evolution and
the kilonova light curves obtained by radiative-transfer simulations. Finally,
we discuss the implication of this paper in Section 5. Throughout this paper,
$c$ denotes the speed of light.

##  2 Method

Merger ejecta of a BNS are expected to be homologously expanding at the time
of kilonova emission ($\gtrsim 0.1\,{\rm d}$). To obtain the ejecta profile in
the homologously expanding phase, we follow the same procedures as in the
previous work (Kawaguchi et al., 2021, 2022); adopting the outflow data
obtained by NR simulations as the inner boundary condition (Fujibayashi et
al., 2023), the hydrodynamics evolution of merger ejecta is calculated by
employing an axisymmetric relativistic hydrodynamics code developed in
Kawaguchi et al. (2021, 2022). In the following, to distinguish between the
present simulation and NR simulation, we refer to the present hydrodynamics
simulations as the HD simulations.

In the hydrodynamics code, relativistic hydrodynamics equations in the
spherical coordinates are solved taking into account the effect of fixed-
background gravity of a non-rotating BH metric in the isotropic coordinates.
Radioactive-decay heating of heavy elements is also taken into account by
referring to the nucleosynthesis results computed for each ejecta fluid
element in the NR simulation (see Fujibayashi et al. (2023) for the details).
We employ the ideal-gas EoS with the adiabatic index of $\Gamma=4/3$. For the
HD simulations, the uniform grid spacing with $N_{\theta}$ grid points is
prepared for the polar angle $\theta$, while for the radial direction, the
following non-uniform grid structure is employed; the $j$-th radial grid point
is given by

$\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})=(2048,256)$, and $r_{\rm in}$ and $r_{\rm ext}$ are
initially set to be $8,000\,{\rm km}$ and $10^{3}\,r_{\rm in}$, respectively.
We employ the same time origin for the HD simulations as in the NR simulations
for the post-merger evolution.

To import the outflow data from the NR simulations of Fujibayashi et al.
(2023) to the present HD simulations, the time-sequential hydrodynamics
property of the outflow is extracted at $r=r_{\rm in}$ in the NR simulations,
and is used as the boundary condition at the inner radius, $r=r_{\rm in}$, of
the HD simulations. The NR simulation data are run out at $t>5$ s, and after
then, the HD simulation is continued by setting a very small floor-value,
which is negligible for the ejecta dynamics, 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}$. We note that the total
mass lost by removing the inner radial grids is always much smaller ($\lesssim
10^{-4}\,M_{\odot}$) than the post-merger ejecta mass.

The light curves of kilonovae are calculated using a wavelength-dependent
radiative transfer 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 element 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. The impact of this assumption will be
discussed in Appendix A.

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) are 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, line data are
replaced with those calibrated with the atomic data of VALD and NIST database
(Kramida et al., 2021).

The radiative transfer 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}$. The spatial distributions of the heating rate
and element abundances are determined by the table obtained by the
nucleosynthesis calculations referring to the injected time and angle of the
fluid elements. Note that the element abundances at $t=1\,{\rm d}$ are used
during the entire time evolution in the radiative transfer 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).

##  3 Model

Table 1: Key model parameters. The columns describe the model name, the EoS adopted, the masses of the NSs, type of the MNS evolution, the ejecta mass evaluated in the NR simulations ($M^{\rm NR}_{\rm eje}$, $M^{\rm NR}_{\rm dyn}$, and $M^{\rm NR}_{\rm post}$ denote the total, dynamical, and post-merger masses, respectively; see Fujibayashi et al. 2020c; Shibata et al. 2021b; Fujibayashi et al. 2023), and the ejecta mass evaluated in the HD simulations at $t=0.1\,{\rm d}$, $M^{\rm HD}_{\rm eje}$, respectively. “short-lived”, “long-lived”, and “long-lived with strong dynamo” denote the cases for which the remnant MNS collapses to a BH within 20 ms, survives for $\gg$ 1 s, and survives for $\gg$ 1 s with significant magnetic dynamo effects, respectively. The values for $M^{\rm NR}_{\rm eje}$ are calculated by integrating the mass flux at the sphere with radius 8,000 km over time in 2D NR simulations. We then subtract from $M^{\rm NR}_{\rm eje}$ the mass of dynamical ejecta $M^{\rm NR}_{\rm dyn}$, which is evaluated in the corresponding 3D NR simulations with the Bernoulli criterion to obtain the contribution of the post-merger ejecta $M^{\rm NR}_{\rm post}$. Model | EoS | $(m_{1}\,[M_{\odot}],m_{2}\,[M_{\odot}])$ | MNS evolution | $M^{\rm NR}_{\rm eje}(M^{\rm NR}_{\rm dyn},M^{\rm NR}_{\rm post})\,[10^{-2}M_{\odot}]$ | $M^{\rm HD}_{\rm eje}\,[10^{-2}M_{\odot}]$
---|---|---|---|---|---
SFHo-135135 | SFHo | $(1.35,1.35)$ | short-lived | $1.0\,(0.73,\,0.25)$ | $1.0$
SFHo-130140 | SFHo | $(1.30,1.40)$ | short-lived | $1.0\,(0.48,\,0.50)$ | $0.9$
SFHo-125145 | SFHo | $(1.25,1.45)$ | short-lived | $1.2\,(0.64,\,0.60)$ | $1.1$
SFHo-120150 | SFHo | $(1.20,1.50)$ | short-lived | $1.6\,(0.45,\,1.1)$ | $1.5$
SFHo-125155 | SFHo | $(1.25,1.55)$ | short-lived | $1.5\,(0.95,\,0.55)$ | $1.4$
DD2-135135 | DD2 | $(1.35,1.35)$ | long-lived | $7.6\,(0.15,\,7.5)$ | $6.5$
MNS75a | DD2 | $(1.35,1.35)$ | long-lived with strong dynamo | $9.4\,(0.15,\,9.3)$ | $8.4$

In this work, we employ the NR outflow profiles obtained in Fujibayashi et al.
(2023) as the input for the HD simulations. The key quantities of each model
are summarized in Table 1. The first four models listed in Table 1 are BNSs
with the total gravitational mass (at the infinite separation) of
$2.7\,M_{\odot}$ but with various mass ratios in the range of 0.8–1.0. We also
study an unequal mass BNS with a larger total gravitational mass
($2.8\,M_{\odot}$), which we refer to as SFHo-125155. The SFHo EoS (Steiner et
al., 2013) supplemented by the Timmes (Helmholtz) EoS (Timmes & Swesty, 2000)
for the low density part is employed. For all the models employing the SFHo
EoS, a remnant MNS is formed after the merger, but it collapses to a BH within
$\approx 20\,{\rm ms}$. We note that these mass ranges of the BNSs with a
short-lived remnant broadly cover the range of the mass estimation obtained by
the gravitational-wave data analysis of GW170817 (Abbott et al., 2017b, 2019).
The BNS models which result in the formation of an MNS surviving for a long
time $(>1\,{\rm s};$ Fujibayashi et al. 2020c; Shibata et al. 2021b) are also
shown in Table 1 for comparison purposes (see also Kawaguchi et al. 2021,
2022). The NR simulations are performed by a general-relativistic viscous
neutrino-radiation hydrodynamics code with the dimensionless alpha viscous
parameter of $\alpha=0.04$ (Fujibayashi et al., 2020c, 2023) except for MNS75a
in which general-relativistic neutrino-radiation resistive-
magnetohydrodynamics code is employed to take the magnetic dynamo effects into
account (Shibata et al., 2021b).

The ejecta mass evaluated in the NR simulations is also listed in Table 1. The
total ejecta mass increases as the mass ratio of the BNS deviates from unity
due to the increase in the torus mass, and hence, the ejecta mass of the post-
merger component. Broadly speaking, the mass of the dynamical ejecta tends to
decrease as the binary becomes more asymmetric (but not so monotonically).
This reflects the fact that, for an asymmetric binary, the tidal-interaction-
driven component dominates the dynamical ejecta rather than the collisional
shock-driven component, of which the launching mechanism is more efficient in
mass ejection than the former. The total ejecta mass of the BNS merger for
which the remnant MNS collapses to a BH in a short time is an order of
magnitude smaller than that for the BNS which results in the formation of an
MNS surviving for a long time $(>1\,{\rm s};$ Fujibayashi et al. 2020c;
Shibata et al. 2021b). Note that for the latter case, the total ejecta mass is
dominated by the post-merger ejecta.

##  4 Results

###  4.1 Ejecta profiles

For all the models, we find that the total internal energy of ejecta is
smaller by $\approx 4$ order of magnitudes than the total kinetic energy at
$t=0.1\,{\rm d}$ and that the mass-averaged deviation of the velocity field
from that in the later homologous expanding phase ($v^{r}=r/t$ with $v^{r}$
being the radial velocity) is as small as $10^{-3}$ at $t=0.1\,{\rm d}$. This
shows that the homologous expansion is well achieved for $t\geq 0.1\,{\rm d}$.

The total mass in the computational domain measured at $t=0.1\,{\rm d}$,
$M^{\rm HD}_{\rm eje}$, is listed in Table 1. Note that the matter is in the
homologously expanding phase at $t=0.1\,{\rm d}$, and hence, $M^{\rm HD}_{\rm
eje}$ can be regarded as the total ejecta mass. It is found that $M^{\rm
HD}_{\rm eje}$ is slightly smaller than $M^{\rm NR}_{\rm eje}$ for some of the
models. This is a consequence of the fact that a fraction of the matter falls
back across the inner boundary as the pressure support from the inner boundary
vanishes when the outflow data run out. While a fraction of the matter can
actually experience such fall-back due to the deceleration by the pressure
from the precedingly ejected matter, our treatment of suddenly vanishing
pressure support on the inner boundary at the run-out time of NR data may
artificially increase the mass of the fall-back matter. Nevertheless, as found
in our previous studies (Kawaguchi et al., 2021, 2022), the contribution of
such marginally unbound matter to the kilonova emission is minor because it
has only low velocity and has only a small contribution to the emission due to
the long diffusion time scale.

Figure 1: Rest-mass density profiles at $t=0.1\,{\rm d}$ obtained by the HD
simulations. The top-left, top-right, bottom-left, and bottom-right panels
display the results for models SFHo-135135, SFHo-130140, SFHo-125145, and
SFHo-120150, respectively. The gray curves in each panel denote the contour
lines of $10^{-15}$, $10^{-14}$, $10^{-13}$, $10^{-12}$, $10^{-11}$,
$10^{-10}$, and $10^{-9}\,{\rm g/cm^{3}}$ from outside.

Figure 2: The same as Fig. 1 but for the electron fraction, $Y_{e}$. The value
of $Y_{e}$ is evaluated when the temperature of the fluid element decreases to
$T=5\,{\rm GK}$. Note that only the region of which the rest-mass density at
$t=0.1\,{\rm d}$ is higher than $10^{-14}\,{\rm g/cm^{3}}$ is shown.

Figure 3: The same as Figures 1 and 2 but for SFHo-125155.

First, we focus on the BNS models of which the total mass is $2.7\,M_{\odot}$
to see the effect of the binary mass ratio. Fig. 1 shows the rest-mass density
profiles at $t=0.1\,{\rm d}$ obtained by the HD simulations for models
SFHo-135135, SFHo-130140, SFHo-125145, and SFHo-120150. The dynamical ejecta
component located at $x/ct\gtrsim 0.05$ or $z/ct\gtrsim 0.15$ exhibits a
broadly spherical morphology in the rest-mass density structure. On the other
hand, the post-merger ejecta component, which is present in $x/ct\lesssim
0.05$ and $z/ct\lesssim 0.15$, exhibits a mildly prolate shape (see Fig. 2 for
a clearer distinction between the dynamical and post-merger ejecta
components). These characteristics of the density profile are in broad
agreement with the ejecta profile obtained in our previous studies (Kawaguchi
et al., 2021, 2022), in which BNSs result in long-lived MNSs (with the
lifetime of $>1\,{\rm s}$).

Taking a closer look, the dynamical ejecta show a relatively more prolate
shape for an equal-mass BNS (SFHo-135135), while relatively more oblate shapes
are seen for unequal mass cases (SFHo-125145 and SFHo-120150). This reflects
the fact that the tidally driven component which spreads preferentially toward
the equatorial direction dominates in the dynamical ejecta for an asymmetric
binary over the collisional-shock-driven component which spreads in a more
spherical manner.

Fig. 2 shows the electron fraction ($Y_{e}$) profiles at $t=0.1\,{\rm d}$ for
models SFHo-135135, SFHo-130140, SFHo-125145, and SFHo-120150. Here, the value
of $Y_{e}$ is evaluated when the temperature of the fluid element decreases to
$T=5\,{\rm GK}$ ($=5\times 10^{9}$ K). A clear boundary-like feature starting
from $x/ct\approx 0.05$ on the equatorial plane to $z/ct\lesssim 0.15$ along
the polar axis is seen for all the models. This corresponds to the boundary
between the dynamical and post-merger ejecta components. The dynamical ejecta
has a clear angular dependence in the $Y_{e}$ profile. With $\theta$ being the
angle measured from the polar axis, the value of $Y_{e}$ of the dynamical
ejecta is higher than $0.3$ for $\theta\lesssim 45^{\circ}$–$60^{\circ}$,
while it is lower than $0.3$ for $\theta\gtrsim 45^{\circ}$–$60^{\circ}$. This
clearly reflects the difference in the mass ejection mechanism; the former is
shock-heating-driven and the latter is tidally driven. The dynamical ejecta
for unequal mass BNSs have relatively more extended distribution and lower
$Y_{e}$ values along the equatorial direction than those for the equal-mass
case. This also reflects the fact that the tidally driven component dominates
the dynamical ejecta and the ejecta experience a relatively small rise in
temperature resulting from the shock heating for the unequal mass cases. On
the other hand, the post-merger ejecta has only weak angular dependence in the
$Y_{e}$ value, which is always $\gtrsim 0.3$. These profiles of $Y_{e}$ are
also in broad agreement with the previous results of BNS mergers that result
in long-lived remnant MNSs (Kawaguchi et al., 2021, 2022) and the results of
BNS mergers in which the remnant survives for a moderately long time (0.1–1 s)
(Just et al., 2023).

Fig. 3 shows the rest-mass density and electron fraction profiles for model
SFHo-125155. The qualitative features of the rest-mass density and $Y_{e}$
profiles for this model are the same as those for other models with the total
mass of 2.7 $M_{\odot}$, but the oblate shape and low-$Y_{e}$ value region are
more pronounced than those for the models shown in Figs. 1 and 2. This
reflects the fact that SFHo-125155 has the largest dynamical ejecta mass
dominated by the tidally driven component as the consequence of the large
asymmetry in the NS masses.

Figs. 1–3 illustrate that the profiles of the rest-mass density and electron
fraction depend sensitively on the total mass and mass ratio of the binaries.
In the following we will show that the light curve and the spectral evolution
depend on these differences, although the type of the remnant (either a short-
lived or long-lived neutron star is formed as a remnant) has more impact on
the brightness of the kilonova light curve.

###  4.2 Kilonova light curves

Figure 4: (Left panel) The bolometric light curves for all the models
considered in this paper. The solid and dashed curves denote the total and
isotropically equivalent bolometric luminosities (the latter measured from the
polar direction, $0^{\circ}\leq\theta\leq 20^{\circ}$), respectively. The
isotropically equivalent bolometric luminosity observed in AT2017gfo is shown
by the filled circles adopting the data in Waxman et al. (2018) with the
distance of 40 Mpc. (Right panel) the ratios of the bolometric fluxes measured
from the polar (the dashed curves; $0^{\circ}\leq\theta\leq 20^{\circ}$) and
equatorial directions (the dotted curves; $86^{\circ}\leq\theta\leq
90^{\circ}$) to those of spherical average.

The left panel of Fig. 4 shows the results of the bolometric light curves
obtained by radiative-transfer simulations. The solid and dashed curves
denote, respectively, the total and isotropically equivalent bolometric
luminosities (the latter measured from the polar direction,
$0^{\circ}\leq\theta\leq 20^{\circ}$). For all the models, the bolometric
light curves show approximately flat features with the luminosity of $\sim
10^{41}\,{\rm erg/s}$ for $0.3\,{\rm d}\leq t\leq$3–5 d, and decline rapidly
after 3–5 d. As the ejecta mass increases, the epoch at which the bolometric
light curve starts rapidly declining is delayed, and the luminosity after the
decline becomes larger. This reflects the larger total optical depth and
deposition energy for larger ejecta mass models.

The right panel of Fig. 4 shows the ratios of the bolometric fluxes measured
from the polar ($0^{\circ}\leq\theta\leq 20^{\circ}$) and equatorial
directions ($86^{\circ}\leq\theta\leq 90^{\circ}$) to those of spherical
average. The isotropically equivalent luminosities measured from the polar and
equatorial directions are brighter and fainter by a factor of $\approx 2$,
respectively, at $t\sim 1\,{\rm d}$ due to the preferential diffusion of
photons in the presence of optically thick dynamical ejecta around the
equatorial plane (Kawaguchi et al., 2018; Kawaguchi et al., 2020). However,
such effects become less significant in the late phase ($\gtrsim 10\,{\rm d}$)
as the optical depth of the ejecta decreases due to the expansion. The
viewing-angle dependence of the bolometric light curves is sustained for a
longer time scale as the binary becomes more asymmetric. This reflects the
fact that the tidally driven component of dynamical ejecta has more mass and a
lower value of $Y_{e}$ for more asymmetric binaries, resulting in more opaque
ejecta.

None of the model light curves in the left panel of Fig. 4 can explain the
observed brightness of the kilonova associated with GW170817 (AT2017gfo). The
bolometric light curves are always below the observational data from 0.5 d to
17 d except for the last two data points in the plot. This is the case even if
the enhancement of the brightness due to geometrical effects is taken into
account (see the dashed curves in the left panel of Fig. 4, which denote the
light curves measured from the polar direction). This is primarily due to the
smallness of the ejecta mass, which leads to insufficient total radioactive
deposition energy to explain the observation of AT2017gfo. Our results
indicate that a BNS for which a remnant MNS collapses to a BH in a short time
($t\lesssim 20\,{\rm ms}$) is unlikely to be the progenitor of GW170817. We
note that our light curves are fainter than the results of Just et al. (2023),
which considers the cases that a remnant MNS survives for a relatively longer
time scale before it collapses to a BH (at $t=0.1$–1 s after the onset of the
merger). This simply reflects the fact that the total ejecta mass is smaller
for our present models.

Figure 5: gzK-band light curves. The top, middle, and bottom panels denote the
light curves observed from $0^{\circ}\leq\theta\leq 20^{\circ}$,
$41^{\circ}\leq\theta\leq 46^{\circ}$, and $86^{\circ}\leq\theta\leq
90^{\circ}$, respectively. The purple, green, and red curves denote the g, z,
and K-band light curves, respectively. The data points denote the observation
data of AT2017gfo taken from Villar et al. (2017) with the distance of 40 Mpc.

Fig. 5 shows the gzK-band light curves for all the models of the short-lived
cases listed in Table 1. The obtained light curves show the broadly similar
properties to those obtained by the observation of AT2017gfo as well as the
previous studies for a kilonova with multiple ejecta components (e.g., Kasen
et al., 2015; Wollaeger et al., 2018; Kawaguchi et al., 2018; Bulla, 2019);
the optical emission lasts for a short time scale ($\sim 1\,{\rm d}$), and the
near-infrared (NIR) emission lasts for a longer time scale ($\sim 10\,{\rm
d}$). The emission becomes faint as the viewing angle measured from the axis
of symmetry increases. This primarily reflects the spatial dependence of
element abundances (see Figs. 2 and 3). The viewing-angle dependence is more
pronounced for the emission in the optical wavelength (i.e., in the g-band)
due to the so-called lanthanide-curtain effects in the presence of low-$Y_{e}$
dynamical ejecta around the equatorial plane (Kasen et al., 2015; Wollaeger et
al., 2018; Kawaguchi et al., 2020; Bulla, 2019; Zhu et al., 2020; Darbha &
Kasen, 2020; Korobkin et al., 2021).

Interestingly, the peak magnitudes in the NIR wavelengths (i.e., in the
K-band) do not significantly differ among the models regardless of the
difference in the ejecta mass. However, the time scale for the emission to
sustain the brightness close to the peak becomes shorten as the total ejecta
mass decreases. The light curves in the optical wavelengths observed from the
polar direction also show similar shapes among the models except for the most
asymmetric BNSs (SFHo-120150 and SFHo-125155) for which the g-band light
curves are fainter by $\geq 1$ mag than those for the other models. We find
that the strong suppression of the optical emission for the most asymmetric
BNS models is due to the fact that the polar regions are more polluted by the
lanthanide elements. The difference in the brightness of the optical emission
observed from the equatorial direction among the models simply reflects the
difference in the dynamical ejecta mass (see Table 1).

Our result implies that an earlier follow-up observation than in
GW170817/AT2017gfo is needed to observe the kilonova emission in the optical
band for the short-lived BNS formation. For example, for the hypothetical
distance of 200 Mpc, the g-band emission can only be detected by the
observation within 0.5–1 d with the sensitivity deeper than 22 mag, which
requires telescopes larger than 2 m-classes (Nissanke et al., 2013). Also,
such a detection can be achieved only for the case that the event is face-on,
but we should note that it could be hidden by the GRB afterglow emission. In
the z band, the emission lasts for a longer time scale, but yet the
observation within 1 d is needed with 2 m and 4 m-class telescopes,
respectively, to find kilonovae for the case of $\theta\leq 45^{\circ}$. The
NIR follow-up observation by a telescope larger than 4-m classes, such as
VISTA (Ackley et al., 2020), can detect the kilonova emission up to 5 d after
the onset of the merger with the hypothetical distance of 200 Mpc and 100 Mpc
for face-on and edge-on events, respectively. However, since the field of view
of an NIR telescope is not as large as that of the optical one (Sutherland et
al., 2015), the improvement in the source localization by the gravitational-
wave observation is crucial.

###  4.3 Comparison with different BNS models

Figure 6: Comparison of the gzK-band light curves among the models in which
remnant MNSs survive for a short time scale (the dashed curves; SFHo-125145)
and for a long time scale (the solid curves; DD2-135135, Fujibayashi et al.
(2020c); Kawaguchi et al. (2022)), and for the case in which significant
magnetic dynamo effects are present in a long-surviving remnant MNS (the
dotted curves; MNS75a, Shibata et al. (2021b); Kawaguchi et al. (2022)). The
g, z, and K-band light curves are shown in the top, middle, and bottom panels,
respectively.

Fig. 6 compares the gzK-band kilonova light curves for the BNS models for
which the remnant MNS survives for a short time scale (SFHo-125145) and for a
long time scale (DD2-135135, Fujibayashi et al. (2020c); Kawaguchi et al.
(2022)), and for the case that significant magnetic dynamo effects are
hypothetically present in a long-surviving remnant MNS (MNS75a, Shibata et al.
2021b; Kawaguchi et al. 2022). The time scale for the emission to rapidly
decline is much shorter for the model with a short-lived remnant MNS than
those for the models associated with the formation of a long-lived MNS simply
because the ejecta mass for the short-lived MSN models is smaller by a factor
of 5–10 than that for the latter cases. The brightness at the peak is also
high for the case with a long-lived MNS, and the difference is more
significant in a shorter wavelength.

As already mentioned, none of the merger models that result in a short-lived
remnant MNS can explain the peak kilonova brightness of AT2017gfo observed in
the gz-band, nor the brightness in the K-band in the late phase ($\gtrsim
5\,{\rm d}$). This is likely to be the case even if we consider a possible
enhancement in the optical-band emission due to the modification in the
ionization states by the non-LTE effects (see Appendix A). On the other hand,
the kilonova model of a BNS that results in a long-surviving MNS (DD2-135135)
reproduces the peak brightness in the optical wavelengths as well as the
brightness and declining time scale in the NIR wavelengths, although a
deviation from the observation is present in the optical wavelengths in the
late phase ($t\gtrsim 2$ d)111Taking the non-LTE effects on the ionization
populations into account may solve the tension; see Kawaguchi et al. (2021,
2022) for the discussion.. This suggests that the formation of a short-lived
remnant MNS is unlikely the case for GW170817 and the formation of an MNS
which survives for a longer time scale ($\gtrsim 0.1\,{\rm s}$) is more likely
from the viewpoint of kilonova light curves. However, for the case that the
significant magnetic dynamo effects are present in the long-surviving remnant
MNS (MNS75a), the kilonova emission will be significantly brighter than the
observed data (see the light curves of MNS75a in Fig. 6). This suggests that
the remnant MNS of GW170817 should have not survived for too long time (i.e.,
over the time scale of the dynamo magnetic-field amplification) if the
magnetic dynamo effect played a significant role in the post-merger phase (see
also the discussion below for the viewpoint of the nucleosynthesis yields).

###  4.4 Approximate scaling law of kilonova light curves

Figure 7: The time (top) and AB absolute magnitude (bottom) at which the
decline power of the z-band magnitude, $dM_{\it z}/d{\rm log}_{10}t$, reaches
to $2.5$ as functions of the viewing angle. The solid, dashed, and dash-dotted
curves denote the cases in which the remnant MNS survives for a long period
($t\gg 1\,{\rm s}$; DD2-125 and DD2-135 in Fujibayashi et al. 2020c; Kawaguchi
et al. 2021, 2022), the remnant MNS collapses to a BH in a short time
($t\lesssim 20\,{\rm ms}$; see Table 1), and the magnetic dynamo effects in
the long-lived MNS are considered (MNS70a, MNS75a, and MNS80 in Shibata et al.
2021b; Kawaguchi et al. 2022). Figure 8: The gzK-band light curves for various
models and viewing angles for which the time and magnitude are scaled by those
at which the decline power of the z-band magnitude, $dM_{\it z}/d{\rm
log}_{10}t$, reaches $2.5$. The solid, dashed, and dash-dotted curves denote
the long-lived, short-lived, and long-lived dynamo cases, respectively, as in
Fig. 7. The light curves observed from $0^{\circ}\leq\theta\leq 20^{\circ}$,
$28^{\circ}\leq\theta\leq 35^{\circ}$, $59^{\circ}\leq\theta\leq 64^{\circ}$,
and $86^{\circ}\leq\theta\leq 90^{\circ}$ are considered. The scaled
observational data of AT2017gfo in the gzK-band taken from Villar et al.
(2017) are also plotted by circles with error bars.

While the peak brightness and the time scale of the emission differ among
different BNS models and setups, Fig. 6 implies that the shapes of the light
curves as well as their relative brightness among different wavelengths share
similar behaviour among the models. To examine this idea, we compare the gzK-
band light curves for various models and viewing angles with the time and
magnitude of each light curve being scaled by those at a certain reference
time.

For this purpose, we chose the reference time for each light curve to be the
decline time of the z-band emission, $t_{\rm z,dec}$, defined as the time at
which the decline power of the z-band magnitude, $dM_{\it z}/d{\rm
log}_{10}t$, reaches $2.5$. Fig. 7 shows the reference time and z-band
magnitude as functions of the viewing angle for various kilonova models
(Kawaguchi et al., 2021, 2022). The reference time and magnitude largely vary
among the models and viewing angles. As expected from Fig. 6, the reference
time and magnitude tend to be earlier and fainter, respectively, for the
short-lived cases than the long-lived cases. The viewing-angle dependence is
more pronounced for the short-lived cases, which reflects the fact that the
dynamical component has a larger fraction in the total ejecta compared to the
long-lived cases.

Fig. 8 compares the gzK-band light curves for various models and viewing
angles, which are scaled with the reference time and $z$-band magnitude for
each case. The g-band light curves show a large diversity among the models
even after the scaling, for which we find no clear trend among the models and
viewing-angles. On the other hand, although the reference time and magnitude
largely vary among the models and viewing angles, the K-band light curves show
relatively a less diversity after the scaling. In particular, the value of the
K-band magnitude is always within $\approx 1$ mag relative to the value of the
reference z-band magnitude for $0.6\lesssim t/t_{\rm z,dec}\lesssim 4$. We
find that this is also the case for the H band. Hence, this suggests that the
HK-band follow-up observation should be at least $1$ mag deeper than the value
of the z-band reference magnitude and earlier than 4 times the reference time.

Once the kilonova candidate is found and the decline time is determined by the
z-band observation in a few days after the event, this approximate scaling law
can be used as a guideline for the NIR follow-up observation by letting us
know how rapid and how deep the observation should be. For example, let us
suppose the case for which an EM candidate is found in the z band and $dM_{\it
z}/d{\rm log}_{10}t$ reaches $2.5$ with the z-band magnitude being 20 mag at
1.5 d after the gravitational-wave trigger. Then, our approximate scaling-law
suggests that the follow-up observation deeper than 21 mag within 6 d is at
least needed not to miss the peak brightness of the HK-band counterparts.

Notably, the K-band emission tends to decline within $t/t_{\rm z,dec}\approx
5$–$10$ for the cases with a long-lived remnant MNS, while the K-band
magnitude for the cases with a short-lived remnant MNS tends to keep the value
close to the peak until a larger value of $t/t_{\rm z,dec}$. The observational
data of AT2017gfo in the gzK-band scaled in the same way tend to follow the
trend of the cases with a long-lived remnant MNS, which also supports our
hypothesis that the remnant MNS for GW170817 did not collapse to a BH within a
short time ($<20\,{\rm ms}$).

##  5 Discussions

We found that the kilonova light curves of a BNS of which the remnant MNS
survives for a short time are too faint and last for a too short duration to
explain the brightness of the optical and NIR observation of
GW170817/AT2017gfo. This is primarily due to the smallness of ejecta mass.
Instead, kilonova models of a BNS which results in a long-surviving MNS
(DD2-135135) are more consistent with the observation. This indicates that the
remnant MNS of GW170817 might not have collapsed within a short time
($\lesssim 20\,{\rm ms})$ but survived for a longer time ($\gtrsim 0.1\,{\rm
s}$). On the other hand, our previous study (Kawaguchi et al., 2022) indicated
that, if the dynamo effects play a significant role for an efficient
amplification of magnetic fields in a long-lived remnant MNS, the kilonova as
well as the synchrotron emission stemming from the interaction between the
ejecta fast tail and inter-stellar medium becomes too bright to be consistent
with the EM observations associated with GW170817 (see also the discussion in
Sarin et al. (2022)). Hence, the remnant MNS should have collapsed to a BH
within the dynamo time scale of the magnetic-field growth, or the dynamo
effect in the post-merger phase was subdominant.

Figure 9: Isotropic equivalent ejecta mass for various models and latitudinal
angles. The solid, dashed, and dash-dotted curves denote the long-lived,
short-lived, and long-lived dynamo cases, respectively, as in Fig. 7. The
purple, green, and orange curves denote the results for $\theta=0^{\circ}$,
$30^{\circ}$, and $90^{\circ}$, respectively.

We find that the mass distribution of the ejecta in the polar region for the
long-lived case is also compatible with the required property of the fast blue
component, for which the origin is often discussed to be mysterious (e.g.,
Kasliwal et al., 2017; Cowperthwaite et al., 2017; Kasen et al., 2017; Villar
et al., 2017; Waxman et al., 2018; Kawaguchi et al., 2018; Kawaguchi et al.,
2020; Bulla, 2019; Almualla et al., 2021; Kedia et al., 2023; Bulla, 2023).
Fig. 9 shows the isotropic equivalent ejecta mass, $M^{\rm iso}_{\rm
eje}(v^{r},\theta)$, for various models and latitudinal angles, which is
defined by

$\displaystyle M^{\rm iso}_{\rm
eje}(v^{r},\theta)=4\pi\int_{>v^{r}}\rho(r,\theta)r^{2}dr,$ (2)

where $\rho$ denotes the rest-mass density. For the case of long-lived MNS
formation (DD2-135135), the polar value of $M^{\rm iso}_{\rm eje}$ for
$v^{r}\gtrsim 0.2\,c$ is larger than $10^{-2}\,M_{\odot}$. This matches the
property of the ejecta which is required to explain the luminosity and photo-
spheric velocity of the blue component in AT2017gfo (see also Just et al. 2023
for similar findings). Such a polar ejecta component is originated from the
dynamical ejecta component and the post-merger ejecta component of which the
velocity is enhanced by neutrino-radiation from the MNS. While the spectral
analysis with the non-LTE effects being taken into account is needed for a
more quantitative argument, our finding suggests that the photo-spheric
velocity of the blue component can be naturally explained by the setup
obtained by NR simulations.

Fig. 9 suggests that the presence of the diversity in the evolution of photos-
spheric velocities reflects the different types of the MNS evolution. For the
case of short-lived MNS formation (SFHo-125145), the value of $M^{\rm
iso}_{\rm eje}$ only reaches $10^{-2}\,M_{\odot}$ for $v_{r}<0.05\,c$, simply
reflecting the smallness of the ejecta mass. This suggests that the photo-
spheric velocity of the short-lived case is $\lesssim 0.05\,c$ for $t\gtrsim
1$ d. On the other hand, the result of MNS75a shows that an applicable amount
of ejecta is distributed in the very high velocity components. This is due to
the acceleration of ejecta in the presence of significant magnetic dynamo
effects in the long-lived MNS, and the photo-spheric velocity of $>0.8\,c$ is
expected be observed in the early phase of emission for such a case.

Figure 10: Comparison of nucleosynthesis yields among the cases in which
remnant MNSs survive for a short time (red; SFHo-125145) and for a long time
(blue; DD2-135135, Fujibayashi et al. 2020c), and for the case in which
significant magnetic dynamo effects are present in a long-surviving remnant
MNS (olive; MNS75a, Shibata et al. 2021b). The $r$-process residuals to the
solar system abundances (Lodders et al., 2009) are also shown by gray curves,
which are scaled to match the abundance of 153Eu for SFHo-125145 as well as
that for DD2-135135.

As described above, the BNS that results in a long-lived MNS is more likely
the case for GW170817 than the BNS that results in a short-lived remnant MNS
from the viewpoint of kilonova light curves. However, the calculated
nucleosynthesis yields for such long-lived MNS cases (DD2-135135 and MNS75a in
Fig. 10) exhibit overproduction of the nuclei between the first and second
$r$-process abundance peaks ($A\sim 80$–130) when compared to the solar
$r$-process abundances (see also Fujibayashi et al. 2020c; Shibata et al.
2021b for the details, and Just et al. 2023 for similar results). This fact
suggests that such long-lived MNSs should not be the major outcomes of BNSs
that merge in a Hubble time if the dominant sources of $r$-process elements
are BNS mergers. This implies that GW170817 may not be a typical type of BNS
mergers in the universe.

However, we should note that the total nucleosynthesis yields can be sensitive
to the setups and physical ingredients of the numerical simulation. A latest
work suggests that a more self-consistent magnetohydrodynamics treatment of
angular momentum transfer could result in more production of elements heavier
than the first $r$-process peak in the post-merger ejecta (Kiuchi et al.,
2022a). Hence, there may still exist a room that both the observation of
GW170817 and the robustness of the solar abundance pattern (Cowan et al.,
2021) can be explained by some configuration of a BNS, while we should remind
that the presence of an MNS which survives for a long time scale ($t>1\,{\rm
s}$) with significant dynamo effects is unlikely the case of GW170817 as
discussed above. For example, a BNS which results in a remnant MNS with
significant dynamo effects but collapses to a BH at $O(0.1)$ s can be a
plausible model for interpreting GW170817 from this point of view.

For the BNS resulting in a short-lived MNS, the kilonova emission lasts over a
time scale appreciably shorter than that of GW170817/AT2017gfo, in particular
for the optical band. This implies that for detecting kilonovae of this type,
we need observation earlier than that for AT2017gfo. This is in particular the
case for a large value of $\theta$. It is also likely that the optical light
curves could be more easily hidden by the afterglow light curves of GRBs for
the small value of $\theta$. Hence, the NIR light curves may be the primary
target of the observation in the simultaneous detection of a GRB.

Figure 11: Comparison between the optical and NIR observation in GRB130603B
and various kilonova models. The r\- and H-band light curves in the observer
frame are calculated by employing the redshift value of the source ($z=0.356$,
Thone et al., 2013; Cucchiara et al., 2013). The solid, dashed, and dash-
dotted curves denote the long-lived, short-lived, and long-lived dynamo cases,
respectively, as in Fig. 7. The gray dashed lines denote the GRB afterglow
light curves. The observational data points for the $r$\- and $H$-band
magnitudes (circles) in the GRB130603B observation and the GRB afterglow model
light curves are taken from Tanvir et al. (2013).

In fact, the comparison of our model light curves with the observation of
GRB130603B (Berger et al., 2013; Tanvir et al., 2013), with which a plausible
kilonova candidate is associated, indicates that the r-band emission for the
case of short-lived MNS formation (SFHo-125145) is likely hidden by the
afterglow emission (see Fig. 11). The brightness in the H band for the case of
short-lived MNS formation is also at most only comparable to that of the
afterglow emission. Hence, the progenitor of GRB130603B was unlikely to be a
BNS which results in the formation of a short-lived remnant MNS assuming the
excess in the H band is due to the kilonova emission. This also indicates that
GRB-associated kilonovae from BNSs leading to short-lived MNS formation could
be missed by being entirely hidden by the afterglows, which should result in a
number of simultaneous detection of gravitational waves with short GRBs but
lack of kilonovae in future. Indeed a statistical study shows that there are a
substantial fraction of previous short GRBs that are not associated with
kilonovae (Troja, 2023).

As the brightness of AT2017gfo is known to be broadly comparable with the
optical and NIR counterparts of GRB130603B (Rossi et al., 2020), the kilonova
model light curves for the cases of long-lived MNS formation (DD2-135135 and
MNS75a) are also consistent with the observation of GRB130603B; while the
r-band emission is hidden by the afterglow emission, the H band emission for
the long-lived cases is brighter than the afterglow emission, and is
consistent with the observed excess. This suggests that the progenitor of
GRB130603B is likely to be a BNS which results in the formation of a MNS that
survives more than $\sim 10\,{\rm ms}$.

Figure 12: Sr mass density profiles for SFHo-135135 and DD2-135135
(Fujibayashi et al., 2020c; Kawaguchi et al., 2022) at $t=1\,{\rm d}$.

In Watson et al. (2019); Domoto et al. (2021); Domoto et al. (2022);
Gillanders et al. (2022), spectral features observed in the data of AT2017gfo
are interpreted as the p-Cygni profiles by Sr (note, however, Perego et al.
(2022); Tarumi et al. (2023) suggested that the spectral features could be
also well interpreted by the absorption lines by He if non-LTE effects are
considered). Recently, Sneppen et al. (2023) performed a more detailed
analysis for those spectral features, and show that the Sr distribution of the
ejecta should have nearly spherical morphology. Fig. 12 shows the Sr mass
density profiles at $t=1\,{\rm d}$ for SFHo-135135 and DD2-135135 (Fujibayashi
et al., 2020c; Kawaguchi et al., 2022). The Sr distribution with the velocity
larger than $0.15\,c$ approximately exhibits a spherical morphology for
SFHo-135135. On the other hand, the Sr distributions for DD2-135135 as well as
the low-velocity part ($<0.15\,c$) for SFHo-135135 show mildly prolate shapes.
These aspherical features, which are in broad agreement with the results of
Just et al. (2023), are inconsistent with the implication of Sneppen et al.
(2023). Detailed quantitative spectral analysis taking into account various
uncertainties is nevertheless needed to clarify how severe the current tension
from the observational implication is, which we leave it for a future task.

## Acknowledgements

KK thanks Masaomi Tanaka and Eli Waxman for the valuable discussions. We also
thank Kenta Hotokezaka for helpful discussions. Numerical computation was
performed on Yukawa21 at Yukawa Institute for Theoretical Physics, Kyoto
University and the Sakura, Cobra, Raven clusters at Max Planck Computing and
Data Facility. The simulations were performed on Fugaku provided by RIKEN
through the HPCI System Research Project (Project ID: hp220174, hp230084), 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 (JP20H00158, JP21K13912, JP23H04900, 22KJ0317, 23H01772) of
JSPS/MEXT.

## References

  * Abbott et al. (2017a) Abbott B., et al., 2017a, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2017b) Abbott B. P., et al., 2017b, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2017c) Abbott B. P., et al., 2017c, Astrophys. J. Lett., 848, L13
  * Abbott et al. (2017d) Abbott B. P., et al., 2017d, Astrophys. J. Lett., 851, L16
  * Abbott et al. (2019) Abbott B. P., et al., 2019, Phys. Rev. X, 9, 011001
  * 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
  * 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
  * Bauswein et al. (2013) Bauswein A., Goriely S., Janka H. T., 2013, Astrophys. J., 773, 78
  * Berger (2014) Berger E., 2014, Ann. Rev. Astron. Astrophys., 52, 43
  * Berger et al. (2013) Berger E., Fong W., Chornock R., 2013, Astrophys. J., 774, L23
  * Bernuzzi et al. (2020) Bernuzzi S., et al., 2020, Mon. Not. Roy. Astron. Soc., 497, 1488
  * Bovard et al. (2017) Bovard L., Martin D., Guercilena F., Arcones A., Rezzolla L., Korobkin O., 2017, Phys. Rev., D96, 124005
  * Bulla (2019) Bulla M., 2019, Mon. Not. Roy. Astron. Soc., 489, 5037
  * Bulla (2023) Bulla M., 2023, Mon. Not. Roy. Astron. Soc., 520, 2558
  * Bulla et al. (2021) Bulla M., et al., 2021, Mon. Not. Roy. Astron. Soc., 501, 1891
  * Christie et al. (2019) Christie I. M., Lalakos A., Tchekhovskoy A., Fernández R., Foucart F., Quataert E., Kasen D., 2019, Mon. Not. Roy. Astron. Soc., 490, 4811
  * Ciolfi & Kalinani (2020) Ciolfi R., Kalinani J. V., 2020, Astrophys. J., 900, L35
  * Collins et al. (2023) Collins C. E., Bauswein A., Sim S. A., Vijayan V., Martínez-Pinedo G., Just O., Shingles L. J., Kromer M., 2023, Mon. Not. Roy. Astron. Soc., 521, 1858
  * 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
  * Cowperthwaite et al. (2017) Cowperthwaite P. S., et al., 2017, Astrophys. J., 848, L17
  * Cucchiara et al. (2013) Cucchiara A., Prochaska J. X., Perley D. A., Cenko S. B., Werk J., Cao Y., Bloom J. S., Cobb B. E., 2013, Astrophys. J., 777, 94
  * Curtis et al. (2022) Curtis S., Mösta P., Wu Z., Radice D., Roberts L., Ricigliano G., Perego A., 2022, Mon. Not. Roy. Astron. Soc., 518, 5313
  * Curtis et al. (2023) Curtis S., Bosch P., Mösta P., Radice D., Bernuzzi S., Perego A., Haas R., Schnetter E., 2023, arXiv
  * Darbha & Kasen (2020) Darbha S., Kasen D., 2020, arXiv
  * Dessart et al. (2009) Dessart L., Ott C., Burrows A., Rosswog S., Livne E., 2009, Astrophys. J., 690, 1681
  * Dietrich et al. (2017) Dietrich T., Ujevic M., Tichy W., Bernuzzi S., Bruegmann B., 2017, Phys. Rev., D95, 024029
  * Domoto et al. (2021) Domoto N., Tanaka M., Wanajo S., Kawaguchi K., 2021, Astrophys. J., 913, 26
  * 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
  * 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
  * Fernández et al. (2020) Fernández R., Foucart F., Lippuner J., 2020, Mon. Not. Roy. Astron. Soc., 497, 3221
  * Fernández et al. (2019) Fernández R., Tchekhovskoy A., Quataert E., Foucart F., Kasen D., 2019, Mon. Not. Roy. Astron. Soc., 482, 3373
  * Foucart et al. (2016) Foucart F., et al., 2016, Phys. Rev., D93, 044019
  * Foucart et al. (2020) Foucart F., Duez M. D., Hebert F., Kidder L. E., Pfeiffer H. P., Scheel M. A., 2020, Astrophys. J. Lett., 902, L27
  * Foucart et al. (2021) Foucart F., Moesta P., Ramirez T., Wright A. J., Darbha S., Kasen D., 2021, arXiv
  * Foucart et al. (2022) Foucart F., Duez M. D., Haas R., Kidder L. E., Pfeiffer H. P., Scheel M. A., Spira-Savett E., 2022, arXiv
  * 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
  * Fujibayashi et al. (2018) Fujibayashi S., Kiuchi K., Nishimura N., Sekiguchi Y., Shibata M., 2018, Astrophys. J., 860, 64
  * 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., Shibata M., Wanajo S., Kiuchi K., Kyutoku K., Sekiguchi Y., 2020b, Phys. Rev. D, 102, 123014
  * Fujibayashi et al. (2020c) Fujibayashi S., Wanajo S., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2020c, 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
  * Gillanders et al. (2022) Gillanders J. H., Smartt S. J., Sim S. A., Bauswein A., Goriely S., 2022, Mon. Not. Roy. Astron. Soc., 515, 631
  * Grossman et al. (2014) Grossman D., Korobkin O., Rosswog S., Piran T., 2014, Mon. Not. Roy. Astron. Soc., 439, 757
  * 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. (2013) Hotokezaka K., Kiuchi K., Kyutoku K., Okawa H., Sekiguchi Y.-i., Shibata M., Taniguchi K., 2013, Phys. Rev., D87, 024001
  * 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
  * Just et al. (2015) Just O., Bauswein A., Pulpillo R. A., Goriely S., Janka H. T., 2015, Mon. Not. Roy. Astron. Soc., 448, 541
  * 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., Vijayan V., Xiong Z., Bauswein A., Goriely S., Guilet J., Janka H.-T., Martínez-Pinedo G., 2023, arXiv
  * 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
  * Kasen et al. (2017) Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Nature
  * Kasliwal et al. (2017) Kasliwal M. M., et al., 2017, Science, 358, 1559
  * 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
  * Kedia et al. (2023) Kedia A., et al., 2023, Phys. Rev. Res., 5, 013168
  * Kiuchi et al. (2018) Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2018, Phys. Rev. D, 97, 124039
  * Kiuchi et al. (2022a) Kiuchi K., Fujibayashi S., Hayashi K., Kyutoku K., Sekiguchi Y., Shibata M., 2022a, arXiv
  * Kiuchi et al. (2022b) Kiuchi K., Held L. E., Sekiguchi Y., Shibata M., 2022b, Phys. Rev. D, 106, 124041
  * 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
  * Kupka et al. (1999) Kupka F., Piskunov N., Ryabchikova T. A., Stempels H. C., Weiss W. W., 1999, A&AS, 138, 119
  * Kurucz & Bell (1995) Kurucz R. L., Bell B., 1995, Atomic line list
  * Lattimer & Schramm (1974) Lattimer J. M., Schramm D. N., 1974, Astrophys. J., 192, L145
  * Li & Paczynski (1998) Li L.-X., Paczynski B., 1998, Astrophys. J., 507, L59
  * Lippuner et al. (2017) Lippuner J., Fernández R., Roberts L. F., Foucart F., Kasen D., Metzger B. D., Ott C. D., 2017, Mon. Not. Roy. Astron. Soc., 472, 904
  * Lodders et al. (2009) Lodders K., Palme H., Gail H. P., 2009, Landolt Börnstein, 4B, 712
  * Margalit & Metzger (2017) Margalit B., Metzger B. D., 2017, Astrophys. J., 850, L19
  * Margalit & Piran (2020) Margalit B., Piran T., 2020, MNRAS, 495, 4981
  * Metzger & Fernández (2014) Metzger B. D., Fernández R., 2014, Mon. Not. Roy. Astron. Soc., 441, 3444
  * Metzger et al. (2010) Metzger B. D., et al., 2010, Mon. Not. Roy. Astron. Soc., 406, 2650
  * Miller et al. (2019) Miller J. M., et al., 2019, Phys. Rev. D, 100, 023008
  * Mösta et al. (2020) Mösta P., Radice D., Haas R., Schnetter E., Bernuzzi S., 2020, Astrophys. J. Lett., 901, L37
  * Nakar (2007) Nakar E., 2007, Phys. Rept., 442, 166
  * Nakar & Piran (2011) Nakar E., Piran T., 2011, Nature, 478, 82
  * 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
  * Nedora et al. (2021) Nedora V., et al., 2021, Astrophys. J., 906, 98
  * Neuweiler et al. (2023) Neuweiler A., Dietrich T., Bulla M., Chaurasia S. V., Rosswog S., Ujevic M., 2023, Phys. Rev. D, 107, 023016
  * Nissanke et al. (2013) Nissanke S., Kasliwal M., Georgieva A., 2013, Astrophys. J., 767, 124
  * Paczynski (1991) Paczynski B., 1991, Acta Astron., 41, 257
  * Perego et al. (2014) Perego A., Rosswog S., Cabezón R. M., Korobkin O., Kaeppeli R., Arcones A., Liebendorfer M., 2014, Mon. Not. Roy. Astron. Soc., 443, 3134
  * Perego et al. (2019) Perego A., Bernuzzi S., Radice D., 2019, Eur. Phys. J. A, 55, 124
  * Perego et al. (2022) Perego A., et al., 2022, ApJ, 925, 22
  * 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
  * Price & Rosswog (2006) Price D., Rosswog S., 2006, Science, 312, 719
  * Radice et al. (2016) Radice D., Galeazzi F., Lippuner J., Roberts L. F., Ott C. D., Rezzolla L., 2016, Mon. Not. Roy. Astron. Soc., 460, 3255
  * Rezzolla et al. (2018) Rezzolla L., Most E. R., Weih L. R., 2018, Astrophys. J. Lett., 852, L25
  * Rossi et al. (2020) Rossi A., et al., 2020, Mon. Not. Roy. Astron. Soc., 493, 3379
  * Rosswog et al. (1999) Rosswog S., Liebendoerfer M., Thielemann F. K., Davies M. B., Benz W., Piran T., 1999, Astron. Astrophys., 341, 499
  * Rosswog et al. (2014) Rosswog S., Korobkin O., Arcones A., Thielemann F., Piran T., 2014, Mon. Not. Roy. Astron. Soc., 439, 744
  * Ruffert et al. (2001) Ruffert M., Ruffert H. T., Janka H. T., 2001, Astron. Astrophys., 380, 544
  * Ruiz et al. (2018) Ruiz M., Shapiro S. L., Tsokaros A., 2018, Phys. Rev., D98, 123017
  * 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
  * Sarin et al. (2022) Sarin N., Omand C. M. B., Margalit B., Jones D. I., 2022, MNRAS, 516, 4949
  * Sekiguchi et al. (2015) Sekiguchi Y., Kiuchi K., Kyutoku K., Shibata M., 2015, Phys. Rev., D91, 064059
  * Sekiguchi et al. (2016) Sekiguchi Y., Kiuchi K., Kyutoku K., Shibata M., Taniguchi K., 2016, Phys. Rev., D93, 124046
  * Shibata & Hotokezaka (2019) Shibata M., Hotokezaka K., 2019, Ann. Rev. Nucl. Part. Sci., 69, 41
  * Shibata et al. (2017) Shibata M., Fujibayashi S., Hotokezaka K., Kiuchi K., Kyutoku K., Sekiguchi Y., Tanaka M., 2017, Phys. Rev., D96, 123012
  * Shibata et al. (2019) Shibata M., Zhou E., Kiuchi K., Fujibayashi S., 2019, Phys. Rev. D, 100, 023015
  * Shibata et al. (2021a) Shibata M., Fujibayashi S., Sekiguchi Y., 2021a, Phys. Rev. D, 103, 043022
  * Shibata et al. (2021b) Shibata M., Fujibayashi S., Sekiguchi Y., 2021b, Phys. Rev. D, 104, 063026
  * Siegel & Metzger (2017) Siegel D. M., Metzger B. D., 2017, Phys. Rev. Lett., 119, 231102
  * Siegel & Metzger (2018) Siegel D. M., Metzger B. D., 2018, Astrophys. J., 858, 52
  * Sneppen et al. (2023) Sneppen A., Watson D., Bauswein A., Just O., Kotak R., Nakar E., Poznanski D., Sim S., 2023, Nature, 614, 436
  * Steiner et al. (2013) Steiner A. W., Hempel M., Fischer T., 2013, Astrophys. J., 774, 17
  * Sutherland et al. (2015) Sutherland W., et al., 2015, A&A, 575, A25
  * Tanaka & Hotokezaka (2013) Tanaka M., Hotokezaka K., 2013, Astrophys. J., 775, 113
  * 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
  * Tarumi et al. (2023) Tarumi Y., Hotokezaka K., Domoto N., Tanaka M., 2023, arXiv
  * Thone et al. (2013) Thone C. C., de Ugarte Postigo A., Gorosabel J., Tanvir N., Fynbo J. P. U., 2013, GRB Coordinates Network, 14744, 1
  * Timmes & Swesty (2000) Timmes F. X., Swesty F. D., 2000, ApJS, 126, 501
  * Troja (2023) Troja E., 2023, Universe, 9, 245
  * Villar et al. (2017) Villar V. A., et al., 2017, Astrophys. J., 851, L21
  * Wanajo et al. (2014) Wanajo S., Sekiguchi Y., Nishimura N., Kiuchi K., Kyutoku K., Shibata M., 2014, Astrophys. J., 789, L39
  * Watson et al. (2019) Watson D., et al., 2019, Nature, 574, 497
  * 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. (2016) Wu M.-R., Fernández R., Martínez-Pinedo G., Metzger B. D., 2016, Mon. Not. Roy. Astron. Soc., 463, 2323
  * 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
  * 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 Estimate of uncertainties due to non-LTE Effects

Figure 13: Comparison of the gK-band light curves observed from
$0^{\circ}\leq\theta\leq 20^{\circ}$ for models in which the presence of both
neutral and first ionized atoms (“w/o I and II”; dotted curves) are
prohibited. The light curves shown in the solid curves (“Default”) are the
results with the default setting and are the same as those in Fig. 5. The left
and right panels show the results of the g and K-band light curves,
respectively. The top and bottom panels show the results for models
SFHo-135135 and SFHo-125145, respectively. The data points denote the
observational data of AT2017gfo taken from Villar et al. (2017) for a
hypothetical distance to the source of 40 Mpc.

In our present radiative-transfer simulations, the LTE condition is assumed to
determine the ionization/excitation populations of atoms. This assumption can
be invalid for the later phase of kilonova emission, at which the ionization
of the atoms caused by radioactive decay becomes more significant than the
recombination of ions (Hotokezaka et al., 2021; Kawaguchi et al., 2021).
Hotokezaka et al. (2021); Pognan et al. (2022) indicate that such non-LTE
effects may suppress the neutral and first ionized ions in the outer part of
the ejecta even in the earlier phases. Such modifications in the
ionization/excitation populations of atoms can have a great impact to the
opacity and resultant light curves.

Because computing the ionization/excitation population is challenging due to
its computational complexity and lack of the atomic data for $r$-process
elements (see Hotokezaka et al. 2021; Pognan et al. 2022), we here provide
qualitative estimates for the impacts of the non-LTE effects to the kilonova
light curves following the same prescription which we applied in our previous
studies (Kawaguchi et al., 2021, 2022); we perform the radiative transfer
simulations with a hypothetical setup in which both neutral and first ionized
atoms are artificially forced to be ionized to the second ionization states.
Note that this prescription is applied to whole ejecta, including high-density
regions for simplicity.

Fig. 13 shows the g and K-band light curves for models SFHo-135135 and
SFHo-125145 obtained with these hypothetical setups. As is also found in
Kawaguchi et al. (2021, 2022), the emission in the optical wavelengths is
enhanced by artificially increasing the ionization degrees. Yet, the
brightness of the g-band emission is not high enough to explain the brightness
of AT2017gfo. The emission in the NIR wavelengths in the late phase becomes
even fainter and more inconsistent with the observation. These results
indicate that the BNS that results in a short-lived remnant MNS is, at least
in our studied range of binary configurations, likely to be different from the
BNS of GW170817.