#  Fallback Accretion Halted by $R$-process Heating in Neutron Star Mergers
and Gamma-Ray Bursts

Wataru Ishizaki Center for Gravitational Physics, Yukawa Institute for
Theoretical Physics, Kyoto University, Kyoto, 606-8502, Japan  Kenta Kiuchi
Center for Gravitational Physics, Yukawa Institute for Theoretical Physics,
Kyoto University, Kyoto, 606-8502, Japan  Max Planck Institute for
Gravitational Physics (Albert Einstein Institute), Am Mühlenberg, Potsdam-
Golm, 14476, Germany  Kunihito Ioka Center for Gravitational Physics, Yukawa
Institute for Theoretical Physics, Kyoto University, Kyoto, 606-8502, Japan
Shinya Wanajo Max Planck Institute for Gravitational Physics (Albert Einstein
Institute), Am Mühlenberg, Potsdam-Golm, 14476, Germany  Interdisciplinary
Theoretical and Mathematical Sciences Program (iTHEMS), RIKEN, Wako, Saitama
351-0198, Japan

###### Abstract

The gravitational wave event GW170817 with a macronova/kilonova shows that a
merger of two neutron stars ejects matter with radioactivity including
$r$-process nucleosynthesis. A part of the ejecta inevitably falls back to the
central object, possibly powering long-lasting activities of a short gamma-ray
burst (sGRB), such as extended and plateau emissions. We investigate the
fallback accretion with the $r$-process heating by performing one-dimensional
hydrodynamic simulations and developing a semi-analytical model. We show that
the usual fallback rate $dM/dt\propto t^{-5/3}$ is halted by the heating
because pressure gradients accelerate ejecta beyond an escape velocity. The
suppression is steeper than Chevalier’s power-law model through Bondi
accretion within a turn-around radius. The characteristic halting timescale is
$\sim 10^{4}$–$10^{8}$ sec for the GW170817-like $r$-process heating, which is
longer than the typical timescale of the long-lasting emission of sGRBs. The
halting timescale is sensitive to the uncertainty of the $r$-process. Future
observation of fallback halting could constrain the $r$-process heating on the
year scale.

hydrodynamics — accretion, accretion discs — nucleosynthesis — gamma-ray
burst: general

##  1 INTRODUCTION

The discovery of the short-duration gamma-ray burst (sGRB), GRB 170817A,
coincided with the detection of GW170817 by the Laser Interferometer
Gravitational-Wave Observatory (LIGO) and the Virgo Consortium (LVC), is a
direct evidence that binary neutron star (BNS) mergers are one of the sources
of sGRBs (Abbott et al., 2017a, b). The simultaneous detection of
macronova/kilonova emission is a strong indication for the ejection of
neutron-rich matter and $r$-process element synthesis in this BNS merger
(Arcavi et al., 2017; Chornock et al., 2017; Coulter et al., 2017;
Cowperthwaite et al., 2017; Drout et al., 2017; Kasen et al., 2017; Kasliwal
et al., 2017; Kilpatrick et al., 2017; McCully et al., 2017; Nicholl et al.,
2017; Shappee et al., 2017; Smartt et al., 2017; Soares-Santos et al., 2017;
Tanaka et al., 2017; Tanvir et al., 2017). The observed properties of the
gravitational waves and electromagnetic counterparts are broadly consistent
with the predictions from a series of theoretical studies including numerical
relativity calculations (Shibata et al., 2017; Cowperthwaite et al., 2017;
Kasen et al., 2017; Kasliwal et al., 2017; Villar et al., 2017; Ioka &
Nakamura, 2018, 2019) and these simultaneous detection has established the
paradigm that sGRBs originate from compact binary mergers (Paczynski, 1986;
Goodman, 1986; Eichler et al., 1989).

Although sGRBs are classified as GRBs with the duration of prompt emission
being less than $2\sec$ (Kouveliotou et al., 1993), their central engines are
thought to remain active for longer time (Burrows et al., 2005; Ioka et al.,
2005). This was suggested from the observations of extended emission which
lasts for about $100\sec$ with a luminosity of $10^{48}$–$10^{49}~{}{\rm
erg}~{}{\rm s}^{-1}$ and possibly from the observations of plateau emission
which lasts for about $10^{4}\sec$ with a luminosity of
$10^{45}$–$10^{46}~{}{\rm erg}~{}{\rm s}^{-1}$ (Barthelmy et al., 2005;
Rowlinson et al., 2013; Kisaka et al., 2017; Kagawa et al., 2019). These long-
lasting emissions are fainter than the prompt emission in luminosity. However,
because of the long duration of these emissions, they are comparable with or
greater than the prompt emission in fluence (Kisaka et al., 2017). Some of the
extended emissions have been observed to darken abruptly, which are not
expected in the sGRB afterglow (Ioka et al., 2005). Therefore, a central
engine activity is very likely to explain these long-lasting emission111Note
that, although such a rapid fade-out has not been reported for the plateau
emission, it may also reflect the late-time activity of the central engine
(see also Matsumoto et al., 2020).. There are two models for the energy source
of the long-lasting emissions, which have been extensively discussed in
previous works. One is the release of rotational or magnetic field energy from
a strongly magnetized massive neutron star that forms after a BNS merger
(Zhang & Mészáros, 2001; Metzger et al., 2008; Bucciantini et al., 2012; Fan
et al., 2013; Murase et al., 2018). Currently, this scenario is not conclusive
because the expected observational features, e.g., late radio emission, have
not been detected (Metzger & Bower, 2014; Horesh et al., 2016).

The other scenario for the energy source of long-lasting activity is the
fallback accretion of ejecta (Rosswog, 2007; Lee & Ramirez-Ruiz, 2007; Rossi &
Begelman, 2009; Kisaka & Ioka, 2015; Kisaka et al., 2017). Numerical
calculations of compact binary mergers show that a part of the ejecta is still
gravitationally bound (e.g., Rosswog et al., 1999; Bauswein et al., 2013;
Kyutoku et al., 2015; Radice et al., 2016; Kiuchi et al., 2017). Kyutoku et
al. (2015) calculated the coalescence of a black hole and a neutron star
(BH–NS) based on numerical relativity and showed that a part of the ejecta
falls back to the merger remnant. While the gravitational energy released by
fallback accretion is large enough to explain the extended emission and the
plateau emission, the simple theory of fallback accretion, which assumes a
zero-temperature fluid, predicts that the mass accretion rate is proportional
to the power of time as $t^{-5/3}$ with no typical timescale (Rees, 1988;
Michel, 1988). However, the observed light curve of the long-lasting emission
of sGRB clearly has a certain timescale, which is not compatible with the
simple theory (Kagawa et al., 2019).

The coincidence of the macronova/kilonova emission with the gravitational wave
source GW170817 indicates that the ejecta of the BNS merger is heated by the
radioactive decay of $r$-process elements. Therefore, the assumption of zero-
temperature fluid is inappropriate (Smartt et al., 2017; Waxman et al., 2018;
Kawaguchi et al., 2018). Metzger et al. (2010) discussed the effect of ejecta
heating by the radioactive decay of $r$-process elements on the mass accretion
using a test-particle model. Desai et al. (2019) performed a more
sophisticated calculation based on the model in Metzger et al. (2010), using
the ejecta profiles obtained from numerical relativity simulations and the
radioactive heating rates obtained from nucleosynthesis calculations. They
showed that the mass accretion stops and resumes after a finite time, the so-
called “gap”, because the marginally bound fluid elements become unbound by
the heating due to the radioactive decay of $r$-process elements. Furthermore,
they argued that the timescale of this resumption is $\mathcal{O}(100)\sec$,
which is in agreement with the timescale of the extended emission.

In the test-particle model of Metzger et al. (2010) and Desai et al. (2019),
they assumed that all the radioactive energy is converted to the kinetic
energy. In reality, the energy from the radioactive heating is converted to
the kinetic energy through the pressure gradient force resulting from the
increased internal energy. It is unclear whether the pressure gradient is
large enough to convert all the internal energy of the fluid element into the
kinetic energy, so that the validity of their assumption is also unclear.
Because their assumption cannot be verified within the framework of the test-
particle model, it is necessary to solve the hydrodynamic equations
incorporating the effects of radioactive heating.

In this paper, we reconsider the effect of the radioactive heating due to
decaying $r$-process elements on the fallback accretion in BNS mergers by
numerically solving one-dimensional fluid equations. We also construct a semi-
analytical model that reproduces the hydrodynamical simulation results and
explores it over a large parameter space including the realistic $r$-process
heating. This paper is organized as follows. In Section 2, we describe the
method of the numerical calculation of the fluid equations and show the
results. In Section 3 and Section 4, the semi-analytical models of the
accretion flow are developed for a constant heating rate and a heating rate of
a broken power-law form, respectively. In Section 5, the semi-analytical model
is applied to explore the parameters relevant to the radioactive heating. In
Section 6, we summarize this work, and discuss the scope of application of the
spherically symmetric modeling used in this study and future prospects.

##  2 Hydrodynamical Simulation of fallback accretion

###  2.1 Numerical method

In order to investigate the effect of the radioactive heating due to
$r$-process nuclei on the fallback accretion in the BNS merger, we have
performed long-term one-dimensional hydrodynamic simulations of the matter
ejected during the merger. The ejecta profiles of the velocity and the mass
density are derived from the numerical relativity simulations performed by
Kiuchi et al. (2017). The hydrodynamical equations for spherically symmetric
and purely radial flow, including the heating and the point source gravity,
are written as follows:

$\frac{\partial\rho}{\partial t}+\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\rho v\right)=0,$ (1) $\frac{\partial v}{\partial
t}+v\frac{\partial v}{\partial r}=-\frac{1}{\rho}\frac{\partial P}{\partial
r}-\frac{GM}{r^{2}},$ (2) $\frac{\partial}{\partial t}\left(\frac{1}{2}\rho
v^{2}+\epsilon_{\rm int}\right)+\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left[r^{2}v\left(\frac{1}{2}\rho v^{2}+\epsilon_{\rm
int}+P\right)\right]\\\ =-\rho v\frac{GM}{r^{2}}+Q_{\rm heat},$ (3)

where $\rho$ is the mass density, $v$ is the radial bulk velocity of the
fluid, $P$ is the pressure, $\epsilon_{\rm int}$ is the internal energy, $M$
is the mass of the central object, and $Q_{\rm heat}$ is the radioactive
heating rate per unit volume per unit time.222The radioactive heating due to
the nuclear decay is the conversion of the mass deficit of the nucleus into
the internal energy. Therefore, the sum of the rest mass energy and the
internal energy should be conserved, and it is not strictly correct to simply
add the external heating term to the right-hand side of Equation (3). Here, we
ignore the term $-Q_{\rm heat}/c^{2}$, which should be in the right-hand side
of Equation (1). In the decay of a nucleus, the ratio of the mass deficit to
the nucleon mass is about $\mathcal{O}({\rm MeV})/\mathcal{O}({\rm GeV})$.
Since the heating due to this mass deficit is comparable to the internal
energy, the magnitude of the term $-Q_{\rm heat}/c^{2}$ relative to the left
hand side of Equation (1) is of the magnitude $\mathcal{O}(v^{2}/c^{2})$. In
this paper, we ignore the $\mathcal{O}(v^{2}/c^{2})$ term because we discuss
the motion of ejecta in the Newtonian limit. For a more formal formulation,
see Uchida et al. (2017). The equation of state is assumed to follow the
$\Gamma$-law,

$P=\left(\Gamma-1\right)\epsilon_{\rm int},$ (4)

where $\Gamma$ is the adiabatic index. In Equation (4), we adopt the radiation
dominant case of $\Gamma=4/3$. Here we neglect the radiative loss.
Furthermore, we also neglect the self-gravity of the ejecta, because the
ejecta mass is much smaller than the central mass.

We solve Equations (1)–(4) by using a one-dimensional hydrodynamics code with
the Newtonian gravity. The advection term of the hydrodynamic equations is
solved by using the HLL method (e.g., Del Zanna & Bucciantini, 2002) and the
3rd-order MUSCL reconstruction (see Balsara, 2017, for a review). In the range
of the parameters used in this paper, the system can be well described in a
non-relativistic regime. The calculation is performed by dividing the radius
from $40\ {\rm km}$ to $80,000\ {\rm km}$ into a uniform grid of 16,384 cells
(about $\sim 4.9\ {\rm km}$ per grid). The boundary conditions imposed on the
inner and outer boundaries are such that the radial differential coefficient
is set to zero for all physical quantities. Note that, because the inner
boundary is far inside from the sonic radius, the inner boundary condition
corresponds to a free-flow condition to the central object. In order to check
the convergence, we perform the calculation with a spatial resolution that is
twice as fine; as a result, the maximum relative error in the mass accretion
rate is found to be about 0.2% (the result is shown in Figure 3).

###  2.2 Initial profile of merger ejecta

Figure 1:  Radial profiles of the radial velocity (left axis) and the mass
density (right axis). The thick lines represent the model curves at $t=0$ used
in the calculations (see text). The thin lines with different styles show the
result of Kiuchi et al. (2017) at $\approx 10$ ms after the merger for the
zenith angles given in the legend. Note that only a profile of the dynamical
ejecta is shown without a subsequent component (such as viscosity-driven
wind). The bold yellow dashed line indicates the escape velocity. The region
inside $\sim 490\ {\rm km}$ is gravitationally bound.

