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

In order to investigate the effect of 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 heating and the point-source gravity, are written as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\rho v\right)=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial v}{\partial t}+v\displaystyle \frac{\partial v}{\partial r}=-\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial P}{\partial r}-\displaystyle \frac{{GM}}{{r}^{2}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{1}{2}\rho {v}^{2}+{\epsilon }_{\mathrm{int}}\right)+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left[{r}^{2}v\left(\displaystyle \frac{1}{2}\rho {v}^{2}+{\epsilon }_{\mathrm{int}}+P\right)\right]\\ \quad =\ -\rho v\displaystyle \frac{{GM}}{{r}^{2}}+{Q}_{\mathrm{heat}},\end{array}\end{eqnarray} \tag{ 3 }$$

where ρ is the mass density, v is the radial bulk velocity of the fluid, P is the pressure, _int is the internal energy, M is the mass of the central object, and Q_heat is the radioactive heating rate per unit volume per unit time. ⁵ The equation of state is assumed to follow the Γ-law,

$$\begin{eqnarray}&&P=\left({\rm{\Gamma }}-1\right){\epsilon }_{\mathrm{int}},\end{eqnarray} \tag{ 4 }$$

where Γ is the adiabatic index. In Equation (4), we adopt the radiation dominant case of Γ = 4/3. Here, we neglect 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 Newtonian gravity. The advection term of the hydrodynamic equations is solved by using the HLL method (e.g., Del Zanna & Bucciantini 2002) and the third-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 nonrelativistic regime. The calculation is performed by dividing the radius from 40 km to 80,000 km into a uniform grid of 16,384 cells (about ∼4.9 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).

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_⊙. The equation of state (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_⊙ is 12.3 km.

Figure 1 shows the radial profiles of radial velocity and mass density. The thin lines in Figure 1 show the result of Kiuchi et al. (2017) at ≈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_⊙. It can be confirmed that the fluid elements located at r ≲ 490 km are gravitationally bound. Around the radius r ∼ 490 km, which is important for the calculation of mass accretion rates, the radial dependence of velocity and density is found to be weakly dependent on latitude. This is one justification of our treatment of spherically symmetric modeling.

Zoom In Zoom Out Reset image size

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 velocity is proportional to the radius (out to 3000 km):

$$\begin{eqnarray}\begin{array}{l}v(t=0,r)=\left\{\begin{array}{ll}0.26\,c\left(\tfrac{r}{1000\ \mathrm{km}}\right) & \left(r\lt 3000\ \mathrm{km}\right)\\ 0 & \left(r\gt 3000\ \mathrm{km}\right)\end{array}\right.,\end{array}\end{eqnarray} \tag{ 5 }$$

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; ≈10 ms after the merger). The radius of r ∼ 3000 km corresponds to that which the outermost ejecta reaches with nearly the speed of light during about 10 ms. For the mass density, we adopt a broken power-law model, which has a break at 600 km, namely

$$\begin{eqnarray}\rho (t=0,r)=\left\{\begin{array}{ll}{\rho }_{0}{\left(\tfrac{r}{600\ \mathrm{km}}\right)}^{-2.4} & \left(r\lt 600\ \mathrm{km}\right)\\ {\rho }_{0}{\left(\tfrac{r}{600\ \mathrm{km}}\right)}^{-7.5} & \left(r\gt 600\ \mathrm{km}\right)\end{array}\right.,\end{eqnarray} \tag{ 6 }$$

where ρ₀ = 2.0 × 10⁶ g cm⁻³ is the value at r = 600 km. We assume P = 10⁻⁵ ρ c² 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)\,{\rm{s}}$, 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. (2021) 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 ∼10 s, the velocity profile of the wind approximately matches that of a homologous expansion in Equation (5). The density distribution is also expected to be close to ρ ∝ r⁻² such as that of a steady-state supersonic flow, which is qualitatively similar to the radial dependence of the density distribution expressed in Equation (6) for r < 600 km. However, for the viscosity-driven wind, to take into account the difference in total ejecta mass, it will be necessary to increase the value of mass density ρ₀ by one order of magnitude compared with the value of the dynamical ejecta. Here, we presume that the disk wind is modeled by Equations (5) and (6) with a larger value of density ρ₀ 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.

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 β-decay, α-decay, and fission of r-process nuclei produced in the dynamical ejecta. Each thin curve shows the heating rate in units of MeV nuc⁻¹ s⁻¹ as a function of time (since the merger) for the Lagrangian tracer particle of the ejecta with a given initial Y_e . The color indicates the value of Y_e from 0.09 (purple) to 0.44 (yellow) with an interval of ΔY_e = 0.01. The heating rate averaged over the ejecta mass is also shown by the red thick curve.

Zoom In Zoom Out Reset image size

As can be seen in Figure 2, the radioactive heating rates exhibit two phases: one in which the value is approximately constant over time (the constant phase) and the other in which the value decays with time (the decay phase). It can be seen that the duration of the constant phase varies and tends to be longer for larger values of Y_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:

$$\begin{eqnarray}{Q}_{\mathrm{heat}}=\left\{\begin{array}{ll}\rho {\dot{q}}_{0} & (t\lt {t}_{\mathrm{heat}})\\ 0 & (t\gt {t}_{\mathrm{heat}})\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

where ${\dot{q}}_{0}$ is the radioactive energy (except for that in neutrinos) released per unit mass per unit time during the constant phase ⁶ and t_heat is the duration of the constant phase. The value of Y_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_e in the dynamical ejecta of BNS mergers are about Y_e ∼ 0.1–0.4. As mentioned in Section 2.2, the dominant component of ejecta can be the late-time viscosity-driven wind. This component is expected to have a large Y_e compared to that of the dynamical ejecta, with Y_e ∼ 0.3–0.4, as shown in previous work, e.g., Fujibayashi et al. (2020). Considering the uncertainties in Y_e over the different components of ejecta, we adopt the simple model described by Equation (7) and treat ${\dot{q}}_{0}$ and t_heat as parameters. Furthermore, in Section 3, we model the accretion flows semi-analytically. For this purpose, a simplified treatment in Equation (7) is convenient. The modeling of the case with heating rates including the decay phase is given in Section 4.

Here, we neglect a possible effect of heating due to the jet that 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 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).

The thick red line in Figure 3 presents the time evolution of the mass accretion rate calculated with ${\dot{q}}_{0}=3\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}$ and t_heat = 5 s in Equation (7). Here, in order to clearly illustrate the effect of radioactive heating on the mass accretion rate, we have adopted the value of t_heat longer than those shown in Figure 2 for the relevant ${\dot{q}}_{0}$. We assume that all of the radioactive energy (except for that in neutrinos) is converted into the internal energy 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 ∝ t^−5/3 (where the initial ∼0.05 s is affected by the initial conditions). Here, the mass accretion rate is evaluated using the mass flux at r = 650 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 s corresponds to the termination of the heating t = t_heat. After t = t_heat, the mass accretion rate continues roughly in proportion to t^−5/3.

