"Extended Emission" from Fallback Accretion onto Merger Remnants

As mentioned in the Introduction, what distinguishes our work from similar studies in the literature is the realistic modeling of BNS mergers so as to have not only an accurate description of the matter ejection on the nonlinear physics involved in the merger but also an accurate dependence of the fallback dynamics on the properties of the binaries. To this scope, we have considered four representative binaries described by two different and temperature-dependent EOSs. The first one is the SFHO EOS (Hempel & Schaffner-Bielich 2010), which is based on relativistic mean field calculations and predicts a maximum nonrotating Tolman–Oppenheimer–Volkoff (TOV) mass of M_TOV ≃ 2.06 M_⊙; the second EOS is the TNTYST (Togashi et al. 2017), which is based on the variational many-body theory and instead allows for higher masses, with M_TOV ≃ 2.23 M_⊙. For all four binaries, the total gravitational mass at infinite separation exceeds the expected maximum mass for a uniformly rotating remnant with their respective EOSs assuming a ratio ${M}_{\max }/{M}_{\mathrm{TOV}}\simeq 1.2\mbox{--}1.25$ (Breu & Rezzolla 2016; Musolino et al. 2024). On the other hand, all the binaries have the same chirp mass corresponding to the estimated value for the GW170917 event (Abbott et al. 2017b), and considering M₂ (M₁) to be the gravitational mass of the secondary (primary) star in the binary, we consider two different mass ratios q:=M₂/M₁ ≤ 1, namely, q = 1.00 and q = 0.75, to model either an equal or an unequal-mass merger. Overall, this choice of EOS and mass ratio allows us to explore and contrast the scenarios modeled by the simulation in which the merger remnant remains a hypermassive neutron star (HMNS) from the one in which it collapses to a black hole a few millisecond after the merger. Details on the properties of the four binaries considered and of their outcome can be found in Table 1.

As mentioned earlier, the evolution of the four binaries is carried out making use of the state-of-the-art numerical 3+1 code suite developed in Frankfurt, which consists of the FIL code for the higher-order finite-difference solution of the GRMHD equations using high-resolution shock-capturing methods (Toro 1999; Rezzolla & Zanotti 2013) and a fourth-order accurate conservative finite-difference scheme (Del Zanna et al. 2007). The evolution of the spacetime is instead carried out with the Antelope spacetime solver (Most et al. 2019), which solves the constraint damping formulation of the Z4 formulation of the Einstein equations (Bernuzzi & Hilditch 2010; Alic et al. 2012). The code suite evolves the full set of equations in conjunction with the EinsteinToolkit (Loeffler et al. 2012; Zlochower et al. 2022), exploiting the Carpet box-in-box adaptive mesh refinement (AMR) driver in Cartesian coordinates (Schnetter et al. 2006),

In addition, to accurately capture the temperature and composition of the ejected and bound material, the effect of weak interactions needs to be properly included. We do so by employing the recently developed code FIL-M1 (Musolino & Rezzolla 2024), which employs a moment-based neutrino-transport scheme. More specifically, FIL-M1 evolves the first two moments of the energy-integrated Boltzmann equations for neutrinos using standard, second-order, finite-volume, high-resolution shock-capturing methods. Thanks to the moment-based radiation transport scheme, we are able to include all the leading-order weak reactions in our simulations, thereby accounting for neutrino heating and cooling and momentum transfer, as well as the composition changes resulting from weak interactions.

All our simulations are performed on a grid spanning Ω = [−1500, 1500]³ km with seven fixed levels of refinement and the grid spacing on the finest level being dx = 310 m. Moreover, we include the effect of magnetic fields in all our simulations, initializing it as a purely poloidal magnetic field confined inside the stars, with a maximum initial strength of ≃10^16.5 G for the equal-mass systems SFHO-q1.0 and TNTYST-q1.0, ≃10^16.7 G for SFHO-q.75, and ≃10¹⁵ G for TNTYST-q.75 (see Chabanov et al. 2023 for a recent discussion of why magnetic fields confined to the crust are not likely). This value, despite being much larger than what would be expected for old neutron stars at the end of their inspiral, has a corresponding magnetic energy that is much smaller than the binding energy of the system and hence does not significantly affect the dynamics of the system before merger (see, e.g., Giacomazzo et al. 2009). In addition, and as is customary in this type of simulation, starting with an unrealistically large magnetic field allows us to better capture realistic magnetic field values after the BNS merger despite our resolution being insufficient for correctly capturing the dynamical magnetic field amplification due to the Kelvin–Helmholtz instability.

Clearly, a very important aspect of our research consists of properly distinguishing matter that is gravitationally bound from matter that is not and hence will not fall back onto the merger remnant. To this scope, and as is customary in these simulations, we construct a matter "detector" consisting of a spherical surface at a distance R_d = 300 M_⊙ from the center of coordinates and record the properties of the matter flowing out of such a 2-sphere. In particular, we use the covariant time component of the material's four-velocity u_t to determine the matter that is gravitationally bound (i.e., u_t > −1) or unbound (i.e., u_t < −1). We favor this criterion, which is also referred to as the "geodesic criterion," over the other common alternative represented by the so-called "Bernoulli criterion," hu_t > −1 (Rezzolla & Zanotti 2013), where h is the specific enthalpy, because the pure kinematic quantity u_t relates directly to the subsequent dynamics of the bound material in the gravitational field of the remnant that we will use to determine the fallback time below. We note that, in principle, hydrodynamical processes can convert internal energy to kinetic energy, thus possibly unbinding material that was initially deemed bound through the u_t criterion; similarly, matter that is considered unbound can undergo shock heating and become bound (see Bovard et al. 2017 for a discussion). However, as we will show in Section 2.3 by varying the position of the matter detector, the conversion of internal energy into kinetic energy does not occur in practice in our simulations, likely because of the small rest-mass density of the matter once it reaches the detector.

Because the dynamics of the bound matter outside of the detection sphere cannot be followed in detail (even if bound, the matter past the detector at R_d travels outward to distances that are orders of magnitude larger than R_d ), we model it in terms of geodesic motion in the gravitational field of the merger remnant. This is a reasonable assumption given that the rest-mass densities past the detector are comparatively very small (i.e., ρ ≲ 10⁹ g cm⁻³), and hydrodynamical effects can be neglected. More importantly, and as we show in Section 2.4, the kinetic energy of the material detected decreases over the course of the postmerger evolution for all the systems studied. Thus, the material ejected at later times does not have the opportunity to collide with matter ejected earlier (potentially leading to shocks), so that we do not expect any interaction between the trajectories of the bound material.

