# Shocks Power Tidal Disruption Events

Taeho Ryu Max Planck Institute for Astrophysics, Karl-Schwarzschild-Strasse 1,
85748 Garching, Germany  Physics and Astronomy Department, Johns Hopkins
University, Baltimore, MD 21218, USA  Julian Krolik  Physics and Astronomy
Department, Johns Hopkins University, Baltimore, MD 21218, USA  Tsvi Piran
Racah Institute of Physics, Hebrew University, Jerusalem 91904, Israel  Scott
C. Noble  Gravitational Astrophysics Laboratory, NASA Goddard Space Flight
Center, Greenbelt, MD 20771, USA  Mark Avara  Institute of Astronomy,
University of Cambridge, Madingley Road, Cambridgeshire, England

###### Abstract

Accretion of debris seems to be the natural mechanism to power the radiation
emitted during a tidal disruption event (TDE), in which a supermassive black
hole tears apart a star. However, this requires the prompt formation of a
compact accretion disk. Here, using a fully relativistic global simulation for
the long-term evolution of debris in a TDE with realistic initial conditions,
we show that at most a tiny fraction of the bound mass enters such a disk on
the timescale of observed flares. To “circularize” most of the bound mass
entails an increase in the binding energy of that mass by a factor $\sim 30$;
we find at most an order unity change. Our simulation suggests it would take a
time scale comparable to a few tens of the characteristic mass fallback time
to dissipate enough energy for “circularization”. Instead, the bound debris
forms an extended eccentric accretion flow with eccentricity $\simeq 0.4-0.5$
by $\sim 2$ fallback times. Although the energy dissipated in shocks in this
large-scale flow is much smaller than the “circularization” energy, it matches
the observed radiated energy very well. Nonetheless, the impact of shocks is
not strong enough to unbind initially bound debris into an outflow.

black hole physics $-$ gravitation $-$ hydrodynamics $-$ galaxies:nuclei $-$
stars: stellar dynamics

††journal: ApJ††software:  python (Van Rossum & Drake, 2009); matplotlib
(Hunter, 2007); MESA (Paxton et al., 2011); Harm3d (Noble et al., 2009).

##  1 Introduction

A tidal disruption event (TDE) takes place when a star that has wondered into
the vicinity of a supermassive black hole (SMBH) is torn apart by the SMBH’s
gravitational field. About half of the stellar mass is left unbound and is
ejected to infinity while the other (bound) half returns to the vicinity of
the black hole. The end result is a flare of optical/UV, X-ray, and at times,
radio emission.