We model the initial profile of dynamical ejecta based on the results of the
numerical relativity simulation of the BNS merger performed in Kiuchi et al.
(2017). The mass for each NS is 1.35 $M_{\odot}$. The EOS of the NS matter is
described by the two segments piece-wise polytropic EOS, which is referred to
as the H model in Kiuchi et al. (2017). In this particular model, the radius
of the NS with $1.35M_{\odot}$ is 12.3 km.

Figure 1 shows the radial profiles of the radial 3-velocity and the mass
density. The thin lines in Figure 1 show the result of Kiuchi et al. (2017) at
$\approx 10$ ms after the merger for each zenith angle given in the legend.
The dashed yellow line in the figure shows the escape velocity from the
gravitational field of the central object with mass $M=2.7M_{\odot}$. It can
be confirmed that the fluid elements located at $r\lesssim 490\ {\rm km}$ are
gravitationally bound. Around the radius $r\sim 490\ {\rm km}$, which is
important for the calculation of mass accretion rates, the radial dependence
of velocity and density is found to be weekly dependent on latitude. This is
one justification of our treatment of spherically symmetric modelling.

In this paper, the initial velocity and density profiles shown in Figure 1 are
modeled as described below. Given that the expansion of the dynamical ejecta
can be regarded as a homologous expansion, we adopt a model in which the
radial 3-velocity is proportional to the radius (out to $3000\ {\rm km}$):