A most important property of the bound matter is the time it will take to fall back onto the remnant. We could calculate this "fallback timescale" t_FB in terms of the geodesic equations for the matter collected at the detector (and indeed we have done so), but a much simpler and computationally less expensive estimate can be obtained using simple and analytic Newtonian expressions. In particular, matter with four-velocity u_t and specific angular momentum ℓ will have an orbit in the remnant's gravitational field with mass M_rem with specific energy

$$\begin{eqnarray}&&\epsilon :=-(1+{u}_{t})\lt 0,\end{eqnarray} \tag{ 1 }$$

orbital eccentricity e

$$\begin{eqnarray}&&e:=\sqrt{1+\displaystyle \frac{2\epsilon {{\ell }}^{2}}{{M}_{\mathrm{rem}}^{2}}},\end{eqnarray} \tag{ 2 }$$

and orbit semimajor axis a

$$\begin{eqnarray}&&a:=-\displaystyle \frac{{M}_{\mathrm{rem}}}{2\epsilon }.\end{eqnarray} \tag{ 3 }$$

In general, in its motion from the detector surface away from the merger remnant and back, the ejected matter does not trace a full orbit, and the duration of its freefall orbit is given by

$$\begin{eqnarray}&&{t}_{\mathrm{FB}}:=\displaystyle \frac{{a}^{2}{(1-{e}^{2})}^{2}}{{\ell }}{\int }_{{\theta }_{s}}^{2\pi -{\theta }_{s}}\displaystyle \frac{1}{{(1+e\cos \theta )}^{2}}d\theta ,\end{eqnarray} \tag{ 4 }$$

where the θ is the angle traced by the (planar) orbit and θ_s can be determined from a, e , and R_d (see, e.g., Rosswog 2007). This is the explicit solution to the Newtonian equation of motion in the gravity field of the central object.

However, we have the particular aim here of describing fallback flows on timescales of tens of seconds, much larger than the dynamical timescale around the central object. In this case, the duration of the fallback orbit is almost that of the full Newtonian orbit, whose period is

$$\begin{eqnarray}&&{t}_{\mathrm{FB}}=\displaystyle \frac{\pi }{\sqrt{2}}{M}_{\mathrm{rem}}| \epsilon {| }^{-3/2},\end{eqnarray} \tag{ 5 }$$

which is the well-known result of Newton's two-body problem.

Equations (4) and (5) ultimately represent an approximation of the actual fallback time that would be computed for the geodesic orbit of a massive test particle having the same energy and angular momentum moving in a Kerr spacetime with mass M_rem and dimensionless spin χ_rem. In Appendix A.2, we report a comparison between the fallback times computed with the two approaches. However, we can anticipate here that for long-period orbits, i.e., orbits with t_FB ≳ 1 s, the differences between the Kerr geodetic estimate and the Newtonian estimate are always ≲0.2%, which is not surprising given that the orbits are mostly in a weak-field region where r/M_rem ≳ 200. Such an error is well below other uncertainties in our modeling, and given that Equation (4) comes with a considerable computational gain, it has been set as the method of choice in our analysis. Note that, by construction, we cannot monitor bound material with orbits that are totally contained within the detector 2-sphere. Stated differently, the fallback times we compute are necessarily larger than

$$\begin{eqnarray}&&{t}_{\mathrm{FB},\ \min }\gtrsim 0.1\,\sqrt{{\left(\displaystyle \frac{{R}_{d}}{300\,{M}_{\odot }}\right)}^{3}\left(\displaystyle \frac{2.7\,{M}_{\odot }}{{M}_{\mathrm{rem}}}\right)}\quad {\rm{s}},\end{eqnarray} \tag{ 6 }$$

which is much smaller than the typical timescales considered in our analysis but is useful to bear in mind when considering the discussion in Section 2.3.

What still needs to be determined is the mass M_rem of the remnant, which we obtain with two different methods providing comparable results. More specifically, in the case where a black hole forms within the simulation timescale, we determine the black hole mass using the quasi-local integral on the apparent horizon surface as determined by the vanishing of a null geodesic congruence expansion (Thornburg 1996; Dreyer et al. 2003). In the cases where the remnant is a metastable HMNS, we rely on the gravitational field measured at the detector surface. In particular, in a 3+1 decomposition of spacetime (see, e.g., Gourgoulhon 2007; Rezzolla & Zanotti 2013), the covariant tt metric component is given by

$$\begin{eqnarray}&&{g}_{{tt}}=-{\alpha }^{2}+{\beta }_{i}{\beta }^{i},\end{eqnarray} \tag{ 7 }$$

where α and βⁱ are the lapse function and shift vector, respectively. On the other hand, in the weak-field approximation of general relativity, such a metric function at the detector surface ⁴ can be related to the Newtonian gravitational potential Φ as

$$\begin{eqnarray}&&{g}_{{tt}}\sim -1-2{\rm{\Phi }}=\displaystyle \frac{2{M}_{\mathrm{rem}}}{{R}_{d}}-1,\end{eqnarray} \tag{ 8 }$$

so that

$$\begin{eqnarray}&&{M}_{\mathrm{rem}}\simeq \displaystyle \frac{{R}_{d}}{2}\left({g}_{{tt}}+1\right).\end{eqnarray} \tag{ 9 }$$

Equation (9) thus provides an independent measurement of the mass relevant for the fallback dynamics, including, in the case of black hole formation, the eventual accretion disk, though it is negligible in the cases studied here. The gravitational mass measured via Equation (9) in this method agrees almost exactly with the Arnowitt–Deser–Misner (ADM) mass in the simulation domain. By applying these two methods, we find their results to agree to within 2% for the two binaries where a black hole forms (see Table 1). For convenience, since it can be used equally well in the case of a black hole or HMNS, we will use the remnant mass computed via Equation (9) for all of the binary systems.

As a concluding remark, we note that to avoid oversampling the polar region of our spherical detectors while only having a few cells near the equator, as would be the case when using spherical coordinates to represent the discrete grid on the 2-sphere, we developed a code where the surface of the detector is discretized following the HEALpix scheme (Górski et al. 2005). ⁵ In this coordinate system, all the cells have the same surface area, and a lower overall resolution is sufficient to capture the relevant characteristics of the ejected matter as it crosses the detector.