During the last decades, TDEs were transformed from a theoretical prediction
(Hills, 1988; Rees, 1988) to an observational reality (Gezari, 2021). With the
use of numerous detectors, ranging from ROSAT (and now eROSITA) in X-rays, to
GALEX in the UV, and multiple telescopes in the optical band (including
systematic surveys like SDSS (Blanton et al., 2017), (Blanton et al., 2017),
ASSAS-N (Shappee et al., 2014), Pan-STARRS (Kaiser et al., 2002), and, more
recently, ZTF (Bellm et al., 2019), more than a hundred TDEs have now been
detected. Upcoming observations with the Rubin Observatory will provide an
overwhelming amount of data in the near future.

Although TDEs are of great interest in their own right, their properties offer
a wide range of opportunities to learn about other astrophysical questions.
They can reveal quiescent SMBHs and possibly permit inference of their masses.
They present otherwise unobtainable information about non-steady accretion
onto SMBHs and the conditions for jet-launching. In addition, understanding
the rates of TDEs would reveal valuable information concerning the stellar
dynamics in galaxy cores.

Early theoretical predictions pointed out that the bound debris returns to the
BH at a rate $\propto t^{-5/3}$ (Rees, 1988; Phinney, 1989). It was then
speculated that, having returned, matter quickly forms a disk whose outer
radius is comparable in size to the pericenter of the star’s original
trajectory. With such a small scale, the inflow time through the disk should
be so short that the light emitted by the disk would follow the matter
fallback, i.e., likewise $\propto t^{-5/3}$. With a characteristic temperature
$\sim 10^{5}-10^{6}$ K, the emitted light would be in the FUV/EUV or soft
X-ray band, and the peak luminosity would be much larger than Eddington.

However, the observed luminosity rarely reaches the Eddington luminosity for
the expected SMBH masses. In addition, when more TDEs were discovered in the
optical, it was realized that typical temperatures are a few $10^{4}$ K
(Gezari, 2021), implying that the radiating area is much larger than that of a
small disk whose radial scale is similar to the star’s pericenter. Moreover,
the total energy radiated during the flare period is generally two orders of
magnitude smaller than the energy expected if half the stellar mass had been
efficiently accreted onto a SMBH. In fact, it is an order of magnitude smaller
than the “circularization” energy that would have been emitted during the
initial formation of the small disk envisioned. That the radiated energy is so
small is sometimes called the “inverse energy crisis” (Piran et al., 2015;
Svirski et al., 2017).

A possible explanation for the low luminosity and low temperatures observed is
that the energy produced by the accretion disk is reprocessed by a radiation-
driven wind ejected from the disk itself (Strubbe & Quataert, 2009). If this
wind carries a significant kinetic energy, it would also resolve the “inverse
energy crisis” (Metzger & Stone, 2016) (see, however Matsumoto & Piran 2021).
An alternative possibility is that matter does not circularize quickly and
most of the observed emission arises from self-intersection shocks at the
apocenter (Shiokawa et al., 2015; Piran et al., 2015; Krolik et al., 2016).
These shocks, which are expected to be strong at and shortly after the time of
maximum mass return, take place at $\simeq O(10^{3})r_{\rm g}$ from the BH
($r_{\rm g}\equiv GM_{\bullet}/c^{2}$ for BH mass $M_{\bullet}$), are
consistent with the observed luminosity, temperature, line width and total
energy generated (Ryu et al., 2020a).

The long-term fate of bound debris is less clear. Because it is created on
highly eccentric orbits, with only a small further diminution in angular
momentum it may be able to fall directly into the black hole, releasing very
little energy (Svirski et al., 2017). Alternatively, when there has been time
for the magnetorotational instability (MRI) to build strong
magnetohydrodynamic (MHD) turbulence and for the gas to lose energy to
radiation, the debris may accrete gradually, while radiating efficiently
(Shiokawa et al., 2015).

Although it may be possible by observational means to determine whether one of
these two different scenarios occurs, numerical simulations of the disruption
and subsequent accretion process may provide an alternative way to resolve
this issue. However, such simulations are hindered by the extreme contrast in
length and timescales involved. Adding fully relativistic features that are
critical for some of the physics ingredients also poses technical challenges.

Aiming to clarify the many questions regarding the evolution of the bound
debris’ energy, angular momentum, and location, in this work we present a
fully relativistic numerical simulation of a complete tidal disruption of a
realistic $3M_{\odot}$ star by a $10^{5}M_{\odot}$ SMBH in which we follow the
system long enough to see the majority of the bound mass return to the black
hole. Several of these features have never previously appeared in a global TDE
simulation: fully relativistic treatment, both in hydrodynamics and stellar
self-gravity; a main-sequence internal structure for the star’s initial state;
and its long duration. Our simulation scheme does not include time-dependent
radiation transfer; instead, we make the approximation (for our parameters,
well-justified for most of the debris mass) that the radiation pressure is the
value achieved in LTE and does not diffuse relative to the gas.

The structure of the work is as follows: we define the physical problem and
present our numerical scheme in §2. We discuss the results in §3 and their
implications in §4. In this last section we also compare our results to
previous work. We conclude in §5, where we also discuss the possible
observational implications of this work.

##  2 The Calculation

###  2.1 Physics Problem

To provide a context for our choice of numerical methods and parameters, we
begin with a brief discussion of our problem’s overall structure. We begin
with a $3\;M_{\odot}$ star (radius $R_{\star}=2.4\;R_{\odot}$) whose internal
structure is taken from a MESA (Paxton et al., 2011) evolution to middle-age
on the main-sequence, the age at which the hydrogen mass fraction in its core
has fallen to 0.5. The star (initially placed at $r\simeq 900~{}r_{\rm g}$
from the black hole) approaches a SMBH of $M_{\bullet}=10^{5}\;M_{\odot}$ on a
parabolic orbit with a pericenter distance of $r_{\rm p}\simeq 110~{}r_{\rm
g}$, just close enough to be fully disrupted. Although the nominal tidal
radius $r_{\rm t}=R_{\star}(M_{\bullet}/M_{\star})^{1/3}$ is $\simeq
370~{}r_{\rm g}$ for this black hole mass and stellar mass, the critical
distance within which a star is fully disrupted is given by
$\Psi(M_{\bullet},M_{\star})r_{\rm t}$, where the order-unity factor $\Psi$
encodes general relativistic corrections through its $M_{\bullet}$ dependence
and stellar structure corrections through its $M_{\star}$ dependence Ryu et
al. (2020c).

As the star passes through pericenter, it begins to be torn apart; the process
is complete by the time its center-of-mass reaches $\sim 7500r_{\rm g}$.
During the debris’ first orbit, its trajectory is ballistic, with specific
orbital energy $E$ and specific angular momentum $J$ very close to that of the
star’s center-of-mass. Immediately after the disruption, the distribution of
mass with energy $dM/dE$ is roughly a square-wave and is approximately
symmetric with respect to $E=0$. Although order-of-magnitude arguments suggest
that the half-width of the energy square-wave $\Delta\varepsilon\simeq
GM_{\bullet}R_{\star}/r_{\rm t}^{2}$, Ryu et al. (2020c) showed that this
estimate should be multiplied by an order-unity correction factor
$\Xi(M_{\bullet},M_{\star})$, which depends on both $M_{\bullet}$ and
$M_{\star}$ as it accounts for both relativistic and stellar internal
structure effects. For our parameters, $\Xi\approx 1.64$, so that the energy
half-width is

$\displaystyle\Delta E\simeq 1.4\times 10^{-4}c^{2}$
$\displaystyle\left(\frac{\Xi}{1.64}\right)\left(\frac{M_{\bullet}}{10^{5}\;M_{\odot}}\right)^{1/3}$
$\displaystyle\times\left(\frac{M_{\star}}{3\;M_{\odot}}\right)^{2/3}\left(\frac{R_{\star}}{2.4R_{\odot}}\right)^{-1}\
.$ (1)

The semimajor axis $a$ and the eccentricity $e$ of the bound debris with
$E=-\Delta E$ are

$\displaystyle a$ $\displaystyle=\frac{GM_{\bullet}}{2E}\simeq 3600r_{\rm
g}\left(\frac{\Xi}{1.64}\right)^{-1}$
$\displaystyle\left(\frac{M_{\bullet}}{10^{5}\;M_{\odot}}\right)^{-1/3}\left(\frac{M_{\star}}{1\;M_{\odot}}\right)^{-2/3}\left(\frac{R_{\star}}{2.4\;R_{\odot}}\right),$
(2)

and

$\displaystyle e\simeq
1-0.07\left(\frac{\Xi}{1.64}\right)\left(\frac{M_{\bullet}}{10^{5}\;M_{\odot}}\right)^{-1/3}\left(\frac{M_{\star}}{3\;M_{\odot}}\right)^{1/3},$
(3)

respectively. The apocenter distance of the debris with $E=-\Delta E$ is then
$(1+e)a\simeq 7000r_{\rm g}$.

The bound debris, i.e., debris with $E<0$, must return to the vicinity of the
SMBH. The maximal mass return rate occurs when the debris with $E\simeq-\Delta
E$ returns. The fallback time scale of the material with $E\simeq-\Delta E$ is
very nearly its orbital period,

$\displaystyle t_{0}$
$\displaystyle=\frac{\uppi}{\sqrt{2}}\frac{GM_{\bullet}}{\Delta{E}^{3/2}}\simeq
7.6~{}{\rm days}~{}\left(\frac{\Xi}{1.64}\right)^{-3/2}$
$\displaystyle\left(\frac{M_{\bullet}}{10^{5}\;M_{\odot}}\right)^{1/2}\left(\frac{M_{\star}}{3\;M_{\odot}}\right)^{-1}\left(\frac{R_{\star}}{2.4\;R_{\odot}}\right)^{3/2}.$
(4)

The maximal fallback rate is then

$\displaystyle\dot{M}_{0}$ $\displaystyle\simeq\frac{M_{\star}}{3t_{\rm
0}}\simeq 0.13\;M_{\odot}\,{\rm days}^{-1}\left(\frac{\Xi}{1.64}\right)^{3/2}$
$\displaystyle\times\left(\frac{M_{\bullet}}{10^{5}\;M_{\odot}}\right)^{-1/2}\left(\frac{M_{\star}}{3\;M_{\odot}}\right)^{2}\left(\frac{R_{\star}}{2.4\;R_{\odot}}\right)^{-3/2}.$
(5)

The different scales involved are why the problem is a difficult one for
numerical simulation. Fluid travels through regions where the characteristic
scale on which gravity changes runs from $\sim r_{\rm g}$ to $\sim
10^{4}r_{\rm g}$, while the structure of the star varies on a scale that is a
fraction of $R_{\star}\sim 15r_{\rm g}$. Similarly, the characteristic
dynamical timescale ranges from $\sim r_{\rm g}/c$ to $\sim 10^{6}r_{\rm
g}/c$.

Figure 1: Successive moments in a full tidal disruption event simulated using
two domains: The density distribution of the star before (1) and after (2) the
pericenter passage and the debris after the pericenter passage (3, 4) in the
equatorial plane. The yellow dot in the big left panel indicates the BH. The
cyan box in the four right panels depicts the boundary of Domain1. Once the
star is disrupted, we increase the size of the box (from 2 to 3) to keep the
debris inside Domain1 where the self-gravity is calculated as much as
possible. Panels are not drawn to scale.

###  2.2 Code

To treat this problem, we perform a fully relativistic global hydrodynamics
simulation using a software package comprising three core codes, PatchworkMHD
(Avara et al., in prepartion), Harm3d (Noble et al., 2009), and a relativistic
self-gravity solver (Ryu et al., 2020d).

We use the intrinsically conservative general relativistic magneto-
hydrodynamics (GRMHD) code Harm3d (Noble et al., 2009) to solve the equations
of relativistic pure hydrodynamics. This code employs the Lax–Friedrichs
numerical flux formula and uses a parabolic interpolation method (Colella &
Woodward, 1984) with a monotonized central-differenced slope limiter. Because
of its robust algorithm, it has been used for studying a wide variety of SMBH
accretion problems (e.g., Noble et al., 2009, 2012; Shiokawa et al., 2015; Ryu
et al., 2020c, d). We take radiation pressure into account by setting the
pressure $p=\rho kT/\bar{m}+aT^{4}/3$ when the internal energy density
$u=(3/2)\rho kT/\bar{m}+aT^{4}$. Here $\rho$ is the proper rest-mass density,
$\bar{m}$ is the mass per particle, and $T$ is the temperature. We can then
define an equation of state with an “effective adiabatic index” (Shiokawa et
al., 2015) that varies between $\gamma=4/3$ and $5/3$ depending on the ratio
of the gas pressure and radiation pressure.

Our self-gravity solver is described in detail in Ryu et al. (2020d). In
brief, it constructs a metric valid in the star’s center-of-mass frame by
superposing the potential found by a Poisson solver operating in a tetrad
frame defined at the center-of-mass on top of the metric for the BH’s Kerr
spacetime described in the center-of-mass frame. Because this is a compact
free-fall frame, the spacetime for the star and its close surroundings is
always very nearly flat.

The large contrasts in length and time scales noted earlier make such a
computation prohibitively expensive if the entire range of scales is resolved
in a single domain. To overcome this difficulty, we introduce multiple domains
with PatchworkMHD(Avara et al., in preparation, first application in Avara et
al. 2023). PatchworkMHD is an extension of the original purely hydrodynamic
Patchwork code (Shiokawa et al., 2018) that both enables MHD and introduces a
number of algorithmic optimizations. Both versions create a multiframe /
multiphysics / multiscale infrastructure in which independent programs
simultaneously evolve individual patches with their own velocities, internal
coordinate systems, grids, and physics repertories. The evolution of each
patch is parallelized in terms of subdomains; where different patches overlap,
the infrastructure coordinates boundary data exchange between the relevant
processors in the different patches. A further extension of PatchworkMHD,
specialized to problems involving numerical relativity, is described in Bowen
et al. (2020).

###  2.3 Domain setup

During the first part of a TDE, a star travels through extremely rarefied gas
as it traverses a nearly parabolic trajectory around the black hole. This
situation is a prime example of the contrasts motivating our use of a
multipatch system: the characteristic lengthscales inside the star are much
smaller than those in the surrounding gas; in addition, self-gravity is
important inside and near the star, but not in the remainder of the volume.
Consequently, during this portion of the event we employ PatchworkMHD and
apply it to two patches: one covering the star, the other covering the
remainder of the SMBH’s neighborhood. Once the star is fully disrupted, there
is no further need of the star patch, and it is removed, its content
interpolated onto a single remaining patch, which evolves the entire region
around the SMBH.

####  2.3.1 The star’s pericenter passage—two-patch simulation

This first stage begins with the initial approach of the star to pericenter
and ends when the star’s center-of-mass trajectory reaches a distance from the
SMBH $\simeq 20r_{\rm t}$. It ends here because the star is then fully
disrupted. During this stage, a rectangular-solid Cartesian domain denoted
Domain1 covers the volume around the star; it is completely embedded in a
larger spherical-coordinate domain denoted Domain2 that ultimately covers a
large spherical region centered on the SMBH.

Following the methods described in Ryu et al. (2020d), in Domain1 we follow
the hydrodynamics of gas with self-gravity in a frame that follows the star’s
center of mass. Initially, the domain is a cube with edge-length $5R_{\star}$,
and the cell size in each dimension is $\simeq R_{\star}/25$. The orientation
of this box relative to the black hole is rotated during the disruption in
order to follow the direction of the bulk of the tidal flow, which is easily
predicted from Kepler’s laws. As the debris expands, the cube is adaptively
extended with constant cell-size to keep the debris inside the domain for as
long as possible.

As the debris expands, a fraction of the debris crosses smoothly from Domain1
to Domain2, where we continue to evolve it under the gravity of the SMBH
without any gas self-gravity. For Domain2, we adopt modified spherical
coordinates in Schwarzschild spacetime. That is, if the code’s three spatial
coordinates are $(x_{1},x_{2},x_{3})$, they can represent spherical
coordinates $(r,\theta,\phi)$ through the relations

$\displaystyle r$ $\displaystyle=e^{x_{1}},$ $\displaystyle\theta$
$\displaystyle=0.5\uppi[1.0+h_{1}x_{2}+(1-h_{1}-\frac{2.0\theta_{0}}{\uppi}){x_{2}}^{h_{2}}],$
$\displaystyle\phi$ $\displaystyle=x_{3}.$ (6)

where $h_{1}$ and $h_{2}$ are tuning parameters that determine the vertical
coordinate structure ($h_{1}\simeq 0.03$ and $h_{2}\simeq 9$), and
$\theta_{0}$ is the angle from the polar axis to the $\theta$-boundaries. In
this coordinate system, the radial grid cells have a constant ratio of cell-
dimension to radius, and the $\theta$ grid cells are concentrated near the
mid-plane. This coordinate system is suitable for modeling systems that
involve a wide range of radial scales and contain a disk-like structure near
the mid-plane.

To be computationally efficient, instead of fixing the size and the tuning
parameters throughout the simulation, we flexibly adjust the size and the
resolution of Domain2 so that it is large enough to contain the entire debris,
but we do not waste cells on regions where there is no debris. This strategy
reduces the computational cost by a large factor. To ensure proper resolution
during the period of flexible domain size, we require at least 15-20 cells per
scale height in all three dimensions. In addition, at all times during the
two-patch evolution, we keep the cell sizes in the overlapping regions of both
domains comparable.

At its largest extent, Domain2 runs from $r_{\rm min}=40r_{\rm g}$ to $r_{\rm
max}=18000r_{\rm g}$. The maximum radius is chosen to be greater than the
apocenter of debris that would return to the black hole within a time
$4t_{0}$; our simulation ran for $3t_{0}$. The minimum radius was chosen by
balancing two opposing goals: minimizing the mass lost through the inner
radial boundary while maximizing the time step so as to limit computational
expense. Similarly, when Domain2 has its largest volume, the polar angle
extends from $\theta_{0}=2^{\circ}$ to $\pi-\theta_{0}=178^{\circ}$. The
azimuthal angle $\phi$ covers a full $2\pi$ when Domain2 is largest.

We adopt outflow boundary conditions at all boundaries of Domain2. All the
primitive variables are extrapolated to the ghost cells with zero gradients.
However, to ensure outflow, if the extrapolated normal component of the fluid
velocity in the ghost cells is directed inward, it is set to zero. When this
domain is maximally extended azimuthally, we provide boundary conditions to
the processor domains having surfaces at $\phi=0$ and $\phi=2\pi$ by matching
those with the same radial and polar angle locations.

We show in Figure 1 the density distribution of the star before the pericenter
passage (top-left of the four small panels), during the passage (top-right)
and the debris after the passage (bottom) in the equatorial plane in the two-
domain simulation. The cyan box demarcates the boundary of Domain1. These
figures demonstrate how, as the star is disrupted, the stellar debris crosses
smoothly the inter-patch boundary between Domain1 and Domain2.

####  2.3.2 Single-domain simulation

When only a tiny fraction of the star’s original mass remains within Domain1,
self-gravity becomes irrelevant (measured in the star center-of-mass frame,
the acceleration throughout Domain 1 due to self-gravity is $\lesssim 10\%$ of
the tidal acceleration), and gradients in gas properties are no longer
connected with $R_{\star}$. We therefore remove Domain1 and continue the
simulation using only Domain2. In this stage, it takes its maximum extent.

Shocks caused by the stream-stream intersection dissipate the orbital energy
into thermal energy, which results in the vertical expansion of streams. The
vertical coordinate system used for the early evolution (Equation 2.3.1) is
not suitable to resolve gas at a large height as the resolution becomes
increasingly crude as $z$ increases. To better resolve the gas at high $z$, we
reduce the concentration of $\theta$ cells to the midplane while leaving the
$r$ and $\phi$ grids untouched. To do this, we redefine $\theta(x_{2})$:

$\displaystyle\theta$
$\displaystyle=\alpha(\tanh[b(x_{2}-a)]+\tanh[b(x_{2}+a)])+0.5\uppi,$ (7)

where $\alpha=-(0.5\uppi-\theta_{0})/[\tanh(b(-0.5-a))+\tanh(b(-0.5+a))]$.
Here, $a$ and $b$ are a set of tuning parameters that determine the vertical
structure. At this stage, we fix the domain extent of $r$ and $\phi$. But we
keep adjusting the domain extent and the resolution of the vertical structure
flexibly, whenever it is necessary, to ensure sufficiently high vertical
resolution (at least 20 cells per vertical scale heights) by properly choosing
the cell number (within $80-120$), $\theta_{0}$ ($2-15^{\circ}$) and the
tuning parameters ($a\simeq 0.32-0.34$ and $b\sim 9.8$).

Figure 2: The energy distribution of debris at $t\simeq 0.4t_{0}$. The
vertical dashed line ($E/\Delta\epsilon=\Xi\simeq-1.64$) shows the
characteristic energy at which fallback rate peaks (see Figure 3).  Figure 3:
The mass fallback rate $\dot{M}$, using the energy distribution shown in
Figure 2. The rate and time are normalized by $\dot{M}_{0}=M_{\star}/3t_{0}$
and $t_{0}$, respectively. The horizontal dashed line indicates the Eddington
limit with the radiative efficiencty of 0.01 and the diagonal dashed line a
power-law of $t^{-5/3}$.

Figure 4: The density distribution around the black hole excision (yellow dot)
in the equatorial plane at four different times, $t/t_{0}=0.5$, $1$, $2$ and
3. The extent of the inset is $500r_{\rm g}$.

##  3 Results

Figure 5: Accumulated mass normalized by $M_{\star}$ as a function of the
distance from the SMBH at $t/t_{0}\simeq 0.5$, 1, 2, and 3. The dashed black
vertical line indicates the pericenter distance of the original stellar orbit.
The mass expelled through the radial inner boundary at $r=40r_{\rm g}$ is
included as if the accreted mass is confined within $r=40r_{\rm g}$.

###  3.1 Overview

The star becomes strongly distorted as it passes through the pericenter (see
Figure 1) and then falls apart entirely as it travels farther away from the
black hole. When its entire mass has been dispersed, the orbital energy
distribution of the debris $dM/dE$ is not confined within the order of
magnitude estimate of the energy distribution’s width $\Delta\epsilon$; it is
about twice as wide and has extended wings (Figure 2). The mass of bound gas
($E<0$) is $\simeq 0.48M_{\star}$, slightly smaller than that of the unbound
gas ($\simeq 0.52M_{\star}$). The energy at which the peak mass return rate
(Figure 3) occurs is $E\simeq-1.64\Delta\epsilon$, or $\Xi=1.64$. If we assume
the debris follows ballistic orbits, the energy distribution presented in
Figure 2 may be translated into a mass fallback rate as a function of time
(Figure 3). The mass return rate peaks earlier by a factor of
$|\Xi|^{-1.5}\simeq 0.5$ than the traditional order of magnitude estimate
predicts, and the maximal fallback rate is greater by a factor $\simeq 2.5$.

After the bound debris goes out through its orbital apocenter and returns to
the region close to the BH, it undergoes multiple shocks (initially near
pericenter and apocenter), and finally forms an eccentric flow, as illustrated
in Figure 4 (see Section 3.3 for more details)111The density fluctuations
visible in newly-returning matter are likely numerical artifacts; they are
erased by the first shock the gas encounters and have no subsequent influence.
Near pericenter, strengthening vertical gravity and orbital convergence
compress the returning debris, creating a “nozzle” shock (as predicted by
Evans & Kochanek 1989) visible at $t\gtrsim 0.5t_{0}$. Over time, the shocked
gas becomes hotter and thicker. On the way to apocenter, the gas cools
adiabatically222We ignore here the effect of recombination as this energy is
negligible compared with the orbital energy: even at apocenter, the ratio of
gas kinetic energy to recombination energy is $\sim(GM_{*}/R_{*})/E_{\rm
recomb}\sim O(100)$.. Near apocenter, the previously-shocked outgoing debris
collides with fresh incoming debris, creating another shock (the “apocenter”
shock). Previously-shocked gas close to the orbital plane is deflected inward
toward the BH, while the portion farther from the plane is deflected above and
below the incoming stream (see Section 3.4). These deflections broaden the
angular momentum distribution. A small part of the debris loses enough angular
momentum that it acquires a pericenter smaller than the star’s pericenter.
Other gas gains angular momentum, which results in the nozzle shock-front
gradually extending to larger and larger radii.

At $t\simeq(1-2)t_{0}$, the debris in the apocenter region undergoes a
dramatic transition in shape, from well-defined incoming and outgoing streams
to an extended eccentric accretion flow. By the time the return rate of the
newly incoming debris declines, the mass that had arrived earlier becomes
large enough to significantly disrupt the newly incoming debris’ orbit. The
space inside the apocenter region is then quickly filled with gas. The outcome
is an extended eccentric accretion flow ($e\simeq 0.4-0.5$ at $t\simeq
3t_{0}$), most of whose mass resides at radii $\sim 10^{3}-10^{4}r_{\rm g}$
(see Figure 5). At $t\simeq 3t_{0}$, only $2-3\times 10^{-2}\;M_{\odot}$ can
be found inside $2r_{\rm p}$. By contrast, from $t\approx t_{0}$ onward,
nearly all the thermal energy is found at small radii, a condition that has
consequences for the time-dependence of escaping radiation. Throughout the
volume occupied by the debris, radiation pressure dominates: it is generally
larger than gas pressure by a factor $\sim 10^{3}-10^{4}$.

Figure 6: The azimuthally integrated density distribution at $t/t_{0}=0.5$, 1,
2, and 3.

In other words, circularization is not prompt: the flow retains significant
eccentricity, and the great majority of the gas remains at a distance
$10^{3}-10^{4}r_{\rm g}\gg r_{\rm p}$ even after several characteristic
timescales. That this is so can also be seen from another point of view. At
$t\simeq 3t_{0}$, the total amount of dissipated energy is only 10% of the
energy, $E_{\rm circ}\equiv GM_{\bullet}/4r_{\rm p}$, required for the debris
to fully “circularize” into a compact disk on the commonly-expected radial
scale of $2r_{\rm p}$. Extrapolating this slow energy dissipation rate to late
times suggests that true “circularization” would take a few tens of $t_{0}$
(see Figure 11). As we will discuss in detail in Section 4.1, the total
dissipation rate is roughly constant with time from $t\simeq 0.5t_{0}$ until
the end of our simulation at $t\simeq 3t_{0}$. Thus, there is no runaway
increase in dissipation of the sort suggested by Steinberg & Stone (2022).

The shocked gas expands outward quasi-symmetrically (Figure 6). Because the
intrinsic binding energy of the debris is much smaller than $E_{\rm circ}$,
the dissipated energy is large enough to be comparable to the specific orbital
energy. As a result, the expanding material is marginally bound. The radial
expansion speed of the gas near the photosphere ($r\simeq 7000-10000r_{\rm
g}$, Figure 7) at $t\simeq 3t_{0}$ is $0.005-0.01c\simeq 1500-3000$ km/s; the
associated specific energy is $\lesssim 10^{-4}c^{2}$, comparable to the
intrinsic energy scale, $\Delta E$.

Figure 7: The location of the thermalization photosphere (magenta curves) as
seen by a distant observer plotted over the density distribution at $\phi=0$
at $t/t_{0}=1$, $2$, $3$ (top panels) and at $\phi=90^{\circ}$, $180^{\circ}$
and $270^{\circ}$ at $t/t_{0}=3$ (bottom panels). We define the thermalization
optical depth as $\sqrt{\tau_{\rm t}\tau_{\rm ff}}$, for $\tau_{\rm T}$
($\tau_{\rm ff}$) the Thomson (absorption) optical depths integrated radially
inwards from the outer boundary.

Although we do not incorporate radiation transfer into the simulation, we
estimate the luminosity in post-processing (see § 3.6). The bolometric
luminosity rises to $(6-7)\times 10^{43}$erg/s in $t\simeq t_{0}$; its
effective temperature on the photosphere is generally close to $\sim 2\times
10^{5}$ K.

Figure 8: The density (left panels) and temperature (right panels)
distribution near the black hole (yellow dot) at four different times. Note
the stationary spiral shock structure that appears around $2t_{0}$.

Figure 9: The temperature distribution (left) at $t=t_{0}$ in the equatorial
plane. The right panel shows a zoom-in of the region demarcated by the green
box in the left panel near where two streams collide. The arrows indicate the
motion of gas. Notice the different color scales in the left and right panels.

Figure 10: The mass-weighted average of the specific orbital energy of debris
(top-left), angular momentum (top-right), eccentricity (bottom-left) and
aspect ratio (bottom-right) for given $r$ at $t/t_{0}=0.5$, 1, 2 and 3. The
dashed grey vertical line in all panels indicate the pericenter distance of
the original stellar orbit. The dotted grey horizontal line in some panels
show different quantities: $E=\Delta\epsilon$ (top-left), initial angular
momentum (top-right), and $R_{\star}/r_{\rm p}$ (bottom-right). The diagonal
black line in the top-left panel shows a power-law of $r^{-1}$.

###  3.2 Formation of shocks

It is convenient to divide the multiple shocks by approximate location:
pericenter or apocenter. The compression and heating of the gas near the
pericenter are depicted in detail in Figure 8. When the incoming stream is
narrow and well-defined (the two upper panels at $t\lesssim t_{0}$), the
nozzle shock structure can be described in terms of two components. As the
different portions of the returning stream converge, adiabatic compression
raises the temperature at the center of the stream. The shock itself runs more
or less radially across the stream, extending both inward and outward from the
stream center. However, at later times (beginning at $t\gtrsim t_{0}$), the
structure becomes more complex. The matter that has been deflected onto lower
angular momentum orbits circulates in the region inside $r_{\rm p}$ and
develops a pair of nearly stationary spiral shocks. The shock closest to the
position of the nozzle shock stretches progressively farther outward, reaching
$\sim 2000r_{\rm g}$ by $t\simeq 2t_{0}$ (see Figure 4). However, while the
shock extends to greater radii, it also loses strength. A similar progressive
widening and weakening of the nozzle shock were found by Shiokawa et al.
(2015).

Outgoing previously-returned matter intersects the path of newly-arriving
matter in the apocenter region because a combination of apsidal rotation due
to the finite duration of the disruption and relativistic apsidal precession
cause earlier and later stream orbits to be misaligned (Shiokawa et al.,
2015). When the apocenter shock first forms, it is relatively close to the
black hole ($\gtrsim 1000r_{\rm g}$) because the very first debris to return
has orbital energy more negative than $-\Delta E$. As the mass-return rate
rises, its orbital energy also increases, so the debris apocenter moves
outward. However, even at $t\simeq t_{0}$, when the shock is located at
$r\simeq 6000r_{\rm g}$, it is found closer to the BH than the apocenter
distance corresponding to $E=-\Delta E$ because the outgoing stream has lost
orbital energy to dissipation in the nozzle shock. At still later times, the
apocenter shock moves further inward as the mean energy of the previously
shocked matter decreases further.

Some of the outgoing material, upon collision with the incoming stream, is
deflected both horizontally and vertically. In the left panel of Figure 9, we
show the temperature distribution in the equatorial plane at $t\simeq t_{0}$.
In the right panel, one can see a clear boundary between the incoming stream
(with temperature $T\simeq 10^{5}$ K) and the outgoing stream (with $T\gtrsim
3\times 10^{5}$ K). At this boundary, the outgoing stream is deflected towards
the SMBH.

The two panels of Figure 9 also reveal that, just as found in Shiokawa et al.
(2015), the apocenter shock splits into two (called shocks 2a and 2b by
Shiokawa et al.). Shock 2a occurs where the outgoing stream encounters the
incoming stream, while shock 2b (visible in these figures as the surface on
which the outgoing gas temperature rises from $\simeq 3\times 10^{5}$ K to
$\simeq 1\times 10^{6}$ K) occurs where one portion of the outgoing gas
catches up with another portion that has been decelerated by gravity.

###  3.3 Debris orbit evolution

Undergoing multiple shocks, the debris becomes less eccentric and thicker both
in the orbital plane and vertically. To be more quantitative, in Figure 10, we
present the spherically mass-weighted averages of $E/E_{\rm circ}$, angular
momentum $l$, eccentricity $e$ and the aspect ratio of debris $h/r$ as
functions of $r$ at $t/t_{0}=0.5$, 1, 2 and 3. At $t=0.5t_{0}$, only the
debris with $|E|\gtrsim 3\Delta\epsilon$ has returned or is returning, which
means only a few $10^{-2}\;M_{\odot}$ is within $4000r_{\rm g}$, the apocenter
distance for $|E|=3\Delta\epsilon$. At this point, the majority of the mass
within $4000r_{\rm g}$ is unshocked incoming debris; only the earlier-
returned, outgoing debris has been shocked once near the pericenter.
Consequently, all the mass-weighted quantities primarily reflect the
properties of the newly incoming gas. The debris at $r_{\rm p}\lesssim
r<4000r_{\rm g}$ has $|E|\simeq 3\Delta\epsilon\simeq 0.1E_{\rm circ}$ and
angular momentum very close to that of the original stellar orbit, giving an
eccentricity almost unity: $1-e\lesssim 0.06$. Similarly, the aspect ratio of
the debris is $\simeq 0.2$, comparable to the value implied by
$R_{\star}/r_{\rm p}\simeq 0.1$. At $0.5\lesssim t/t_{0}\lesssim 1$, the
shocks begin to dissipate more energy. This results in a small fraction of the
gas dropping to orbital energy $\sim-E_{\rm circ}$ and orbiting within $r_{\rm
p}$. However, the global orbital properties of the debris ($r\gtrsim r_{\rm
p}$) remain largely unchanged. By the second half of the simulation, enough
orbital energy has been dissipated to reduce the orbital period of much of the
gas by a factor order unity (i.e., as previously noted, the dissipated energy
becomes comparable to the orbital energy). The result is a structure close to
virialization: the mean orbital energy as a function of radius is $\approx-
GM_{\rm BH}/2r$.

A similar evolution occurs in the radial distribution of angular momentum $l$.
Until $t\simeq t_{0}$, the specific angular momentum for nearly all the gas is
essentially the same as that of the star. At later times, however, the shocks
have greatly broadened the angular momentum distribution, and the gas has
sorted itself with higher angular momentum material located at larger radii.
It is, in fact, this broadening of the angular momentum distribution that
leads to the outward radial extension of the nozzle shock noted in the
previous subsection. The combined changes in $E$ and $l$ lead to a decrease in
eccentricity from $\gtrsim 0.8-0.9$ at $t\simeq t_{0}$ to $0.4-0.5$ at
$t\simeq 2t_{0}$. The aspect ratio changes rather little with time over the
entire inner region, $r\lesssim 2000r_{\rm g}$: it rises only from $\simeq
0.2$ to $\simeq 0.5$. At larger radii, it rises from $\sim 10^{-2}$
immediately after the disruption to $\simeq 0.5$; like the inner region, this
is accomplished by $t\simeq 2t_{0}$.

As mentioned in § 3.1, progress toward circularization can be thought of in
terms of the rate at which orbital energy is converted into thermal energy.
Evaluating this rate in units of $E_{\rm circ}$ and $t_{0}$ gives a measure of
the circularization “efficiency”:

$\displaystyle\eta\equiv\frac{|dU_{\rm gas+rad}/dt|}{M_{\rm gas}E_{\rm
circ}/t_{0}}.$ (8)

Here $U_{\rm gas+rad}$ is the total thermal (gas + radiation) energy,
including thermal energy that has been carried out of the simulation domain by
gas flows; from $t\simeq 0.5t_{0}$ onward, its rate of change is roughly
constant at $\simeq 1.4\times 10^{44}$ erg s-1, integrating to total thermal
energy $\simeq 2.2\times 10^{50}$ erg at the end of the simulation. $M_{\rm
gas}$ is the total gas mass in the domain at $r<10^{4}r_{\rm g}$.
Conveniently, the mass within this region stays very nearly constant over the
duration of the simulation (see the top panel of Figure 5), so the total rate
of thermal energy creation is very close to directly proportional to $\eta$.
This “circularization efficiency” as a function of time is shown in Figure
11333Note that this definition of $\eta$ differs from the one used by
Steinberg & Stone (2022) who compare the instantaneous energy dissipation rate
to the instantaneous fall back rate, $\eta^{\prime}=|dU_{\rm
gas+rad}/dt|(\dot{M}_{\rm fb}E_{\rm circ})^{-1}$. . $\eta$ rises very rapidly
during the time from $t\simeq 0.5t_{0}$ to $t\simeq t_{0}$, quickly reaching a
maximum $\sim 0.03$. However, from $t\simeq t_{0}$ to the end of the
simulation, its value remains very nearly flat. If it were to stay at that
level until circularization was complete, the process would take $\simeq
30t_{0}$.

Figure 11: Circularization efficiency $\eta$, defined in Equation 8, as a
function of time. The red line shows the time average of $\eta$, estimated as
$\eta(t)=\int_{0}^{t}dt^{\prime}\eta(t^{\prime})/\int_{0}^{t}dt^{\prime}$.

###  3.4 Radial motion

The outer bound of the region in which bound debris is found expands quasi-
spherically due to the combined effects of the radiation pressure gradient
built by shock heating and deflection caused by stream-stream collisions.
Figure 6 illustrates the azimuthally integrated density at four different
times, $t/t_{0}=0.5$, 1, 2 and 3. At $t=0.5t_{0}$, the outgoing gas that had
been heated by the nozzle shock forms a vertically thick density structure
within $3000r_{\rm g}$. Later on, at $t>0.5t_{0}$, the apocenter shock
contributes further to the expansion of the gas. At $t\gtrsim 2t_{0}$, most of
the mass that had been pushed outward starts to fall back towards the SMBH.

However, in the outermost $\lesssim 1\%$ of the flow, i.e., radii
$7000-10000r_{\rm g}$, the gas continues moving outward with a speed of
$0.005-0.01c\simeq 1500-3000$ km/s even at late times. To test whether this is
an incipient wind, we define unbound gas by the total energy criterion
$E+E_{\rm therm}>0$, where $E$ is the orbital (kinetic and gravitational)
energy and $E_{\rm therm}$ is the thermal energy. We find almost no matter
that has been made unbound after the initial disruption at any time during our
simulation (i.e., up to $3t_{0}$). Figure 12 depicts $E+E_{\rm therm}$ at
$\phi=0^{\circ}$ (near the nozzle shock) and $180^{\circ}$ (near the apocenter
shock) at $t\simeq 3t_{0}$. As already noted, although essentially all the
mass is bound, its specific binding energy is small. It is therefore not
straightforward to predict the final fate of the expanding envelope based on
the energy distribution measured at a specific time: energy is readily
transferred from one part of the system to another, or from one form of energy
to another. However, because the fallback rate declines beyond $3t_{0}$, the
major energy source at later times would be effectively the interactions of
gas in the accretion flow that has formed. Hence, the energy distribution in
the outer envelope is unlikely to evolve much over time. In addition, because
we allow for no radiative losses, the thermal energy content measured in the
simulation data, particularly in the outer layers, is an upper bound to the
actual value. This material’s most likely long-term evolution is, therefore, a
gradual deceleration followed by an eventual fallback.

Moving gas carries energy. Defining the mechanical luminosity by

$\displaystyle L_{\rm mech}(r)=\int_{\Omega}d\Omega r^{2}\rho v^{r}(E+E_{\rm
therm}),$ (9)

we find (as shown in Fig. 13) that the net $L_{\rm mech}$ integrated over
spherical shells is nearly always positive for $t\gtrsim t_{0}$ and is super-
Eddington. The predominantly negative slope in $L_{\rm mech}(r)$ at $r\lesssim
10^{3}r_{\rm g}$ and $r\gtrsim 5000r_{\rm g}$ indicates that these regions are
gaining energy, while the relatively constant mechanical luminosity at
$10^{3}r_{\rm g}\lesssim r\lesssim 5\times 10^{3}r_{\rm g}$ shows relatively
insignificant energy exchange in that range of radii.

In interpreting this radial flow of mechanical luminosity, it is important to
note that it is due to a mix of outwardly-moving unbound matter and inwardly-
moving bound matter. Both signs of radial velocity are represented on almost
every spherical shell; in fact, the mass-weighted mean radial velocity is
generally inward with magnitude $\sim 300-1000$ km s-1. Thus, the regions
gaining energy do so in large part, but not exclusively, by losing strongly
bound mass.

Figure 12: The distribution of the total energy $E+E_{\rm therm}$, normalized
by the local gravitational potential $E_{\rm grav}=-GM_{\bullet}/r$, at
$\phi=0^{\circ}$ (left, near nozzle shock) and $180^{\circ}$ (right, near
apocenter shock) at $t=3t_{0}$.  Figure 13: Radial total energy
($E+E_{\rm{therm}}$) flux integrated over a spherical shell at given $r$,
normalized by the Eddington luminosity $L_{\rm Edd}$, at four different times.
Positive values mean net positive (or negative) energy carried by gas moving
radially outward (or inward).  Figure 14: The rate at which gas is expelled
through the inner radial boundary at $r=40r_{\rm g}$, divided by
$\dot{M}_{0}$, as a function of $t/t_{0}$.

###  3.5 Matter loss through inner radial boundary

Figure 14 shows the rate at which mass falls through the inner radial boundary
at $40r_{\rm g}$ and leaves the computational domain. This rate rises rapidly
from $t=0$ to $t=t_{0}$, and from then until $t=3t_{0}$ fluctuates about a
nearly-constant mean value $\simeq 0.035\dot{M}_{0}$. The total lost mass up
to $t\simeq 3t_{0}$ is $\lesssim 0.025M_{\star}$. This is a factor $\sim 3$
smaller than the rate, at $t=3t_{0}$, of mass-loss through a similarly-placed
inner cut-out in the simulation of Shiokawa et al. (2015).

Although the simulation provides no information about what happens to this gas
after it passes within $r=40r_{\rm g}$, we can make certain informed
speculations. By computing its mass-weighted mean specific energy and angular
momentum, we find that its mean eccentricity is $\simeq 0.7-0.8$ and does not
evolve over time. Its mean pericenter and apocenter distances are $10-15r_{\rm
g}$ and $80-120r_{\rm g}$, respectively. These values imply that the mass lost
through the boundary would not accrete onto the SMBH immediately. In fact, if
it follows such an orbit, it should re-emerge from the inner cutout. If it did
so without suffering any dissipation, it would erase the positive energy flux
we find at the inner boundary, which is due to bound matter leaving the
computational domain. On the other hand, if it suffered the maximum amount of
dissipation consistent with an unchanging angular momentum and settled onto
circular orbits inside $40r_{g}$, it might release as much as $\simeq 10^{51}$
erg, more than enough to unbind the rest of the bound debris, whose binding
energy is only $\sim 3\times 10^{50}$ erg. However, this estimate of
dissipation is an upper bound, and likely a very loose one, with the actual
dissipation far smaller. The location of the shocks suffered by this gas would
be much nearer its apocenter at $\sim 100r_{g}$ than its pericenter, reducing
the kinetic energy available for dissipation by an order of magnitude or more.
Oblique shock geometry, as seen in the directly-simulated shocks, sharply
diminishes how much kinetic energy is dissipated per shock passage. In
addition, if matter does fall onto weakly eccentric orbits with semimajor axes
$<40r_{\rm g}$, it will block further inflow, thereby decreasing the net flow
across the $r=40r_{\rm g}$ surface. Whatever accumulation of gas occurs within
$\lesssim 100r_{g}$ is also unlikely to have much effect on the bulk of the
bound debris because, as visible in Figure 8, the radial range at which the
bulk of the returning gas passes through the nozzle shock moves steadily
outward over time, reaching $\gtrsim 500r_{g}$ by $t=3t_{0}$.

Figure 15: (Left) The photosphere temperature distribution $dL/dT$ at four
different times $t/t_{0}=0.5$, 1, 2, and 3, and (right) the angular
distribution of the temperature at $t=3t_{0}$.

###  3.6 Radiation

To infer the bolometric luminosity of this event, we post-process the
simulation results. We first identify the thermalization photosphere with the
surface where $\sqrt{\tau_{\rm T}\tau_{\rm ff}}\simeq 1$. Here, $\tau_{\rm T}$
($\tau_{\rm ff}$) is the Thomson (absorption) optical depth integrated
radially inward from the outer $r$ boundary. The absorption cross section is
calculated using an OPAL table for solar metallicity. The upper panels of
Figure 7 show that the photosphere expands quasi-spherically, which is
expected from the radial expansion of the outer debris. At $t=t_{0}$, the
photosphere is quasi-spherical located at $r\simeq 4000-5000r_{\rm g}$. It
expands to $9000-10000r_{\rm g}$ at $t=2t_{0}$ and to $\simeq 12000r_{\rm g}$
at $t=3t_{0}$. Given the eccentric orbit of debris, the radius of the
photosphere depends on $\phi$ at a quantitative level, even while the its
overall shape is qualitatively round. To demonstrate the $\phi$-dependence, we
show in the bottom panels of Figure 7 the density distribution and the
photosphere at four different azimuthal angles, $\phi=0$, $90^{\circ}$,
$180^{\circ}$ and $270^{\circ}$ at $t=3t_{0}$.

We then estimate the cooling time $t_{\rm cool}$ at all locations inside the
photosphere as

$\displaystyle t_{\rm cool}(r)=\frac{h_{\rho}\tau(r)}{c}{(1+u_{\rm gas}/u_{\rm
rad})},$ (10)

where $h_{\rho}$ is the first-moment density scale height of the gas along a
radial path, $\tau$ is the optical depth (radially integrated) to $r$, and
$u_{\rm gas}/u_{\rm rad}$ is the ratio of the local internal energy density to
the radiation energy (a ratio that is often only slightly greater than unity).
We then estimate the luminosity by integrating the energy escape rate over the
volume within the photosphere, but including only those locations for which
$t_{\rm cool}$ is smaller than the elapsed time in the simulation. This
condition accounts for the fact that in order to leave the debris by time $t$,
the cooling time from the light’s point of origin must be less than $t$. The
resulting expression is

$\displaystyle L=\int_{0}^{2\pi}\int_{\theta_{\rm c}}^{\pi-\theta_{\rm
c}}\int_{r=r(t_{\rm cool}<t)}^{r=r(\tau=1)}\frac{aT^{4}}{t_{\rm
cool}}r^{2}\sin\theta drd\theta d\phi,$ (11)

where $a$ is the radiation constant.

The local effective temperature at each individual cell near the photosphere
is then calculated as $T_{\rm ph}=[(dL/dA)/\sigma]^{1/4}$, where $A$ is the
surface area of the photosphere and $\sigma$ is the Stefan–Boltzmann constant.
We estimate that the peak luminosity is $\simeq 10^{44}$ erg/s $\simeq
10L_{\rm Edd}$, which occurs at $t\simeq t_{0}$. This is roughly the mean rate
of thermal energy creation during the simulation. The photospheric temperature
distribution at $t=0.5t_{0}$ can be described as nearly flat within the range
$\times 10^{4}K\lesssim T\lesssim 2\times 10^{5}K$, as shown in the left panel
of Figure 15. At $t\gtrsim t_{0}$, the distribution becomes narrower: the
distribution at $t\gtrsim 2t_{0}$ has a single peak at $T\simeq(5-6)\times
10^{4}$ K. In the right panel of Figure 15 we show the photospheric
temperature distribution as a function of observer direction at $t\simeq
3t_{0}$. The temperature is $5-6\times 10^{4}$ K over almost the entire
photosphere except for a noticeably low-temperature spot at $\phi\simeq\pi$
and $\theta=0.5\pi$, corresponding to the low-$T$ incoming stream.

##  4 Discussion

###  4.1 Circularization - fast or slow?

The pace of “circularization” has long played a central role in understanding
how TDE flares are powered. If it is rapid, i.e., takes place over a time
$\lesssim t_{0}$, the debris joins a small ($r\lesssim r_{\rm p}$) accretion
disk as soon as it first returns. In addition, accretion takes place on a
timescale short compared to $t_{0}$ because the orbital period on this scale
is shorter than $t_{0}$ by a factor $(M_{\star}/M_{\bullet})^{1/2}\sim
10^{-3}$. Even after waiting $\sim 10$ orbital periods for MRI turbulence to
saturate and then consuming many more orbital periods to flow inward by
magnetic stresses, the total inflow is still short compared to $t_{0}$. The
dissipation rate at the time of peak mass-return would then be strongly super-
Eddington.

The result of our simulation, however, is that “circularization” is actually
very slow. We find that the returning debris forms a large cloud that
stretches all the way from the pericenter of the original stellar orbit to the
apocenter of the most-bound debris, a dynamic range $\sim 30-100$. Throughout
the first $3t_{0}$ after disruption, only a small fraction of the debris
resides within the pericenter. The mass-weighted mean eccentricity falls from
$\simeq 0.9$ to $\simeq 0.4-0.6$ by $t\gtrsim 2t_{0}$, but doesn’t decrease
further from that time to until at least $3t_{0}$. Thus, by this late time the
debris has neither achieved a circular orbit nor been compressed within $r\sim
r_{\rm p}$.

Such slow circularization is consistent with the low energy dissipation rate.
The thermal energy within the system is very small compared with
$E_{\mathrm{circ}}$, the energy that must be removed from the bound debris’
orbital energy in order to fully circularize it. Similarly, the
circularization efficiency parameter $\eta$ suggests that several dozen
$t_{0}$ are required in order to dissipate $E_{\mathrm{circ}}$ of energy (see
Figure 11). Thus, we may conclude that little circularization is accomplished
during the time in which most of the debris mass returns to the BH.

Our conclusions in this matter agree with earlier findings of Shiokawa et al.
(2015), who used a somewhat cruder computational scheme and less realistic
conditions (these authors considered a disruption of a white dwarf of
$0.64\;M_{\odot}$ by a BH of $500\;M_{\odot}$). On the other hand, they differ
with those of Steinberg & Stone (2022), who analyzed the “circularization
efficiency” in terms of the heating rate per returning mass rather than our
definition, the heating rate per returned mass over a time $t_{0}$. On the
basis of tracking this definition of circularization efficiency up to
$t=t_{0}$, they argued that it was growing exponentially on a timescale $\sim
t_{0}$, so that full circularization might be achieved quickly. Interestingly,
our definition of efficiency also grows rapidly with time during the first
$t_{0}$; in this respect we agree with Steinberg & Stone (2022). However, we
also find that it flattens out shortly after $t_{0}$. Thus, one possible
explanation of the contrast in our conclusions about the magnitude of energy
dissipation is simply that our simulation ran longer than theirs when measured
in $t_{0}$ units. It is also possible that some of the difference in the
results could be attributed to differences in our physical assumptions.
Steinberg & Stone (2022) used a spherical harmonic oscillator potential at
$r<30~{}r_{\rm g}$ (private conversation with Elad Steinberg) and a Paczynski-
Wiita potential at larger radii, whereas we used a Schwarzschild spacetime
with a cut-out at $40r_{\rm g}$; they described radiation transport by a flux-
limited diffusion scheme, whereas we included radiation only as a contribution
(often the dominant one) to the pressure. On balance, though, because the
gravity descriptions used are not very different on the relevant lengthscales
and the long cooling times in the system severely limit radiative diffusion,
these contrasts are unlikely to explain this disagreement. Lastly, it is
possible that the difference in parameters (our $M_{\rm BH}=10^{5}M_{\odot}$
and $M_{*}=3M_{\odot}$ vs. their $M_{\rm BH}=10^{6}M_{\odot}$ and
$M_{*}=1M_{\odot}$) may also play a role. Further simulations will be
necessary in order to test this possibility.

###  4.2 Energy Dissipation: Shocks vs. Accretion

The physical assumptions in our simulation restrict the creation of thermal
energy to two mechanisms: shocks and compressive work done within the fluid.
There is no energy release due to classical accretion because our equations
contain neither MHD turbulence nor phenomenological viscosity. Nonetheless, we
have demonstrated that shocks and compression can, without these other
processes, generate enough energy during a few $t_{0}$ to power the observed
luminosity of TDEs. We estimate a photon luminosity during this period of
$\sim 10^{44}$ erg s-1, and all of this energy was generated by shocks and
compressive work. As discussed in the previous subsection, we have
demonstrated the absence of orbital energy loss that is a prerequisite for
forming a classical accretion flow.

###  4.3 Outflow

A third interesting finding is that we do not find a significant unbound
outflow emerging from the bound debris. Very nearly all the bound material
that has returned to the vicinity of the SMBH remains bound by the end of our
simulations. Although we do see outward motion, its slow speed indicates that
the material remains bound (see Figure 12). It should therefore eventually
slow down and fall back.

This result places an even stronger upper bound on the dissipated energy than
the earlier result that there was too little dissipation to circularize the
matter, as the specific energy needed to unbind the debris is significantly
smaller than that needed to circularize it around the original pericenter.
Whereas the circularization energy is $\sim 3\times 10^{51}$ erg, the binding
energy is only $\sim 3\times 10^{50}$ erg. That almost no initially bound
debris is rendered unbound is consistent with observational limits on outflows
from both radio and optical TDEs (Matsumoto & Piran, 2021).

This conclusion, which is contrary to a number of predictions (e.g., Jiang et
al., 2016; Bonnerot et al., 2021; Huang et al., 2023), also casts some doubt
on the possibility (Metzger & Stone, 2016) that the kinetic energy of an
outflow is the solution to the “inverse energy crisis” mentioned earlier. When
the source of heating is shocks, we find negligible transport of energy to
infinity associated with outflows. Interestingly, although Steinberg & Stone
(2022) do find an unbound outflow, its mechanical luminosity is only $\simeq
1.6\times 10^{42}$ erg s-1 if the outflow velocity they quote, 7500 km s-1, is
its velocity at infinity. This is such a small fraction of the heating rate
that even this sort of wind does not play a significant role in the energy
budget. Moreover, even if all the mass lost through our inner boundary were
quickly accreted in a radiatively efficient manner, as we have already
estimated, the associated heat produced would be only a factor of 4 – 5
greater than the thermal energy generated by shocks in the first $3t_{0}$
after the disruption. In this sense, we have also placed a strong limit on the
ability of a wind dependent upon accretion energy to carry away a large
quantity of energy.

###  4.4 The ultimate fate of the bound debris

Our simulation ends at $3t_{0}$ with nearly all the bound debris $\sim
10^{3}-10^{4}r_{\rm g}$ from the BH, spread over a large eccentric cloud. The
question naturally arises: what happens next? Extrapolating from their
qualitatively similar results, Shiokawa et al. (2015) suggested that, after
the usual $\sim 10$ orbital-period time necessary for saturation of MHD
turbulence driven by the magnetorotational instability (MRI), the gas would
accrete in more or less the fashion of circular accretion disks.

Since that work, it has been shown (Chan et al., 2018, 2022) that, indeed, the
MRI is a genuine exponentially-growing instability in eccentric disks and, in
its nonlinear development, creates internal magnetic stresses comparable to
those seen in circular disks. However, its outward transport of angular
momentum may, in the context of eccentric disks, cause the innermost matter to
grow in eccentricity while outer matter, the recipient of the angular momentum
removed from the inner matter, becomes more circular (Chan et al., 2022).

If this is the generic result of MRI-driven turbulent stresses in an eccentric
disk, accretion might be radiatively inefficient, as matter can plunge
directly into the BH if it has sufficiently small angular momentum (Svirski et
al., 2017). The condition for this to happen is for the angular momentum
transport to be accompanied by very little orbital energy loss. It is then
possible for fluid elements of very low angular momentum to fall ballistically
into the SMBH after having radiated only a small amount of energy. In this
case, the system will dim rapidly after the thermal energy created by shocks
has diffused out in radiation.

However, it remains to be determined whether this is, in fact, the situation
in TDE eccentric accretion flows. If, instead, the work done by torques
associated with angular momentum transport is substantial, a compact, more
nearly circular, accretion disk eventually forms. This disk will then behave
much more like a conventional accretion flow, radiating soft X-rays until most
of the disk mass has been consumed.

If the energy lost per unit accreted mass comes anywhere near the $\sim
0.1c^{2}$ of radiatively-efficient accretion, the total energy radiated over
this prolonged accretion phase could be quite large: $0.1M_{\star}c^{2}\approx
10^{53}$ erg. However, sufficiently long accretion timescales might keep the
luminosity relatively low. There is some observational evidence for such
radiation on multi-year timescales, both in X-rays (e.g., Jonker et al., 2020;
Kajava et al., 2020) and UV (e.g., van Velzen et al., 2021; Hammerstein et
al., 2023). In these long-term observations, the luminosity declines gradually
enough ($\propto t^{-1}$) to make the total energy radiated logarithmically
divergent.

A related question is posed by the matter that passed through our inner radial
boundary. To the extent that some portion of it does dissipate enough energy
to achieve a near-ISCO orbit, there is the possibility of significant energy
release in excess of what was seen in our simulation. In fact, in order to
generate soft X-ray luminosities comparable to those often seen ($\sim
10^{44}$ erg s-1 at peak), all that is required is a mass accretion rate $\sim
3\times 10^{-3}M_{\star}/t_{0}$. Thus, if $\sim 0.3$ of the matter passing
through our inner boundary were able to accrete onto the BH, it might be able
to account for the X-ray luminosity sometimes seen, given an optically thin
path to infinity. For the parameters of our simulation, there appears to be
little or no solid angle through which such a path exists (see Figure 7), but,
as shown by Ryu et al. (2020b), the ratio $t_{\rm cool}/t_{0}$ falls to
$\lesssim O(1)$ when $M_{\bullet}\gtrsim 10^{6}M_{\odot}$. Consequently,
radiative cooling might make the flow geometrically thinner for larger
$M_{\bullet}$ events, permitting X-rays emitted near the center to emerge
during the time of the optical/UV flare. Alternatively, for those cases that,
like our simulation, have relatively long cooling times, X-ray emission may
become visible only after a significant delay relative to the optical/UV
light, a delay that has been observed in several TDEs (Gezari et al., 2017;
Kajava et al., 2020; Hinkle et al., 2021; Goodwin et al., 2022).

###  4.5 Comparison with Ryu et al. (2020b)

Ryu et al. (2020b) introduced a parameter-inference method TDEmass for
$M_{\bullet}$ and $M_{\star}$ built on the assumption that optical TDEs are
powered by the energy dissipated by the apocenter shock. In this method, one
assumes that the peak luminosity and temperature occur at $t\simeq 1.5t_{0}$
when the most-bound debris collide with the incoming stream at the apocenter.
Using our numerical results to determine the two parameters of TDEmass
(setting $c_{1}$, the ratio of the photospheric radius to the apocenter
distance, to 1.2 and the solid angle of the photosphere to $4\pi$), we find
that the luminosity and temperature at the peak of the bolometric lightcurve
would be $3\times 10^{44}$ erg/s and 70000 K (see Equations. 1, 2, 6 and 9 of
Ryu et al., 2020b).

These values can be compared with the estimates derived from our cooling time
method, $L\approx 10^{44}$ erg s-1 and $T\approx 60000$ K, measured at
$t\simeq 1.5t_{0}$. The contrast in luminosity may be a consequence of an
assumption made in the method of Ryu et al. (2020b): that the heating due to
shocks is radiated promptly. Although this is a reasonable approximation for
$M_{\bullet}\gtrsim 10^{6}M_{\odot}$, our simulation has shown that when
$M_{\bullet}$ is as small as $\sim 10^{5}M_{\odot}$, cooling is significantly
retarded (in fact, Ryu et al. 2020b pointed out that $t_{\rm
cool}/t_{0}\propto\Xi^{5/2}M_{\bullet}^{-7/6}M_{\star}^{4/9}$). Although our
simulation suggests that this method may require some refinement in the range
of small SMBH masses, overall, whether with or without the corrections
suggested by the detailed numerical simulation, the peak luminosity is in the
range of optical/UV bright TDEs. The temperature estimated from the simulation
is larger by only a factor of 1.2, which is reasonable given the approximate
treatment of the radiation in our scheme.

##  5 Conclusions

Following the energy often provides a well-marked path toward understanding
the major elements of a physical event. It is especially useful for TDEs
because one might define their central question as “How does matter whose
initial specific orbital energy is $\sim 10^{-4}c^{2}$ dissipate enough energy
to both power the observed radiation and then, in the long-run, fall into the
black hole?”

This question can be made more specific by pointing out certain milestones in
energy. In a typical TDE flare, $\sim 3\times 10^{50}$ erg is radiated during
its brightest period, although in a number of cases an order of magnitude more
is radiated over multiple year timescales (e.g., van Velzen et al., 2021;
Hammerstein et al., 2023). The immediately post-disruption binding energy of
the bound gas in the simulation described here is very similar to this number,
$2\times 10^{50}$ erg. The energy required to circularize all the bound gas is
1.5 dex larger, $7.5\times 10^{51}$ erg. Lastly, the energy that might be
liberated through conventional relativistic accretion of all the bound
material is $\sim 3\times 10^{53}$ erg.

Comparing the results of our simulation—$\sim 1.5\times 10^{50}$ erg radiated
over a time $3t_{0}$ long and final gas binding energy less than a factor of 2
greater than in the initial state ($3\times 10^{50}$ erg)—to these milestones
points to a number of strong implications.

First, and most importantly, the radiation we estimate as arising from our
simulation is very close to the typical radiated energy during the brightest
portion of the flare. In other words, the hydrodynamics we have computed, in
which shocks dissipate orbital energy into heat, succeed in matching the most
important quantity describing TDE flares.

Second, over this period the binding energy of the debris does not change
appreciably. It immediately follows from the virial theorem that the scale of
the region occupied by the debris likewise does not change appreciably. The
only modification that might be made to this conclusion is that radiation
losses would increase the binding energy by a factor of $\sim 1.5-2$. The area
of the photosphere is determined by the scale of the region containing the
bound debris.

Third, swift “circularization”, that is, confinement of the bound debris to a
circular disk with outer radius $\sim r_{\rm p}$, does not happen. This
process requires the bulk of the debris to increase its binding energy by a
factor $\sim 30$; this did not happen.

Fourth, radiatively efficient accretion of most of the debris mass onto the
black hole certainly did not happen. If this had occurred, the mass remaining
on the grid would be substantially smaller, and the energy released would have
rendered the remaining mass strongly unbound, as it corresponds to a total
dissipated energy $\sim 10^{3}\times$ larger than seen.

Lastly, we have also found that, contrary to some expectations, essentially no
debris gas that was bound immediately after the disruption was rendered
unbound by shock dynamics.

## Acknowledgements

We are thankful to the anonymous referee for constructive comments and
suggestions. We thank Elad Steinberg and Nick Stone for helpful conversations.
We also thank Suvi Gezari for informing us about delayed X-ray flares in TDEs.
This research project was conducted using computational resources (and/or
scientific computing services) at both the Texas Advanced Computing Center and
the Max-Planck Computing & Data Facility. At TACC, we used Frontera under
allocations PHY-20010 and AST-20021. In Germany, the simulations were
performed on the national supercomputer Hawk at the High Performance Computing
Center Stuttgart (HLRS) under the grant number 44232. TP is supported by ERC
grant MultiJets. JK is partially supported by NSF grants AST-2009260 and
PHY-2110339.

## References

  * Avara et al. (2023) Avara, M. J., Krolik, J. H., Campanelli, M., et al. 2023, arXiv e-prints, arXiv:2305.18538, doi: 10.48550/arXiv.2305.18538
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002, doi: 10.1088/1538-3873/aaecbe
  * Blanton et al. (2017) Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28, doi: 10.3847/1538-3881/aa7567
  * Bonnerot et al. (2021) Bonnerot, C., Lu, W., & Hopkins, P. F. 2021, M.N.R.A.S., 504, 4885, doi: 10.1093/mnras/stab398
  * Bowen et al. (2020) Bowen, D. B., Avara, M., Mewes, V., et al. 2020, arXiv e-prints, arXiv:2002.00088.  https://arxiv.org/abs/2002.00088
  * Chan et al. (2018) Chan, C.-H., Krolik, J. H., & Piran, T. 2018, ApJ, 856, 12, doi: 10.3847/1538-4357/aab15c
  * Chan et al. (2022) Chan, C.-H., Piran, T., & Krolik, J. H. 2022, ApJ, 933, 81, doi: 10.3847/1538-4357/ac68f3
  * Colella & Woodward (1984) Colella, P., & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174, doi: 10.1016/0021-9991(84)90143-8
  * Evans & Kochanek (1989) Evans, C. R., & Kochanek, C. S. 1989, ApJL, 346, L13, doi: 10.1086/185567
  * Gezari (2021) Gezari, S. 2021, arXiv e-prints, arXiv:2104.14580.  https://arxiv.org/abs/2104.14580
  * Gezari et al. (2017) Gezari, S., Cenko, S. B., & Arcavi, I. 2017, ApJL, 851, L47, doi: 10.3847/2041-8213/aaa0c2
  * Goodwin et al. (2022) Goodwin, A. J., van Velzen, S., Miller-Jones, J. C. A., et al. 2022, M.N.R.A.S., 511, 5328, doi: 10.1093/mnras/stac333
  * Hammerstein et al. (2023) Hammerstein, E., van Velzen, S., Gezari, S., et al. 2023, ApJ, 942, 9, doi: 10.3847/1538-4357/aca283
  * Hills (1988) Hills, J. G. 1988, Nat., 331, 687, doi: 10.1038/331687a0
  * Hinkle et al. (2021) Hinkle, J. T., Holoien, T. W. S., Auchettl, K., et al. 2021, M.N.R.A.S., 500, 1673, doi: 10.1093/mnras/staa3170
  * Huang et al. (2023) Huang, X., Davis, S. W., & Jiang, Y.-f. 2023, arXiv e-prints, arXiv:2303.17443, doi: 10.48550/arXiv.2303.17443
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jiang et al. (2016) Jiang, Y.-F., Guillochon, J., & Loeb, A. 2016, ApJ, 830, 125, doi: 10.3847/0004-637X/830/2/125
  * Jonker et al. (2020) Jonker, P. G., Stone, N. C., Generozov, A., van Velzen, S., & Metzger, B. 2020, ApJ, 889, 166, doi: 10.3847/1538-4357/ab659c
  * Kaiser et al. (2002) Kaiser, N., Aussel, H., Burke, B. E., et al. 2002, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4836, Survey and Other Telescope Technologies and Discoveries, ed. J. A. Tyson & S. Wolff, 154–164, doi: 10.1117/12.457365
  * Kajava et al. (2020) Kajava, J. J. E., Giustini, M., Saxton, R. D., & Miniutti, G. 2020, A&A, 639, A100, doi: 10.1051/0004-6361/202038165
  * Krolik et al. (2016) Krolik, J., Piran, T., Svirski, G., & Cheng, R. M. 2016, ApJ, 827, 127, doi: 10.3847/0004-637X/827/2/127
  * Matsumoto & Piran (2021) Matsumoto, T., & Piran, T. 2021, M.N.R.A.S., 507, 4196, doi: 10.1093/mnras/stab2418
  * Metzger & Stone (2016) Metzger, B. D., & Stone, N. C. 2016, M.N.R.A.S., 461, 948, doi: 10.1093/mnras/stw1394
  * Noble et al. (2009) Noble, S. C., Krolik, J. H., & Hawley, J. F. 2009, ApJ, 692, 411, doi: 10.1088/0004-637X/692/1/411
  * Noble et al. (2012) Noble, S. C., Mundim, B. C., Nakano, H., et al. 2012, ApJ, 755, 51, doi: 10.1088/0004-637X/755/1/51
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJ Supp., 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Phinney (1989) Phinney, E. S. 1989, in The Center of the Galaxy, ed. M. Morris, Vol. 136, 543
  * Piran et al. (2015) Piran, T., Svirski, G., Krolik, J., Cheng, R. M., & Shiokawa, H. 2015, ApJ, 806, 164, doi: 10.1088/0004-637X/806/2/164
  * Rees (1988) Rees, M. J. 1988, Nat., 333, 523, doi: 10.1038/333523a0
  * Ryu et al. (2020a) Ryu, T., Krolik, J., & Piran, T. 2020a, ApJ, 904, 73, doi: 10.3847/1538-4357/abbf4d
  * Ryu et al. (2020b) —. 2020b, ApJ, 904, 73, doi: 10.3847/1538-4357/abbf4d
  * Ryu et al. (2020c) Ryu, T., Krolik, J., Piran, T., & Noble, S. C. 2020c, ApJ, 904, 98, doi: 10.3847/1538-4357/abb3cf
  * Ryu et al. (2020d) —. 2020d, ApJ, 904, 99, doi: 10.3847/1538-4357/abb3cd
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48, doi: 10.1088/0004-637X/788/1/48
  * Shiokawa et al. (2018) Shiokawa, H., Cheng, R. M., Noble, S. C., & Krolik, J. H. 2018, ApJ, 861, 15, doi: 10.3847/1538-4357/aac2dd
  * Shiokawa et al. (2015) Shiokawa, H., Krolik, J. H., Cheng, R. M., Piran, T., & Noble, S. C. 2015, ApJ, 804, 85, doi: 10.1088/0004-637X/804/2/85
  * Steinberg & Stone (2022) Steinberg, E., & Stone, N. C. 2022, arXiv e-prints, arXiv:2206.10641, doi: 10.48550/arXiv.2206.10641
  * Strubbe & Quataert (2009) Strubbe, L. E., & Quataert, E. 2009, M.N.R.A.S., 400, 2070, doi: 10.1111/j.1365-2966.2009.15599.x
  * Svirski et al. (2017) Svirski, G., Piran, T., & Krolik, J. 2017, M.N.R.A.S., 467, 1426, doi: 10.1093/mnras/stx117
  * Van Rossum & Drake (2009) Van Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: CreateSpace)
  * van Velzen et al. (2021) van Velzen, S., Gezari, S., Hammerstein, E., et al. 2021, ApJ, 908, 4, doi: 10.3847/1538-4357/abc258