$\displaystyle v(t=0,r)=\left\\{\begin{array}[]{ll}\displaystyle
0.26\,c\left(\frac{r}{1000\ {\rm km}}\right)&\left(r<3000\ {\rm km}\right)\\\
\displaystyle 0&\left(r>3000\ {\rm km}\right)\end{array}\right.,$ (7)

where $t=0$ is set at the beginning of the fluid calculations in this study
and corresponds to the time slice of the numerical relativity simulation in
Kiuchi et al. (2017) ($\approx 10$ ms after the merger). The radius of $r\sim
3000\ {\rm km}$ corresponds to that the outermost ejecta with nearly the speed
of light reach during about $10~{}{\rm ms}$. For the mass density, we adopt a
broken power-law model, which has a break at $600\ {\rm km}$, namely

$\rho(t=0,r)=\left\\{\begin{array}[]{ll}\displaystyle\rho_{0}\left(\frac{r}{600\
{\rm km}}\right)^{-2.4}&\left(r<600\ {\rm km}\right)\vspace{1mm}\\\
\displaystyle\rho_{0}\left(\frac{r}{600\ {\rm km}}\right)^{-7.5}&\left(r>600\
{\rm km}\right)\end{array}\right.,$ (8)

where $\rho_{0}=2.0\times 10^{6}~{}{\rm g~{}cm^{-3}}$ is the value at $r=600\
{\rm km}$. We assume $P=10^{-5}\rho c^{2}$ as the initial pressure
distribution in the ejecta, because the initial internal energy is
sufficiently low and does not affect the motion of the ejecta.

Fujibayashi et al. (2020) have pointed out that the viscosity-driven wind can
be launched with the timescale of $\mathcal{O}(1)\sec$, which dominates the
total mass of the ejecta from a BNS merger (see also Fernández & Metzger,
2013; Just et al., 2015; Fernández et al., 2019). Kawaguchi et al. (2020) have
studied further the long term temporal evolution of the viscosity-driven wind
based on the results of Fujibayashi et al. (2020). According to these studies,
after a time $\sim 10\sec$, the velocity profile of the wind approximately
matches that of a homologous expansion in Equation (7). The density
distribution is also expected to be close to $\rho\propto r^{-2}$ as that of a
steady-state supersonic flow, which is qualitatively similar to the radial
dependence of the density distribution expressed in Equation (8) for $r<600$
km. However, for the viscously driven wind, to take into account the
difference in total ejecta mass, it will be necessary to increase the value of
mass density $\rho_{0}$ by one order of magnitude compared with the value of
the dynamical ejecta. Here we presume that the disk wind is modeled by
Equations (7) and (8) with a larger value of density $\rho_{0}$ than that for
the dynamical ejecta case. Note that, at the time of interest for considering
the effect of radioactive heating, the mass accretion rate is determined by
the profile of the marginally bound ejecta rather than that of the overall
ejecta. Therefore, the detailed modeling of the mass density profile is
probably not very important.

###  2.3 Radioactive heating

Figure 2: Temporal evolution of the radioactive heating rates (excluding the
energies in neutrinos) adopted from Wanajo et al. (2014). The thick line
represents the mass-averaged heating rate. For the fluid elements of various
electron fractions $Y_{\rm e}$, the respective heating rate is plotted as a
thin line with the color indicated by the color bar.

Figure 2 plots the radioactive heating rates adopted from the results of
nucleosynthesis calculations based on the numerical model of a BNS merger
(Wanajo et al., 2014). The heating is due to $\beta$-decay, $\alpha$-decay,
and fission of $r$-process nuclei produced in the dynamical ejecta. Each thin
curve shows the heating rate in units of MeV nuc-1 s-1 as a function of time
(since the merger) for the Lagrangian tracer-particle of the ejecta with a
given initial $Y_{\mathrm{e}}$. The color indicates the value of
$Y_{\mathrm{e}}$ from 0.09 (purple) to 0.44 (yellow) with an interval of
$\Delta Y_{\mathrm{e}}=0.01$. The heating rate averaged over the ejecta mass
is also shown by the red thick curve.

As can be seen in Figure 2, the radioactive heating rates exhibit two phases:
one in which the value is approximately constant over time (constant phase)
and the other in which the value decays with time (decay phase). It can be
seen that the duration of the constant phase has a diversity and tends to be
longer for larger values of $Y_{\rm e}$. The decay phase can be well described
by a power law. Since the power-law indices are generally smaller than $-1$,
the most amount of radioactive energy is added during the constant phase. As a
first step, we ignore the decay phase (see Section 4 for the case with the
decay phase) and take the duration and heating rate of the constant phase as
parameters and model the heating rate per unit volume (see Equation (3)) as
follows:

$Q_{\rm heat}=\left\\{\begin{array}[]{ll}\rho\dot{q}_{0}&(t<t_{\rm heat})\\\
0&(t>t_{\rm heat})\end{array}\right.,$ (9)

where $\dot{q}_{0}$ is the radioactive energy (except for that in neutrinos)
released per unit mass per unit time during the constant phase333 For the
convenience of comparisons with the results by Metzger et al. (2010) and Desai
et al. (2019), we use the unit $~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$
($\sim 10^{18}$ erg g-1 s-1) for the heating rate $\dot{q}_{0}$.  and $t_{\rm
heat}$ is the duration of the constant phase. The value of $Y_{\rm e}$ can
take a variety of values depending on the outflow mechanism of ejecta as well
as the time of ejection. The values of $Y_{\rm e}$ in the dynamical ejecta of
BNS mergers are about $Y_{\rm e}\sim 0.1$–$0.4$. As mentioned in Section 2.2,
the dominant component of ejecta can be the late-time viscously driven wind.
This component is expected to have a large $Y_{\rm e}$ compared to that of the
dynamical ejecta, with $Y_{\mathrm{e}}\sim 0.3$–$0.4$, as shown in previous
work, e.g., Fujibayashi et al. (2020). Considering the uncertainties in
$Y_{\rm e}$ over the different components of ejecta, we adopt the simple model
described by Equation (9) and treat $\dot{q}_{0}$ and $t_{\rm heat}$ as
parameters. Furthermore, in Section 3, we model the accretion flows semi-
analytically. For this purpose, a simplified treatment in Equation (9) is
convenient. The modeling of the case with the heating rates including the
decay phase is given in Section 4.

Here, we neglect a possible effect of heating due to the jet which interacts
with the preceding ejecta by making a hot cocoon. A part of the energy of the
jet is converted into the internal energy of the cocoon, so that we can regard
the jet as an additional heating source for the ejecta. However, since this
occurs only for a narrow solid angle about the jet axis, the mass of the
heated ejecta (i.e., cocoon) is expected to be small relative to the total
ejecta mass (e.g., Hamidani & Ioka, 2021).

###  2.4 Results

Figure 3:  Temporal evolution of the mass accretion rate evaluated at $r=650\
{\rm km}$, where $t=0$ corresponds to the beginning of the calculation. The
thick red curve presents the result for the heating model with
$\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$ and $t_{\rm
heat}=5\sec$. The orange curve is the result for the same model but with the
twice finer spatial resolution. The relative error between the two different
resolutions of calculations is at most 0.2%. The gray curve represents the
result for the model without radioactive heating. The green dashed curve
represents the mass accretion rate calculated based on the test-particle model
(see Appendix A for details).

The thick red line in Figure 3 presents the time evolution of the mass
accretion rate calculated with $\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm
nuc^{-1}}~{}{\rm s^{-1}}$ and $t_{\rm heat}=5~{}\sec$ in Equation (9). Here,
in order to clearly illustrate the effect of radioactive heating on the mass
accretion rate, we have adopted the value of $t_{\rm heat}$ longer than those
shown in Figure 2 for the relevant $\dot{q}_{0}$. We assume that all the
radioactive energies (except for those in neutrinos) are converted to the
internal energies of the fluid elements. The thin gray line shows the fallback
rate in the absence of radioactive heating, which can be well described by
$dM/dt\propto t^{-5/3}$ (where the initial $\sim 0.05\sec$ is affected by the
initial conditions). Here, the mass accretion rate is evaluated using the mass
flux at $r=650\ {\rm km}$. It can be seen that the mass accretion rate is
substantially suppressed by the effect of radioactive heating. The break in
the red curve at $t=5~{}\sec$ corresponds to the termination of the heating
$t=t_{\rm heat}$. After $t=t_{\rm heat}$, the mass accretion rate continues
roughly in proportion to $t^{-5/3}$.

The dashed green line shows the mass accretion rate calculated using the same
method as the test-particle model of Metzger et al. (2010). The detail of the
method is described in Appendix A. Although the heating rate is the same as
that in the fluid model (red line), the time-dependence of the mass accretion
rate for the test-particle model differs from that for the fluid model. In the
fluid model, the mass accretion rate does not show a sharp cutoff as observed
in the test-particle model at $t\sim 5$ sec but slowly decreases taking a few
times longer duration. This is due to the difference in conversion of the
internal energy to the kinetic energy between the test-particle model and the
fluid model. In the test-particle model of Metzger et al. (2010), all
radioactive energy injected to a fluid element is assumed to be converted into
the kinetic energy. On the other hand, in the fluid model, the radioactive
heating first increases the internal energy (or pressure) of a fluid element,
and then is converted to the kinetic energy through the pressure gradient.
Actually, in the fluid model, a part of the internal energy from radioactive
heating is not converted to the kinetic energy and hence falls to the central
object with the fluid element. As the rate of conversion to kinetic energy is
higher, the more amount of ejecta tends to be unbound. As a result, the mass
accretion rate decreases more slowly in the fluid model than that in the test-
particle model.

Figure 4:  Dependence of the mass accretion rate on the heating parameters.
The gray and red thick curves are the same as those in Figure 3. The blue
curves represent the case for $\dot{q}_{0}=30~{}{\rm MeV}~{}{\rm
nuc^{-1}}~{}{\rm s^{-1}}$. The thick, dashed, and dot-dashed curves indicate
the results for the different lengths of $t_{\rm heat}$ as shown in the
legend. The black curves show the mass accretion rates calculated by the semi-
analytical model (see Section 3 for detail). The semi-analytical model is not
valid during the times indicated by the black dotted curves.

The dependence of the mass accretion rate on the heating parameters is shown
in Figure 4. The solid red curve is the same as the red curve in Figure 3. The
dot-dashed and thin-dashed curves show the results for $t_{\rm heat}=1\sec$
and $t_{\rm heat}=10\sec$, respectively. As can be seen from the figure, the
longer the radioactive energy injects, the longer the suppression of the mass
accretion rate continues. The blue curves are those for $\dot{q}_{0}=30~{}{\rm
MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$, where the solid and dashed curves
represent the cases for $t_{\rm heat}=1~{}\sec$ and $t_{\rm heat}=4~{}\sec$,
respectively. This value of $\dot{q}_{0}$ is a factor of a few greater than
that reached by radioactive heating (see Figure 2); we take this value for a
possible case with additional energy sources such as shock heating (e.g., due
to the subsequent viscosity-driven wind) or strong magnetic field. It can be
seen that the mass accretion rate is suppressed on a shorter timescale than
that for $\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$ because
$\dot{q}_{0}$ is larger. The semi-analytical modeling described in Section 3
explains these behaviors.

Figure 5:  Mass accretion rate using the mass-averaged heating rate of the
dynamical ejecta calculated by Wanajo et al. (2014) (red bold curve). The gray
curve shows the mass accretion rate without radioactive heating. At $10\sec$,
the mass accretion rate with heating is suppressed by about 70% compared to
that without heating.

Finally, we present the mass accretion rate calculated using a more realistic
heating rate rather than using a constant approximation in Equation (9).
Figure 5 shows the temporal evolution of the mass accretion rate calculated
using the mass-averaged heating rate of the dynamical ejecta in Wanajo et al.
(2014) (see Figure 2). At $t=10\sec$, the mass accretion rate with radioactive
heating becomes about 30% of that without heating. Although Figure 5 only
shows the results of numerical calculations until about $10\sec$, as will be
discussed in Section 5, the mass accretion rate is expected to be suppressed
to less than $\sim 10$% of that without heating after
$\mathcal{O}(10^{4})$–$\mathcal{O}(10^{8})$ sec due to the continuous heating
in the decay phase.

##  3 semi-analytical modeling of fallback accretion

###  3.1 Characteristic scales of the accretion flow

For better understanding of the numerical results obtained in Section 2, we
develop a semi-analytical modeling. We first introduce the characteristic
scales for the basic Equations (1)–(4). For this purpose, it is convenient to
rewrite Equation (3) into that in terms of $P/\rho$. By using Equations (1),
(2), and (4), the energy conservation law (3) can be written as follows:

$\frac{\partial}{\partial
t}\left(\frac{P}{\rho}\right)+v\frac{\partial}{\partial
r}\left(\frac{P}{\rho}\right)\\\
=(\Gamma-1)\left[\dot{q}_{0}-\frac{P}{\rho}\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}v\right)\right].$ (10)

The characteristic parameters of the fluid equations (1), (2), and (10) are
$GM$ and $\dot{q}_{0}$. Hence, the system has the characteristic scales of
length and time. From a dimensional analysis, these can be expressed in the
following forms:

$r_{c}=\left[\frac{(GM)^{3}}{\dot{q}_{0}^{2}}\right]^{1/5}\\\ \sim 3540~{}{\rm
km}~{}\left(\frac{M}{2.7~{}M_{\odot}}\right)^{3/5}\left(\frac{\dot{q}_{0}}{3~{}{\rm
MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}}\right)^{-2/5},$ (11)
$t_{c}=\left[\frac{(GM)^{2}}{\dot{q}_{0}^{3}}\right]^{1/5}\\\ \sim 0.35~{}{\rm
s}~{}\left(\frac{M}{2.7~{}M_{\odot}}\right)^{2/5}\left(\frac{\dot{q}_{0}}{3~{}{\rm
MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}}\right)^{-3/5}.$ (12)

Combining these quantities, a characteristic scale of velocity is also
obtained as

$v_{c}=\left(GM\dot{q}_{0}\right)^{1/5}\\\ \sim
0.033c~{}\left(\frac{M}{2.7~{}M_{\odot}}\right)^{1/5}\left(\frac{\dot{q}_{0}}{3~{}{\rm
MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}}\right)^{1/5}.$ (13)

Then, we introduce the dimensionless length $\xi\equiv r/r_{c}$, time
$\chi\equiv t/t_{c}$, density $\phi\equiv\rho/\rho_{c}$, velocity $u\equiv
v/v_{c}$, and pressure $\theta\equiv P/\rho
v_{c}^{2}=P/\phi\rho_{c}v_{c}^{2}$, where $\rho_{c}$ is an arbitrary constant
with the dimension of density. Rewriting Equations (1), (2) and (10) by using
these dimensionless variables, we obtain the following dimensionless
equations:

$\frac{\partial\phi}{\partial\chi}+\frac{1}{\xi^{2}}\frac{\partial}{\partial\xi}\left(\xi^{2}\phi
u\right)=0,$ (14) $\frac{\partial u}{\partial\chi}+u\frac{\partial
u}{\partial\xi}=-\frac{1}{\phi}\frac{\partial\left(\phi\theta\right)}{\partial\xi}-\frac{1}{\xi^{2}},$
(15)
$\frac{\partial\theta}{\partial\chi}+v\frac{\partial\theta}{\partial\xi}=(\Gamma-1)\left[1-\theta\frac{1}{\xi^{2}}\frac{\partial}{\partial\xi}\left(\xi^{2}u\right)\right].$
(16)

Because these non-dimensional equations have the same form for any $\rho_{c}$,
there is no characteristic scale of ejecta mass (equivalently mass density) in
this system (i.e., Equations (1), (2), and (10)). As mentioned in Section 2.2,
we presume that the profile of the viscosity-driven wind can be modeled by
enhancing the density $\rho$ (see Equation (8)). The invariance to the density
scale ensures that the subsequent semi-analytical model can be used for both
the viscosity-driven wind and the dynamical ejecta.

Figure 6:  Temporal evolution of turn-around radii (where $v=0$). The vertical
and horizontal axes are normalized by the characteristic scales in Equations
(11) and (12), respectively. The red and blue curves show the results for
$\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$ and
$\dot{q}_{0}=30~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$, respectively.
The thick and dashed curves differ in the length of $t_{\rm heat}/t_{c}$. The
calculation indicated by the red dashed curve is terminated at $t=32\,t_{c}$.
Note that the small difference between the red and the blue curves is due to
the fact that the initial conditions are not scaled.

Equations (1), (2), and (10) can be normalized by Equations (11)–(13) to
eliminate the parameters $GM$ and $\dot{q}_{0}$. Therefore, the accretion flow
follows the same equations under the characteristic scales (11)–(13) up to
$t=t_{\rm heat}$. Figure 6 shows the time evolution of the turn-around radius
$r_{\rm turn}$, at which the velocity becomes $v=0$. The vertical axis is
normalized by $r_{c}$ and the horizontal axis by $t_{c}$. As can be seen from
Figure 6, the difference between these curves is due solely to the difference
in $t_{\rm heat}/t_{c}$. This clearly shows the effectiveness of the scaling.
It is also seen that from $t=6\,t_{c}$ to $t_{\rm heat}$, the turn-around
radius $r_{\rm turn}$ is approximately constant $r_{\rm turn}\sim r_{c}$ over
time. Using the chain rule with the differential coefficients of velocity, we
can obtain the time evolution of $r_{\rm turn}$ as follows:

$\frac{dr_{\rm turn}(t)}{dt}=-\frac{\left(\partial v/\partial
t\right)_{r=r_{\rm turn}}}{\left(\partial v/\partial r\right)_{r=r_{\rm
turn}}}.$ (17)

At the turn-around radius $r=r_{\rm turn}$, the velocity is $v=0$ by
definition, so that the Lagrangian time derivative becomes $Dv/Dt=\partial
v/\partial t$. The fact that $r_{\rm turn}$ is approximately constant over
time, $dr_{\rm turn}(t)/dt\sim 0$ means $Dv/Dt\sim 0$; i.e., at the turn-
around radius, the gravity and pressure gradient forces are almost balanced
(see Equation (2)). This indicates that the ejecta outside $r_{\rm turn}$ does
not fall back.

###  3.2 The hydrodynamical structure inside $r_{\rm turn}$

Figure 7:  The top panel shows the time evolution of the ejecta mass inside
the turn-around radius $r_{\rm turn}$ (excluding that within a sufficiently
small radius taken to be $r_{\rm count}=650\ {\rm km}$). The horizontal axis
is normalized by the characteristic timescale in Equation (12). The solid red
curve shows the result obtained from the numerical fluid calculation for
$\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$ and
$t_{\mathrm{heat}}=5\sec$. The blue dashed curve indicates the enclosed mass
inside $r_{\rm turn}$ evaluated with $f_{0}=1$ in Equation (18), where the
values of $r_{\rm turn}$ and $\rho_{\rm turn}$ are adopted from the numerical
fluid calculations. The bottom panel shows the value of $f_{0}$ required to
reproduce the results of the fluid calculations.

In the following, we discuss the phase after $r_{\rm turn}$ becomes constant,
i.e., $t>6~{}t_{c}$. When the position of $r=r_{\rm turn}$ is constant, the
forces are balanced (see Equation (17)) and the velocity is $0$; thus the mass
does not fall back from the radius beyond $r_{\rm turn}$. Therefore, the
decreasing rate of the ejecta mass between a sufficiently small radius (taken
to be $r_{\rm count}=650\ {\rm km}$ in the numerical calculation) and $r_{\rm
turn}$, which has not been accreted yet to the central object, is equal to the
mass accretion rate to the central object. The ejecta mass contained within
$r=r_{\rm turn}$ can be written as

${M}_{\rm turn}=\frac{4\pi}{3}f_{0}r_{\text{turn }}^{3}\rho_{\text{turn }},$
(18)

where the subscript “turn” means the value at $r=r_{\rm turn}$ and $f_{0}$ is
an $\mathcal{O}(1)$ constant that corrects the difference originating from the
mass density distribution. In Figure 7, the enclosed mass within $r=r_{\rm
turn}$ evaluated with Equation (18) is compared to that of the numerical fluid
calculation (excluding that within $r_{\rm count}$) for $\dot{q}_{0}=3~{}{\rm
MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$ and $t_{\mathrm{heat}}=5\sec$. The
bottom panel of Figure 7 shows the value of $f_{0}$ in order for $M_{\rm
turn}$ to match the numerical value. Generally, there are three terms
contributing to the time derivative of $M_{\rm turn}$ as follows:

$\frac{\dot{M}_{\rm turn}}{M_{\rm
turn}}=\frac{\dot{f}_{0}}{f_{0}}+3\frac{\dot{r}_{\rm turn}}{r_{\rm
turn}}+\frac{\dot{\rho}_{\rm turn}}{\rho_{\rm turn}}.$ (19)

Evaluating the value of each term from the fluid calculations, we find that,
between $t\sim 6\,t_{c}$ and $t=t_{\rm heat}$, $r_{\rm turn}$ and $f_{0}$ vary
on longer timescales than the mass density $\rho_{\rm turn}$. We therefore at
first take $r_{\rm turn}$ and $f_{0}$ to be constants and assume that the time
variation of $M_{\rm turn}$ is due only to $\rho_{\rm turn}$.

Figure 8:  The infall velocity $|v|$ (blue line) and the sound speed (red
line) of the accretion flow for $\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm
nuc^{-1}}~{}{\rm s^{-1}}$ and $t_{\mathrm{heat}}=5\sec$ at $t=3\sec\sim
8.5\,t_{c}>6\,t_{c}$. The vertical and horizontal axes are normalized by the
characteristic scales in Equation (13) and Equation (12), respectively. The
black dotted and dashed curves represent the inflow velocity and the sound
velocity calculated from the Bondi accretion flow, respectively. The green
long-dashed curve represents the escape velocity at each radius. The vertical
and horizontal thin lines mark the sonic radius $r_{\rm sonic}$ and the sonic
velocity $v_{\mathrm{sonic}}$, respectively.  Figure 9:  The top panel shows
the time evolution of the mass accretion rate. The horizontal axis is
normalized by the characteristic timescale in Equation (12). The solid red
curve presents the mass accretion rate calculated using the mass flux at
$r=650\ {\rm km}\sim 0.18r_{c}$ in the numerical fluid calculation for
$\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$ and
$t_{\mathrm{heat}}=5\sec$. The blue dashed curve is evaluated by Equation (20)
with $f_{1}=1$. Here, $\rho_{\rm turn}$ and $a_{\rm turn}$ are the values
obtained from the numerical fluid calculation. The bottom panel shows the
values of $f_{1}$ required to reproduce the result of the numerical fluid
calculation.

The hydrodynamical structure at the time $t=8.5~{}t_{c}>6~{}t_{c}$, where the
turn-around radius $r_{\rm turn}$ has become constant, is shown in Figure 8.
In the region where $|v|\gtrsim|v_{c}|$, the flow and sound velocities can be
approximated well by the Bondi accretion flows. The accretion rate of the
Bondi solution can be written as follows using the physical quantities at
$r=r_{\rm turn}$ (Bondi, 1952):

$\dot{M}_{\rm Bondi}=2\sqrt{2}\pi f_{1}(GM)^{2}\rho_{\text{turn
}}/a_{\text{turn }}^{3},$ (20)

where $a=\sqrt{\Gamma P/\rho}$ is the sound speed, $a_{\rm turn}=a(r_{\rm
turn})$ and $f_{1}$ is an $\mathcal{O}(1)$ constant that corrects the
difference owing to the fact that $r_{\rm turn}$ is not infinite. Originally,
Equation (20) is a relation between the gas that is sufficiently distant and
stationary, but here we evaluate this value at $r_{\rm turn}$. The fluid
element at the turn-around radius $r_{\rm turn}$ is almost at rest (see Figure
6 and Section 3.1). Although the pressure gradient is comparable to the
gravity, the gravitational potential is rather smaller than the internal
energy of the fluid, so that it can be approximated in this way. In order to
check the validity of this assumption, in Figure 9, the mass accretion rate
obtained by the numerical fluid calculation is compared to that evaluated
using Equation (20). After $t=4\,t_{c}$, it can be seen that the mass
accretion rate is well approximated by Equation (20). Equating the mass
accretion rate in Equation (20) to that in Equation (19) (in the absolute
values), we obtain the equation for the time evolution of the mass density
$\rho$ at $r=r_{\rm turn}$ as follows:

$\frac{\dot{\rho}_{\rm turn}}{\rho_{\rm
turn}}=-\frac{3\sqrt{2}}{2}\frac{f_{1}}{f_{0}}\frac{1}{t_{c}}\left(\frac{a_{\rm
turn}}{v_{c}}\right)^{-3}.$ (21)

Here, $f_{0}$, $f_{1}$ and $r_{\rm turn}\sim r_{c}$ are assumed to be constant
and Equations (11)–(13) are used. In order to solve Equation (21), a model of
the time evolution of the sound speed $a_{\rm turn}$ at $r=r_{\rm turn}$ is
needed.

###  3.3 Sound speed at $r_{\rm turn}$

Figure 10:  Time evolution of the square of the sound speed at $r=r_{\rm
turn}$. The vertical and horizontal axes are normalized by the characteristic
scales in Equation (13) and Equation (12), respectively. Note that both the
axes are given in linear scales. The red solid and dashed curves show the
results of the numerical fluid calculations with $\dot{q}_{0}=3~{}{\rm
MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$ for $t_{\rm heat}=5\sec$ and $t_{\rm
heat}=10\sec$, respectively. The black dotted line, which is a linear function
of time, is a model curve of the sound speed given by Equation (22). It can be
seen that the model approximates the results of the fluid calculations well at
$t\gtrsim 6~{}t_{c}$.

In this section, we model the temporal evolution of the sound speed at the
turn-around radius $r_{\rm turn}$. Figure 10 shows the time evolution of the
square of $a_{\rm turn}$. As can be seen from the figure, the values of
$a_{\rm turn}^{2}$ in the numerical fluid calculations are proportional to
time after $t\gtrsim 6\,t_{c}$ up to $t=t_{\rm heat}$. This time dependence is
expected from the characteristic scale $\dot{q}_{0}t$, which has the dimension
of the square of the velocity. We assume the following equation as a model for
the temporal evolution of the sound speed at $r=r_{\rm turn}$:

$a_{\rm turn}=\sqrt{\frac{\Gamma P_{\rm turn}}{\rho_{\rm
turn}}}=\beta_{0}v_{c}\left(\frac{t}{t_{c}}\right)^{1/2},$ (22)

where the value of $\beta_{0}$ includes the effect of adiabatic cooling in the
accretion flow (that will be $\sqrt{\Gamma\left(\Gamma-1\right)}\sim 0.67$ in
the absence of adiabatic cooling). Note that, in this paper, we consider the
case where the total heating per nucleon is $\mathcal{O}({\rm MeV})$, so that
this value will never reach the speed of light. By solving Equation (22) for
$P_{\rm turn}$ and by using Equations (4), (12), and (13), we obtain the value
of the internal energy, $\epsilon_{\rm int}\sim 0.42\,\rho\dot{q}_{0}t$, which
indicates that about half of the added radioactive energy is converted to
internal energy and the rest to kinetic energy. The test-particle model of
Metzger et al. (2010) and Desai et al. (2019) assumes that all of the
radioactive energy is converted into kinetic energy; however as we have shown
here, such a 100% conversion efficiency from internal energy to kinetic energy
is not the case.

###  3.4 Halting of the fall-back accretion by $r$-process

The mass accretion rate $\dot{M}_{\rm turn}$ at $r=r_{\rm turn}$, can be
calculated if the time evolution of $\rho_{\rm turn}$ is obtained with
Equations (21) and (22). Using Equations (18), (21), and (22), we obtain the
following relation:

$\frac{d\ln\dot{M}_{\rm turn}}{d\ln
t}=-\frac{3}{2}\left[1+\sqrt{2}\frac{f_{1}}{f_{0}}\frac{1}{\beta_{0}^{3}}\left(\frac{t}{t_{c}}\right)^{-1/2}\right].$
(23)

By integrating this over time, we obtain

$\frac{\dot{M}_{\rm
turn}}{\dot{M}_{0}}=\left(\frac{t}{t_{0}}\right)^{-3/2}\exp\left[-\frac{3\sqrt{2}}{\beta_{0}^{3}}\frac{f_{1}}{f_{0}}\sqrt{\frac{t_{c}}{t_{0}}}\left(1-\sqrt{\frac{t_{0}}{t}}\right)\right],$
(24)

where $\dot{M}_{0}$ is the initial condition for $\dot{M}_{\rm turn}$ at a
given time $t_{0}$. For the evaluation of $f_{0}$, $f_{1}$, and $\beta_{0}$,
we adopt the values at $t_{0}=6\,t_{c}$, at which the turn-around radius
becomes approximately constant (see Figure 6), and use $f_{1}/f_{0}=0.5$ and
$\beta_{0}=0.43$. The resultant mass accretion rates using Equation (24) are
shown in Figure 4 by black curves. The dotted curves indicate those for
$t<6\,t_{c}$, i.e., the range of time when the assumptions necessary to derive
equation (24) are not valid. It is clear from Figure 4 that the mass accretion
rate can be well approximated by Equation (24) after $t=6\,t_{c}$, which
decreases rapidly along the theoretical curve. Although the model is
calibrated based on the result for $\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm
nuc^{-1}}~{}{\rm s^{-1}}$, it can be seen that the theoretical curve also
explains well those for $\dot{q}_{0}=30~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm
s^{-1}}$.