In Figure 1, we present snapshots at representative times of the properties of the material crossing the detector sphere for the two SFHO binaries. For each column, the top hemispheres report the instantaneous rest-mass flux per unit solid angle $d{\dot{M}}_{* }/d{\rm{\Omega }}$, while the bottom hemispheres show the values of u_t , with red (blue) shading marking matter that is gravitationally bound (unbound). The color bar for u_t was chosen to easily distinguish the bound material in red from the unbound material in blue. Above each panel, we report the times of the snapshots relative to the merger time, which, as is customary, is defined as the coordinate time at which the amplitude of the ℓ = m = 2 component of the gravitational-wave strain reaches its first maximum. In Figure 2, we present the same data for the two TNTYST binaries.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the properties of the ejected material in BNS mergers have been discussed in a number of papers (see, e.g., Sekiguchi et al. 2015; Bovard et al. 2017; Most et al. 2021; Camilletti et al. 2022; Papenfort et al. 2022), Figures 1 and 2 show rather clearly a feature that has not been remarked on sufficiently so far. More specifically, and irrespective of the mass ratio of the binary, the ejected material is initially mostly unbound (left column), and then the fraction of bound ejecta progressively increases (middle column) until a stage is reached where the ejecta is predominantly bound (right column). In addition, the geometrical distribution of the ejected matter is such that the front of the bound material progresses from the equatorial regions up to the polar directions, ending up by filling the whole sphere. Finally, while the early unbound ejecta is strongly anisotropic as a result of the tidal tails, the bound ejecta is quite isotropic. By comparing the upper and lower panels for each column, we find that the mass asymmetry in the binary increases the anisotropy in the polar direction of the early unbound ejecta (right column), as expected from its origin from the tidal tails, retaining such a property also at later times (right column). A more careful analysis of the polar distribution of the bound ejecta (not shown in Figures 1 and 2) reveals that the anisotropies are smaller than 15% for t − t_mer ≳ 20 ms; i.e., for the bound matter, the mass flux differs by less than 15% across the radial directions at these times. Therefore, the bound ejecta can be considered as essentially isotropic at late times.

By contrasting Figures 1 and 2, we can appreciate the impact of the EOS and in particular the role played by the stiffness (the SFHO EOS is stiffer than the TNTYST). More specifically, when considering a softer EOS in Figure 2, we find that the transition from unbound to bound is still present, however to a lesser extent, where significant unbound flux is found even at later times in the high-latitude regions. Similarly, the distribution of bound material at intermediate and late times (middle and right columns) is less isotropic, with a clear presence of unbound material at very high latitudes as a result of a collimated outflow in the case of the TNTYST EOS. The origin of the slightly different behavior is to be found in the different nature of the remnant, which is a black hole in the case of the SFHO EOS, while it is an HMNS in the case of the TNTYST EOS that is metastable at least over the time window of the simulations (see Table 1). Under these conditions, and as shown in a number of recent works (see, e.g., Fujibayashi et al. 2017, 2023), the HMNS is responsible for a neutrino-driven, collimated, and mildly relativistic matter outflow, which is clearly visible in the blue polar regions in the middle and right columns of Figure 2. Clearly, a similar behavior can be found when analyzing the neutrino fluxes in the two sets of binaries. Overall, what Figures 1 and 2 help appreciate rather transparently is that the ejection of unbound material is active as long as a stable merger remnant is present and is confined in the high-latitude region. However, such unbound outflow is rapidly shut down as soon as a black hole is formed, and when this happens, the ejected matter is mostly bound and with an essentially isotropic distribution. These considerations apply in the same manner to equal- and unequal-mass binaries and quite independently of the EOS, at least when considering the different nature of the remnant.

Figure 3 provides a more quantitative history of the mass ejection at the detector of our four binaries as measured, distinguishing the total (black solid line) rest-mass flux from the flux of bound (red solid line) and unbound material (blue solid line). The gray dashed lines report the times of the three snapshots presented in Figures 1 and 2, thus allowing a clear identification over the ejection history and highlighting that they correspond respectively to the maxima of the total and bound rest-mass fluxes and of the late-time behavior. Furthermore, in the case of the SFHO binaries (top row), vertical black dotted lines mark the time of black hole formation.

Zoom In Zoom Out Reset image size

Figure 3 also allows one to appreciate in great detail the general behaviors discussed above and relative to Figures 1 and 2. More precisely, it is possible to note that for the two SFHO binaries (top row), a sharp transition from unbound to bound ejection appears at about 10 ms postmerger, thus pointing out that the episode of unbound ejection is rather short-lived and restricted to the phase of dynamical ejection. Soon after the remnant collapses to a black hole at t − t_mer ∼ 4 ms, in fact, the rest-mass flux drops significantly, leaving only a small disk with mass ratio M_disk/M_rem ≃ 0.02 that will be responsible for the subsequent secular mass ejection. Indeed, this is a very small disk, and this is most likely due to the fact that in these binaries, the black hole is produced very rapidly after merger, and the angular momentum redistribution in the merger remnant has not taken place yet, hence most of the rest-mass is on unstable circular orbits when the black hole forms and will be rapidly accreted. We expect that binaries with the same SFHO EOS but slightly smaller chirp mass would lead to a longer equilibrium of the HMNS and hence to larger disks once the black hole is formed.

Note that in our two SFHO binaries, the disk produced is very small due to the prompt collapse of the remnant. We conjecture that this small mass inhibits the ejection of unbound material from the disk, which would be due to MHD processes or neutrino-driven processes (which in any case would be very subdominant; see comparison in Gill et al. 2019) in a larger disk that could sustain a higher pressure, as seen, e.g., by Fahlman & Fernández (2018) and Fujibayashi et al. (2018). The low mass of the disk in our SFHO cases only allows for bound ejecta, as shown in Figure 3. Note also that the transition from unbound to bound ejecta suggests that the fallback material will always remain below the unbound matter, with a somewhat precise geometrical separation. Indeed, inspecting the velocities of the bound and unbound components and extrapolating in time leads us to conclude that the two components will remain distinct in the subsequent evolution. This behavior has important consequences that we will discuss in more detail in Section 4.

Interestingly, this interplay between the mass of the disk and the ratio of the bound to unbound ejected material may suggest an anticorrelation between the energy in the prompt GRB emission and that in the extended emission. Indeed, one would expect more massive disks to be able to sustain a powerful jet for a longer time (before the disk either depletes or enters a magnetically arrested state, leading in both cases to a decrease in the jet power; Gottlieb et al. 2023) and thus to a larger energy release in the prompt emission. On the other hand, a less massive disk would lead to a comparatively larger amount of bound material being ejected and therefore to a larger energy reservoir for the extended emission when compared to the prompt GRB signal.