Zoom In Zoom Out Reset image size

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 the Appendix. 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 ∼ 5 s but slowly decreases, lasting a few times longer duration. This is due to the difference in conversion of internal energy to 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. By contrast, 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, a larger 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.

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 dotted–dashed and thin dashed curves show the results for t_heat = 1 s and t_heat = 10 s, 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\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}$, where the solid and dashed curves represent the cases for t_heat = 1 s and t_heat = 4 s, 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\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}$ because ${\dot{q}}_{0}$ is larger. The semi-analytical modeling described in Section 3 explains these behaviors.

Zoom In Zoom Out Reset image size

Finally, we present the mass accretion rate calculated using a more realistic heating rate rather than using a constant approximation in Equation (7). 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 s, 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 s, as will be discussed in Section 5, the mass accretion rate is expected to be suppressed to less than ∼10% of that without heating after ${ \mathcal O }({10}^{4})$–${ \mathcal O }({10}^{8})$ s due to continuous heating in the decay phase.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{P}{\rho }\right)+v\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{P}{\rho }\right)\\ =\ ({\rm{\Gamma }}-1)\left[{\dot{q}}_{0}-\displaystyle \frac{P}{\rho }\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}v\right)\right].\end{array}\end{eqnarray} \tag{ 8 }$$

The characteristic parameters of the fluid Equations (1), (2), and (8) 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:

$$\begin{eqnarray}&&\begin{array}{l}{r}_{c}={\left[\displaystyle \frac{{({GM})}^{3}}{{\dot{q}}_{0}^{2}}\right]}^{1/5}\\ \sim \ 3540\,\mathrm{km}\,{\left(\displaystyle \frac{M}{2.7\,{M}_{\odot }}\right)}^{3/5}{\left(\displaystyle \frac{{\dot{q}}_{0}}{3\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}}\right)}^{-2/5},\end{array}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\begin{array}{l}{t}_{c}={\left[\displaystyle \frac{{({GM})}^{2}}{{\dot{q}}_{0}^{3}}\right]}^{1/5}\\ \sim \ 0.35\,{\rm{s}}\,{\left(\displaystyle \frac{M}{2.7\,{M}_{\odot }}\right)}^{2/5}{\left(\displaystyle \frac{{\dot{q}}_{0}}{3\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}}\right)}^{-3/5}.\end{array}\end{eqnarray} \tag{ 10 }$$

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

$$\begin{eqnarray}&&\begin{array}{l}{v}_{c}={\left({GM}{\dot{q}}_{0}\right)}^{1/5}\\ \sim \ 0.033c\,{\left(\displaystyle \frac{M}{2.7\,{M}_{\odot }}\right)}^{1/5}{\left(\displaystyle \frac{{\dot{q}}_{0}}{3\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}}\right)}^{1/5}.\end{array}\end{eqnarray} \tag{ 11 }$$

Then, we introduce the dimensionless length ξ ≡ r/r_c , time χ ≡t/t_c , density ϕ ≡ ρ/ρ_c , velocity u ≡ v/v_c , and pressure $\theta \equiv P/\rho {v}_{c}^{2}=P/\phi {\rho }_{c}{v}_{c}^{2}$, where ρ_c is an arbitrary constant with the dimension of density. Rewriting Equations (1), (2), and (8) using these dimensionless variables, we obtain the following dimensionless equations:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \phi }{\partial \chi }+\displaystyle \frac{1}{{\xi }^{2}}\displaystyle \frac{\partial }{\partial \xi }\left({\xi }^{2}\phi u\right)=0,\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial u}{\partial \chi }+u\displaystyle \frac{\partial u}{\partial \xi }=-\displaystyle \frac{1}{\phi }\displaystyle \frac{\partial \left(\phi \theta \right)}{\partial \xi }-\displaystyle \frac{1}{{\xi }^{2}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \theta }{\partial \chi }+v\displaystyle \frac{\partial \theta }{\partial \xi }=({\rm{\Gamma }}-1)\left[1-\theta \displaystyle \frac{1}{{\xi }^{2}}\displaystyle \frac{\partial }{\partial \xi }\left({\xi }^{2}u\right)\right].\end{eqnarray} \tag{ 14 }$$

Because these nondimensional equations have the same form for any ρ_c , there is no characteristic scale of ejecta mass (equivalently mass density) in this system (i.e., Equations (1), (2), and (8)). As mentioned in Section 2.2, we presume that the profile of the viscosity-driven wind can be modeled by enhancing the density ρ (see Equation (6)). 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.

Equations (1), (2), and (8) can be normalized by Equations (9)–(11) to eliminate the parameters GM and ${\dot{q}}_{0}$. Therefore, the accretion flow follows the same equations under the characteristic scales (9)–(11) up to t = t_heat. Figure 6 shows the time evolution of the turn-around radius r_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_heat/t_c . This clearly shows the effectiveness of the scaling. It is also seen that from t = 6 t_c to t_heat, the turn-around radius r_turn is approximately constant r_turn ∼ r_c over time. Using the chain rule with the differential coefficients of velocity, we can obtain the time evolution of r_turn as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{{dr}}_{\mathrm{turn}}(t)}{{dt}}=-\displaystyle \frac{{\left(\partial v/\partial t\right)}_{r={r}_{\mathrm{turn}}}}{{\left(\partial v/\partial r\right)}_{r={r}_{\mathrm{turn}}}}.\end{eqnarray} \tag{ 15 }$$