$\dot{M}_{0}$ is a free parameter that cannot be determined due to the lack of
a typical scale of ejecta mass in this system. For the cases in which only
$\dot{q}_{0}$ is different but the ejecta profile is the same, we can derive a
scaling law of $\dot{M}_{0}$ for various $\dot{q}_{0}$. For $t\ll 6\,t_{c}$,
the mass accretion rate exhibits almost the same temporal evolution as in
$\dot{q}_{0}=0$ case (see the gray line in Figure 4), so that we can
approximate the mass accretion rate as that without radioactive heating.
Because the ejecta mass follows the same dimensionless equation under the
normalized scales given by Equations (11)–(13), the ratio of mass accretion
rates between with and without heating, which is a dimensionless quantity, has
the same time evolution as a function of the normalized time $t/t_{c}$ for
various $\dot{q}_{0}$. If we choose the reference point $t_{0}$ in time units
of $t_{c}$ (as we chose the reference point of $\dot{M}_{0}$ as
$t_{0}=6\,t_{c}$), this ratio at $t=t_{0}$ is uniquely determined. Without
heating, the mass accretion rate has the time evolution proportional to
$t^{-5/3}$ so that the scaling law for $\dot{M}_{0}$ with respect to
$\dot{q}_{0}$ is given as follows:

$\dot{M}_{0}\propto t_{c}^{-5/3}\propto\left(GM\right)^{-2/3}\dot{q}_{0}.$
(25)

In the setting of this paper, at $t=6\,t_{c}$, the value of $\dot{M}_{0}$ is
about 13% of that without heating, which is derived from the numerical result
with $\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$. In the case
of $\dot{q}_{0}=30~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$ in Figure 4,
we adopt $\dot{M}_{0}$ calculated using Equation (25).

As represented by Equation (24) and Figure 4, the mass accretion rate is
suppressed compared to the case without radioactive heating. This can be
understood from the following two effects by considering the mass accretion
rate $\dot{M}\propto a_{\rm B}r_{\rm B}^{2}\rho_{\rm B}$ evaluated at the
sonic point $r=r_{\rm B}$. The first effect is that $r_{\rm turn}$ is nearly
constant over a certain period of time. As mentioned in Section 3.1, this
corresponds to the fact that no mass can flow into the inside of the sphere of
$r_{\rm turn}$ from the outside. Thus the mass density $\rho_{\rm turn}$ at
the turn-around radius will necessarily decrease with accretion. Furthermore,
the flow outside the sonic radius $r_{\rm B}$ is almost incompressible because
it is a subsonic flow, so that the mass density $\rho_{\rm B}$ at the sonic
point is comparable to $\rho_{\rm turn}$. Therefore, $\rho_{\rm B}$ is also a
decreasing function of time. The other effect is the increase in the speed of
sound with time. As can be seen in Figure 8, the speed of sound at the sonic
point of the accretion flow is about the same order of magnitude as that at
the turn-around radius. Since the speed of sound increases in proportion to
$t^{1/2}$ (see Equation (22)), the radius of the sonic point $r_{\rm
B}=GM/2a_{\rm B}^{2}$ (Bondi, 1952) shrinks as $r_{\rm B}\propto t^{-1}$.
Because the mass accretion rate is $\dot{M}\propto a_{\rm B}r_{\rm
B}^{2}\rho_{\rm B}\propto\rho_{\rm B}t^{-3/2}$, the above two effects both
work to reduce the accretion rate.

Chevalier (1989) calculated the accretion rate with radioactive heating,
mainly due to 56Ni, in the context of fallback accretion to the proto-neutron
star in a supernova explosion. He showed analytically that the fallback
accretion rate declines sharply as $\dot{M}\propto t^{-9/2}$ well within the
half-life of 56Ni, i.e., for a period of time when the heating rate per unit
mass is approximately constant. Although the solution we obtained is
exponential rather than the power-law of time (see Equation (24)), the result
in Chevalier (1989) is qualitatively similar to ours in terms of the
suppression of the mass accretion rate in the presence of heating. The
difference is that the model in Chevalier (1989) considered only a self-
similar expansion of ejecta and assumed the mass density decreasing as
$t^{-3}$. In fact, provided that $\rho_{B}\propto t^{-3}$, we obtain
$\dot{M}\propto\rho_{\rm B}t^{-3/2}\propto t^{-9/2}$ by using the same
argument as described earlier. The difference arises because we also consider
the reduction of mass density $\rho_{\rm B}$ due to accretion.

Figure 11:  The parameter regions for radioactive heating, in which the mass
accretion is sufficiently suppressed. The solid blue and red lines represent
the boundaries of the suppression for the case of a constant heating rate (see
Equation (26)) and the case of a broken-power-law heating rate as expressed in
Equation (27) ($q_{0}=3.0~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$,
$t_{\rm heat}=1.0~{}\sec$, and $p=1.4$), respectively, for a central mass of
$2.7M_{\odot}$. The thick black dash lines are shown to guide the dependence
of the halting condition $\dot{q}\propto t_{\rm heat}^{-5/3}$ (see Equations
(26) and (27)). In the upper-right regions above these lines, the mass
accretion is well suppressed by the effect of radioactive heating. The star-
shaped symbols indicate the locations for the parameters used in the numerical
calculations of this paper, for which the halting has been observed. The thin
blue dashed line shows the boundary for the “cutoff condition” by Desai et al.
(2019) (see Appendix for details). The thick red line is the locus of the
mass-averaged heating rate (excluding that in neutrinos) calculated by Wanajo
et al. (2014) (see Figure 2). The thick red dotted line is that for the model
in which the total radioactive energy per nucleon is $1~{}{\rm MeV}~{}{\rm
nuc^{-1}}$, the constant phase lasts for $1\sec$, and the heating rate
decreases proportionally to $t^{-1.4}$ in the decay phase.

If $t_{\rm heat}$ is shorter than $6~{}t_{c}$, the stagnation of $r_{\rm
turn}$ and the increase of the sound speed will not occur for the constant
heating model. Thus, $t\gtrsim 6~{}t_{c}$ is the necessary condition for these
two processes to work, namely,

$t_{\rm heat}>Kt_{c}\\\ \sim 2.1\sec~{}K_{\rm
6}\left(\frac{M}{2.7~{}M_{\odot}}\right)^{2/5}\left(\frac{\dot{q}_{0}}{3~{}{\rm
MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}}\right)^{-3/5},$ (26)

where $K$ is the time (with respect to $t_{c}$) that takes for $r_{\rm turn}$
to become constant and $K=6\,K_{6}$ (see Figure 6). After $t\sim 6~{}t_{c}$,
the analytic solution (24) becomes valid, and the mass accretion rate rapidly
decreases. Therefore, we call this time, $t_{\rm halt}\sim 6~{}t_{c}$,
“halting time” and this suppression of mass accretion “halting”. The halting
time corresponds to the timescale in which the mass accretion rate decreases
to $\sim 13$% of its original value (see the description below Equation (25)).
For given $\dot{q}_{0}$, $t_{\rm heat}$, and $M$, the halting condition is
determined by Equation (26). Figure 11 shows whether or not the halting occurs
in the $\dot{q}_{0}t_{\rm heat}$–$t_{\rm heat}$ space. Here, because we are
considering only the constant heating phase (see Equation (9), and see also
Equation (27) for the case with a power-law decay phase), the vertical axis in
Figure 11 indicates the total radioactive energy injected into the fluid. As
can be seen from Figure 11, even if the total radioactive energy is the same,
the halting is more likely to occur as $t_{\rm heat}$ increases. This is
because the ejecta which is accreted at the later phase has less binding
energy, and therefore, the accretion is more easily disturbed by the later
injection of radioactive energy.

Desai et al. (2019) also argued that the mass accretion stops at a finite time
if there is sufficient heating. We compare our model with that of Desai et al.
(2019) as shown by the thin dashed lines in Figure 11. The method for
calculating their theoretical curve is summarized in Appendix A. We consider
that the “halting time” we obtained corresponds to the “cutoff” of the mass
accretion claimed by Desai et al. (2019). As can be seen from the figure,
their results and ours have the same dependence on the variables, but the
total heating energy per nucleon $\dot{q}t_{\mathrm{heat}}$ to halt the
fallback accretion differ by about an order of magnitude for the constant
heating case. In other words, compared to their test-particle model, our fluid
model requires about $10$ times larger heating rate to halt the mass accretion
for given $t_{\rm heat}$. This difference is likely due to the fact that in
the test-particle model, all of the radioactive energy is converted into the
kinetic energy, whereas in the fluid model, this energy conversion is
incomplete remaining internal energy.

##  4 Semi-analytical model with power-law decaying heating rate

As we find in Figure 2, there is actually not only a constant phase in the
heating profile but also a decay phase. The heating profile can be reasonably
approximated by