When considering instead the evolution of the ejected material for the TNTYST binaries (bottom row in Figure 3), it is possible to recognize the same features discussed for the SFHO binaries, with the important difference that the ejection of unbound matter does not drop rapidly since no black hole is formed in this case during the simulation. As a result, the HMNS is able to eject, especially in the polar region, a neutrino-driven, collimated, and mildly relativistic outflow of unbound matter. Such an unbound flow, however, is essentially always smaller than that of the bound material, which represents the dominant contribution to the matter crossing the detector. Note also that the bound flow does not exhibit a monotonically decreasing evolution with time and has episodes in which its rate increases, although by a factor of a few only (see, e.g., t − t_mer ∼ 35 ms for the TNTYST-q1.0 binary and t − t_mer ∼ 35 ms for the TNTYST-q.75 binary). Given the very complex dynamics of the HMNS, it is difficult to disentangle the origin of these "rebrightening" episodes, which we attribute to the material being ejected after being shock-heated by the bound material that in the meantime has fallen back onto the HMNS. As we will show later, the fallback time of this late-time bound ejecta is very short and hence will have no impact on the EM signatures of the binary. The situation would be rather different if this rebrightening material had longer fallback times, as we discuss in Appendix A.1.

Finally, Table 1 reports not only the total masses ejected but also the corresponding components that are in bound and unbound matter. Note that quite independently of the EOS, unequal-mass binaries tend to have larger values of ejected mass; on the other hand, the softer TNTYST EOS favors bound ejection by almost a factor of 2 in the equal-mass case. We have also estimated that, quite generically, about 10⁻⁶ M_⊙ of bound ejecta has a fallback time larger than 100 s; this mass increases by a factor of about 5 if we consider a fallback time larger than 10 s.

We next concentrate on illustrating the general characteristics of the significant amount of fallback inflow and discuss how it can lead to long-term EM emission. We start by presenting in Figure 4 the distribution of the orbital specific energy := −(1 + u₁) in the fallback matter for our four binaries. Since = 0 corresponds to marginally bound matter, we obviously consider the distributions of bound matter only for < 0, the with lower bounds ( ≃ −0.0045) being set by the tightest orbits we can measure at the detector (see the discussion around Equation (6)).

Zoom In Zoom Out Reset image size

Remarkably, these distributions are essentially flat near the ≲ 0 bound, and they remain nearly constant over the entire range of specific energies for the SFHO binaries (blue and green solid lines). The same can be said for the binaries with the TNTYST EOS (red and purple solid lines in Figure 4). The fact that these distributions are essentially flat has important implications since, as already discussed early on in the study of fallback of debris from TDEs (see, e.g., Rees 1988; Evans & Kochanek 1989; Phinney 1989), in this case, it is not difficult to show that the fallback accretion rate ${\dot{m}}_{\mathrm{FB}}$ must follow a power law in time; namely, using Equation (5), it is trivial to obtain that

$$\begin{eqnarray}&&{\dot{m}}_{\mathrm{FB}}:=\displaystyle \frac{{{dm}}_{\mathrm{FB}}}{{{dt}}_{\mathrm{FB}}}=\displaystyle \frac{{{dm}}_{\mathrm{FB}}}{d\epsilon }\displaystyle \frac{d\epsilon }{{{dt}}_{\mathrm{FB}}}\propto \displaystyle \frac{{{dm}}_{\mathrm{FB}}}{d\epsilon }\times {t}_{\mathrm{FB}}^{-5/3}\propto {t}_{\mathrm{FB}}^{-5/3}.\end{eqnarray} \tag{ 10 }$$

The importance of the result of Equation (10) stems not only from the very simple hypotheses employed to obtain it but also because it has been indirectly confirmed by some late-time observations of TDEs (see, e.g., Komossa 2015 for a review).

Interestingly, the power-law behavior predicted by Equation (10) is exhibited also by our fallback material, which obviously includes, in addition to the tidal ejecta, the bound material coming from the neutrino-driven and magnetically driven winds from the merger remnants. More specifically, using the distribution of fallback times, we can compute the fallback accretion rate as a function of time and report the results in Figure 5, where we also plot a t^−5/3 power-law decay to guide the eye. As anticipated, for all the cases considered, and hence rather independently of the EOS or the mass ratio, the fallback accretion proceeds with a power-law decay ${\dot{m}}_{\mathrm{FB}}\propto {t}^{-5/3}$ (see black dashed line). This result, which was already known for neutron star–black hole binaries, either in Newtonian calculations (e.g., SPH Rosswog 2007) or in full general relativity (Chawla et al. 2010; Kyutoku et al. 2015), is now also confirmed from rather generic conditions of BNS mergers. Figure 5 also highlights that the power-law result for the fallback accretion rate depends only very weakly on the specific-energy distributions. As shown in Figure 4, in fact, the specific-energy distributions of our TNTYST binaries are only approximately constant, showing instead an excess (deficit) at lower (higher) specific energies. Yet, this does not affect the ∼t^−5/3 behavior of the fallback accretion rate of the TNTYST binaries. We believe this is because material with low specific energy will affect only the onset of the fallback process; however, the long-term features of the fallback accretion are essentially determined by the matter that is only marginally bound and hence with ≲ 0. Because all of the distributions in Figure 5 are essentially constant at ≲ 0, it is not surprising that all of our binaries reproduce the t^−5/3 behavior irrespective of the EOS, mass ratio, or nature of the merger remnant. Note also from Figure 5 that the measured accretion rates in the range t_FB ∼ 10²–10³ s are of the order ${\dot{m}}_{\mathrm{FB}}\sim {10}^{-9}\mbox{--}{10}^{-7}\,{M}_{\odot }\,{{\rm{s}}}^{-1}$. As a result, using a crude conversion of accretion rate to luminosity, i.e., ${L}_{\mathrm{acc}}\sim 0.1\times {\dot{m}}_{\mathrm{FB}}{c}^{2}$, it is straightforward to obtain a very rough estimate of the luminosities produced by this inflow, namely, L ≃ 10⁴⁴ − 10⁴⁶ erg s⁻¹. In Section 3, we will confirm and refine this result using a more detailed EM model.

Zoom In Zoom Out Reset image size