At the turn-around radius r = r_turn, the velocity is v = 0 by definition, so that the Lagrangian time derivative becomes Dv/Dt = ∂v/∂t. The fact that r_turn is approximately constant over time, dr_turn(t)/dt ∼ 0 means Dv/Dt ∼ 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_turn does not fall back.

Zoom In Zoom Out Reset image size

In the following, we discuss the phase after r_turn becomes constant, i.e., t > 6 t_c . When the position of r = r_turn is constant, the forces are balanced (see Equation (15)) and the velocity is 0; thus the mass does not fall back from the radius beyond r_turn. Therefore, the decreasing rate of the ejecta mass between a sufficiently small radius (taken to be r_count = 650 km in the numerical calculation) and r_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_turn can be written as

$$\begin{eqnarray}&&{M}_{\mathrm{turn}}=\displaystyle \frac{4\pi }{3}{f}_{0}{r}_{\mathrm{turn}\ }^{3}{\rho }_{\mathrm{turn}\ },\end{eqnarray} \tag{ 16 }$$

where the subscript "turn" means the value at r = r_turn and f₀ 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_turn evaluated with Equation (16) is compared to that of the numerical fluid calculation (excluding that within r_count) for ${\dot{q}}_{0}=3\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}$ and t_heat = 5 s. The bottom panel of Figure 7 shows the value of f₀ in order for M_turn to match the numerical value. Generally, there are three terms contributing to the time derivative of M_turn as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{M}}_{\mathrm{turn}}}{{M}_{\mathrm{turn}}}=\displaystyle \frac{{\dot{f}}_{0}}{{f}_{0}}+3\displaystyle \frac{{\dot{r}}_{\mathrm{turn}}}{{r}_{\mathrm{turn}}}+\displaystyle \frac{{\dot{\rho }}_{\mathrm{turn}}}{{\rho }_{\mathrm{turn}}}.\end{eqnarray} \tag{ 17 }$$

Evaluating the value of each term from the fluid calculations, we find that, between t ∼ 6 t_c and t = t_heat, r_turn and f₀ vary on longer timescales than the mass density ρ_turn. We therefore at first take r_turn and f₀ to be constants and assume that the time variation of M_turn is due only to ρ_turn.

Zoom In Zoom Out Reset image size

The hydrodynamical structure at the time t = 8.5 t_c > 6 t_c , where the turn-around radius r_turn has become constant, is shown in Figure 8. In the region where ∣v∣ ≳ ∣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_turn (Bondi 1952):

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{Bondi}}=2\sqrt{2}\pi {f}_{1}{({GM})}^{2}{\rho }_{\mathrm{turn}}/{a}_{\mathrm{turn}\ }^{3},\end{eqnarray} \tag{ 18 }$$

where $a=\sqrt{{\rm{\Gamma }}P/\rho }$ is the sound speed, a_turn = a(r_turn), and f₁ is an ${ \mathcal O }(1)$ constant that corrects for the difference owing to the fact that r_turn is not infinite. Originally, Equation (18) is a relation between the gas that is sufficiently distant and stationary, but here we evaluate this value at r_turn. The fluid element at the turn-around radius r_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 (18). After t = 4 t_c , it can be seen that the mass accretion rate is well approximated by Equation (18). Equating the mass accretion rate in Equation (18) to that in Equation (17); (in the absolute values), we obtain the equation for the time evolution of the mass density ρ at r = r_turn as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{\rho }}_{\mathrm{turn}}}{{\rho }_{\mathrm{turn}}}=-\displaystyle \frac{3\sqrt{2}}{2}\displaystyle \frac{{f}_{1}}{{f}_{0}}\displaystyle \frac{1}{{t}_{c}}{\left(\displaystyle \frac{{a}_{\mathrm{turn}}}{{v}_{c}}\right)}^{-3}.\end{eqnarray} \tag{ 19 }$$