$\dot{q}(t)=\left\\{\begin{array}[]{ll}\dot{q_{0}}&(t\leq t_{\rm heat,0})\\\
\dot{q_{0}}\left(\frac{t}{t_{\rm heat,0}}\right)^{-p}&(t>t_{\rm
heat,0})\end{array}\right.,$ (27)

where $t_{\rm heat,0}$ is the duration of the constant phase and $p>1$. The
total radioactive energy $Q_{\rm tot}$ can be written as

$Q_{\rm tot}=\frac{p}{p-1}\dot{q_{0}}t_{\rm heat,0}.$ (28)

If the halting occurs within the constant phase of radioactive heating in
Equations (27), that is $t_{\rm halt}\leq t_{\rm heat,0}$, the halting time
$t_{\rm halt}=Kt_{c}$ (cf. Equation (26)) can be written as follows:

$t_{\mathrm{halt}}=\left[K^{5}(GM)^{2}\left(\frac{p}{p-1}\right)^{3}Q_{\rm
tot}^{-3}t_{\mathrm{heat},0}^{3}\right]^{\frac{1}{5}}\\\ \sim 265\sec
K_{6}\left(\frac{M}{2.7M_{\odot}}\right)^{2/5}\left(\frac{Q_{\rm
tot}}{1~{}{\rm MeV}~{}{\rm nuc^{-1}}}\right)^{-3/5}\\\
\times\left(\frac{t_{\mathrm{heat},0}}{300\sec}\right)^{3/5}\ \ (p=1.4).$ (29)

Even if the mass accretion does not halt during the constant phase, the
halting will eventually occur. Here, we construct a semi-analytical model of
the accretion flow with the heating that decays with the power-law of time. In
the decay phase, the typical scales given by Equations (11) and (12) become
time-varying as follows:

$\displaystyle
r_{c}(t)=\left[\frac{(GM)^{3}}{\dot{q}(t)^{2}}\right]^{1/5}=\left[\frac{(GM)^{3}}{\dot{q}_{0}^{2}}\right]^{1/5}\left(\frac{t}{t_{\text{heat,0}}}\right)^{2p/5},$
(30)
$t_{c}(t)=\left[\frac{(GM)^{2}}{\dot{q}(t)^{3}}\right]^{1/5}=\left[\frac{(GM)^{2}}{\dot{q}_{0}^{3}}\right]^{1/5}\left(\frac{t}{t_{\text{heat,0}}}\right)^{3p/5}.$
(31)

Introducing the dimension-less variables $\xi=r/r_{c}(t)$ and
$\chi=t/t_{c}(t)$, the fluid equations (1)–(3) can be normalized as follows:

$\left(1-\frac{3p}{5}\right)\frac{\partial\ln\phi}{\partial\ln\chi}+\left(V-\frac{2p}{5}\right)\frac{\partial\ln\phi}{\partial\ln\xi}\\\
+\frac{\partial V}{\partial\ln\xi}+3V=0,$ (32)
$\left(1-\frac{3p}{5}\right)\frac{\partial
V}{\partial\ln\chi}+\left(V-\frac{2p}{5}\right)\frac{\partial
V}{\partial\ln\xi}\\\ +\frac{\partial
Z}{\partial\ln\xi}+Z\frac{\partial\ln\phi}{\partial\ln\xi}+2Z+V(V-1)=-\frac{\chi^{2}}{\xi^{3}},$
(33) $\left(V-\frac{2p}{5}\right)\frac{\partial
Z}{\partial\ln\xi}+\left(1-\frac{3p}{5}\right)\frac{\partial
Z}{\partial\ln\chi}\\\ +(\Gamma-1)Z\frac{\partial
V}{\partial\ln\xi}+[(3\Gamma-1)V-2]Z=(\Gamma-1)\frac{\chi^{3}}{\xi^{2}},$ (34)

where $\phi=\rho/\rho_{c}$ is the normalized mass density, $\rho_{c}$ is
arbitrary constant with dimension of mass density, $V=vt/r$ is the normalized
radial velocity, and $Z=(P/\rho)(t/r)^{2}$ is the normalized pressure. Note
that the normalized equations (32)–(34) do not explicitly include $t_{\rm
heat,0}$. As we will see below, $t_{\rm heat,0}$ is relevant to the evolution
of ejecta only as an initial condition in Equations (32)–(34).

The temporal evolution of the accretion flow under the heating rate (27) is as
follows. Up to time $t_{\rm heat,0}$, as seen in Section 3, the accretion flow
evolves according to Equations (14)–(16) with the normalization of Equations
(11) and (12). If $t_{\rm heat,0}$ is longer than $\sim 6~{}t_{c}$ (see
Section 3.4), the halting occurs during the constant phase, and the halting
time is expressed by Equation (29). On the other hand, if $t_{\rm
heat,0}\lesssim 6~{}t_{c}$, the ejecta evolves according to Equations
(14)–(16) and (32)–(34), respectively, before and after

$\chi_{\rm heat}=t_{\rm heat,0}/t_{c}(t_{\rm heat,0}),$ (35)

at which the heating rate switches from the constant phase to the power-law
decaying phase. As seen in Section 3, during the constant phase, since the
temporal evolution of the ejecta with various parameters are identical under
normalized variables, almost independently of the initial conditions.
Therefore, even if model parameters are different but $\chi_{\rm heat}$ and
$p$ are the same, these accretion flows follows the same temporal evolution
according to Equations (32)–(34) from the same initial conditions (i.e., the
states at the end of the constant phase) with normalized variables.

In the case of the halting in the decay phase, unlike the case in the constant
phase, the normalized halting time $K$ takes a different value from $6$ (see
Section 3), where $K$ can be written as

$K=t_{\rm halt}/t_{c}(t_{\rm halt}).$ (36)

As seen in Section 3.4, the halting time $t_{\rm halt}$ is the time it takes
for the mass accretion rate with the heating to be suppressed to about $13\%$
of that without heating. As mentioned earlier, the normalized halting time $K$
is basically a function of $\chi_{\rm heat}$ and $p$ only. Also, if $t_{\rm
heat,0}$ (or $\chi_{\rm heat}$) is long enough (i.e., $\chi_{\rm heat}>6$),
this will match the constant model, and $K=6$. We investigate the dependence
of $K$ on $\chi_{\rm heat}$ and $p$ based on the numerical calculation. As a
heating rate profile, we take various values of $t_{\rm heat,0}$ and $p$ in
Equation (27). We choose $\dot{q}_{0}$ and $t_{\rm heat,0}$ for $p=1.2$,
$1.3$, and $1.4$, respectively, with the condition
$\dot{q}_{0}=2.0\left(t_{\rm heat,0}/1.0~{}\sec\right)^{-p}~{}{\rm MeV}~{}{\rm
nuc^{-1}}~{}{\rm s^{-1}}$ for various $t_{\rm heat,0}$. The ejecta profiles
(see Figure 1) and the calculation method are the same as in Section 2.

Figure 12 shows the dependence of the normalized halting time $K$ on the
normalized duration $\chi_{\rm heat}$ of the constant phase. The red, blue,
and green dots represent the results of the numerical calculation for $p=1.2$,
$1.3$, and $1.4$, respectively. We fit these results with the following
monotonically increasing functional form:

$K^{\alpha}=K_{0}^{\alpha}+\left[1-\left(\frac{K_{0}}{6}\right)^{\alpha}\right]\chi_{\rm
heat}^{\alpha},$ (37)

where

$K_{0}=A_{K}p+B_{K},~{}~{}~{}\alpha=C_{K}p+D_{K}.$ (38)

Here we fix the value of $K=6$ for $\chi_{\rm heat}=6$ to recover the result
of the model with the constant heating exactly. The resultant parameters are
$A_{K}=-1.68$, $B_{K}=4.75$, $C_{K}=-3.14$, and $D_{K}=-1.80$. The resultant
fitting functions are shown as solid curves in Figure 12. Using the
radioactive heating rate in the BNS merger as shown in Figure 2, we find that
for that heating rate, $K=2.6$ (with $\chi_{\rm heat}=0$ and $p=1.3$) is
appropriate for almost all $Y_{\rm e}$ cases. Once we obtain the value of $K$,
we can calculate the halting time in the decay phase by solving the equation,

$t_{\rm halt}=K\left[\frac{\left(GM\right)^{2}}{\dot{q}(t_{\rm
halt})^{3}}\right]^{1/5},$ (39)

and the halting time is obtained as follows,

$t_{\mathrm{halt}}=\left[K^{5}(GM)^{2}\left(\frac{p}{p-1}\right)^{3}Q_{\rm
tot}^{-3}t_{\mathrm{heat},0}^{3(1-p)}\right]^{\frac{1}{5-3p}}\\\ \sim
3.9\times 10^{3}\sec
K_{2.6}^{6.25}\left(\frac{M}{2.7M_{\odot}}\right)^{2.5}\left(\frac{Q_{\rm
tot}}{1~{}{\rm MeV}~{}{\rm nuc^{-1}}}\right)^{-3.75}\\\
\times\left(\frac{t_{\mathrm{heat},0}}{1\sec}\right)^{-1.5}\ \ (p=1.4).$ (40)

Here we adopt $K=2.6$ corresponding to $\chi_{\rm heat}\sim 0$, which is a
good approximation for the case of realistic heating rates shown in Figure 2.
This formula holds if the halting does not occur during the constant phase,
and thus the above expression is only valid if $t_{\rm halt}$ is longer than
$t_{\rm heat,0}$. Note that, as can be seen from Figure 11, the power-law
index of time in the decay phase of the heating rate $\dot{q}(t)$ must be
shallower than $t^{-5/3}$ for halting to occur.

Figure 12:  Dependence of the normalized halting time $K$ on the normalized
duration $\chi_{\rm heat}$ of the constant phase and the decay index $p$ of
the heating profile. Here we adopt the mass of the central object
$M=2.7~{}M_{\odot}$, and the parameters of the heating rate with the condition
$\dot{q}_{0}=2.0\left(t_{\rm heat,0}/1.0~{}\sec\right)^{-p}~{}{\rm MeV}~{}{\rm
nuc^{-1}}~{}{\rm s^{-1}}$ for various $t_{\rm heat,0}$. The red, blue, and
green points show the results of numerical calculations with $p=1.2$, $1.3$,
and $1.4$, respectively. The solid curves show the fitting functions in
Equations (37) and (38) for various $p$. The black dash line represents the
line of $K=\chi_{\rm heat}$. For $\chi_{\rm heat}=6$, the condition $K=6$
shown in Section 3 is recovered because the mass accretion halts within the
constant phase.

Let us evaluate the mass accretion rate in the decay phase of the heating
$\dot{q}(t)\propto t^{-p}$ at $t>t_{\rm halt}$. Here, we assume that $p<5/3$
so that the halting occurs during the decay phase. From Equation (30), the
turn-around radius $r_{\rm turn}$ is expected to evolve in time as $r_{\rm
turn}\propto\dot{q}(t)^{-2/5}\propto t^{2p/5}$. Furthermore, from Equation
(22), the sound speed at turn-around radius is expected to be $a_{\rm
turn}\propto(\dot{q}(t)t)^{1/2}\propto t^{(1-p)/2}$. We have checked that
these time dependencies are approximately consistent with those obtained by
numerical calculations, except for the constant factor (e.g., $r_{\rm
turn}/r_{\rm c}(t)\sim 0.8$ rather than unity). By solving Equation (18) for
$\rho_{\rm turn}$ and substituting it into Equation (20), we obtain the
differential equation for the mass $M_{\rm turn}$ within the turn-around
radius,

$\frac{d\ln M_{\rm turn}}{dt}=\frac{3\sqrt{2}}{2}f_{\rm
PL}\frac{f_{1}}{f_{0}}\frac{\left(GM\right)^{2}}{a_{\rm turn}^{3}r_{\rm
turn}^{3}},$ (41)

where $r_{\rm turn}$ depends on time. We find that $f_{0}$, $f_{1}$,
$\beta_{0}$, and $r_{\rm turn}/r_{c}$ in the power-law heating model slightly
deviate from those in the constant heating model shown in Section 3, and
depend on the normalized duration of the constant phase $\chi_{\rm heat}$ and
the decay index of the heating rate $p$. Here, we introduce the time-constant
factor $f_{\rm PL}=f_{\rm PL}(\chi_{\rm heat},p)$ to adjust these differences.
By integrating this equation over time and differentiating the obtained
solution with $t$, we get

$\frac{\dot{M}_{\rm turn}}{\dot{M}_{0}}=\left(\frac{t}{t_{\rm
halt}}\right)^{(3p-15)/10}\\\ \times\exp\left[-\frac{5}{5-3p}f_{\rm
PL}\frac{3\sqrt{2}}{\beta_{0}^{3}}\frac{f_{1}}{f_{0}}\frac{1}{\sqrt{K}}\left(1-\left(\frac{t_{\rm
halt}}{t}\right)^{\frac{5-3p}{10}}\right)\right].$ (42)

It can be easily checked that Equation (24) is recovered when $p=0$ and
$f_{\rm PL}=1$. From the calculation results used for determining the
$\chi_{\rm heat}$ dependence of $K$, we can also obtain the fitting function
of $f_{\rm PL}$. We fit the calculation results of $f_{\rm PL}$ with the
following functional form:

$f_{\rm PL}=\left(S_{\rm PL}p+T_{\rm PL}\right)\chi_{\rm
heat}^{2}+\left(U_{\rm PL}p+V_{\rm PL}\right).$ (43)

The resultant parameters are $S_{\rm PL}=-1.22\times 10^{-2}$, $T_{\rm
PL}=2.65\times 10^{-2}$, $U_{\rm PL}=-0.324$, and $V_{\rm PL}=0.661$.

Figure 13:  Mass accretion rate calculated using the heating rate profile of
broken power-law (Equation (27)) with $q_{0}=3~{}{\rm MeV}~{}{\rm
nuc^{-1}}~{}{\rm s^{-1}}$, $t_{\rm heat}=1.0~{}\sec$, and $p=1.4$ (red curve).
The gray curve shows the mass accretion rate without heating. The black curve
is the mass accretion rate calculated by the semi-analytical model of Equation
(42). The dotted line is before the halting time, and the Equation (42) is not
valid.

In Figure 13, the numerical calculation result with the parameters
$q_{0}=3~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$, $t_{\rm
heat,0}=1.0~{}\sec$, and $p=1.4$ (shown in the red curve) is compared with the
mass accretion rate calculated from Equation (42) (shown in the black curve).
Note that the numerical results shown in Figure 13 are for parameters
independent of the calculations used to calibrate Equations (38) and (43). As
can be seen from the figure, the semi-analytical model reproduces the
numerical results well. Furthermore, generally, with the parameters used in
our calculations (i.e., $f_{\rm PL}\sim 0.3$, $f_{1}/f_{0}\sim 0.5$ and
$\beta_{0}\sim 0.43$), the mass accretion rate in Equation (42) decreases more
sharply than $t^{-5/3}$ at $t=t_{\rm halt}$, so that the mass accretion is
expected to be strongly suppressed after the halting time even in the decay
phase.

##  5 Application to BNS mergers

In this section, we consider the halting process with the realistic heating
rate in BNS mergers as shown in Figure 2. We can estimate the halting time
from the intersection of the heating rate curve and the line above which
halting occurs (red solid line in Figure 11 for $M=2.7M_{\odot}$ and $K\sim
2.6$). As shown by the thick red curve in Figure 11, when using the mass-
averaged heating rate of the nucleosynthesis calculations in Figure 2, we can
see that the halting occurs by the heating only after a time $\sim 3.7\times
10^{5}$ sec.

Figure 14:  Dependence of the halting time on the radioactive heating for the
fallback accretion with the central object mass $M=2.7M_{\odot}$. The circles
indicate the halting times calculated by using the nucleosynthesis results for
individual initial $Y_{\mathrm{e}}$ values in Wanajo et al. (2014). The star
represents the results for the mass-averaged heating rate (see Figure 2). The
red and green lines represent the results of the calculations with $p=1.3$ and
$p=1.4$ power-law index in the decay phase of radioactive heating (see
Equation (27)), respectively, for $t_{\rm heat}=0.1$ sec (solid) and $t_{\rm
heat}=1.0$ sec (dotted). Note that, in all cases, the halting occurs in the
decay phase. Shaded areas are the ranges of typical observed timescales for
the extended emission and the plateau emission of sGRBs. We can see that the
accretion ceases, long after these emissions decay.

Figure 14 shows the dependence of the halting time on the parameters of the
heating, i.e., $t_{\rm heat}$ and $Q_{\rm tot}$. In addition, in order to
investigate the dependence on $Y_{\rm e}$, we also show the halting times
calculated using the heating rates in Figure 2. Reflecting the tendency of
larger $Q_{\rm tot}$ for smaller $Y_{\rm e}$, the halting time becomes shorter
for smaller $Y_{\rm e}\lesssim 0.3$. For larger $Y_{\rm e}$, the dependence of
the halting time on $Y_{\rm e}$ is not monotonic. This is because as $Y_{\rm
e}$ increases, the length of $t_{\rm heat}$ increases. The longer $t_{\rm
heat}$ leads to the greater heating rate in the decay phase, so that the
halting time becomes shorter. We find that the halting by $r$-process does not
occur within a timescale of $\sim 10^{4}$ sec, although the radioactive
heating may decelerate the accretion flow to some extent.

The radioactive heating of ejecta is the thermalization process of non-thermal
particles due to ionization losses of charged particles and repeated
scattering and absorption of gamma-rays. It has been pointed out that, after
about some of ten days, the timescale of thermalization becomes longer than
the expansion timescale of ejecta and thus the thermalization becomes
inefficient (Hotokezaka et al., 2016; Barnes et al., 2016; Waxman et al.,
2018; Kasen & Barnes, 2019; Hotokezaka & Nakar, 2020). Because the halting
time we obtained is 10 days or even longer, we expect that the thermalization
is actually insufficient. However, all the charged particles associated with
radioactive decay, which are efficiently trapped inside the ejecta by the
magnetic field, contribute to pressure being independent of their
thermalization. Since the halting is essentially due to an increase in
pressure rather than an increase in temperature, the thermalization efficiency
is less important for charged particles. On the other hand, gamma-rays (about
a half of $\beta$-decay energy is emitted in the form of gamma-rays) can
escape from the ejecta. Thus, gamma-rays make little contribution to pressure
after a sufficient time and the heating rate becomes effectively small. As can
be seen from Equation (40), when the heating rate becomes about twice as
small, the halting time becomes about an order of magnitude longer. Therefore,
the halting time shown in Figure 14 should be considered as the lower limit.
Note that, if the thermalization is insufficient, the radiation efficiency is
small, so that the radiative cooling is negligible.

The uncertainty in nuclear physics may also affect the estimated value of the
halting time. According to Barnes et al. (2020), there is a systematic
variation of about one order of magnitude in the radioactive heating rate at
$\mathcal{O}(1)$–$\mathcal{O}(10)$ days depending on the adopted nuclear
ingredients. Considering this uncertainty in the heating rate and the
inefficient thermalization, the halting time for a mass-averaged heating rate
is around $\mathcal{O}(10^{4})$–$\mathcal{O}(10^{8})\sec$. Note that,
according to Zhu et al. (2021), the uncertainty in the heating rate after
$\mathcal{O}(10)~{}{\rm days}$ become larger, and they suggested the
uncertainty range is about three orders of magnitude at
$\mathcal{O}(10^{8})~{}\sec$, allowing for a possibility of the halting time
larger than $\mathcal{O}(10^{8})~{}{\rm sec}$.

Typical observed timescales for the extended emission and the plateau emission
are about $\mathcal{O}(10^{3})\sec$ and $\mathcal{O}(10^{4})\sec$,
respectively (Rowlinson et al., 2013; Kagawa et al., 2019). As can be seen in
Figure 14, the estimated halting time is much longer than these timescales,
when using the nucleosynthesis results (Wanajo et al. (2014); see Figure 2).
If the energy source for these emissions were purely the fallback accretion in
BNS mergers (Rosswog, 2007; Rossi & Begelman, 2009; Kisaka et al., 2017), the
radioactive heating from decaying $r$-process nuclei appears insufficient to
disturb the energy supply from the accretion flow. Some additional heating
sources which inhibit accretion or other mechanisms, for example, the time-
varying radiation efficiency of the accreting matter (e.g., Kawanaka et al.,
2013), may be required to explain the characteristic timescales of the
extended and plateau emissions.

##  6 Conclusions

The discovery of GRB 170817A associated with GW170817 established that a BNS
merger is a source of sGRBs. However, the origin of the late-time emission in
sGRBs, namely, extended emission and plateau emission, is still unknown. We
have investigated the fallback accretion model of these long-lasting emission
of sGRBs (e.g., Rosswog, 2007; Lee & Ramirez-Ruiz, 2007; Rossi & Begelman,
2009; Kisaka & Ioka, 2015; Kisaka et al., 2017). While the canonical fallback
accretion rate of $t^{-5/3}$ has no typical timescales, Metzger et al. (2010)
and Desai et al. (2019) discussed that the effect of radioactive heating
results in a timescale of $\mathcal{O}(10)$–$\mathcal{O}(100)\sec$ for mass
accretion by using a test-particle model. We have revisited the effect of the
radioactive heating due to decaying $r$-process nuclei on the fallback
accretion by using a hydrodynamic model rather than a test-particle model. We
have shown that the timescale for the suppression of mass accretion becomes an
order of magnitude longer than that in the test-particle model. Furthermore,
we have found no temporal gap (i.e., halt and revival) of mass accretion,
being opposed to the results by Metzger et al. (2010) and Desai et al. (2019).
Their model assumes that all of the radioactive energy is promptly converted
to the kinetic energy. In addition, they also assume that a fluid element does
not fall back once it becomes unbound. However, our fluid calculations have
revealed that these assumptions are inappropriate.

We have developed a semi-analytical model for the temporal evolution of mass
accretion (see Equation (24) and (42)), which reproduces the numerical
results. The fallback accretion has characteristic length and time scales that
depend on the mass of the central object and the radioactive heating rate (see
Equations (11) and (12)). Normalizing the hydrodynamical equations with these
scales, we have obtained the scale-free equations (see Equations (14)–(16)).
Semi-analytical modeling of these normalized equations has allowed us to
investigate a wide parameter range of accretion flow. While the radioactive
heating with a constant rate continues, the radius at which the accretion flow
stagnates (turn-around radius) becomes nearly a constant value, being
approximately equal to the characteristic length scale. We have found that the
accretion flow inside the turn-around radius can be well approximated by the
Bondi accretion flow, and that the mass accretion rate is well reproduced by
the Bondi accretion rate evaluated at the turn-around radius. Furthermore, we
have derived the conditions on the heating rate for the substantial
suppression of mass accretion (see Equation (26)). We have extended this
condition to more general heating profiles that decay with time (see Figure
11). For the case where the heating rate can be written as a combination of
constant and decay phases (see Equation (27)), we have found that as long as
the heating rate decays more slowly than $t^{-5/3}$ in the decay phase, the
halting will occur after a sufficient amount of time (see Equation (40)). For
typical BNS mergers (GRB 170817-like events), the halting timescale for the
suppression of mass accretion is found to be
$\mathcal{O}(10^{4})$–$\mathcal{O}(10^{8})$ sec, which is, however, much
longer than the timescales in the late-time activity of sGRBs
$\mathcal{O}(10^{2})$–$\mathcal{O}(10^{4})$ sec (Rowlinson et al., 2013;
Kisaka et al., 2017; Kagawa et al., 2019). The observations of
macronovae/kilonovae associated with sGRBs suggest that the amount and
distribution of $r$-process product differ from event to event (Gompertz et
al., 2018; Ascenzi et al., 2019). For events such as GW170817, where the
macronova/kilonova light curve can be observed in detail, the halting time
will be determined by modeling the heating rate in the same way. Even if the
detailed light curves can not be obtained, by estimating the abundance
distribution (such as the lanthanide fraction) from color evolution, we may be
able to obtain the halting time from $Y_{\rm e}$ as seen in Figure 14.
Besides, our model will be applicable not only to BNS mergers but also to the
fallback accretion of proto-neutron stars in supernova explosions with
r-process nucleosynthesis (e.g., Nishimura et al., 2006, 2015; Mösta et al.,
2015).

Our results imply the existence of different mechanisms or different sources
of heating, which can stop the late-time activity of sGRBs. For example, the
shock heating by the interaction between the viscously driven wind and the
accretion flow may occur. In order to examine such a mechanism, multi-
dimensional hydrodynamical simulations with the effects of radioactive heating
will be necessary, in which both outflow (such as late-time viscously driven
wind) and inflow (such as the fallback accretion of early dynamical ejecta)
exist (see Kawaguchi et al., 2020, for a recent development). Note that, for a
system with only inflow, as in our calculation, the multidimensionality has a
minor effect. In the time evolution of the mass accretion rate, the ejecta
profile near the boundary between the gravitationally bound and unbound states
is essential. As seen in Figure 1, the radial dependence of velocity and
density around the radius $r\sim 490\ {\rm km}$ (the boundary of the bound
ejecta) is independent of latitude, which justifies the calculation with the
spherically symmetric model. Alternatively, magnetic reconnection or other
magnetic field dissipation processes may play a role in heating in the ejecta.
Instead of invoking other heating sources, the time scales of extended and
plateau emission may be explained by considering a mechanism in which the
conversion rate from gravitational energy to radiation decreases rapidly. For
instance, there may be a rapid change in radiation efficiency due to the state
transition of the accretion disk as the accretion rate decreases over time
(e.g., Kawanaka et al., 2013). We leave these issues for our future work.

The halting time is sensitive to the uncertainty of the radioactive heating
rate in the $r$-process elements, which ranges from $10^{4}$–$10^{8}$ sec for
one order of magnitude ambiguity in the heating rate (Barnes et al., 2020).
Furthermore, it has been suggested that the uncertainty becomes larger in the
later stages ($\mathcal{O}(1)$–$\mathcal{O}(10)\ {\rm yr}$) (Zhu et al.,
2021). This indicates that, if we can obtain the halting time for a
macronova/kilonova event, we may be able to constrain the physical conditions
for the $r$-process as well as the relevant nuclear ingredients. One possible
observational sign is the X-ray excess in the yearly-scale light curve of
GW170817 (e.g., Hajela et al., 2019; Balasubramanian et al., 2021; Hajela et
al., 2021), which we are currently investigating (Ishizaki et al., 2021).

## Acknowledgments

We thank the anonymous referee for fruitful comments. We are grateful to
Masaru Shibata, Kazuya Takahashi, Hamidani Hamid, Tomoki Wada, Koutarou
Kyutoku, Sho Fujibayashi, Kyohei Kawaguchi, Hiroki Nagakura, Shota Kisaka,
Kazumi Kashiyama, and Shuta Tanaka for fruitful discussion and valuable
comments. We thank the participants and the organizers of the workshops with
the identification number YITP-T-19-04, YKIS2019 and YITP-T-20-19 for their
generous support and helpful comments. This work is supported by Grants-in-Aid
for Scientific Research No. 21J01450 (WI), 18H01213 (KK), 20H01901, 20H01904,
20H00158, 18H01213, 18H01215, 17H06357, 17H06362, 17H06131 (KI) from the
Ministry of Education, Culture, Sports, Science and Technology (MEXT) of
Japan.

##  Appendix A Test-particle model

###  A.1 Cold case

First, we describe the fallback accretion in a system with negligible
pressure. The velocity of a fluid particle evolves over time according to the
equation of motion:

$\frac{dv}{dt}=-\frac{GM}{r^{2}},$ (A1)

where $M$ is the mass of the central object. The first integral of Equation
(A1) gives the dynamical energy per mass of the fluid element, which is
written as

$E_{0}=\frac{1}{2}v^{2}-\frac{GM}{r}.$ (A2)

Let us introduce dimensionless variables $x=r/r_{s}$ and $\beta=v/c$, where
$r_{s}=2GM/c^{2}$. The dimensionless time is also defined as $\tau=t/t_{s}$,
where $t_{s}=r_{s}/c$. The dimensionless energy per mass $\lambda$ is
determined by

$\lambda\equiv\frac{1}{x_{0}}-\beta_{0}^{2}=\frac{1}{x}-\beta^{2}=-2E_{0}/c^{2},$
(A3)

where the subscript $0$ indicates the values in the initial state.

The turn-around time, i.e., the time it takes for the fluid particle to change
the direction of motion, can be written as:

$\tau_{f}=\frac{\beta_{0}}{\lambda\left(\beta_{0}^{2}+\lambda\right)}+\tan^{-1}\left(\frac{\beta_{0}}{\sqrt{\lambda}}\right)\lambda^{-3/2}.$
(A4)

Note that $\tau_{f}$ is a function of the initial velocity and energy. When
$\lambda\ll\beta_{0}^{2}<1$, $\tau_{f}$ can be written as:

$\tau_{f}\sim\frac{\pi}{2}\lambda^{-3/2}.$ (A5)

As can be seen from Equation (A5), the turn-around time of a marginally bound
fluid particle ($\lambda\sim 0$) hardly depends on the initial velocity.

Let us find the mass per unit time, which falls back through the sphere of
$r=r_{\rm fin}$. Since the dynamical energy $E_{0}$ is conserved, the time it
takes to return to $r_{\rm fin}$ from the point at which $v=0$ coincides with
$\tau_{f}\left(r_{\rm fin},E_{0}(r_{0})\right)$. Therefore, the time $t_{\rm
fb}$ required for the fluid particle launched at a velocity $v_{0}$ from a
radius $r_{0}$ to fall back to $r_{\rm fin}$ can be written as follows:

$t_{\rm
fb}(r_{0},E_{0}(r_{0}))=t_{s}\left[\tau_{f}(r_{0},E_{0}(r_{0}))+\tau_{f}\left(r_{\rm
fin},E_{0}(r_{0})\right)\right].$ (A6)

Once the fallback time is determined as a function of $r_{0}$, the mass
accretion rate $\dot{M}$ is calculated as follows:

$\dot{M}\left(t_{\mathrm{fb}}\left(r_{0},E_{0}\left(r_{0}\right)\right)\right)\equiv\frac{dM}{dt_{\mathrm{fb}}}=4\pi
r_{0}^{2}\rho_{0}\left(r_{0}\right)\left(\frac{dt_{\mathrm{fb}}}{dr_{0}}\right)_{r=r_{0}}^{-1},$
(A7)

where $\rho(r_{0})$ is the mass density in the initial state.

###  A.2 Test-particle model for the $r$-process halting

According to Metzger et al. (2010) and Desai et al. (2019), we calculate the
fallback time with radioactive heating per unit time, $\dot{q}$. Assuming that
all the radioactive energy is converted to the kinetic energy, the total
energy of the particle at the turn-around time can be estimated as follows:

$E_{f}(r_{0})=E_{0}(r_{0})+\int_{t_{\text{start}}}^{t_{s}\tau_{f}\left(r_{0},E_{f}\right)}\dot{q}(t)\,dt,$
(A8)

where the subscript $f$ represents the values in the final state, i.e., the
values at the turn-around radius. Here, in order to introduce the effect that
the turn-around time becomes longer as the energy of the particle increases,
the value of $\tau_{f}$ at the upper end of the integration is evaluated by
using $E_{f}$. In fact, the internal energy injected to the fluid element is
converted into the kinetic energy via the pressure gradient forces. Since it
is difficult to deal with this process in the test-particle model, we adopt
Equation (A8), which is the same prescription in Metzger et al. (2010) and
Desai et al. (2019). Furthermore, according to Desai et al. (2019), after the
turn-around time (or, equivalently, fluid particles with $v<0$), we neglect
the effect of radioactive heating. Therefore, the fallback time is written as:

$t_{\text{fb}}\left(r_{0}\right)=t_{s}\left[\tau_{f}\left(r_{0},E_{f}\left(r_{0}\right)\right)+\tau_{f}\left(r_{\text{fin}},E_{f}\left(r_{0}\right)\right)\right].$
(A9)

Using Equations (A7) and (A9), the mass accretion rate $\dot{M}$ when
radioactive heating is effective can be obtained.

###  A.3 Halting condition for the test-particle model

Figure 15:  Left (blue line) and right (red curves) hand sides (with negative
signs) of the algebraic Equation (A10) as a function of $-E_{f}$ in the test-
particle model. The parameters $\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm
nuc^{-1}}~{}{\rm s^{-1}}$ and $t_{\rm heat}=0.6\sec$ are adopted, which
correspond to a case such that mass accretion stops and resumes, i.e., making
a “gap” (for details see text below Equation (A17) or Metzger et al., 2010).
The vertical dotted line indicates the value of $E_{f}$ such that the turn-
around time becomes equal to $t_{\rm heat}$. The different types of red curves
represent differences in the initial positions of the fluid particles, with
the upper lines corresponding to the inner initial positions, i.e., closer to
the central object. The thick red line indicates the case of the initial
radius $r_{0,c}$ being the boundary at which the solution series vary
discontinuously.