The robustness of the result presented in Figure 5 is also underlined by the two panels in Figure 6, where we show the same distributions as in Figure 4 for the SFHO binaries but when computed at detectors that are at either larger (i.e., 400 M_⊙) or smaller (i.e., 250 M_⊙) radii. In these cases, the lower bounds in the specific-energy distributions obviously change, increasing (decreasing) as the detector radius is decreased (increased) simply because of the differences in the tightest orbits that the detectors can capture. Clearly, the overall distributions do not show significant changes with the location of the detector, and, more importantly, they are essentially identical at the highest specific energies, i.e., at ≃ 0. As a result, following the logic discussed above, the corresponding fallback accretion rates all yield the same ∼t^−5/3 power-law behavior (not shown here). Furthermore, since the specific-energy distributions all agree around ≃ 0, we can conclude that the conversion of internal energy into kinetic energy is not unbinding material significantly between the different radii, hence justifying our use of u_t over hu_t to determine the bound material.

Zoom In Zoom Out Reset image size

Finally, in Figure 7, we present the distributions in time and fallback time of the matter accreted onto the remnant for the four binaries considered in our analysis. Note that in the first 10 ms after merger, all binaries exhibit large ejections of bound matter and with a broad distribution of fallback times, i.e., from 10⁻¹ s up to 10⁴ s, although most of the bound mass has t_FB ≲ 10² s. As the evolution proceeds, however, differences brought in by the EOS and the mass ratio start to emerge. In particular, in the case of the stiffer SFHO binaries and equal-mass binaries, the bound ejecta are essentially suppressed after ≲20 ms (see also top left panel of Figure 3), while they continue to be present for the unequal-mass case but with a much reduced variance in fallback times, so that t_FB ≲ 10² s after about 20 ms after merger for SFHO-q.75. On the other hand, in the case of the softer TNTYST binaries, large amounts of bound ejecta are produced for much longer timescales (i.e., t − t_mer ≲ 40 ms) and with a broad distribution of fallback times. The main difference in this case is to be found in the mass ratio, with the ejection being essentially shut down for t − t_mer ≲ 40 ms in the case of the equal-mass binary while being persistent for the unequal-mass TNTYST-q.75 binary. Interestingly, for the latter, the distribution of fallback times reduces significantly at late times (i.e., t_FB ≲ 10 s for t − t_mer ≳ 50 ms), most likely because of the cooling of the HMNS, which reduces the specific energy of the ejecta in the neutrino and magnetically driven winds (see Equation (5)). Nevertheless, one could expect that by running the simulation further, provided that the HMNS remains stable against gravitational collapse, a secondary, collimated, and magnetically driven wind would emerge, as shown by several recent works (Kiuchi et al. 2024; Combi & Siegel 2023). The fact that this wind does not appear in our simulation time might depend on several factors. First, the physical mechanism at the basis of this ejection, while argued to be likely tied to an α − Ω dynamo, is still poorly understood. Second, the resolution of our simulations is too low to consistently resolve the magnetorotational instability in the outer layer of the massive star, as well as in the high-density and high-temperature regions of the disk, therefore delaying the onset of the dynamo and resulting wind.

Zoom In Zoom Out Reset image size

In Section 2, we showed that significant fallback accretion of material can occur on long timescales, i.e., for t_FB ≳ 100 s. The fallback rates reported in Figure 5 show that this accretion is largely super-Eddington, such that the conversion of the mass accretion rate ${\dot{m}}_{\mathrm{FB}}$ to the accretion luminosity L_acc is not straightforward and requires a more precise modeling than the crude one we have employed in the previous section. Furthermore, the fallback material, just like the unbound material, is subject to nucleosynthetic processes that contribute to the generation of EM radiation that needs to be taken properly into account. However, because the dynamics of the outflow–inflow cannot be followed on such timescales by our simulations, a semianalytical model is needed in order to determine such a dynamics, from which a model for the EM emission can be derived.

As indicated by the basic principles of orbital dynamics (see Equations (2) and (5)) and by the inspection of the numerical data from the simulations, material with long fallback times, e.g., t_FB ≳ 100 s, has very eccentric orbits. Furthermore, as noted in Section 2.3, the ejection of bound material is largely isotropic, such that the fallback dynamics can also be considered independent of θ and ϕ. Given these considerations, it is reasonable to simplify our semianalytic model by considering the inflow–outflow dynamics to be purely radial. In this case, the properties of the fallback inflow reduce to the rest-mass density and to the radial velocity fields ρ and v that are functions only of time and radius (r, t). However, since the electron fraction of the material does depend on the latitude, we will account for this in the determination of the fallback radiation in Section 3.2.

To determine the dynamics of the infalling material, we consider it to be in radial freefall. In addition, we consider all the fallback matter to have been ejected at the same time, which marks t = 0 in our study of the fallback dynamics. This assumption corresponds to neglecting the duration of the ejection episode in the merger, which is Δt_ej ≲ 100 ms (see Figure 7), which is valid for the fallback timescales considered. A fluid element from freefalling from an initial radius R_i with zero velocity to a radius R_f will have a Newtonian freefall time given by

$$\begin{eqnarray}&&{t}_{\mathrm{FF}}({R}_{i},{R}_{f})=\sqrt{\displaystyle \frac{{R}_{i}^{3}}{2{GM}}}\left(\sqrt{\displaystyle \frac{{R}_{f}}{{R}_{i}}\left(1+\displaystyle \frac{{R}_{f}}{{R}_{i}}\right)}+\arccos \sqrt{\displaystyle \frac{{R}_{f}}{{R}_{i}}}\right),\end{eqnarray} \tag{ 11 }$$

and its velocity at the final radius R_f is simply

$$\begin{eqnarray}&&{v}_{\mathrm{FF}}({R}_{i},{R}_{f})=\sqrt{2{GM}\left(\displaystyle \frac{1}{{R}_{f}}-\displaystyle \frac{1}{{R}_{i}}\right)}.\end{eqnarray} \tag{ 12 }$$

Consider now a radial segment of matter of size d L₀ initially at radius R_i . After freefalling to radius R_f < R_i , this radial segment will have spread out to a size dL_f such that

$$\begin{eqnarray}&&{{dL}}_{f}-{{dL}}_{0}={v}_{\mathrm{FF}}({R}_{i},{R}_{f})\displaystyle \frac{\partial {t}_{\mathrm{FF}}}{\partial {R}_{i}}({R}_{i},{R}_{f}){{dL}}_{0},\end{eqnarray} \tag{ 13 }$$

since the trailing matter will take time t_FF(R_i , R_f ) − t_FF(R_i − dL₀, R_f ) to catch up to the leading matter, with velocity v_FF(R_i , R_f ) (note that we here assume that the velocity vanishes at both edges simultaneously). By conservation of matter in the freefall, the densities ρ_i and ρ_f before and after the freefall are related by

$$\begin{eqnarray}&&{\rho }_{i}{R}_{i}^{2}{{dL}}_{0}={\rho }_{f}{R}_{f}^{2}{{dL}}_{f},\end{eqnarray} \tag{ 14 }$$