Here, f₀, f₁, and r_turn ∼ r_c are assumed to be constant and Equations (9)–(11) are used. In order to solve Equation (19), a model of the time evolution of the sound speed a_turn at r = r_turn is needed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we model the temporal evolution of the sound speed at the turn-around radius r_turn. Figure 10 shows the time evolution of the square of a_turn. As can be seen from the figure, the values of ${a}_{\mathrm{turn}}^{2}$ in the numerical fluid calculations are proportional to time after t ≳ 6 t_c up to t = t_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_turn:

$$\begin{eqnarray}&&{a}_{\mathrm{turn}}=\sqrt{\tfrac{{\rm{\Gamma }}{P}_{\mathrm{turn}}}{{\rho }_{\mathrm{turn}}}}={\beta }_{0}{v}_{c}{\left(\displaystyle \frac{t}{{t}_{c}}\right)}^{1/2},\end{eqnarray} \tag{ 20 }$$

where the value of β₀ includes the effect of adiabatic cooling in the accretion flow (that will be $\sqrt{{\rm{\Gamma }}\left({\rm{\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 }(\mathrm{MeV})$, so that this value will never reach the speed of light. By solving Equation (20) for P_turn and by using Equations (4), (10), and (11), we obtain the value of the internal energy, ${\epsilon }_{\mathrm{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.

Zoom In Zoom Out Reset image size

The mass accretion rate ${\dot{M}}_{\mathrm{turn}}$ at r = r_turn, can be calculated if the time evolution of ρ_turn is obtained with Equations (19) and (20). Using Equations (16), (19), and (20), we obtain the following relation:

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}{\dot{M}}_{\mathrm{turn}}}{d\mathrm{ln}t}=-\displaystyle \frac{3}{2}\left[1+\sqrt{2}\displaystyle \frac{{f}_{1}}{{f}_{0}}\displaystyle \frac{1}{{\beta }_{0}^{3}}{\left(\displaystyle \frac{t}{{t}_{c}}\right)}^{-1/2}\right].\end{eqnarray} \tag{ 21 }$$

By integrating this over time, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{M}}_{\mathrm{turn}}}{{\dot{M}}_{0}}={\left(\displaystyle \frac{t}{{t}_{0}}\right)}^{-3/2}\exp \left[-\displaystyle \frac{3\sqrt{2}}{{\beta }_{0}^{3}}\displaystyle \frac{{f}_{1}}{{f}_{0}}\sqrt{\displaystyle \frac{{t}_{c}}{{t}_{0}}}\left(1-\sqrt{\displaystyle \frac{{t}_{0}}{t}}\right)\right],\end{eqnarray} \tag{ 22 }$$

where ${\dot{M}}_{0}$ is the initial condition for ${\dot{M}}_{\mathrm{turn}}$ at a given time t₀. For the evaluation of f₀, f₁, and β₀, we adopt the values at t₀ = 6 t_c , at which the turn-around radius becomes approximately constant (see Figure 6), and use f₁/f₀ = 0.5 and β₀ = 0.43. The resultant mass accretion rates using Equation (22) 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 (22) are not valid. It is clear from Figure 4 that the mass accretion rate can be well approximated by Equation (22) 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\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}$, it can be seen that the theoretical curve also explains well those for ${\dot{q}}_{0}=30\,\mathrm{MeV}\,{\mathrm{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 ≪ 6 t_c , the mass accretion rate exhibits almost the same temporal evolution as in the ${\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 (9)–(11), the ratio of mass accretion rates 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₀ in time units of t_c (as we chose the reference point of ${\dot{M}}_{0}$ as t₀ = 6 t_c ), this ratio at t = t₀ 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:

$$\begin{eqnarray}&&{\dot{M}}_{0}\propto {t}_{c}^{-5/3}\propto {\left({GM}\right)}^{-2/3}{\dot{q}}_{0}.\end{eqnarray} \tag{ 23 }$$

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\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}$. In the case of ${\dot{q}}_{0}=30\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}$ in Figure 4, we adopt ${\dot{M}}_{0}$ calculated using Equation (23).

As represented by Equation (22) 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_B. The first effect is that r_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 inside of the sphere of r_turn from the outside. Thus the mass density ρ_turn at the turn-around radius will necessarily decrease with accretion. Furthermore, the flow outside the sonic radius r_B is almost incompressible because it is a subsonic flow, so that the mass density ρ_B at the sonic point is comparable to ρ_turn. Therefore, ρ_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 (20)), the radius of the sonic point ${r}_{{\rm{B}}}={GM}/2{a}_{{\rm{B}}}^{2}$ (Bondi 1952) shrinks as r_B ∝ t⁻¹. 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 ⁵⁶Ni, 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 ⁵⁶Ni, 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 a power-law function of time (see Equation (22)), 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 decreases as t⁻³. In fact, provided that ρ_B ∝ t⁻³, 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 ρ_B due to accretion.

If t_heat is shorter than 6 t_c , the stagnation of r_turn and the increase of the sound speed will not occur for the constant heating model. Thus, t ≳ 6 t_c is the necessary condition for these two processes to work, namely,

$$\begin{eqnarray}&&\begin{array}{l}{t}_{\mathrm{heat}}\gt {{Kt}}_{c}\\ \sim \ 2.1\,{\rm{s}}{K}_{6}{\left(\displaystyle \frac{M}{2.7\,{M}_{\odot }}\right)}^{2/5}{\left(\displaystyle \frac{{\dot{q}}_{0}}{3\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}}\right)}^{-3/5},\end{array}\end{eqnarray} \tag{ 24 }$$