Let us analytically evaluate the test-particle model for the case in which
$\dot{q}$ is constant with time. Considering only marginally bound ejecta, we
evaluate Equation (A9) by using Equation (A5). In this case, Equation (A8) can
be written as an algebraic equation for $E_{f}$ as follows:

$E_{f}(r_{0})=\left\\{\begin{array}[]{ll}\displaystyle
E_{0}(r_{0})+\dot{q}_{0}t_{s}\frac{\pi}{2}\left(-\frac{2E_{f}}{c^{2}}\right)^{-3/2}&(t_{s}\tau_{f}<t_{\rm
heat})\vspace{1mm}\\\ \displaystyle E_{0}(r_{0})+\dot{q}_{0}t_{\rm
heat}&(t_{s}\tau_{f}\geq t_{\rm heat})\end{array}\right.,$ (A10)

where we assume $t_{\rm start}\ll t_{s}\tau_{f}$. This can be calculated by
finding the intersection of the two curves represented by the right-hand side
(red) and the left-hand side (blue) as shown in Figure 15. If there are
multiple intersections, the solution with the smallest $E_{f}$ (i.e., the
solution with the shortest turn-around time) is the physical solution.

The series of solutions discontinuously change at which Equation (A10) has a
double root. Since $E_{f,c}$, which corresponds to the point such that the
blue line comes into contact with the red thick curve in Figure 15, is the
root of the derivative of Equation (A10), it can be written as follows:

$E_{f,c}=-\frac{1}{2}\left({3\pi\dot{q}_{0}GM}\right)^{2/5}.$ (A11)

The corresponding fallback time can be written as

$t_{\rm fb,c}\sim
2t_{s}\tau_{f}\left(r_{0,c},E_{f,c}\right)\sim\left(\frac{32}{27}\pi^{2}\right)^{1/5}t_{c},$
(A12)

or numerically,

$t_{\rm fb,c}\sim 1.64t_{c}\sim 0.58\sec\left(\frac{Q_{\rm tot}}{3~{}{\rm
MeV}~{}{\rm
nuc^{-1}}}\right)^{-3/5}\left(\frac{M}{2.7~{}M_{\odot}}\right)^{2/5}\left(\frac{t_{\rm
heat}}{1\sec}\right)^{3/5}.$ (A13)

Here we used $Q_{\rm tot}\sim\dot{q}_{0}t_{\rm heat}$, assuming that the
offset $t_{\rm start}$ of the calculation start time is sufficiently shorter
than $t_{\rm heat}$. The corresponding radius $r_{0,c}$ is determined from

$E_{0}(r_{0,c})=-\frac{5}{6}\left(3\pi\dot{q}_{0}GM\right)^{2/5}.$ (A14)

The right-hand side of this equation corresponding to the solution
($r_{0,c},E_{f,c}$) is shown as the thick red curve in Figure 15. As can be
seen from the figure, the velocities of fluid particles released from the
radius smaller than $r_{0,c}$ become $v=0$ before reaching $t=t_{\rm heat}$
and then the particles start infalling. On the other hand, the fluid particles
released from the radius greater than $r_{0,c}$, which corresponds to those
represented by the red curves below the thick red curve, continue to be heated
until $t=t_{\rm heat}$. In order for a fluid particle released from the radius
greater than $r_{0,c}$ to have a bound solution (i.e., $E_{f}<0$), the
following condition is required:

$E_{0}(r_{0,c})+\dot{q}_{0}t_{\rm heat}<0.$ (A15)

Rewriting this condition in terms of $t_{\rm heat}$ and $Q_{\rm tot}$, we
obtain

$t_{\rm heat}<\sqrt{\frac{3125}{3456}}\pi
t_{s}\left(\frac{Q_{\mathrm{tot}}}{c^{2}}\right)^{-3/2}.$ (A16)

If this condition is satisfied, even a particle released from the outside of
the sphere of $r_{0,c}$ by an infinitesimal distance (see the red curve for
$r_{0}=449.74\ {\rm km}$) has a finite fallback time $t_{\rm fb,r}$ longer
than $t_{\rm fb,c}$, namely,

$t_{\rm
fb,r}=\frac{3}{4}\left(\frac{\pi^{4}}{3}\right)^{1/10}\left(\frac{5}{4}-\frac{t_{\rm
heat}}{t_{\rm fb,c}}\right)^{-3/2}t_{c}.$ (A17)

Evaluating the value of $t_{\rm fb,r}$ for $\dot{M}=2.7M_{\odot}$,
$\dot{q}_{0}=3~{}{\rm MeV}~{}{\rm nuc^{-1}}~{}{\rm s^{-1}}$, and $t_{\rm
heat}=0.6\sec$ gives $t_{\rm fb,r}\sim 3.95\sec$. Further, outwardly released
fluid particles (see the red curve of Figure 15 for $r_{0}=456.03\ {\rm km}$)
have a longer fallback time than $t_{\rm fb,r}$ and thus the mass accretion
continues. This is exactly the “gap”, the suspension of mass accretion between
$t=t_{\rm fb,c}$ and $t_{\rm fb,r}$, which has been shown in Metzger et al.
(2010). On the other hand, if $t_{\rm heat}$ is sufficiently long such that
Equation (A16) is not satisfied, the mass accretion halts and never resumes.
This is what has been demonstrated as a “cutoff” case in Metzger et al.
(2010). In fact, we find a cutoff at the time calculated from Equation (A13)
for the test-particle model shown in Figure 3.

Equation (A16) is only a necessary condition for which a gap of mass accretion
appears. For this case, there must be a double root $r_{0,c}$, in other words,
$t_{\rm heat}$ must be sufficiently long enough for mass accretion to stop
once. This can be given by $t_{s}\tau_{f}\left(r_{0,c},E_{f,c}\right)>t_{\rm
heat}$. This can be also written as a condition for $t_{\rm heat}$ and $Q_{\rm
tot}$:

$\frac{1}{6\sqrt{3}}\pi t_{s}\left(\frac{Q_{\text{tot
}}}{c^{2}}\right)^{-3/2}<t_{\text{heat }}.$ (A18)

A gap of mass accretion appears in the test-particle model if both Equations
(A16) and (A18) are satisfied. Metzger et al. (2010) classified the
qualitative behavior of mass accretion by introducing a parameter $\eta\equiv
t_{\rm heat}/t_{\rm fb,c}$. Using this parameter, we have

$\frac{1}{2}<\eta<\frac{5}{4}.$ (A19)

The lower and upper limits of the inequality represent the conditions under
which mass accretion stops and resumes, respectively. The upper bound of
$1.25$ was also obtained in Desai et al. (2019), which confirms our
theoretical framework being equivalent to theirs. The cutoff condition
represented by the dashed lines for the test-particle model in Figure 11 is
determined such that the lower bound of $\eta$ becomes $1/2$.

## References

  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, ApJ, 848, L13, doi: 10.3847/2041-8213/aa920c
  * Abbott et al. (2017b) —. 2017b, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Arcavi et al. (2017) Arcavi, I., Hosseinzadeh, G., Howell, D. A., et al. 2017, Nature, 551, 64, doi: 10.1038/nature24291
  * Ascenzi et al. (2019) Ascenzi, S., Coughlin, M. W., Dietrich, T., et al. 2019, MNRAS, 486, 672, doi: 10.1093/mnras/stz891
  * Balasubramanian et al. (2021) Balasubramanian, A., Corsi, A., Mooley, K. P., et al. 2021, arXiv e-prints, arXiv:2103.04821.  https://arxiv.org/abs/2103.04821
  * Balsara (2017) Balsara, D. S. 2017, Living Reviews in Computational Astrophysics, 3, 2, doi: 10.1007/s41115-017-0002-8
  * Barnes et al. (2016) Barnes, J., Kasen, D., Wu, M.-R., & Martínez-Pinedo, G. 2016, ApJ, 829, 110, doi: 10.3847/0004-637X/829/2/110
  * Barnes et al. (2020) Barnes, J., Zhu, Y. L., Lund, K. A., et al. 2020, arXiv e-prints, arXiv:2010.11182.  https://arxiv.org/abs/2010.11182
  * Barthelmy et al. (2005) Barthelmy, S. D., Barbier, L. M., Cummings, J. R., et al. 2005, Space Sci. Rev., 120, 143, doi: 10.1007/s11214-005-5096-3
  * Bauswein et al. (2013) Bauswein, A., Goriely, S., & Janka, H. T. 2013, ApJ, 773, 78, doi: 10.1088/0004-637X/773/1/78
  * Bondi (1952) Bondi, H. 1952, MNRAS, 112, 195, doi: 10.1093/mnras/112.2.195
  * Bucciantini et al. (2012) Bucciantini, N., Metzger, B. D., Thompson, T. A., & Quataert, E. 2012, MNRAS, 419, 1537, doi: 10.1111/j.1365-2966.2011.19810.x
  * Burrows et al. (2005) Burrows, D. N., Romano, P., Falcone, A., et al. 2005, Science, 309, 1833, doi: 10.1126/science.1116168
  * Chevalier (1989) Chevalier, R. A. 1989, ApJ, 346, 847, doi: 10.1086/168066
  * Chornock et al. (2017) Chornock, R., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L19, doi: 10.3847/2041-8213/aa905c
  * Coulter et al. (2017) Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556, doi: 10.1126/science.aap9811
  * Cowperthwaite et al. (2017) Cowperthwaite, P. S., Berger, E., Villar, V. A., et al. 2017, ApJ, 848, L17, doi: 10.3847/2041-8213/aa8fc7
  * Del Zanna & Bucciantini (2002) Del Zanna, L., & Bucciantini, N. 2002, A&A, 390, 1177, doi: 10.1051/0004-6361:20020776
  * Desai et al. (2019) Desai, D., Metzger, B. D., & Foucart, F. 2019, MNRAS, 485, 4404, doi: 10.1093/mnras/stz644
  * Drout et al. (2017) Drout, M. R., Piro, A. L., Shappee, B. J., et al. 2017, Science, 358, 1570, doi: 10.1126/science.aaq0049
  * Eichler et al. (1989) Eichler, D., Livio, M., Piran, T., & Schramm, D. N. 1989, Nature, 340, 126, doi: 10.1038/340126a0
  * Fan et al. (2013) Fan, Y.-Z., Yu, Y.-W., Xu, D., et al. 2013, ApJ, 779, L25, doi: 10.1088/2041-8205/779/2/L25
  * Fernández & Metzger (2013) Fernández, R., & Metzger, B. D. 2013, MNRAS, 435, 502, doi: 10.1093/mnras/stt1312
  * Fernández et al. (2019) Fernández, R., Tchekhovskoy, A., Quataert, E., Foucart, F., & Kasen, D. 2019, MNRAS, 482, 3373, doi: 10.1093/mnras/sty2932
  * Fujibayashi et al. (2020) Fujibayashi, S., Wanajo, S., Kiuchi, K., et al. 2020, ApJ, 901, 122, doi: 10.3847/1538-4357/abafc2
  * Gompertz et al. (2018) Gompertz, B. P., Levan, A. J., Tanvir, N. R., et al. 2018, ApJ, 860, 62, doi: 10.3847/1538-4357/aac206
  * Goodman (1986) Goodman, J. 1986, ApJ, 308, L47, doi: 10.1086/184741
  * Hajela et al. (2019) Hajela, A., Margutti, R., Alexander, K. D., et al. 2019, ApJ, 886, L17, doi: 10.3847/2041-8213/ab5226
  * Hajela et al. (2021) Hajela, A., Margutti, R., Bright, J. S., et al. 2021, arXiv e-prints, arXiv:2104.02070.  https://arxiv.org/abs/2104.02070
  * Hamidani & Ioka (2021) Hamidani, H., & Ioka, K. 2021, MNRAS, 500, 627, doi: 10.1093/mnras/staa3276
  * Horesh et al. (2016) Horesh, A., Hotokezaka, K., Piran, T., Nakar, E., & Hancock, P. 2016, ApJ, 819, L22, doi: 10.3847/2041-8205/819/2/L22
  * Hotokezaka & Nakar (2020) Hotokezaka, K., & Nakar, E. 2020, ApJ, 891, 152, doi: 10.3847/1538-4357/ab6a98
  * Hotokezaka et al. (2016) Hotokezaka, K., Wanajo, S., Tanaka, M., et al. 2016, MNRAS, 459, 35, doi: 10.1093/mnras/stw404
  * Ioka et al. (2005) Ioka, K., Kobayashi, S., & Zhang, B. 2005, ApJ, 631, 429, doi: 10.1086/432567
  * Ioka & Nakamura (2018) Ioka, K., & Nakamura, T. 2018, Progress of Theoretical and Experimental Physics, 2018, 043E02, doi: 10.1093/ptep/pty036
  * Ioka & Nakamura (2019) —. 2019, MNRAS, 487, 4884, doi: 10.1093/mnras/stz1650
  * Ishizaki et al. (2021) Ishizaki, W., Ioka, K., & Kiuchi, K. 2021, ApJ, 916, L13, doi: 10.3847/2041-8213/ac1120
  * Just et al. (2015) Just, O., Bauswein, A., Ardevol Pulpillo, R., Goriely, S., & Janka, H. T. 2015, MNRAS, 448, 541, doi: 10.1093/mnras/stv009
  * Kagawa et al. (2019) Kagawa, Y., Yonetoku, D., Sawano, T., et al. 2019, ApJ, 877, 147, doi: 10.3847/1538-4357/ab1bd6
  * Kasen & Barnes (2019) Kasen, D., & Barnes, J. 2019, ApJ, 876, 128, doi: 10.3847/1538-4357/ab06c2
  * Kasen et al. (2017) Kasen, D., Metzger, B., Barnes, J., Quataert, E., & Ramirez-Ruiz, E. 2017, Nature, 551, 80, doi: 10.1038/nature24453
  * Kasliwal et al. (2017) Kasliwal, M. M., Nakar, E., Singer, L. P., et al. 2017, Science, 358, 1559, doi: 10.1126/science.aap9455
  * Kawaguchi et al. (2020) Kawaguchi, K., Fujibayashi, S., Shibata, M., Tanaka, M., & Wanajo, S. 2020, arXiv e-prints, arXiv:2012.14711.  https://arxiv.org/abs/2012.14711
  * Kawaguchi et al. (2018) Kawaguchi, K., Shibata, M., & Tanaka, M. 2018, ApJ, 865, L21, doi: 10.3847/2041-8213/aade02
  * Kawanaka et al. (2013) Kawanaka, N., Piran, T., & Krolik, J. H. 2013, ApJ, 766, 31, doi: 10.1088/0004-637X/766/1/31
  * Kilpatrick et al. (2017) Kilpatrick, C. D., Foley, R. J., Kasen, D., et al. 2017, Science, 358, 1583, doi: 10.1126/science.aaq0073
  * Kisaka & Ioka (2015) Kisaka, S., & Ioka, K. 2015, ApJ, 804, L16, doi: 10.1088/2041-8205/804/1/L16
  * Kisaka et al. (2017) Kisaka, S., Ioka, K., & Sakamoto, T. 2017, ApJ, 846, 142, doi: 10.3847/1538-4357/aa8775
  * Kiuchi et al. (2017) Kiuchi, K., Kawaguchi, K., Kyutoku, K., et al. 2017, Phys. Rev. D, 96, 084060, doi: 10.1103/PhysRevD.96.084060
  * Kouveliotou et al. (1993) Kouveliotou, C., Meegan, C. A., Fishman, G. J., et al. 1993, ApJ, 413, L101, doi: 10.1086/186969
  * Kyutoku et al. (2015) Kyutoku, K., Ioka, K., Okawa, H., Shibata, M., & Taniguchi, K. 2015, Phys. Rev. D, 92, 044028, doi: 10.1103/PhysRevD.92.044028
  * Lee & Ramirez-Ruiz (2007) Lee, W. H., & Ramirez-Ruiz, E. 2007, New Journal of Physics, 9, 17, doi: 10.1088/1367-2630/9/1/017
  * Matsumoto et al. (2020) Matsumoto, T., Kimura, S. S., Murase, K., & Mészáros, P. 2020, MNRAS, 493, 783, doi: 10.1093/mnras/staa305
  * McCully et al. (2017) McCully, C., Hiramatsu, D., Howell, D. A., et al. 2017, ApJ, 848, L32, doi: 10.3847/2041-8213/aa9111
  * Metzger & Bower (2014) Metzger, B. D., & Bower, G. C. 2014, MNRAS, 437, 1821, doi: 10.1093/mnras/stt2010
  * Metzger et al. (2008) Metzger, B. D., Quataert, E., & Thompson, T. A. 2008, MNRAS, 385, 1455, doi: 10.1111/j.1365-2966.2008.12923.x
  * Metzger et al. (2010) Metzger, B. D., Martínez-Pinedo, G., Darbha, S., et al. 2010, MNRAS, 406, 2650, doi: 10.1111/j.1365-2966.2010.16864.x
  * Michel (1988) Michel, F. C. 1988, Nature, 333, 644, doi: 10.1038/333644a0
  * Mösta et al. (2015) Mösta, P., Ott, C. D., Radice, D., et al. 2015, Nature, 528, 376, doi: 10.1038/nature15755
  * Murase et al. (2018) Murase, K., Toomey, M. W., Fang, K., et al. 2018, ApJ, 854, 60, doi: 10.3847/1538-4357/aaa48a
  * Nicholl et al. (2017) Nicholl, M., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L18, doi: 10.3847/2041-8213/aa9029
  * Nishimura et al. (2015) Nishimura, N., Takiwaki, T., & Thielemann, F.-K. 2015, ApJ, 810, 109, doi: 10.1088/0004-637X/810/2/109
  * Nishimura et al. (2006) Nishimura, S., Kotake, K., Hashimoto, M.-a., et al. 2006, ApJ, 642, 410, doi: 10.1086/500786
  * Paczynski (1986) Paczynski, B. 1986, ApJ, 308, L43, doi: 10.1086/184740
  * Radice et al. (2016) Radice, D., Galeazzi, F., Lippuner, J., et al. 2016, MNRAS, 460, 3255, doi: 10.1093/mnras/stw1227
  * Rees (1988) Rees, M. J. 1988, Nature, 333, 523, doi: 10.1038/333523a0
  * Rossi & Begelman (2009) Rossi, E. M., & Begelman, M. C. 2009, MNRAS, 392, 1451, doi: 10.1111/j.1365-2966.2008.14139.x
  * Rosswog (2007) Rosswog, S. 2007, MNRAS, 376, L48, doi: 10.1111/j.1745-3933.2007.00284.x
  * Rosswog et al. (1999) Rosswog, S., Liebendörfer, M., Thielemann, F. K., et al. 1999, A&A, 341, 499.  https://arxiv.org/abs/astro-ph/9811367
  * Rowlinson et al. (2013) Rowlinson, A., O’Brien, P. T., Metzger, B. D., Tanvir, N. R., & Levan, A. J. 2013, MNRAS, 430, 1061, doi: 10.1093/mnras/sts683
  * Shappee et al. (2017) Shappee, B. J., Simon, J. D., Drout, M. R., et al. 2017, Science, 358, 1574, doi: 10.1126/science.aaq0186
  * Shibata et al. (2017) Shibata, M., Fujibayashi, S., Hotokezaka, K., et al. 2017, Phys. Rev. D, 96, 123012, doi: 10.1103/PhysRevD.96.123012
  * Smartt et al. (2017) Smartt, S. J., Chen, T. W., Jerkstrand, A., et al. 2017, Nature, 551, 75, doi: 10.1038/nature24303
  * Soares-Santos et al. (2017) Soares-Santos, M., Holz, D. E., Annis, J., et al. 2017, ApJ, 848, L16, doi: 10.3847/2041-8213/aa9059
  * Tanaka et al. (2017) Tanaka, M., Utsumi, Y., Mazzali, P. A., et al. 2017, PASJ, 69, 102, doi: 10.1093/pasj/psx121
  * Tanvir et al. (2017) Tanvir, N. R., Levan, A. J., González-Fernández, C., et al. 2017, ApJ, 848, L27, doi: 10.3847/2041-8213/aa90b6
  * Uchida et al. (2017) Uchida, H., Shibata, M., Yoshida, T., Sekiguchi, Y., & Umeda, H. 2017, Phys. Rev. D, 96, 083016, doi: 10.1103/PhysRevD.96.083016
  * Villar et al. (2017) Villar, V. A., Guillochon, J., Berger, E., et al. 2017, ApJ, 851, L21, doi: 10.3847/2041-8213/aa9c84
  * Wanajo et al. (2014) Wanajo, S., Sekiguchi, Y., Nishimura, N., et al. 2014, ApJ, 789, L39, doi: 10.1088/2041-8205/789/2/L39
  * Waxman et al. (2018) Waxman, E., Ofek, E. O., Kushnir, D., & Gal-Yam, A. 2018, MNRAS, 481, 3423, doi: 10.1093/mnras/sty2441
  * Zhang & Mészáros (2001) Zhang, B., & Mészáros, P. 2001, ApJ, 552, L35, doi: 10.1086/320255
  * Zhu et al. (2021) Zhu, Y. L., Lund, K. A., Barnes, J., et al. 2021, ApJ, 906, 94, doi: 10.3847/1538-4357/abc69e