or, using Equation (13),

$$\begin{eqnarray}&&{\rho }_{i}={\rho }_{f}{\left(\displaystyle \frac{{R}_{f}}{{R}_{i}}\right)}^{2}\left[1+{v}_{\mathrm{FF}}({R}_{i},{R}_{f})\displaystyle \frac{\partial {t}_{\mathrm{FF}}}{\partial {R}_{i}}({R}_{i},{R}_{f})\right].\end{eqnarray} \tag{ 15 }$$

Hence, after setting R_f = R_d and considering a generic initial radius $\bar{R}$ and time $\bar{t}$, Equation (15) can be written as

$$\begin{eqnarray}&&\rho (\bar{R},\bar{t})={\left(\displaystyle \frac{{R}_{d}}{\bar{R}}\right)}^{2}\rho ({R}_{d},\bar{t}+{t}_{\mathrm{FF}}(\bar{R},{R}_{d}))\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\times \left[1+{v}_{\mathrm{FF}}(\bar{R},{R}_{d})\displaystyle \frac{\partial {t}_{\mathrm{FF}}}{\partial {R}_{i}}(\bar{R},{R}_{d})\right],\end{eqnarray} \tag{ 17 }$$

which allows us to determine the density at an arbitrary radius $\bar{R}\gt {R}_{d}$ and time $\bar{t}\gt 0$ in the flow from its value at the detector R_d . As discussed in Section 2.3, we know that at R_d , the fallback accretion rate is

$$\begin{eqnarray}&&{\dot{m}}_{\mathrm{FB}}({t}_{\mathrm{FB}})=4\pi {R}_{d}^{2}\rho ({R}_{d},{t}_{\mathrm{FB}}){v}_{\mathrm{FF}}({R}_{\max }({t}_{\mathrm{FB}}),{R}_{d}),\end{eqnarray} \tag{ 18 }$$

where ${\dot{m}}_{\mathrm{FB}}({t}_{\mathrm{FB}})\propto {t}_{\mathrm{FB}}^{-5/3}$ (e.g., Figure 5) and ${R}_{\max }({t}_{\mathrm{FB}})$ is the maximum radial coordinate reached by the ejected bound material before falling back onto the remnant. In practice, ${R}_{\max }({t}_{\mathrm{FB}})$ can be computed from Equation (11) after requiring that

$$\begin{eqnarray}&&{t}_{\mathrm{FF}}({R}_{\max }({t}_{\mathrm{FB}}),{R}_{d})=\displaystyle \frac{{t}_{\mathrm{FB}}}{2}.\end{eqnarray} \tag{ 19 }$$

Note that in our semianalytic model, the density of fallback material is 0 for all (r, t) such that

$$\begin{eqnarray}&&r\geqslant {R}_{d}{\left(1+\displaystyle \frac{\sqrt{2{GM}}t}{{R}_{d}}\right)}^{2/3},\end{eqnarray} \tag{ 20 }$$

as the equality in Equation (20) corresponds to the orbit of material with zero energy at infinity. Above the leading edge separating the outgoing and ingoing fallback material (see dashed line in Figure 10), the layer of unbound material will produce the kilonova signal as it moves outward (this is marked with a dashed line in Figure 8). As mentioned in Section 2.3, the velocities of the bound and unbound components of the ejecta are well separated, such that we expect a gap between the flow of fallback material and the unbound outflow. In Section 4, we further discuss the consequences of this geometry on the observation of radiation from the fallback component.

Zoom In Zoom Out Reset image size

To solve for the density in the fallback inflow, we normalize the fallback rate on the detector surface ${\dot{m}}_{\mathrm{FB}}$ using the total mass of bound ejecta (Table 1) and then apply Equations (13), (17), (18), and (19). In this way, it is possible to build a spacetime diagram of the semianalytic model for the infalling matter reporting the world lines (i.e., the trajectories in spacetime) of shells of material with different rest-mass densities (hence with different colors) and that we report in Figure 8 for the representative SFHO-q1.0 binary, which has a total fallback mass of 1.60 × 10⁻³ M_⊙. Note that the spacetime is reported in a double logarithmic scale and hence timelike trajectories appear as spacelike, while the gray line near the vertical axis represents the world line of the detector's surface at R_d = 300 M_⊙, where the fallback accretion flow with ${\dot{m}}_{\mathrm{FB}}\propto {t}^{-5/3}$ is enforced (see Equation (18)). Also shown with white lines are representative fluid lines, i.e., radial freefall orbits with different fallback times, both during the outflow stage (positive slope) and during the inflow stage (negative slope). Clearly, these fluid lines intersect the trajectories of constant rest-mass density shells, since the fallback material is first decompressed in the outflow stage and then recompressed in the inflow part. Also apparent from Figure 8 is that the vast majority of the ejecta has a rest-mass density ρ ≲ 10⁻² g cm⁻³ as of 10 s after merger and that the bound matter with the largest but still negative energy has the longest fallback time. Finally, shown with an orange line in Figure 8 is the photosphere of the EM emission model that we discuss in the next section.

In order to determine the EM radiation arising from the fallback flow discussed in the previous sections, we start by locating the photosphere in the fallback flow. To this effect, we calculate the photospheric radius (photosphere) R_ph(t), defined such that

$$\begin{eqnarray}&&{\int }_{r={R}_{\mathrm{ph}}(t)}^{\infty }\rho (r,t)\kappa {dr}=1,\end{eqnarray} \tag{ 21 }$$

where κ is the opacity of the material. Note that R_ph(t) marks the position of the photosphere within the fallback flow and is different from that of the unbound and outward-moving material that we discuss in Section 4.

In the early stages of the postmerger that we have simulated, the r-process nucleosynthesis is already well underway (see, e.g., Lippuner & Roberts 2015), and the material is thus expected to already be opaque (with κ > 0.1 cm² g⁻¹; Tanaka et al. 2020), though exact opacity calculations at these early times are lacking in the literature. Depending on the initial electron fraction and thermodynamic properties of the material, the exact composition and opacity will change. In addition, the composition of the fallback column at different latitudes above the remnant will differ, according to the differential irradiation by neutrinos that our M1 scheme can capture and we will consider below.