where K is the time (with respect to t_c ) that it takes for r_turn to become constant and K = 6 K₆ (see Figure 6). After t ∼ 6 t_c , the analytic solution (22) becomes valid, and the mass accretion rate rapidly decreases. Therefore, we call this time, t_halt ∼ 6 t_c , the "halting time" and this suppression of mass accretion "halting." The halting time corresponds to the timescale in which the mass accretion rate decreases to ∼13% of its original value (see the description below Equation (23)). For given ${\dot{q}}_{0}$, t_heat, and M, the halting condition is determined by Equation (24). Figure 11 shows whether or not the halting occurs in the ${\dot{q}}_{0}{t}_{\mathrm{heat}}$−t_heat space. Here, because we are considering only the constant heating phase (see Equation (7), and see also Equation (25) 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, halting is more likely to occur as t_heat increases. This is because the ejecta that 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.

Zoom In Zoom Out Reset image size

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 the A.2. 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 differs 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 a 10 times larger heating rate to halt the mass accretion for a given t_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, with some remaining as internal energy.

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:

$$\begin{eqnarray}&&\displaystyle \frac{{dv}}{{dt}}=-\displaystyle \frac{{GM}}{{r}^{2}},\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{E}_{0}=\displaystyle \frac{1}{2}{v}^{2}-\displaystyle \frac{{GM}}{r}.\end{eqnarray} \tag{ A2 }$$

Let us introduce dimensionless variables x = r/r_s and β = v/c, where r_s = 2GM/c². The dimensionless time is also defined as τ = t/t_s , where t_s = r_s /c. The dimensionless energy per mass λ is determined by

$$\begin{eqnarray}&&\lambda \equiv \displaystyle \frac{1}{{x}_{0}}-{\beta }_{0}^{2}=\displaystyle \frac{1}{x}-{\beta }^{2}=-2{E}_{0}/{c}^{2},\end{eqnarray} \tag{ 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:

$$\begin{eqnarray}&&{\tau }_{f}=\displaystyle \frac{{\beta }_{0}}{\lambda \left({\beta }_{0}^{2}+\lambda \right)}+{\tan }^{-1}\left(\displaystyle \frac{{\beta }_{0}}{\sqrt{\lambda }}\right){\lambda }^{-3/2}.\end{eqnarray} \tag{ A4 }$$

Note that τ_f is a function of the initial velocity and energy. When $\lambda \ll {\beta }_{0}^{2}\lt 1$, τ_f can be written as:

$$\begin{eqnarray}&&{\tau }_{f}\sim \displaystyle \frac{\pi }{2}{\lambda }^{-3/2}.\end{eqnarray} \tag{ A5 }$$

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

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

$$\begin{eqnarray}&&{t}_{\mathrm{fb}}({r}_{0},{E}_{0}({r}_{0}))={t}_{s}\left[{\tau }_{f}({r}_{0},{E}_{0}({r}_{0}))+{\tau }_{f}\left({r}_{\mathrm{fin}},{E}_{0}({r}_{0})\right)\right].\end{eqnarray} \tag{ A6 }$$

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

$$\begin{eqnarray}&&\dot{M}\left({t}_{\mathrm{fb}}\left({r}_{0},{E}_{0}\left({r}_{0}\right)\right)\right)\equiv \displaystyle \frac{{dM}}{{{dt}}_{\mathrm{fb}}}=4\pi {r}_{0}^{2}{\rho }_{0}\left({r}_{0}\right){\left(\displaystyle \frac{{{dt}}_{\mathrm{fb}}}{{{dr}}_{0}}\right)}_{r={r}_{0}}^{-1},\end{eqnarray} \tag{ A7 }$$

where ρ(r₀) is the mass density in the initial state.

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 of the radioactive energy is converted to kinetic energy, the total energy of the particle at the turn-around time can be estimated as follows:

$$\begin{eqnarray}&&{E}_{f}({r}_{0})={E}_{0}({r}_{0})+{\int }_{{t}_{\mathrm{start}}}^{{t}_{s}{\tau }_{f}\left({r}_{0},{E}_{f}\right)}\dot{q}(t){dt},\end{eqnarray} \tag{ 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 τ_f at the upper end of the integration is evaluated by using E_f . In fact, the internal energy injected into the fluid element is converted into kinetic energy via 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 as 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:

$$\begin{eqnarray}&&{t}_{\mathrm{fb}}\left({r}_{0}\right)={t}_{s}\left[{\tau }_{f}\left({r}_{0},{E}_{f}\left({r}_{0}\right)\right)+{\tau }_{f}\left({r}_{\mathrm{fin}},{E}_{f}\left({r}_{0}\right)\right)\right].\end{eqnarray} \tag{ A9 }$$

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

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) using Equation (A5). In this case, Equation (A8) can be written as an algebraic equation for E_f as follows:

$$\begin{eqnarray}{E}_{f}({r}_{0})=\left\{\begin{array}{ll}{E}_{0}({r}_{0})+{\dot{q}}_{0}{t}_{s}\displaystyle \frac{\pi }{2}{\left(-\displaystyle \frac{2{E}_{f}}{{c}^{2}}\right)}^{-3/2} & ({t}_{s}{\tau }_{f}\lt {t}_{\mathrm{heat}})\\ {E}_{0}({r}_{0})+{\dot{q}}_{0}{t}_{\mathrm{heat}} & ({t}_{s}{\tau }_{f}\geqslant {t}_{\mathrm{heat}})\end{array}\right.,\end{eqnarray} \tag{ A10 }$$

where we assume t_start ≪ t_s τ_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.

Zoom In Zoom Out Reset image size

The series of solutions discontinuously change, so 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:

$$\begin{eqnarray}&&{E}_{f,c}=-\displaystyle \frac{1}{2}{\left(3\pi {\dot{q}}_{0}{GM}\right)}^{2/5}.\end{eqnarray} \tag{ A11 }$$

The corresponding fallback time can be written as

$$\begin{eqnarray}&&{t}_{\mathrm{fb},{\rm{c}}}\sim 2{t}_{s}{\tau }_{f}\left({r}_{0,c},{E}_{f,c}\right)\sim {\left(\displaystyle \frac{32}{27}{\pi }^{2}\right)}^{1/5}{t}_{c},\end{eqnarray} \tag{ A12 }$$

or numerically,

$$\begin{eqnarray}&&\begin{array}{l}{t}_{\mathrm{fb},{\rm{c}}}\sim 1.64{t}_{c}\sim 0.58\,{\rm{s}}{\left(\displaystyle \frac{{Q}_{\mathrm{tot}}}{3\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}}\right)}^{-3/5}{\left(\displaystyle \frac{M}{2.7\,{M}_{\odot }}\right)}^{2/5}\\ \qquad \times \ {\left(\displaystyle \frac{{t}_{\mathrm{heat}}}{1\,{\rm{s}}}\right)}^{3/5}.\end{array}\end{eqnarray} \tag{ A13 }$$

Here, we used ${Q}_{\mathrm{tot}}\sim {\dot{q}}_{0}{t}_{\mathrm{heat}}$, assuming that the offset t_start of the calculation start time is sufficiently shorter than t_heat. The corresponding radius r_0,c is determined from

$$\begin{eqnarray}&&{E}_{0}({r}_{0,c})=-\displaystyle \frac{5}{6}{\left(3\pi {\dot{q}}_{0}{GM}\right)}^{2/5}.\end{eqnarray} \tag{ 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 a radius smaller than r_0,c become v = 0 before reaching t = t_heat and then the particles start infalling. By contrast, the fluid particles released from a radius greater than r_0,c , which corresponds to the scenarios represented by the red curves below the thick red curve, continue to be heated until t = t_heat. In order for a fluid particle released from a radius greater than r_0,c to have a bound solution (i.e., E_f < 0), the following condition is required:

$$\begin{eqnarray}&&{E}_{0}({r}_{0,c})+{\dot{q}}_{0}{t}_{\mathrm{heat}}\lt 0.\end{eqnarray} \tag{ A15 }$$

Rewriting this condition in terms of t_heat and Q_tot, we obtain

$$\begin{eqnarray}&&{t}_{\mathrm{heat}}\lt \sqrt{\tfrac{3125}{3456}}\pi {t}_{s}{\left(\displaystyle \frac{{Q}_{\mathrm{tot}}}{{c}^{2}}\right)}^{-3/2}.\end{eqnarray} \tag{ 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₀ = 449.74 km) has a finite fallback time t_fb,r longer than t_fb,c, namely,

$$\begin{eqnarray}&&{t}_{\mathrm{fb},{\rm{r}}}=\displaystyle \frac{3}{4}{\left(\displaystyle \frac{{\pi }^{4}}{3}\right)}^{1/10}{\left(\displaystyle \frac{5}{4}-\displaystyle \frac{{t}_{\mathrm{heat}}}{{t}_{\mathrm{fb},{\rm{c}}}}\right)}^{-3/2}{t}_{c}.\end{eqnarray} \tag{ A17 }$$

Evaluating the value of t_fb,r for $\dot{M}=2.7{M}_{\odot }$, ${\dot{q}}_{0}\,=3\,\mathrm{MeV}\,{\mathrm{nuc}}^{-1}\,{{\rm{s}}}^{-1}$, and t_heat = 0.6 s gives t_fb,r ∼ 3.95 s. Further, outwardly released fluid particles (see the red curve of Figure 15 for r₀ = 456.03 km) have a longer fallback time than t_fb,r and thus the mass accretion continues. This is exactly the "gap," the suspension of mass accretion between t = t_fb,c and t_fb,r, which has been shown in Metzger et al. (2010). However, if t_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 in mass accretion appears. For this case, there must be a double root r_0,c ; in other words, t_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)\gt {t}_{\mathrm{heat}}$. This can be also written as a condition for t_heat and Q_tot:

$$\begin{eqnarray}&&\displaystyle \frac{1}{6\sqrt{3}}\pi {t}_{s}{\left(\displaystyle \frac{{Q}_{\mathrm{tot}\ }}{{c}^{2}}\right)}^{-3/2}\lt {t}_{\mathrm{heat}\ }.\end{eqnarray} \tag{ A18 }$$

A gap in 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 η ≡ t_heat/t_fb,c. Using this parameter, we have

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\lt \eta \lt \displaystyle \frac{5}{4}.\end{eqnarray} \tag{ 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 is 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 η becomes 1/2.