Assuming a low opacity of κ = 0.5 cm² g⁻¹, we report with a solid orange line in Figure 8 the trajectory of the photosphere in the spacetime diagram. Clearly, over timescales t_FB ≲ 10⁴ s, the photosphere coincides with the leading edge of the bound ejecta, decoupling from it at later times and inverting its motion only at t_FB ≲ ∼ 10⁵ s. This behavior is found even when considering a value of the opacity as low as κ = 0.1 cm² g⁻¹, with the photosphere moving inward at times no earlier than 5 × 10⁴ s postmerger. As a result, when considering the larger opacity reached through the r-process nucleosynthesis, we conclude that the photosphere will always coincide with the edge of the ejecta, at least on these timescales. Furthermore, the radiation emitted from the ejecta on these timescales and coming from the leading edge of the bound material can be modeled as having a blackbody emission spectrum with a time-varying temperature.

To capture the temperature evolution of this material in time and hence account for r-process nucleosynthesis and expansion cooling, we use the SkyNet nuclear reaction network (Lippuner & Roberts 2017). More specifically, from our GRMHD simulations, we determine the temperature T and the electron fraction Y_e of the material that is least bound, i.e., the material with the largest u_t among the bound matter crossing the detector at R_d . This leading-edge material in the outflow is the source of the radiation of interest, up to the recession of the photosphere. Using our semianalytical model, we extract the density history at the leading edge ρ(R_ph(t), t) and determine the initial nuclear statistical equilibrium at the initial density, temperature, and Y_e using the nuclei libraries for strong reactions, weak reactions, symmetric fission, and spontaneous fission available in SkyNet. Once this information is provided as initial conditions, we use SkyNet's reaction network to obtain the expected temperature history T_ph(t) at the photosphere.

In principle, nucleosynthesis starts from the unbinding of the material from the merger remnant and is active during the propagation of matter up to detector R_d . This can be easily appreciated after considering that the crossing time to the detector is ∼1 ms, which is larger than the neutron-capture timescale

$$\begin{eqnarray}&&{\tau }_{n\,\mathrm{cap}}\gtrsim \displaystyle \frac{1}{v\,{n}_{n}{\sigma }_{n,\mathrm{cap}}}={10}^{-8}\,{\rm{s}},\end{eqnarray} \tag{ 22 }$$

where we used a neutron number density of n_n = ρ/m_n = 10²⁷ cm⁻³ corresponding to a rest-mass density ρ = 10⁴ g cm⁻³ in the early stages of the flow expansion, a thermal velocity $v=\sqrt{2{k}_{{\rm{B}}}T/{m}_{n}}$ with temperature T = 10⁹ K, and a cross section of neutron capture onto protons of σ_n,cap = 10⁻² fm² (Kopecky et al. 1997). In practice, since the detector-crossing time is much smaller than the overall timescale involved in the nucleosynthetic process, ⁶ it is reasonable to ignore this initial contribution, but more refined calculations involving tracer particles (see, e.g., Bovard & Rezzolla 2017) will help validate this assumption.

Knowing R_ph(t) and T_ph(t), we can determine the blackbody luminosity in any spectral band ${ \mathcal B }$ from

$$\begin{eqnarray}&&{L}_{\mathrm{FB}}{(t)=8\pi {R}_{\mathrm{ph}}(t)}^{2}{\int }_{\nu \in { \mathcal B }}{B}_{\nu }({T}_{\mathrm{ph}}(t))d\nu ,\end{eqnarray} \tag{ 23 }$$

where ν is the photon frequency and B_ν is the Planck function, and we have clearly assumed the emission to be isotropic in view of the spherical symmetry of the underlying model.

In our picture, the r-process heats up the photosphere, from which we determine the radiation from the fallback flow. This approach neglects the energy diffusion inside the ejecta, and to judge whether this is reasonable, we recall that the photon diffusion timescale inside the infalling flow can be estimated as (Metzger 2017)

$$\begin{eqnarray}&&{\tau }_{\mathrm{diff}}=\kappa {R}_{\mathrm{ph}}^{2}{\rho }_{\mathrm{ph}},\end{eqnarray} \tag{ 24 }$$

with ρ_ph the rest-mass density at the photosphere. Using an opacity of κ = 1 cm² g⁻¹ and data read off Figure 8, this timescale decreases from τ_diff ∼ 10⁷ s at t − t_mer = 10 s postmerger to τ_diff ∼ 10⁴ s at t − t_mer = 10⁴ s; obviously, the photon diffusion timescale would be even higher with a larger opacity. Overall, these considerations show that τ_diff exceeds the timescales over which we study the radiation (i.e., ≲10⁴ s), such that in our simplified model, we can neglect the thermalization in the outflow and consider the temperature history from the nuclear reaction network output as the photosphere temperature.

In Figure 9, we plot the bolometric luminosity and the luminosity integrated over the Swift/XRT and Swift Burst Alert Telescope (BAT) bands, namely, (0.3–10) keV and (15–350) keV, respectively (Gehrels et al. 2004; Barthelmy et al. 2005b). All the curves in these plots correspond to times before the recession of the photosphere. In the different panels, we consider the fallback inflow deduced from our four GRMHD simulations and with fallback material picked from different latitudes θ in the inflow, as marked in the panels. The initial value of the electron Y_e of that material, also marked in each plot, is determined from our M1 scheme. The dotted line marks the bolometric luminosity from the fallback accretion that one would deduce from a naive conversion from the fallback rate ${L}_{\mathrm{acc}}=\eta {\dot{m}}_{\mathrm{FB}}{c}^{2}$, with an efficiency of η = 10% and ${\dot{m}}_{\mathrm{FB}}$ as in Figure 5 (see discussion at the end of Section 2.4).

Zoom In Zoom Out Reset image size

When inspecting the various panels in Figure 9, a number of considerations follow. First, the duration and luminosity of the fallback radiation are significant and comparable with the observed extended emission in short GRBs (see Sections 1 and 4). Second, through a nontrivial interplay between the photospheric radius and the temperature evolution, the bolometric luminosity displays a power-law decay with a nearly constant slope α ∼ −8/3, which is steeper than the decay of ${\dot{m}}_{\mathrm{FB}}\propto {t}^{-5/3}$. Inspecting the temperature evolution from the nuclear reaction network reveals that the temperature evolves as T ∝ t⁻¹, and a heating rate that also has a power-law behavior as _nuc.(t) ∝ t^−3/2, in agreement with other studies of r-process nucleosynthesis (Wanajo et al. 2014; Lippuner & Roberts 2015).

Because of the temperature decrease with time, the radiation is first mainly output in the gamma-ray bands before dominating the X-ray band when the Planck function's maximum crosses the band. At very early times, i.e., for t < 1 s, the Swift/BAT band is below the blackbody peak temperature T ≳ 100 KeV, such that the Rayleigh regime applies with L_γ ∼ R² T ∝ t^1/3. This near-constant segment of the light curve is outside of the plot in Figure 9, as it is not expected to be visible due to the prompt emission outshining this component. However, an analogous segment is seen in the X-rays, which lasts until tens of seconds, with a near-flat X-ray light curve. Once the blackbody temperature reaches the Swift/BAT band, at t ∼ 1 s, the gamma-ray luminosity follows the bolometric behavior with L_γ ∝^−8/3. This emission lasts until the blackbody peak frequency leaves the Swift/BAT band at nearly 100 s. After that, an exponential cutoff is seen in the gamma rays, leading to the t^−8/3 segment in the X-rays up to 1000 s, followed by the exponential cutoff. In general, the sequence "t^1/3" plus "t^−8/3" plus "exponential cutoff" is present in all bands; however, the earlier segment in the gamma rays and the later segment in the X-rays are expected to be outshone by the prompt and the forward-shock afterglow, respectively.

Since the photospheric radius at early times is determined only by the orbit with u_t = 0 and is therefore essentially the same for all binaries considered, differences in luminosity between binaries and different latitudes in those binaries (Figure 9) are simply due to the initial thermodynamic conditions and Y_e . For example, when comparing the luminosity curves at different latitudes within the same binary (e.g., the two top panels in Figure 9), it is possible to note the effect of a higher initial Y_e : the binary with a lower Y_e experiences a stronger heating rate and a brighter emission. On the other hand, when comparing binaries with similar initial Y_e and different total fallback masses (e.g., the two right panels in Figure 9), it is possible to appreciate the impact of the history of the rest-mass density: the overall luminosity is higher for the higher-mass system.

After the photosphere recedes into the ejecta, our simple emission model no longer applies, as the whole ejecta then shines. The mass of emitting material then follows m_shine ∝ t^−2/3, and, applying a similar temperature evolution T ∝ t⁻¹ to the whole ejecta, one can expect a leveling-out of the luminosity with a shallower decay than in Figure 9. This segment, however, is likely to be unobservable because of the forward-shock afterglow.

The analysis presented here is inevitably limited by the duration of our simulations. While the two SFHO binaries have only a minute mass ejection by the time our simulations end at t ≃ 50 ms (Figure 3), this is not the case for the TNTYST binaries, which do not lead to a black hole, and where the ejection is nonzero at the time the simulations are terminated. Although Figure 7 clearly suggests that the late ejecta has very short fallback timescales, we cannot conclude that ejection on longer timescales than the ones we capture in our simulations would not contribute to longer fallback times. Indeed, if this late-time ejection had longer fallback times, the EM signature would likely be different.

To explore this possibility, we virtually extend the duration of our simulations by also considering ejection of matter at times past the end of our simulations. We do this by replicating the ejection of matter from the TNTYST-q1.0 simulation and relative to a data segment spanning a window in time of 10 ms. This segment is shown with a gray shaded area in the bottom left panel of Figure 7 and refers to an ejection of matter with almost steady-state fallback properties (see also Figure 5). We then replicate this data segment a number of times, i.e., 3, 30, or 300 times, hence mimicking a simulation that would run for ≃0.053, ⁷ 0.3, and 3.0 s, respectively. We then determine the corresponding fallback mass accretion rate and show the results in Figure 11, where different solid colored lines refer to different ejection times. By contrasting Figures 5 and 11, it is clear that the "extra" ejection mimicked here has the effect of introducing a nearly flat accretion segment that lasts for the duration of the extra ejection before connecting to the t^−5/3 falloff for the subsequent times and as described in Section 2.3. This mass accretion stage would lead to a similar phase of constant observed luminosity. Thus, the observational signatures of the bound material are not tied to the timescales of the extended emission but can be prolonged to longer emission components, such as some of the early plateaus in the X-ray afterglows of GRBs (Nousek et al. 2006).

Zoom In Zoom Out Reset image size

As a word of caution, we should remark that when modeling ejection timescales that are much larger than the one considered in the simulations (i.e., ${ \mathcal O }(100\,\mathrm{ms})$), and indeed of the order of seconds as considered in Figure 11, it would be necessary to also account for the interaction of the late-time ejected matter with the infalling matter. The two components could collide, leading to a reverse shock and potentially a new component of the fallback luminosity. The much richer picture that could follow from this additional component is worth exploring in future works.

We next compare our method of determining the fallback time of bound material using a Newtonian description for the fallback orbits (t_FB,N; see Equations (1)–(5)) with a relativistic treatment using the integration of timelike geodesics representing the orbits of massive test particles in a Kerr spacetime.

Taking as a reference the SFHO-q1.0 binary, we perform a quasi-local measurement (Thornburg 1996; Dreyer et al. 2003) of the mass M_BH and dimensionless spin χ_BH of the black hole that forms as a result of the collapse of the HMNS. Since the mass of the torus around the black hole is a few percent at most of the mass of the black hole (Rezzolla et al. 2010), we can ignore it and hence use M_BH and χ_BH to build a Kerr spacetime that we cover with Cartesian Kerr–Schild coordinates having g _KS as a metric tensor (see, e.g., Kerr & Schild 2009). We then collect fluid elements on our detector sphere at R_d = 300 M_⊙ and read off their energy u_t and the local spatial components u_i , which we then rescale for the different coordinate system so as to guarantee that the four-velocity is a unit timelike vector, i.e., that ${({g}_{\alpha \beta })}_{\mathrm{KS}}{u}^{\alpha }{u}^{\beta }=-1$. Using these initial conditions, we integrate the corresponding geodesic equations in Kerr–Schild coordinates until the material reaches the detector surface again, which we set at R_BL = 300 M_⊙, where R_BL is the radial Boyer–Lindquist coordinate (Boyer & Lindquist 1967). In this way, we can determine the corresponding fallback time t_FB as the time coordinate when this terminal condition is met. We apply this procedure to the ∼70,000 fluid elements crossing our 2-sphere over a single time level and find that they span the range in specific energy ∈ [−1 × 10⁻², −3 × 10⁻⁵], which is representative of the material relevant for our study.

Figure 12 reports the relative difference between the general-relativistic and Newtonian fallback times for these particles, i.e., t_{FB, N} and t_{FB, N}, respectively. Naturally, the difference is smaller for longer orbit times that spend more time far away from the black hole, as differences in the gravitational fields are much smaller far from the black hole. Overall, we find that they differ by less than 0.2% on the timescales relevant for our study (i.e., t_FB,N > 1 s), thus fully justifying our use of the Newtonian estimates throughout.

Zoom In Zoom Out Reset image size