# Coupled Disk-Star Evolution in Galactic Nuclei and the Lifetimes of QPE
Sources

Itai Linial Department of Physics and Columbia Astrophysics Laboratory,
Columbia University, New York, NY 10027, USA  Institute for Advanced Study, 1
Einstein Drive, Princeton, NJ 08540, USA  Brian D. Metzger Department of
Physics and Columbia Astrophysics Laboratory, Columbia University, New York,
NY 10027, USA  Center for Computational Astrophysics, Flatiron Institute, 162
5th Ave, New York, NY 10010, USA

###### Abstract

A modest fraction of the stars in galactic nuclei fed towards the central
supermassive black hole (SMBH) approach on low-eccentricity orbits driven by
gravitational-wave radiation (extreme mass ratio inspiral, EMRI). In the
likely event that a gaseous accretion disk is created in the nucleus during
this slow inspiral (e.g., via an independent tidal-disruption event; TDE),
star-disk collisions generate regular short-lived flares consistent with the
observed quasi-periodic eruption (QPE) sources. We present a model for the
coupled star-disk evolution which self-consistently accounts for mass and
thermal energy injected into the disk from stellar collisions and associated
mass ablation. For weak collision/ablation heating, the disk is thermally-
unstable and undergoes limit-cycle oscillations which modulate its properties
and lead to accretion-powered outbursts on timescales of years to decades,
with a time-averaged accretion rate $\sim 0.1\dot{M}_{\rm Edd}$. Stronger
collision/ablation heating acts to stabilize the disk, enabling roughly steady
accretion at the EMRI-stripping rate. In either case, the stellar destruction
time through ablation, and hence the maximum QPE lifetime, is $\sim
10^{2}-10^{3}$ yr, far longer than fall-back accretion after a TDE. The
quiescent accretion disks in QPE sources may at the present epoch be self-
sustaining and fed primarily by EMRI ablation. Indeed, the observed range of
secular variability broadly match those predicted for collision-fed disks.
Changes in the QPE recurrence pattern following such outbursts, similar to
that observed in GSN 069, could arise from temporary misalignment between the
EMRI-fed disk and the SMBH equatorial plane as the former regrows its mass
after a state transition.

Tidal disruption (1696), X-ray transient sources (1852), Supermassive black
holes (1663) , Gravitational waves (678)

##  1 Introduction

Supermassive black holes (SMBH; Salpeter 1964) with millions to billions of
times the mass of our Sun, are ubiquitous in the centers of massive galaxies
(e.g., Kormendy & Richstone 1995). The majority of SMBH growth over cosmic
time results from radiatively-efficient accretion (Soltan, 1982), fueled by
gas funneled towards the nucleus via large-scale galactic processes (e.g.,
Hopkins et al. 2016). Nevertheless, SMBHs can grow via other mechanisms than
those responsible for large-scale active-galactic nuclei (AGN). Most massive
black holes are embedded in dense nuclear star clusters (e.g., Bahcall & Wolf
1976; Gallego-Cano et al. 2018) of comparable or greater mass, the
constituents of which provide additional, albeit “sequestered”, sources of
mass. Stellar capture or tidal disruption events (TDE; Hills 1975; Rees 1988)
offer another SMBH fuel source, one which may even dominate the growth of low-
mass black holes (e.g., Milosavljević et al. 2006; Stone et al. 2017). Most
TDEs involve the destruction of main-sequence stars and are accompanied by
months-long UV/X-ray emission (e.g., Gezari et al. 2006; though most of the
total radiated energy is emitted over years to decades, e.g., van Velzen et
al. 2019; Jonker et al. 2020). The tidal disruption of giant stars can
generate flares lasting up to several decades (e.g., MacLeod et al. 2012), and
may preferentially occur through multiple partial disruptions (e.g., MacLeod
et al. 2013), potentially setting a minimum floor on the accretion rates in
galactic nuclei (e.g., MacLeod et al. 2013; however, see Bonnerot et al.
2016).

While the vast majority of TDE victim stars approach the SMBH on nearly
parabolic orbits from parsec scale distances, a small fraction $\sim 1-10\%$
are predicted to reach the tidal radius due to gravitational-wave radiation on
less eccentric or even quasi-circular orbits. Such stellar “extreme mass-ratio
inspirals” (EMRIs; e.g., Linial & Sari 2017, 2023) have been hypothesized to
produce electromagnetic phenomena as rich or richer than ordinary TDEs, if and
when they begin Roche-lobe overflow onto the SMBH (e.g., Dai & Blandford 2013;
Linial & Sari 2017; Metzger & Stone 2017; King 2020, 2022; Metzger et al.
2022; Krolik & Linial 2022; Linial & Sari 2023; Lu & Quataert 2023). However,
their survival down to such small radial scales $\simeq r_{\rm T}$ is not
assured. In particular, the timescale for a stellar EMRI to inspiral through
gravitational-waves across the range of radii just outside the tidal radius
$r_{\rm T}$ (say, from $4r_{\rm T}$ to $2r_{\rm T}$) is typically millions of
year, far in excess of the typical interval between (ordinary) TDEs in
galactic nuclei, of about once per $10^{4}-10^{5}$ yr (e.g., Stone & Metzger
2016). This suggests that the EMRI will inevitably interact with the transient
gaseous disk formed from a TDE (Linial & Metzger 2023; hereafter LM23). The
generally inclined orbit of the EMRI results in the star periodically passing
through the disk midplane, where it collides with dense gas and has material
stripped from its surface (e.g., Armitage et al. 1996). These interactions can
potentially lead to changes in the star’s orbit (e.g., MacLeod & Lin 2020;
Generozov & Perets 2023; Linial & Quataert 2024) and/or eventually result in
the star’s destruction (LM23).

Repeated star-disk collisions may also be observable directly. “Quasi-periodic
Eruptions” (QPEs; e.g., Miniutti et al. 2019; Giustini et al. 2020; Arcodia et
al. 2021; Chakraborty et al. 2021; Arcodia et al. 2022; Miniutti et al. 2023b;
Webbe & Young 2023; Arcodia et al. 2024a) are a new class of periodic flaring
X-ray sources from low-mass galactic nuclei (e.g., Wevers et al. 2022), with
recurrence times ranging from a few to tens of hours. The short ($\lesssim$
hours-long) eruptions from X-ray QPE systems are characterized by peak
luminosities $\sim 10^{42}$ erg s-1, typically at least comparable to those of
the spectrally softer quiescent luminosity observed between eruptions
(Miniutti et al., 2019; Giustini et al., 2020; Arcodia et al., 2021;
Chakraborty et al., 2021; Arcodia et al., 2022; Miniutti et al., 2023b; Webbe
& Young, 2023). While many models have been proposed for QPEs (Zalamea et al.,
2010; King, 2020; Ingram et al., 2021; Xian et al., 2021; Metzger et al.,
2022; Pan et al., 2022; Zhao et al., 2022; King, 2022; Krolik & Linial, 2022;
Linial & Sari, 2017; Kaur et al., 2023; Linial & Metzger, 2023; Franchini et
al., 2023; Tagawa & Haiman, 2023), only a few account for an important clue: a
regular alternating behavior, observed in at least two of the QPE
sources$-$GSN 069 (Miniutti et al., 2019, 2023b) and eRO-QPE2 (Arcodia et al.,
2021), in which the temporal spacing between consecutive eruptions varies back
and forth by around 10%.

LM23 showed that many of the observed properties of X-ray QPEs can be
reproduced in a scenario in which a stellar EMRI collides twice per orbit with
a gaseous accretion disk (see also Suková et al. 2021; Xian et al. 2021;
Tagawa & Haiman 2023; Franchini et al. 2023; Zhao et al. 2022). The eruptions
are powered by hot shocked disk material which expands from either side of
midplane, akin to dual miniature supernova explosions (e.g., Ivanov et al.
1998). The oscillating long-short recurrence time behavior follows naturally
from the longer time the star spends between collisions on the apocenter side
of the disk. The timing of the eruptions is related to the schedule of disk
passages, which is subject to the precession of the stellar orbit and the
disk, as well as delayed caused by finite light-travel time (e.g., Xian et
al., 2021; Franchini et al., 2023; Chakraborty et al., 2024; Pasham et al.,
2024; Zhou et al., 2024).

One of the outstanding questions in the disk-star-collision scenario for QPEs
is the origin of the accretion flow responsible for the “quiescent” soft
X-rays seen between eruptions. In systems where it is sufficiently luminous to
detect, the quiescent spectra is consistent with thermal emission from the
inner regions (few $R_{\rm g}$ scale) of a radiatively efficient accretion
flow onto the SMBH at a sizable fraction $\gtrsim 0.1$ of the Eddington rate
(Miniutti et al., 2019; Arcodia et al., 2021, 2022; Miniutti et al., 2023b;
Arcodia et al., 2024a).

Long-term evolution of the quiescent luminosity (e.g., over timescales much
longer than the QPE period) offers clues to the origin of the gaseous disks in
these systems. While the quiescent X-ray flux of the QPE source RXJ
1301.9+2747 (Giustini et al., 2020) has been remarkably constant over 20 years
of monitoring, the quiescent luminosities of other QPE sources appear more
variable. Eight years prior to the discovery of QPE activity, the well-studied
source GSN 069 exhibited starting in 2010 an X-ray outburst similar to a TDE
(Shu et al., 2018; Sheng et al., 2021), though exhibiting a longer duration
and slower post-maximum decay than most TDE flares. However, starting at the
end of 2019, and after QPE flares were detected in the system, it exhibited a
second X-ray outburst qualitatively similar to the 2010 brightening, which was
interpreted by Miniutti et al. (2023a) as being another stellar-disruption
event in the same galactic nucleus, or the partial TDE of the same progenitor
that gave rise to the 2010 outburst. Quintin et al. (2023) describe the rise
of a QPE-like eruption (very similar to those seen from eRO-QPE1) in XMM-
Newton, during a gradual decay in X-ray quiescence, several months following
its detection as an optical TDE. Chakraborty et al. (2021) discuss a QPE
candidate XMMSL1 J024916.6-041244 that demonstrated a gradual decay in X-ray
flux reminiscent of GSN 069, with QPE variability seen in 2006, but none
observed in 2021. Of the two new QPE sources recently detected in eROSITA by
Arcodia et al. (2024a), the quiescent flux of eRO-QPE3 declined by an order of
magnitude over a time span of just a few years, while in eRO-QPE4 it showed
evidence of brightening (insofar that the source was not detected in archival
ROSAT or during earlier eROSITA survey epochs taken just a few years prior to
the first detection).

Tentative evidence also exists for some connection between the quiescent
accretion disk and QPE properties. After a year of QPE activity, the eruptions
from GSN 069 become undetected at the onset of the second quiescent-emission
outburst in 2019; however, as the outburst faded, the periodic eruptions
recently reappeared (Miniutti et al., 2023a), albeit initially with different
recurrence properties than before (G. Minutti, private communication). A
secular decrease in the average eruption luminosity from eRO-QPE3 coincided
with a gradual dimming of its quiescent luminosity (Arcodia et al., 2024a).
Chakraborty et al. (2024) and Pasham et al. (2024) reported that the eruptions
from eRO-QPE1 dimmed by a factor of $\sim 10$ over a three year baseline,
possibly with no significant variations in quiescent emission over this
duration.

While most stellar EMRIs will experience a TDE during their gravitational-wave
inspiral (LM23), the complex and diverse evolution of QPE quiescent
emission$-$including both decaying and rising light curves as well as re-
brightenings$-$may call into question whether it is fall-back accretion from
the TDE that continues to feed the SMBH and set the quiescent emission at the
present epoch. Indeed, a second gas source available to feed the disk is mass
stripped from the same stellar EMRI responsible for powering the eruptions.
The mass flux generated from star-disk collisions, for disk and orbital
properties necessary to explain the observed eruptions, was estimated to be
comparable to the black hole accretion rate implied by the quiescent emission
(LM23). While a transient TDE-disk still provides an appealing mechanism to
initiate the overall process, it is natural to ask whether a disk-star
collision system can becoming self-perpetuating even after TDE fall-back has
abated. This longer-lived system would terminate only after the star is
destroyed or if powerful collisions cease, e.g., as a result of changes in the
inclination between the disk and stellar orbit. Lu & Quataert (2023) similarly
proposed gradual stellar mass-loss as the source of quiescent accretion in QPE
systems, albeit through a different mechanism of Roche-lobe overflow regulated
by the effects of ram pressure on the stellar atmosphere from disk collisions;
however, such collisions will strip mass from the star, even for stellar
orbits outside the tidal radius.

In this paper, we develop a model for the coupled long-term evolution of
mutually-inclined star-disk systems, which accounts for mass and thermal
energy added to the disk via collisions and stellar stripping. In Sec. 2 we
describe the model, in which the star is assumed to reside on a quasi-circular
orbit and the disk with which it interacts is modeled as a single-annulus
centered around the collision radius. In Sec. 3 we present our numerical
results, most of which can be understood analytically. Though simplified in
many respects, the model reveals a surprisingly rich behavior, ranging from
limit-cycles driven by the thermal instabilities of radiation-dominated disks,
to approximately steady accretion when heating from heating from stellar-
collisions/ablation is sufficiently powerful to stabilize the disk. In Sec. 4
we discuss implications of our results for the emission properties of observed
QPE sources and their demographics in galactic nuclei. In Sec. 5 we summarize
our findings and conclude.

##  2 Model

A stellar EMRI of mass $M_{\star}=\mathcal{M}_{\star}M_{\odot}$ and radius
$R_{\star}=\mathcal{R}_{\star}R_{\odot}$ orbits the SMBH of mass $M_{\bullet}$
on a quasi-circular orbit of semi-major axis $r_{0}\sim 100R_{g}$ and period
$P_{\rm orb}\simeq 2\pi(r_{0}^{3}/GM_{\bullet})^{1/2}$, where $R_{\rm g}\equiv
GM_{\bullet}/c^{2}$ is the gravitational radius. Collisions between the star
and gaseous disk happen twice per orbit, such that the average interval
between observed eruptions (the QPE period) is given by $P_{\rm QPE}=P_{\rm
orb}/2$ (neglecting orbital eccentricity111Throughout most of the paper we
consider small eccentricities, as is expected for stellar-EMRIs, and as is
interpreted from the sample of known QPEs (e.g., Metzger et al. 2022; Linial &
Sari 2023; Lu & Quataert 2023; Linial & Metzger 2023; however, see King 2020).
Other repeating nuclear transients, such as ASASSN-14ko (Payne et al., 2021),
may also involve star-disk interaction, albeit on highly eccentric orbits
(e.g., Linial & Quataert, 2024).). The collision radius can thus be written:

$r_{0}\approx 1.4\times 10^{13}\,\,{\rm cm}\;\mathcal{P}_{\rm
QPE,4}^{2/3}M_{\bullet,6}^{1/3}\approx 95\,R_{\rm
g}\,\left(\frac{\mathcal{P_{\rm QPE,4}}}{M_{\bullet,6}}\right)^{2/3}\,,$ (1)

where $\mathcal{P}_{\rm QPE,4}\equiv P_{\rm QPE}/(4\,{\rm hr})$ and
$M_{\bullet,6}\equiv M_{\bullet}/(10^{6}M_{\odot})$. The condition that the
star does not overflow its Roche-lobe, $r_{0}\gtrsim 2r_{\rm T}\simeq
2R_{\star}(M_{\bullet}/M_{\star})^{1/3}$, defines a minimum QPE period,

$P_{\rm
QPE,min}\simeq\pi\left(\frac{8R_{\star}^{3}}{GM_{\star}}\right)^{1/2}\simeq
3.9\,{\rm hr}\,\mathcal{R}_{\star}^{3/2}\mathcal{M}_{\star}^{-1/2}.$ (2)

We focus on tight orbits with $P_{\rm orb}\gtrsim P_{\rm QPE,min}$, as these
stars will interact with the gaseous disks created by tidal disruption events
(LM23), the latter expected to form near the circularization radius $\simeq
2r_{\rm T}$ (for properties of the disrupted star broadly similar to the
EMRI).

The star’s orbital angular momentum vector makes an angle $\chi\in[0,\pi/2]$
with the spin-axis of the black hole (a retrograde orbit is not expected after
some period of time, for reasons discussed below), about which it undergoes
nodal precession from the Lense-Thirring torque on a timescale (e.g., LM23):

$\tau_{\rm prec}\simeq\frac{P_{\rm
orb}}{2a_{\bullet}}\left(\frac{r_{0}}{R_{\rm g}}\right)^{3/2}\simeq 1.4\,{\rm
yr}\,\frac{\mathcal{P}_{\rm
QPE,4}^{2}}{M_{\bullet,6}}\left(\frac{a_{\bullet}}{0.3}\right)^{-1},$ (3)

where $a_{\bullet}\leq 1$ is the dimensionless spin parameter of the black
hole. We neglect changes to the inclination angle or semi-major axis of the
star’s orbit, as we will show that for tight orbits, consistent with those
needed to explain observed QPE eruptions, the star is destroyed by collisions
faster than it can undergo appreciable radial migration or have its orbit
ground into the disk midplane (Sec. 4.1).

The angular momentum axis of the gaseous disk is assumed to align with the
black hole spin axis. Other than greatly simplifying the physical picture,
this assumed alignment is justified by two arguments: (1) the Bardeen &
Petterson (1975) effect, through which dissipation of precession-driven warps
in the disk eventually leads to spin-alignment out to hundreds of
gravitational radii, depending on the strength of stresses within the disk
(e.g., Nelson & Papaloizou 2000; Sorathia et al. 2013; Zanazzi & Lai 2019);
(2) any component of angular momentum perpendicular to the black hole spin
axis which results from mass added to the disk by stellar stripping will
nearly cancel out over many stellar nodal precession cycles: we shall show
that the precession time (Eq. (3)) is generally shorter than the timescale
over which the disk’s surface density evolves (see Fig. 1). We expand on this
justification later (Sec. 4.1), and possible observational consequences when
this assumption is violated on the QPE recurrence pattern.

We model the portion of the gaseous disk with which the star interacts as an
annulus of radius $r_{0}$ and radial extent $\sim r_{0}$, with a local surface
density $\Sigma$, local mass $M_{\rm d}=\pi r_{0}^{2}\Sigma$, midplane density
$\rho=\Sigma/h$, midplane temperature $T$, vertical scale-height $h\ll r_{0}$
and orbital speed

$v_{\rm orb}\simeq v_{\rm
k}=\left(\frac{GM_{\bullet}}{r_{0}}\right)^{1/2}\approx 0.1\,{\rm
c}\;\mathcal{P}_{\rm QPE,4}^{-1/3}M_{\bullet,6}^{1/3}.$ (4)

The star loses mass $\Delta M_{\star}$ per collision with the disk,
corresponding to a time-averaged mass-loss rate

$\dot{M}_{\star}=-\frac{\Delta M_{\star}}{P_{\rm QPE}}.$ (5)

Motivated by hydrodynamical simulations of supersonic flows past stars
(Armitage et al., 1996; Liu et al., 2015), we assume that the mass stripped
from the star per midplane passage is proportional to the ram pressure $p_{\rm
ram}=\rho v_{\rm rel}^{2}/2$, viz.

$\Delta M_{\star}=\eta\frac{p_{\rm ram}}{p_{\star}},$ (6)

where $p_{\star}=GM_{\star}^{2}/(4\pi R_{\star}^{4})$ is the interior pressure
of the star and $\eta\sim 10^{-3}$ (Armitage et al. 1996; Liu et al. 2015, Yao
& Quataert, in prep.). The collision speed, corresponding to the relative
motion of the star’s orbit and the rotating disk, can be written $v_{\rm
rel}=\xi_{\rm v}v_{\rm k}$, where $\xi_{\rm v}\in[0,2]$ is a dimensionless
number that depends on the star/disk inclination angle $\chi$. For a quasi-
circular orbit of eccentricity $e\ll 1$

$\xi_{\rm v}\simeq
2\sin{(\chi/2)}\left(1+\frac{3}{4}e\cos{f}\right)+\mathcal{O}(e^{2})\,,$ (7)

where $f$ is collision’s site true anomaly, measured with respect to the
orbit’s argument of periapsis.

Mass-loss from the star occurs sufficiently rapidly relative to its thermal
timescale (typically millions of years), that the star will evolve
adiabatically (Linial & Sari, 2017). For a low-mass star of initial
radius/mass $R_{\star,0}/M_{\star,0}$ whose interior is dominated by gas
pressure, with an effective adiabatic index $\gamma_{\rm ad}=5/3$, the radius
increases upon mass-loss according to:

$R_{\star}=R_{\star,0}\left(\frac{M_{\star}}{M_{\star,0}}\right)^{-1/3}.$ (8)

This expansion of the star will accelerate the mass-loss rate near the end of
its life and potentially induce Roche-lobe overflow onto the black hole,
somewhat hastening its final destruction.

While mass is added to the disk by star-collisions, mass is lost due to
accretion onto the black hole. This occurs at the rate

$\dot{M}_{\rm acc}=\frac{M_{\rm d}}{t_{\rm\nu}}=3\pi\nu_{0}\Sigma,$ (9)

where

$\tau_{\nu}=\frac{1}{3}\frac{r_{0}^{2}}{\nu_{0}}=\frac{1}{3}\frac{r_{0}}{\alpha
v_{\rm k}}\left(\frac{h}{r}\right)^{-2}$ (10)

is the local viscous timescale at radius $r_{0}$ and we have modeled the
effective viscosity using the standard Shakura & Sunyaev (1973) prescription
$\nu_{0}=\alpha r_{0}v_{\rm k}(h/r_{0})^{2},$ where $\alpha=0.01-0.1$ (e.g.,
King et al. 2007; Davis et al. 2010). Even if the gaseous disk is initially
formed at radii $\sim r_{0}$, it will spread to radii larger than $r_{0}$ due
to outward angular momentum redistribution; it might therefore seem
questionable that one can approximate the evolution of the local surface
density at $r_{0}$ using a simple one-zone model. However, by solving the
standard time-dependent diffusion equation for $\Sigma(r,t)$ (e.g., Pringle,
1981), with a mass-source term localized around a fixed radius ($r_{0}$, in
our case due to star-disk collisions) one can show that as long as the source
term evolves slowly compared to the local viscous timescale $\tau_{\nu}$, most
of the injected mass at $r_{0}$ indeed flows to small radii $<r_{0}$, rather
than accompanying the angular momentum carried to large radii (see Metzger et
al. 2012, their Appendix B).

In summary, the disk mass near the collision radius evolves according to:

$\frac{dM_{\rm d}}{dt}=|\dot{M}_{\star}|-\dot{M}_{\rm acc}+\dot{M}_{\rm TDE}.$
(11)

The final term accounts for an external mass source, which is needed to start
the collision-fed disk evolution. The primary example here is fall-back from
the eccentric debris streams of a tidally-disrupted star, which for a full-
disruption obeys (e.g., Phinney 1989)

$\dot{M}_{\rm TDE}=\frac{M_{\rm\star,TDE}}{3t_{\rm TDE}}\left(\frac{t}{t_{\rm
TDE}}\right)^{-5/3},$ (12)

where $t>t_{\rm TDE}$ is measured since the time of disruption,
$M_{\rm\star,TDE}$ is the mass of the disrupted star (distinct from the still-
present quasi-circular EMRI) and $t_{\rm TDE}$ is the fall-back time of the
most tightly bound debris (typically weeks to months for main-sequence stars;
e.g., Stone et al. 2013).

Each midplane passage, the star shocks a quantity of disk material
corresponding to that intercepted by its physical cross section:

$\Delta M_{\rm coll}\simeq 2\pi R_{\star}^{2}\Sigma.$ (13)

where the factor of 2 accounts for the increase in swept-up material mass for
$\chi\approx\pi/4$ (LM23). Although this hot fast expanding material may be
key to generating QPE eruptions as it rises above the disk midplane and
radiates away its energy (LM23; Sec. 4.3), the velocity dispersion it acquires
from the collision $\lesssim v_{\rm k}$ is not sufficient to travel
appreciably away from the collision radius $\sim r_{0}$. We therefore neglect
any such mass-loss term in Eq. (11) insofar that we are modeling an annulus of
radial width $\sim r_{0}$ centered at $r_{0}.$ Even were this not the case, in
practice we typically find that the mass shocked and ejected from the disk by
each collision is much smaller than that stripped from the star, i.e. $\Delta
M_{\rm coll}\ll\Delta M_{\star}$.

The midplane pressure of the disk includes contributions from gas and
radiation,

$p=p_{\rm gas}+p_{\rm rad}=\frac{\rho k_{\rm B}T}{\mu
m_{p}}+\frac{a}{3}T^{4},$ (14)

where $\mu$ is the mean molecular weight (we take $\mu=0.62$ for fully-ionized
solar-composition gas). Vertical hydrostatic balance gives an expression,

$\left(\frac{h}{r}\right)^{2}-\frac{1}{3}\frac{r_{0}}{\Sigma v_{\rm
k}^{2}}aT^{4}\left(\frac{h}{r}\right)-\frac{k_{\rm B}T}{\mu m_{p}v_{\rm
k}^{2}}=0,$ (15)

for the disk aspect ratio in terms of the midplane temperature $T$, yielding

$\frac{h}{r}=\frac{aT^{4}r_{0}}{6\Sigma v_{\rm
k}^{2}}\left(1+\sqrt{1+\frac{k_{\rm B}T}{\mu m_{p}}\left(\frac{6\Sigma v_{\rm
k}}{aT^{4}r_{0}}\right)^{2}}\right)\,.$ (16)

A final equation follows the evolution of the disk’s thermal energy $E_{\rm
th}=\left[(3/2)p_{\rm gas}+3p_{\rm rad}\right]2\pi hr_{0}^{2}$,

$\frac{dE_{\rm th}}{dt}=\dot{E}_{+}-\dot{E}_{-},$ (17)

as a result of various sources of heating $\dot{E}_{+}$ and cooling
$\dot{E}_{-}$. Heating results from both viscous dissipation (“accretion”) and
energy deposited by the stellar collisions and ablation,

$\dot{E}_{+}=\dot{E}_{\rm acc}+\dot{E}_{\rm abl}=\left(\frac{3}{8}\dot{M}_{\rm
acc}+\beta\left(|\dot{M}_{\star}|+\dot{M}_{\rm coll}\right)\right)v_{\rm
k}^{2},$ (18)

where $\dot{M}_{\rm coll}=\Delta M_{\rm coll}/P_{\rm QPE}$ is the time-
averaged rate of mass excavation from the disk due to collisions. In the
collision/ablation heating term, we assume that a fixed fraction of the
specific kinetic energy $\sim v_{\rm k}^{2}$ of the matter stripped
($|\dot{M}_{\star}|$) or intercepted ($\dot{M}_{\rm coll}$) by the star’s
passage through the disk goes into heating the midplane, where the efficiency
of this process is encapsulated in the dimensionless parameter $\beta\lesssim
1$.

For short orbital periods characteristic of most QPE sources we shall find
that $\dot{M}_{\rm coll}\ll|\dot{M}_{\star}|$, i.e. the rate that mass is
stripped from the star far exceeds the rate at which the star intercepts disk
material (this condition is roughly equivalent to the star’s destruction
occurring prior to its orbit undergoing significant drag-induced orbital
evolution). While the detailed evolution of the disk will turn out to be
sensitive to $\beta$, its value is uncertain because it depends on where and
how efficiently the kinetic energy of the ablated stellar mass is incorporated
into the disk midplane prior to being radiated. We note that there is no
a-priori reason that ablation and collision would heat the disk with the same
coupling coefficient $\beta$. However, since typically $\dot{M}_{\rm
coll}\ll|\dot{M}_{\star}|$ in the regimes of interest, this simplification is
not expected to have any significant impact on our results.

Cooling occurs via vertical radiative diffusion and radial advection, viz.

$\dot{E}_{-}=\dot{E}_{\rm rad}+\dot{E}_{\rm adv}=\frac{4}{3}\pi
r_{0}^{2}\frac{\sigma T^{4}}{\tau}+|\dot{M}_{\rm acc}|\frac{P}{\rho}\zeta,$
(19)

where $\tau=\kappa\Sigma$ is the optical depth through the disk midplane, $P$
is the disk’s midplane pressure and the dimensionless parameter $\zeta\propto
d{\rm ln}s/d{\rm ln}r$ depends on the radial entropy gradient and is typically
of order unity (e.g., Di Matteo et al. 2002). We take $\zeta=4$, though our
results are not sensitive to this choice because the disk typically spends
little if any time in the geometrically-thick state where advective cooling is
important. For the opacity law,

$\kappa=\kappa_{\rm T}+\kappa_{\rm K},$ (20)

we include contributions from electron scattering $\kappa_{\rm T}\approx 0.34$
cm2 g-1 and Kramers’ opacity $\kappa_{\rm K}=A_{\rm\kappa}\rho T^{-7/2}$, as
appropriate for temperatures $T\gtrsim 10^{4}$ K, where $A_{\kappa}\approx
5\times 10^{24}\;\rm cm^{5}\,gr^{-2}\,K^{7/2}$.

In thermal equilibrium ($\dot{E}_{+}=\dot{E}_{-}$), Eq. (17) becomes

$\beta v_{\rm k}^{2}\left(|\dot{M}_{\star}|+\dot{M}_{\rm
coll}\right)+\frac{9\pi}{8}\alpha\Sigma r_{0}v_{\rm
k}^{3}\left(\frac{h}{r}\right)^{2}=\frac{4}{3}\pi r_{0}^{2}\frac{\sigma
T^{4}}{\Sigma\kappa}+3\pi\zeta\alpha\Sigma r_{0}v_{\rm
k}^{3}\left(\frac{h}{r}\right)^{4}.$ (21)

To summarize the model, we solve Eqs. (5), (11), (17), supplemented by Eqs.
(15), (20). We resolve the disk’s thermal time (hence, allow for $\dot{E}_{\rm
th}\neq 0$), but assume that hydrostatic equilibrium is always maintained (as
enforced by Eq. (15)). A model is fully specified by the parameters:
$\\{M_{\bullet},M_{\star,0},R_{\star,0},\alpha,\chi,\eta,\beta,\mathcal{P}_{\rm
QPE}\gtrsim\mathcal{P}_{\rm QPE,min}\\}$ in addition to those parameters
$\\{M_{\rm\star,TDE},t_{\rm TDE}\\}$ which characterize the external source-
term $\dot{M}_{\rm TDE}(t)$. Table 1 summarizes the model variables and the
fiducial values adopted in the numerical models we present.

Figure 1: Schematic diagram showing the star-disk system and the long-term effects of Lense-Thirring precession on the stellar orbit. As long as the orbital precession period is rapid compared to the timescale over which the disk’s surface density evolves, the net angular momentum added to the disk by the mass stripped from the star (averaged over many precession cycles) will keep the disk aligned with the black hole spin axis.  Table 1: Model Parameters Symbol | Description | Fiducial Value(s)
---|---|---
$M_{\bullet}$ | SMBH mass | $10^{6}M_{\odot}$
$M_{\star,0}=\mathcal{M}_{\star}M_{\odot}$ | Initial EMRI star mass | $1M_{\odot}$
$R_{\star,0}=\mathcal{R}_{\star}R_{\odot}$ | Initial EMRI star radius | $1R_{\odot}$
$P_{\rm QPE}=P_{\rm orb}/2$ | QPE period | 4 hr
$\chi$ | Star-disk inclination | $\pi/4$
$\beta\propto\dot{E}_{\rm abl}/v_{\rm k}^{2}$ | Collision/ablation heating efficiency | 0…0.25
$\alpha$ | Disk viscosity parameter | 0.1
$\eta\propto\Delta M_{\star}/p_{\rm ram}$ | Collision stripping efficiency | $10^{-3}$
$M_{\rm\star,TDE}$ | Mass of disrupted star in TDE | $1M_{\odot}$
$t_{\rm TDE}$ | Fallback time of TDE debris | 30 d

##  3 Results

We present results for the coupled disk/star system evolution for an initial
EMRI similar to the Sun ($M_{\star,0}=M_{\odot};$ $R_{\star,0}=R_{\odot}$) on
an orbit of half-period $P_{\rm QPE}=4$ hr characteristic of QPE sources.
Assuming that a steady-state can be achieved in which mass accretes inwards
through the disk at the same rate it is stripped from the star (i.e.,
$\dot{M}_{\rm acc}\approx\dot{M}_{\star}$), the ratio of collision/ablation-
and accretion-heating rates (left two terms in Eq. (21)) takes the form:

$\frac{\dot{E}_{\rm abl}}{\dot{E}_{\rm acc}}\approx\frac{8}{3}\beta.$ (22)

We shall find that the system can exhibit qualitatively different evolution
depending on this ratio, i.e. on the value of $\beta$ (even though some
solutions only achieve the steady-state $\dot{M}_{\rm
acc}\approx\dot{M}_{\star}$ in a time-averaged sense).

We first consider the limit of negligible collision/ablation heating
($\beta=0$; Sec. 3.1.1) before considering models with progressively stronger
collision/ablation heating (larger $\beta$; Secs. 3.1.2). While for low values
of $\beta$ the disk is thermally unstable and exhibits limit-cycle behavior,
we find that for sufficiently large $\beta>\beta_{\rm stable}=3/16$
collision/ablation heating can act to stabilize the disk evolution (Sec. 3.2).
Finally, we consider a star on a long-period orbit, for which the disk is gas-
pressure dominated and hence thermally stable even for $\beta=0$ (Sec. 3.3).
The dependence of these different regimes on $\beta$ and $P_{\rm QPE}$ is
discussed in 4.1.

###  3.1 Thermally Unstable Evolution

####  3.1.1 Weak collision/ablation heating

In a standard steady-state accretion flow model (e.g., neglecting
collision/ablation heating), radiation-pressure dominates over gas-pressure in
the disk midplane at small radii (e.g., Linial & Metzger 2024),

$r_{0}\lesssim r_{\rm rad}\approx 5.4\times 10^{13}\,{\rm
cm}\,\alpha_{-1}^{2/21}M_{\bullet,6}^{23/21}\left(\frac{\dot{M}_{\rm
acc}}{0.1\dot{M}_{\rm Edd}}\right)^{16/21},$ (23)

or, correspondingly short QPE periods,

$P_{\rm QPE}\lesssim 30\,{\rm
hr}\,\alpha_{-1}^{1/7}M_{\bullet,6}^{8/7}\left(\frac{\dot{M}_{\rm
acc}}{0.1\dot{M}_{\rm Edd}}\right)^{8/7},$ (24)

where we have normalized the accretion rate to the Eddington value
$\dot{M}_{\rm Edd}\equiv L_{\rm Edd}/(\epsilon c^{2})$, $L_{\rm Edd}\simeq
1.5\times 10^{44}M_{\bullet,6}$ erg s-1 and $\epsilon=0.1$. As most QPE
sources have periods $\sim 2-20$ hr and $\dot{M}\gtrsim 0.1\dot{M}_{\rm Edd}$
(as inferred from their quiescent emission), radiation pressure can be
important in the disk near the collision annulus.

For $r_{0}<r_{\rm rad}$ and $\beta=0$ (i.e., the typical case absent stellar
collisions), Eqs. (15), (21) are known to admit a range of optically-thick
disk solutions, depending on the surface density $\Sigma$ (e.g., Chen et al.
1995). For low $\Sigma$ below a critical threshold value $\Sigma_{\rm min}$
the only solutions are radiatively-cooled ($\dot{E}_{+}\simeq\dot{E}_{\rm
rad}),$ and geometrically-thin $h/r\ll 1.$ The cooler solution is gas
pressure-dominated and thermally-stable, while the hotter solution is
radiation pressure-dominated and thermally unstable (e.g., Pringle & Rees
1972; Hirose et al. 2009; Jiang et al. 2019). On the other hand, for very
large $\Sigma$ above another critical value $\Sigma_{\rm max}$ the only
solution is advectively cooled ($\dot{E}_{+}\simeq\dot{E}_{\rm adv}),$
radiation-dominated, and geometrically-thick $h/r\sim 1$, i.e. even hotter
than the thermally-unstable radiatively-cooled branch (e.g., Abramowicz et al.
1988). At intermediate $\Sigma_{\rm min}<\Sigma<\Sigma_{\rm max}$, all three
solutions are present, though again only the coolest and hottest solutions are
stable (hereafter, the “low” and “high” branch, respectively).

The existence of thermal instabilities in nature, i.e. outside the
idealization of the Shakura & Sunyaev (1973) $\alpha$-model, remains under
debate. AGN observations (e.g., Done et al. 2007) and some global radiation
MHD simulations of radiation-dominated accretion disks do not find evidence
for thermal instabilities (e.g., Jiang et al. 2016, 2019). However, the
proposed physical mechanisms responsible for quenching the instability, such
as the iron opacity bump associated with high-metallicity gas (e.g., Jiang et
al. 2016) or the presence of dynamically-strong magnetic fields (e.g., Dexter
& Begelman 2019; Kaur et al. 2023) amplified and dragged in from larger scales
in the disk (e.g., Jacquemin-Ide et al. 2023), may not exist for the mass
stripped from what is likely a weakly magnetized star. While we therefore take
the presence of thermal instabilities at face value in this section, the
thermally-stable solutions obtained in Sec. 3.2 would be the relevant ones in
their absence.

Figure 2: Disk-star evolution for an initial sun-like star
($\mathcal{M}_{\star}=\mathcal{R}_{\star}=1$), neglecting collision/ablation
heating ($\beta=0$), for fiducial parameters ($M_{\bullet}=10^{6}\,\rm
M_{\odot}$, $P_{\rm QPE}=4$ hr; $\eta=10^{-3}$; $\alpha=0.1$). Top left panel:
Surface density (solid blue) and midplane temperature (dotted red) as a
function of time. Top right panel: $p_{\rm rad}/p_{\rm gas}$,
$\kappa/\kappa_{\rm es}$ and $h/r$ as a function of time. Gas pressure
dominates over radiation pressure and opacity is dominated by Kramer’s
throughout the “low” branch of the limit cycle. Bottom left panel: Stellar and
orbital properties as a function of time ($M_{\star}$, $R_{\star}$ and $P_{\rm
orb}$), from TDE to near depletion of the star. As expected for these
parameters, the stellar ablation rate greatly exceeds the rate of orbital
decay due to hydrodynamical drag. Bottom right panel: Time evolution of
instantaneous and cycle-averaged rates of accretion $\dot{M}_{\rm acc}$,
stellar mass-stripping $|\dot{M}_{\star}|$ and disk mass impacted by the star
$\dot{M}_{\rm coll}$, in comparison to the TDE fall-back rate $\dot{M}_{\rm
TDE}.$

Figure 3: Disk evolution for a sun-like star orbiting a $10^{6}M_{\odot}$
SMBH, in the absence of collision/ablation heating ($\beta=0$) and
$\eta=10^{-3}$. Left panel: Limit-cycle evolution in the $T$ vs. $\Sigma$
plane. The color-map shows the system’s evolution time along the limit cycle,
with its maximal value obtained towards $\Sigma_{\rm max}$ along the lower
thermal equilibrium branch. The black dashed contours correspond to different
values of $h/r$ (Eq. 16), and blue dotted contours correspond to $p_{\rm
rad}/p_{\rm gas}$. The red dotted line is the analytical expression for the
lower thermal equilibrium branch (Eq. 25). Right panel: The fraction of time
spent below $\Sigma$, as a function of $\Sigma$. Bottom panel: QPE eruption
luminosity $L_{\rm QPE}$ and duration $t_{\rm QPE}$ (expressed in ratio to the
QPE period $P_{\rm QPE}$), as a function of time, following the analytic
estimates of Eqs. (51), (50).

As illustrated by Fig. 2, the main feature of the disk evolution for $\beta=0$
is a limit cycle behavior in which $\Sigma$ undergoes periodic see-saw like
oscillations spanning the range $\Sigma_{\rm min}<\Sigma<\Sigma_{\rm max}$. In
each cycle, the disk surface density grows gradually from low
$\Sigma\approx\Sigma_{\rm min}\sim 10$ g cm-2 until reaching
$\Sigma\approx\Sigma_{\rm max}\sim 10^{4}$ g cm-3, before rapidly dropping
back to $\Sigma\approx\Sigma_{\rm min}$, starting the cycle again. The disk
spends most of the time in the gas pressure-dominated low state, in which it
is geometrically-thin $h/r\sim 10^{-3}$ and the viscous timescale
$\tau_{\nu}\propto(h/r)^{-2}$ (Eq. (10)) is effectively infinite. The disk
mass therefore builds up during this phase ($\dot{M}_{\rm
acc}\ll|\dot{M}_{\star}|+\dot{M}_{\rm TDE}$) due to a combination of mass
fall-back from the TDE (early times) and ablation from the orbiting EMRI (late
times). By contrast, once $\Sigma\approx\Sigma_{\rm max}$ the disk transitions
to the high state, which is geometrically thick $h/r\sim 1$ and the viscous
timescale is extremely short, causing rapid accretion ($\dot{M}_{\rm
acc}\gg|\dot{M}_{\star}|+\dot{M}_{\rm TDE})$, and rapidly reducing
$\Sigma\approx\Sigma_{\rm min}.$ The duration of each cycle is shortest at
early times, when the dominant source of mass being added to the disk is fall-
back accretion from the TDE (see Shen & Matzner 2014; Xiang-Gruess et al.
2016). However, throughout most of the disk’s lifetime $t\gtrsim 3$ yr, TDE
fall-back accretion (or whatever process first created the gaseous disk)
becomes subdominant compared to the mass added by stellar stripping,
generating regular cycles every $\approx 2$ yr. This cyclic evolution
continues until the star begins to lose an appreciable fraction of its total
mass and starts to expand near the end of its lifetime, accelerating its
destruction at around $t\approx 70$ yr.222In this particular example, the
EMRI’s evolution would actually terminate prior to complete ablation-induced
destruction, once it overflows its Roche-lobe onto the SMBH; this occurs once
$P_{\rm QPE}<P_{\rm
QPE,min}\propto\mathcal{R}_{\star}^{3/2}\mathcal{M}_{\star}^{-1/2}\propto\mathcal{M}_{\star}^{-1}$
(Eq. (2)) becomes satisfied, i.e. after only moderate fractional stellar mass-
loss (since $P_{\rm QPE}\gtrsim P_{\rm QPE,min}$ in the initial state for the
chosen model parameters).

The main features of the system evolution can be understood analytically.
Equating viscous heating with radiative cooling $\dot{E}_{\rm
rad}\approx\dot{E}_{\rm acc}$ under the assumptions $p_{\rm rad}\ll p_{\rm
gas}$ and $\kappa\simeq\kappa_{\rm K}$ gives the lower gas-pressure dominated
branch (plotted as red-dotted line in the left panel of Fig. 3)

$\left.k_{\rm B}T(\Sigma)\right|_{\rm low}=k_{\rm
B}\left(\frac{27}{8}\frac{A_{\kappa}\alpha}{ac}\left(\frac{v_{\rm
k}}{r_{0}}\right)^{2}\left(\frac{k_{\rm B}}{\mu m_{\rm
p}}\right)^{1/2}\Sigma^{3}\right)^{1/7}\simeq 15\,{\rm
eV}\;\left(\frac{\Sigma}{10^{4}\,\rm
g\,cm^{-2}}\right)^{3/7}\alpha_{-1}^{1/7}\mathcal{P}_{\rm QPE,4}^{-2/7}\,,$
(25)

and a corresponding scale height of

$\left.\frac{h}{r}\right|_{\rm low}\approx\sqrt{\frac{k_{\rm B}T}{\mu m_{\rm
p}v_{\rm k}^{2}}}\approx 1.5\times 10^{-3}\left(\frac{\Sigma}{10^{4}\,\rm
g\,cm^{-2}}\right)^{3/14}\alpha_{-1}^{1/14}\,\mathcal{P}_{\rm
QPE,4}^{4/21}M_{\bullet,6}^{-1/3},$ (26)

where $\alpha_{-1}\equiv\alpha/(0.1).$ The maximal value of $\Sigma$
corresponds to where this lower branch meets the upper radiation-pressure
dominated branch (for which $\kappa\simeq\kappa_{\rm T}$). Equating the
analytic estimates for these two roots of $\dot{E}_{\rm
rad}\approx\dot{E}_{\rm acc}$ show that they converge at a maximal surface
density:

$\Sigma_{\rm
max}\approx\mathcal{C}\cdot\frac{2^{33/12}}{2^{19/12}}\left(\frac{r_{0}/v_{k}}{\alpha^{11}}\right)^{1/12}\left(\frac{c^{11}\mu^{2}m_{\rm
p}^{2}}{A_{\kappa}^{4}k_{\rm B}^{2}\kappa_{\rm
es}^{7}a^{3}}\right)^{1/12}\approx 4\times 10^{4}\;{\rm
g\,cm^{-2}}\;\alpha_{-1}^{-11/12}\mathcal{P}_{\rm QPE,4}^{1/12},$ (27)

where we have included a prefactor $\mathcal{C}\approx 1/3$ to account for the
fact that the transition between branches occurs at a somewhat lower surface
density than where the two analytic branches intersect, resulting in better
consistency with the maximal $\Sigma$ obtained during each cycle in our
numerical calculation (Fig. 2).

When in the low-state, accretion is negligible and the disk’s surface density
grows as $\dot{\Sigma}\approx|\dot{M}_{\star}|/\pi r_{0}^{2}$ (at late times,
when $|\dot{M}_{\star}|\gg\dot{M}_{\rm TDE}$). The time spent at any
$\Sigma_{\rm min}<\Sigma<\Sigma_{\rm max}$ can thus be estimated as

$\displaystyle\tau_{\rm\Sigma}\equiv\frac{\Sigma}{\dot{\Sigma}}=\frac{P_{\rm
QPE}}{\eta}\frac{\pi R_{\star}^{2}\Sigma}{M_{\star}}\frac{p_{\star}}{p_{\rm
ram}}=\frac{P_{\rm
QPE}}{\eta}\left(\frac{r_{0}}{R_{\star}}\right)^{4}\frac{M_{\star}}{M_{\bullet}}\left(\frac{h}{r}\right)\approx
0.8\,{\rm
yr}\,\eta_{-3}^{-1}\left(\frac{h/r}{10^{-3}}\right)\frac{\mathcal{M}_{\star}}{\mathcal{R}_{\star}^{4}}\mathcal{P}_{\rm
QPE,4}^{11/3}M_{\bullet,6}^{1/3},$ (28)

where we have used Eqs. (5), (6) and $\eta_{-3}\equiv\eta/(10^{-3})$. Because
$\tau_{\Sigma}\propto h/r\propto\Sigma^{3/14}$ (Eq. (26)), the disk spends
roughly equal time at all $\Sigma_{\rm min}<\Sigma<\Sigma_{\rm max}$.
Evaluating Eq. (28) at $h/r$ (Eq. (26)) at $\Sigma\approx\Sigma_{\rm max}$
thus gives an estimate of the total cycle duration:

$\tau_{\rm cyc}\approx\tau_{\Sigma_{\rm max}}\approx 1.6\,{\rm
yr}\,\eta_{-3}^{-1}\alpha_{-1}^{-1/8}\frac{\mathcal{M}_{\star}}{\mathcal{R}_{\star}^{4}}\mathcal{P}_{\rm
QPE,4}^{31/8}.$ (29)

This roughly agrees the cycle duration $\approx 2$ yr found in our numerical
calculations (Fig. 2).

Both the mass accretion rate $\dot{M}_{\rm
acc}\propto\Sigma(h/r)^{2}\propto\Sigma^{10/7}$ peaks, and most of the total
time $\tau_{\Sigma}\propto\Sigma^{3/14}$ is spent, near the highest
$\Sigma\approx\Sigma_{\rm max}$ achieved during the cycle. The average
accretion rate in the low-state can therefore be estimated as

$\displaystyle\langle\dot{M}_{\rm low}\rangle\simeq 3\pi\alpha r_{0}v_{\rm
k}\Sigma_{\rm max}\left(\frac{h}{r}\right)_{\Sigma_{\rm max}}^{2}\approx
4\times 10^{-3}\,\dot{M}_{\rm Edd}\alpha_{-1}^{-1/6}\frac{\mathcal{P}_{\rm
QPE,4}^{5/6}}{M_{\bullet,6}}.$ (30)

During this phase, roughly $\left<\dot{M}_{\rm
low}\right>\times\tau_{\Sigma_{\rm max}}\approx 2\times 10^{-4}\,\rm
M_{\odot}$ of mass are accreted onto the SMBH.

However, this is substantially lower than the accretion rate averaged over the
many cycles, because most of the total mass is accreted during the rapid
“flushes” that occur in the brief advective high-state phase. The latter
average, more relevant to what ultimately feeds the black hole, can be
estimated as the typical disk mass built up during the cycle $\sim\pi
r_{0}^{2}\Sigma_{\rm max}$ ($\approx 10^{-2}\,M_{\odot}$ in this example)
divided by the cycle duration: (Eq. (29)),

$\displaystyle\langle\dot{M}_{\rm acc}\rangle\approx\frac{\pi
r_{0}^{2}\Sigma_{\rm max}}{\tau_{\rm cyc}}\approx 0.3\dot{M}_{\rm
Edd}\eta_{-3}\alpha_{-1}^{-19/24}\frac{\mathcal{R}_{\star}^{4}}{\mathcal{M}_{\star}}M_{\bullet,6}^{-1/3}\mathcal{P}_{\rm
QPE,4}^{-59/24}\,,$ (31)

again compatible with the average accretion rate from our numerical
calculation (Fig. 2).

Throughout most of the disk’s evolution the accreted mass originates from the
star (certainly after the initial TDE fallback rate has subsided below
$\langle\dot{M}_{\rm acc}\rangle$). The timescale for the star to be destroyed
can therefore be estimated as:

$\tau_{\rm dest}\sim\frac{3}{10}\frac{M_{\star,0}}{\langle\dot{M}_{\rm
acc}\rangle}\approx\tau_{\rm cyc}\left(\frac{M_{\star}}{\pi
r_{0}^{2}\Sigma_{\rm max}}\right)\approx 40\,{\rm
yr}\,\eta_{-3}^{-1}\alpha_{-1}^{19/24}\frac{\mathcal{M}_{\star}^{2}}{\mathcal{R}_{\star}^{4}}\frac{\mathcal{P}_{\rm
QPE,4}^{59/24}}{M_{\bullet,6}^{2/3}},$ (32)

where the factor $\approx 3/10$ accounts for the fact that $\tau_{\rm
dest}\propto\mathcal{M}_{\star}^{2}/\mathcal{R}_{\star}^{4}\propto\mathcal{M}_{\star}^{10/3}$
decreases as the star sheds mass and expands adiabatically (Eq. (8)), thereby
accelerating the destruction rate. Again, this agrees within a factor of
$\lesssim 2$ of the stellar lifetime $\tau_{\rm dest}\approx 70$ yr found in
our numerical calculation (Fig. 2).

####  3.1.2 Moderate collision/ablation heating

Figure 4: Same as left panel of Fig. 3, but showing evolution in the $T$ vs.
$\Sigma$ plane as the value of the collision/ablation heating parameter
$\beta$ is increased at constant logarithmic intervals between $10^{-3}$ to
$0.25$ (as labeled in color). The thermally-stable solution obtained for
$\beta=0.25>\beta_{\rm stable}=3/16$ is indicated by a red star.

The addition of collision/ablation heating acts to stabilize the disk
evolution, even when radiation pressure dominates. This is because while the
viscous heating $\dot{E}_{\rm acc}\propto(h/r)^{2}\propto p^{2}$ is a
sensitive increasing function of the disk’s midplane pressure $p$, the
collision/ablation heating rate $\dot{E}_{\rm
abl}\propto\dot{M}_{\star}\propto p_{\rm ram}\propto\Sigma/h\propto p^{-1}$
decreases with increasing midplane pressure $p$. In Appendix B we show that
for $\beta>\beta_{\rm stable}\approx 3/16$, collision/ablation heating is
sufficient to stabilize the disk in a steady-state $\dot{M}_{\rm
acc}\approx|\dot{M}_{\star}|$ (we illustrate this regime explicitly in Sec.
3.2).

However, for intermediate efficiencies $0\ll\beta\lesssim\beta_{\rm stable},$
collision/ablation heating is only able to stabilize the disk for a portion of
its limit cycle. This is illustrated by Fig. 4, which shows the limit cycle
evolution in the $\Sigma-T$ space, for models otherwise similar to the
$\beta=0$ case presented in the previous section but now for a range of
different $\beta\in[10^{-3},0.25].$ While the evolution qualitatively
resembles the zero collision/ablation-heating case ($\beta=0$; left panel of
Fig. 2), the ranges of the cycle shrink, e.g., $\Sigma_{\rm min}$ increases
and $\Sigma_{\rm max}$ decreases, with increasing $\beta$. For
$\beta\gtrsim\beta_{\rm stable}=3/16,$ the cycles cease all together and the
solution sets on a unique thermally-stable location (the $\beta=0.25$ case is
shown with an asterisk in Fig. 4), as described in the next section.

As in the $\beta=0$ case, our numerical results for $0\ll\beta<\beta_{\rm
stable}$ can also be understood analytically. In Appendix A, we present
estimates for $\Sigma_{\rm max},\langle\dot{M}_{\rm acc}\rangle$, $\tau_{\rm
dest}$, similar to those derived in Sec. 3.1.1 for $\beta=0$. The key physical
difference is that the maximum surface density in the low-state $\Sigma_{\rm
max}$ is no longer set by when the gas-pressure-dominated branch vanishes, but
instead by the surface density/temperature above which accretion heating
exceeds collision/ablation heating. In other words, while even moderate levels
collision/ablation heating can dominate over accretion heating at low
$\Sigma,T$ and allow the disk to evolve stably at the beginning of the cycle
(despite being radiation-dominated), as $\Sigma,T$ grows, accretion heating
eventually comes to dominate and the disk becomes unstable at $\Sigma_{\rm
max}$.

###  3.2 Thermally Stable Evolution

For sufficiently high $\beta\gtrsim\beta_{\rm stable}=3/16,$
collision/ablation heating is sufficient to stabilize the disk’s thermal
evolution. Our numerical calculation with $\beta=1$, shown in Fig. 5, confirm
this behavior. Although the disk still undergoes limit-cycle behavior at early
times when $\dot{M}_{\rm TDE}\approx\dot{M}_{\rm acc}\gg|\dot{M}_{\star}|,$
once $\dot{M}_{\rm acc}\approx|\dot{M}_{\star}|$ is achieved after
$t\gtrsim{20}$ yr, the accretion rate and surface density settle at roughly
constant values $\dot{M}_{\rm acc}\approx{\ 0.04\dot{M}_{\rm Edd}}$ and
$\Sigma\approx 10^{4}$ g cm-2, respectively. Both are intermediate to those
obtained in the unstable limit cycle evolution (analogous to the solution
shown with a star in Fig. 4). Radiation pressure dominates over gas pressure
in this state ($p_{\rm rad}/p_{\rm gas}\gtrsim 20$) and electron scattering
dominates the opacity, but unlike in the advectively-cooled high state, the
disk is radiatively-cooled and geometrically thin $h/r\sim{10^{-2}}$.

The disk properties in this “collision-supported” state can be derived
analytically. Unlike in the limit-cycle evolution, the accretion rate is
approximately constant once a steady state is reached, the disk scale-height
can be simply estimated by equating the accretion rate $\dot{M}_{\rm acc}$
(Eq. (9)) with the mass-stripping rate $|\dot{M}_{\star}|$ (Eq. (6)). Because
both scale linearly with the disk surface density, this gives an expression
for the disk scale-height:

$\frac{h}{r}=\left(\frac{2\eta\xi_{\rm
v}^{2}}{3\pi\alpha}\frac{R_{\star}^{4}}{r_{0}^{4}}\frac{M_{\bullet}}{M_{\star}}\right)^{1/3}\approx
1.1\times 10^{-2}\left(\frac{\eta_{-3}\xi_{\rm
v}^{2}}{\alpha_{-1}}\right)^{1/3}\frac{\mathcal{R}_{\star}^{4/3}}{\mathcal{M}_{\star}^{1/3}}\mathcal{P}_{\rm
QPE,4}^{-8/9}M_{\bullet,6}^{-1/9}\,.$ (33)

Likewise, thermal equilibrium ($dE_{\rm th}/dt=0;$ Eq. (17)) for a radiation-
supported disk can be written:

$\left(\frac{3}{8}+\beta\right)\dot{M}_{\rm acc}v_{\rm k}^{2}=\frac{L_{\rm
Edd}}{4}\left(\frac{h}{r}\right)=\frac{1}{4}\epsilon\dot{M}_{\rm
Edd}c^{2}\left(\frac{h}{r}\right),$ (34)

giving a second equation for the disk scale-height:

$\frac{h}{r}\simeq\frac{3}{2}\frac{(1+\frac{8}{3}\beta)}{\epsilon}\frac{R_{\rm
g}}{r_{0}}\frac{\dot{M}_{\rm acc}}{\dot{M}_{\rm Edd}},$ (35)

where $\epsilon\equiv L_{\rm Edd}/\dot{M}_{\rm Edd}c^{2}=0.1.$ Equating (33)
with (35) gives the accretion rate:

$\frac{\dot{M}_{\rm acc}}{\dot{M}_{\rm
Edd}}=\frac{2}{3}\frac{\epsilon}{(1+\frac{8}{3}\beta)}\frac{r_{0}}{R_{\rm
g}}\frac{h}{r}=\frac{0.07}{1+\frac{8}{3}\beta}\left(\frac{\eta_{-3}\xi_{\rm
v}^{2}}{\alpha_{-1}}\right)^{1/3}\frac{\mathcal{R}_{\star}^{4/3}}{\mathcal{M}_{\star}^{1/3}}\mathcal{P}_{\rm
QPE,4}^{-2/9}M_{\bullet,6}^{-7/9}.$ (36)

This equilibrium accretion rate is established on the viscous timescale (Eq.
(10)),

$\displaystyle\tau_{\nu}\simeq\frac{1}{3\pi}\frac{P_{\rm
QPE}}{\alpha}\left(\frac{h}{r}\right)^{-2}\approx 4.4\,{\rm
yr}\,\,\alpha_{-1}^{-1/3}\eta_{-3}^{-2/3}\xi_{\rm
v}^{-4/3}\frac{\mathcal{M}_{\star}^{2/3}}{\mathcal{R}_{\star}^{8/3}}\mathcal{P}_{\rm
QPE,4}^{25/9}M_{\bullet,6}^{2/9}\,.$ (37)

The steady-state disk surface density is given by

$\displaystyle\Sigma_{\rm TS}$ $\displaystyle=$
$\displaystyle\frac{\dot{M}_{\rm acc}}{3\pi\alpha r_{0}v_{\rm
k}(h/r)^{2}}\approx 2.4\times 10^{4}\,{\rm
g\,cm^{-2}}\,\frac{1}{1+\frac{8}{3}\beta}(\eta_{-3}\xi_{\rm
v}^{2})^{-1/3}\alpha_{-1}^{-2/3}\frac{\mathcal{M}_{\star}^{1/3}}{\mathcal{R}_{\star}^{4/3}}\frac{\mathcal{P}_{\rm
QPE,4}^{11/9}}{M_{\bullet,6}^{2/9}}.$ (38)

As a consistency check, this is seen to match Eq. (A3) for $\beta=\beta_{\rm
stable}=3/16$.

Accreting at $\dot{M}_{\rm acc}$ (Eq. (36)), the star destruction time is
approximately given by

$\tau_{\rm dest}\approx\frac{9}{16}\frac{M_{\star,0}}{\dot{M}_{\rm
acc}}\approx 320\,{\rm
yr}\,\left(1+\frac{8}{3}\beta\right)\left(\frac{\alpha_{-1}}{\eta_{-3}\xi_{\rm
v}^{2}}\right)^{1/3}\frac{\mathcal{M}_{\star}^{4/3}}{\mathcal{R}_{\star}^{4/3}}\mathcal{P}_{\rm
QPE,4}^{2/9}M_{\bullet,6}^{-2/9}\,.$ (39)

Again, the factor $\approx 9/16$ accounts for the fact that $\tau_{\rm
dest}\propto\mathcal{M}_{\star}^{16/9}$ decreases as the star loses mass.

Figure 5: Same as Fig. 2, but in the presence of substantial
collisions/ablation heating, with $\beta=0.25$. Top left panel: Surface
density (solid blue) and midplane temperature (dotted red) initially fluctuate
when accretion heating due to the TDE fallback is the dominant energy source,
but settle onto a thermally stable evolution throughout most of the system’s
lifetime after collisional/ablation heating from the EMRI comes to dominate.
Dashed blue horizontal line is the analytical estimate of the equilibrium
$\Sigma_{\rm TS}$ (Eq. (38)), evaluated for the initial stellar and orbital
properties. Top right panel: As opposed to the low branch of the $\beta=0$
regime, here $p_{\rm rad}\gg p_{\rm gas}$ and $\kappa\approx\kappa_{\rm es}$
during the thermally stable phase of the system’s evolution. Bottom left
panel: Stellar and orbital properties as a function of time. The star is
completely ablated after $\sim 600\,\rm yr$, in agreement with the analytical
estimate of Eq. (39). Bottom right panel: Instantaneous rates of accretion
$\dot{M}_{\rm acc}$, stellar ablation $|\dot{M}_{\star}|$ and disk mass
impacted by the star $\dot{M}_{\rm coll}$, in comparison to the TDE fallback
rate $\dot{M}_{\rm TDE}$. The transition to a thermally stable evolution
occurs roughly when the TDE fallback rate drops below the equilibrium
accretion rate (Eq. 36).

###  3.3 Thermally Stable Evolution at Long Periods

Thermally stable evolution may be achieved even in the absence of
collision/ablation heating ($\beta\approx 0$), for sufficiently long orbital
periods. At these wider radii, $p_{\rm rad}\lesssim p_{\rm gas}$, and the
system evolves to a radiatively cooling equilibrium, being fed by the mass
collisionally stripped from the star. Fig. 6 summarizes the key timescales of
the star-disk system as a function of $\beta$ and $\mathcal{P}_{\rm QPE}$, for
$R_{\star}=2R_{\odot}$ and otherwise fiducial parameters. For $\beta\lesssim
3/16$, we see that thermal stability is obtained for $P_{\rm orb}\simeq
2P_{\rm QPE}\gtrsim 18\,\rm hr$.

Figure 6: Contours show various key timescales as a function of (half) the
orbital period $P_{\rm QPE}$ and the collisional/ablation heating efficiency
$\beta$. The periods of some known QPE sources are shown along the top axis.
Top left panel: Cycle duration (in the case of a thermally unstable solution).
For $\beta>\beta_{\rm stable}$ or $p_{\rm gas}>p_{\rm rad}$, the system
evolves to a thermal equilibrium, and limit cycles disappear (the boundary
denoted by the yellow line). Top right panel: Stellar destruction time due to
mass-stripping from disk collisions. Bottom left panel: Migration time due to
disk-induced drag. Bottom right panel: Disk viscous time.

##  4 Implications

###  4.1 Secular Evolution of the EMRI Orbit

Our model neglects any changes to the stellar orbit, due to hydrodynamical or
gravitational interactions with the disk, implicitly assuming that they occur
slower than the star’s destruction. We now explore when this is a good
approximation.

The timescale required for the stellar orbit to undergo appreciable radial
migration, or to be ground down into the disk midplane, via gas drag is
approximately given by that required to intercept its own mass in disk
material, $\tau_{\rm mig}$ (e.g., MacLeod & Lin 2020; Generozov & Perets 2023;
Linial & Quataert 2024). In cases when the disk exhibits limit-cycle evolution
and spends most of its evolution near the surface density $\Sigma_{\rm max}$
(Sec. 3.1.1, 3.1.2), the star sweeps up (or rather, collides with) mass at the
rate

$\dot{M}_{\rm coll}\approx\frac{\pi R_{\star}^{2}\Sigma_{\rm
max}f_{\Sigma_{\rm max}}}{P_{\rm QPE}},$ (40)

where the factor $f_{\rm\Sigma_{\rm max}}\sim{\rm ln}\left(\Sigma_{\rm
max}/\Sigma_{\rm min}\right)^{-1}\approx 0.1$ accounts for the fraction of the
limit-cycle spent at $\Sigma\approx\Sigma_{\rm max}$. The migration timescale
is thus given by:

$\displaystyle\tau_{\rm mig}\approx\frac{M_{\star}}{\dot{M}_{\rm
coll}}\approx\frac{r_{0}^{2}}{R_{\star}^{2}}\frac{\tau_{\rm dest}}{\tau_{\rm
cyc}}\frac{P_{\rm QPE}}{f_{\Sigma_{\rm max}}},$ (41)

and hence relative to the destruction time,

$\frac{\tau_{\rm mig}}{\tau_{\rm
dest}}\approx\frac{r_{0}^{2}}{R_{\star}^{2}}\frac{P_{\rm QPE}}{f_{\Sigma_{\rm
max}}\tau_{\rm cyc}}\approx
120\eta_{-3}\alpha_{-1}^{1/8}\frac{\mathcal{R}_{\star}^{2}}{\mathcal{M}_{\star}}\frac{M_{\bullet,6}^{2/3}}{\mathcal{P}_{\rm
QPE,4}^{37/24}},$ (42)

where we have made use of the fact that $\tau_{\rm dest}\approx
M_{\star}/\langle\dot{M}_{\rm acc}\rangle$ and $\langle\dot{M}_{\rm
acc}\rangle\approx\pi r_{0}^{2}\Sigma_{\rm max}/\tau_{\rm cyc}$ and in the
final equality have used Eq. (29) for $\tau_{\rm cyc}$ in the $\beta=0$ limit
(qualitatively similar results are obtained using Eq. (A6) for $\tau_{\rm
cyc}$). The above ratio only grows as the star loses mass, as
$R_{\star}^{2}/M_{\star}\propto M_{\star}^{-5/3}$ for adiabatic expansion
$R_{\star}\propto M_{\star}^{-1/3}$ (Eq. (8)).

In the thermally-stable case ($\beta>3/16$; Sec. 3.2), the surface density of
the disk is approximately constant in time $\Sigma\approx\Sigma_{\rm TS}$ (Eq.
(38)) and the migration to destruction timescale ratio be written (now taking
$f_{\Sigma_{\rm max}}=1$)

$\frac{\tau_{\rm mig}}{\tau_{\rm dest}}\approx\frac{M_{\star}P_{\rm QPE}}{\pi
R_{\star}^{2}\Sigma_{\rm TS}\tau_{\rm dest}}\approx 9(\eta_{-3}\xi_{\rm
v}^{2})^{2/3}\alpha_{-1}^{1/3}\frac{\mathcal{R}_{\star}^{2/3}}{\mathcal{M}_{\star}^{2/3}}\frac{M_{\bullet,6}^{4/9}}{\mathcal{P}_{\rm
QPE,4}^{4/9}},$ (43)

where we have used Eq. (39) for $\tau_{\rm dest}$. As in the thermally
unstable regime, the timescale ratio $\tau_{\rm mig}/\tau_{\rm dest}$
increases as the star loses mass and expands, as the ratio scales as
$(R_{\star}/M_{\star})^{2/3}\propto M_{\star}^{-8/9}$.

Thus, solar-type stars $\mathcal{M}_{\star}\sim\mathcal{R}_{\star}\sim 1$ on
QPE-like orbital periods likely have $\tau_{\rm mig}>\tau_{\rm dest}$ for
$\eta\gtrsim 10^{-4}-10^{-3}$, indicating the star will be destroyed prior to
significant evolution in its orbit and justifying our model’s neglect of
orbital evolution. On the other hand, for smaller $\eta$ or larger orbital
periods ${P}_{\rm orb}\gtrsim$ days, the star may be ground into the disk
prior to being destroyed. Once the star’s angular momentum aligns with that of
the disk, strong collisions between the two will cease and the gaseous disk
will quickly be depleted of mass.

Much more compact stars than the sun, such as M-dwarfs
($\mathcal{R}_{\star}\approx 0.1$) or white dwarfs
($\mathcal{R}_{\star}\approx 0.01$), can also have $\tau_{\rm mig}\gg\tau_{\rm
dest}$. However, in an absolute sense their destruction times $\tau_{\rm
dest}\propto\mathcal{R}_{\star}^{-4}$ (Eq. (32)) can be orders of magnitude
longer than even the gravitational-wave inspiral time, and the gaseous disks
and associated QPE-like eruptions they create would be far less conspicuous.

In addition to hydrodynamical drag, the star and the disk will exert
gravitational torques on each other, potentially leading to precession and
(with dissipation) eventual alignment of their orbital planes. These torques
will act on a characteristic timescale

$\tau_{\rm grav}\approx\frac{\sqrt{GM_{\bullet}r_{0}}}{GM_{\star}/r_{0}}\sim
P_{\rm QPE}\left(\frac{M_{\bullet}}{M_{\star}}\right)\approx 500\,{\rm
yr}\;\frac{M_{\bullet,6}\mathcal{P}_{\rm QPE,4}}{\mathcal{M}_{\star}}\,,$ (44)

where the mass of the disk has been safely neglected relative to that of the
star. Though comparable to the stellar destruction time $\tau_{\rm dest}$, the
timescale for gravitational torques is much slower than other processes
shaping the stellar/disk system such as Lense-Thirring precession (Eq. (3)),
and so we do not expect gravitational torques to significantly alter the
system evolution.

As summarized in Table 2 and Fig. 6, for fiducial parameters and solar-type
stars, we find a typical hierarchy of timescales:

$P_{\rm QPE}\ll\tau_{\rm prec}\ll\tau_{\rm cyc}\ll\tau_{\rm dest}\ll\tau_{\rm
mig}.$ (45)

As mentioned in Sec. 2, the gaseous disk created by stellar mass-stripping
will not align with the stellar orbit as long as the Lense-Thirring precession
time is sufficiently rapid relative to the timescale for the disk to build-up
mass near the minimum of the cycle, $\tau_{\Sigma}(\Sigma_{\rm
min})\sim\tau_{\rm cyc}/10$ (Eq. (28)), a condition that is satisfied for
large QPE periods, low SMBH masses, and/or high SMBH spin (see also Sec. 4.4).
This is because only the component of the angular momentum added by the star
which is parallel to the black hole spin axis will survive the precession-
averaging process (see Fig. 1).

Stated mathematically, if the disk is aligned with the black hole spin along
the $\hat{z}$ direction, then the change in the total angular momentum of the
disk $\vec{J}_{\rm disk}$ due to the addition of the stripped mass carrying
the specific angular momentum of the stellar orbit $\vec{j}_{\star}$ roughly
obeys,

$\frac{d\vec{J}_{\rm
disk}}{dt}\simeq\langle\vec{j}_{\star}\rangle|\dot{M}_{\star}|\underset{\tau_{\rm
prec}\ll\tau_{\rm
cyc}/10}{\simeq}j_{\star,z}|\dot{M}_{\star}|\hat{z}=j_{\star}\cos\chi|\dot{M}_{\star}|\hat{z},$
(46)

since the components of $\vec{j}_{\star}$ perpendicular to the black hole spin
axis (e.g., $\hat{x},\hat{y}$) cancel over multiple precession cycles (where
$\langle...\rangle$ denotes a time-average). The order-unity factor
$\cos\chi<1$ implies that the stripped gas will “circularize” at a radius
$r_{\rm c}\simeq r_{0}\cos^{2}\chi$, slightly smaller than the collision
radius $r_{0}$; however, outwards viscous spreading of the mass deposited near
$r_{\rm c}$ (e.g., Metzger et al. 2012) will result in broadly similar disk
properties as at $r_{0}$.

###  4.2 Eccentric orbits

Stars migrate to orbits of interest (i.e., those resulting in QPE flares, or
generally interacting with TDE disks) primarily through gravitational-wave
circularization from highly eccentric orbits of semi-major axis $a\gg r_{\rm
p}$ (not to be confused with the SMBH spin parameter $a_{\bullet}$).
Throughout the GW-dominated circularization phase, the stellar pericenter
distance remains nearly constant at $r_{\rm p}\approx r_{\rm p,0}$. Thus, no
significant star-disk interaction occurs during the GW inspiral phase, up
until the orbit is mostly circularized, provided that $r_{\rm p,0}\gg r_{\rm
T}$, where $r_{\rm T}$ is the tidal radii of the coincident TDEs progenitors.
In this regime, the remaining eccentricity when TDE-triggered star-disk
interaction first occurs is modest ($e\lesssim 0.5$), and our previous
assumptions and estimates are generally valid.

If, however, $r_{\rm p,0}\gtrsim r_{\rm T}$, star-disk interaction may first
ensue while the orbit is still highly eccentric. This scenario may also
transpire if an AGN disk is present, extending to radii much greater than the
size of a typical TDE disk (e.g., Linial & Quataert, 2024). While the radius
of the collision site between the star and the disk varies due to apsidal
precession, the average (over the precession cycle) star-disk collision occurs
at around $2r_{\rm p}$ (Linial & Quataert, 2024), with collision properties
similar to those of a circular EMRI at radius $r_{0}\approx r_{\rm p}$ and
$P_{\rm QPE}\approx(P_{\rm orb}/2)(r_{\rm p}/a)^{3/2}$. Quantities that depend
on the collision conditions, for example, $\Delta M_{\star}$, are essentially
equivalent to the circular case, up to order unity corrections.

Interestingly, the typical timescale hierarchy diskussed in §4.1 is not
necessarily maintained for highly eccentric orbits. For example, the stellar
destruction time $\tau_{\rm dest}\approx(M_{\star}/\Delta M_{\star})P_{\rm
orb}$ may be longer than the drag-induced circularization timescale

$\tau_{\rm mig}\equiv\tau_{\rm
circ}\approx\frac{a}{|\dot{a}|}\approx\frac{M_{\star}}{\Delta M_{\rm
coll}}\left(\frac{r_{\rm p}}{a}\right)P_{\rm orb}\,.$ (47)

with their ratio scaling as

$\frac{\tau_{\rm circ}}{\tau_{\rm dest}}\approx\alpha^{1/3}\left(\frac{\eta
R_{\star}M_{\bullet}}{r_{\rm p}M_{\star}}\right)^{2/3}\left(\frac{r_{\rm
p}}{a}\right)^{1/2}\,.$ (48)

The reversal of timescale hierarchy as $a\to r_{\rm p}$ suggests that the
stellar destruction and circularization rates become comparable $\tau_{\rm
dest}\approx\tau_{\rm circ}$ at some critical eccentricity, of roughly

$\left.1-e\right|_{\tau_{\rm circ}=\tau_{\rm dest}}\approx
0.2\;\left(\frac{r_{\rm p}}{r_{\rm
T}}\right)^{4/3}\eta_{-3}^{-4/3}\alpha_{-1}^{2/3}M_{\bullet,6}^{-8/9}\mathcal{M}_{\star}^{4/3}\,,$
(49)

suggesting mild residual eccentricity when the stellar ablation exceeds the
circularization rate (as $r_{p}\gtrsim r_{\rm T}$).

Table 2: Key Timescales ($\mathcal{M}_{\star}=\mathcal{R}_{\star}=1$; $\alpha=0.1$; $\eta=10^{-3}$; $\xi_{\rm v}=1$)  Symbol | $\beta$ | Description | Value | Equation
---|---|---|---|---
$P_{\rm QPE}$ | - | QPE period | $4\,{\rm hr}\;\mathcal{P}_{\rm QPE,4}\gtrsim P_{\rm QPE,min}$ | Eq. (2)
$\tau_{\rm prec}$ | - | Nodal precession period | $1.4\,{\rm yr}\;\mathcal{P}_{\rm QPE,4}^{2}M_{\bullet,6}(a_{\bullet}/0.3)^{-1}$ | Eq. (3)
$\tau_{\rm cyc}$ | $\beta=0$ | Limit-cycle period | $7\,{\rm yr}\,\;\mathcal{P}_{\rm QPE,4}^{3.88}$ | Eq. (29)
… | $0\ll\beta<3/16$ | … | $7\,{\rm yr}\,\;(\beta/0.5)^{1/3}\mathcal{P}_{\rm QPE,4}^{2.78}M_{\bullet,6}^{2/9}$ | Eq. (A6)
$\tau_{\nu}$ | $\beta>3/16$ | Viscous time | $4\,{\rm yr}\,\;\mathcal{P}_{\rm QPE,4}^{2.78}M_{\bullet,6}^{2/9}$ | Eq. (37)
$\tau_{\rm dest}$ | $\beta=0$ | EMRI destruction time | $100\,{\rm yr}\,\;\mathcal{P}_{\rm QPE,4}^{2.46}M_{\bullet,6}^{-2/3}$ | Eq. (32)
… | $0\ll\beta<3/16$ | … | $310\,{\rm yr}\,(\beta/0.1)^{2/3}\mathcal{P}_{\rm QPE,4}^{2/9}M_{\bullet,6}^{-2/9}$ | Eq. (A8)
… | $\beta>3/16$ | … | $300\,{\rm yr}\,\;(1+\frac{8}{3}\beta)\mathcal{P}_{\rm QPE,4}^{2/9}M_{\bullet,6}^{-2/9}$ | Eq. (39)
$\tau_{\rm grav}$ | - | Gravitational torquing | $500\,{\rm yr}\;\mathcal{P}_{\rm QPE,4}M_{\bullet,6}$ | Eq. (44)
$\tau_{\rm mig}$ | $\beta=0$ | Drag-induced orbital migration | $6\times 10^{3}\,{\rm yr}\,\mathcal{P}_{\rm QPE,4}^{0.92}$ | Eq. (42)
… | $0\ll\beta<3/16$ | … | $6\times 10^{4}\,{\rm yr}\,(\beta/0.1)^{2/3}\mathcal{P}_{\rm QPE,4}^{-1.32}M_{\bullet,6}^{4/9}$ | …
… | $\beta>3/16$ | … | $2.7\times 10^{3}\,{\rm yr}\,\;(1+\frac{8}{3}\beta)\mathcal{P}_{\rm QPE,4}^{-2/9}M_{\bullet,6}^{2/9}$ | Eq. (43)

###  4.3 Secular Evolution of QPE Quiescent Emission and Eruptions

The time-averaged accretion rates predicted by our model $\dot{M}_{\rm
acc}\gtrsim 0.1\dot{M}_{\rm Edd}$ broadly agree with those of observed QPE
quiescent sources. However, depending on the disk/collision radius and the
(uncertain) value of $\beta$, our calculations reveal that these rates can be
achieved through a wide diversity of time-evolution behavior$-$ranging from
limit cycles to quasi-steady accretion. It is furthermore not realistic to
expect the one-zone model we have presented to capture the entire radially-
dependent disk structure, which determines for example how faithfully the
mass-inflow rate near the collision annulus reflects that reaching the central
SMBH and setting the quiescent X-ray emission. Nevertheless, qualitative
features of the black-hole accretion-rate evolution, such as characteristic
timescales and the potential for both rising and falling light curves, are
likely to be more robust.

For thermally unstable evolution ($\beta<3/16$, small $r_{0}<r_{\rm rad}$),
significant changes in the disk accretion rate occur on the limit-cycle
duration $\tau_{\rm cyc}$ (Eq. (29), (A6)). The latter are typically years to
decades, depending most sensitively on the mass-stripping efficiency $\eta$,
the stellar radius, and the QPE period. Even for thermally-stable accretion
($\beta>3/16$, large $r_{0}>r_{\rm rad}$), gradual long-term variability in
the mass accretion rate may still occur (e.g. due to viscous instabilities;
Lightman & Eardley 1974), on the viscous timescale (Eq. (37)) which is also
years to decades. As with the average accretion rates, these diversity of
timescales match those observed in QPE sources. The quiescent light curves in
GSN 069 (Miniutti et al., 2019, 2023a), XMMSL1 J024916.6-041244 (Chakraborty
et al., 2021), eRO-QPE3, eRO-QPE4 (Arcodia et al., 2024a) exhibit significant
(factor of $\gtrsim\rm few$) changes on timescales from years to decades,
including both rising and decaying emission (or phases of both in GSN 069;
Miniutti et al. 2023b), broadly consistent with our thermally-unstable
accretion solutions. The peak in the mass-accretion rate during the final
“flushing” stage of each cycle evolves even faster in our models (on a
timescale as short as weeks). In contrast, the apparently constant quiescent
X-ray flux of RXJ 1301.9+2747 (Giustini et al., 2020) over 20 years, and of
eRO-QPE2 for $\sim 3.5$ years (Arcodia et al., in preparation), may point to
collisionally-stabilized disk evolution in this system.

Insofar that QPE eruptions are generated by the hot expanding debris clouds
created by disk-star collisions (LM23), the eruption properties themselves
would be expected to co-evolve along with the quiescent disk. This is indeed
supported by observations which show secular changes in the QPE properties
over multiple epochs in parallel with changes in the quiescent emission
(Miniutti et al., 2023a; Chakraborty et al., 2024; Arcodia et al., 2024a;
Pasham et al., 2024). LM23 estimate the duration and luminosities of the
eruptions in terms of the disk properties near the collision annulus:

$t_{\rm QPE}\approx\left(\frac{\Delta M_{\rm coll}\kappa_{\rm T}}{4\pi cv_{\rm
k}}\right)^{1/2}\approx 0.13\,{\rm
hr}\,(1+\frac{8}{3}\beta)^{-1/2}(\eta_{-3}\xi_{\rm
v}^{2})^{-1/6}\alpha_{-1}^{-1/3}\mathcal{M}_{\star}^{1/6}\mathcal{R}_{\star}^{1/3}\mathcal{P}_{\rm
QPE,4}^{7/9}M_{\bullet,6}^{-5/18};$ (50) $L_{\rm QPE}\approx\frac{E_{\rm
ej}(R_{\star}^{2}h)^{1/3}}{3v_{\rm k}t_{\rm QPE}^{2}}\approx\frac{1}{6}L_{\rm
Edd}\frac{(R_{\star}^{2}h)^{1/3}}{r_{0}}\approx 1.6\times 10^{41}\,{\rm
erg\,s^{-1}}\left(\frac{\eta_{-3}\xi_{\rm
v}^{2}}{\alpha_{-1}}\right)^{1/9}\frac{\mathcal{R}_{\star}^{10/9}}{\mathcal{M}_{\star}^{1/9}}\mathcal{P}_{\rm
QPE,4}^{-20/27}M_{\bullet,6}^{20/27}.$ (51)

Here, $\Delta M_{\rm coll}\simeq 2\pi R_{\star}^{2}\Sigma$ (Eq. (13)) and
$E_{\rm ej}\simeq\Delta M_{\rm coll}v_{\rm k}^{2}/2$ are the mass and thermal
energy, respectively, of disk material shocked by the star in a roughly head-
on collision ($\chi\sim 1$), which expands above and below the midplane,
powering the observed eruptions. In the final equality we have made use of our
analytic estimates for the scale-height and surface density of thermally-
stable disks (Eqs. (33),(38)).

The bottom panel of Fig. 3 shows how these eruption properties evolve in time
for the fiducial $\beta=0$ model. The eruption luminosities $L_{\rm QPE}$ and
duty cycles $t_{\rm QPE}/P_{\rm QPE}$ vary by almost an order of magnitude
over a single cycle, illustrating how changes in the quiescent disk properties
(locally, at $r_{0}$) could well manifest through secular evolution of the
observed eruption properties. The cycle time-averaged values $\langle L_{\rm
QPE}\rangle\sim 10^{41}$ erg s-1, $\langle t_{\rm QPE}/P_{\rm QPE}\rangle\sim
10^{-2}$ are each roughly an order of magnitude lower than those of QPE
eruptions ($L_{\rm QPE}\sim 10^{42}$ erg s-1; $t_{\rm QPE}/P_{\rm QPE}\sim
0.1$). This may point to a larger effective stellar radius (for purposes of
interacting with the disk) or a more massive SMBH than assumed in the fiducial
model, or inaccuracies in the analytic estimates.

Several effects may lead to an increase in the effective star-disk cross
section producing the QPE flares. The star’s outer layers may become inflated
due to additional heating that may occur from the fraction of the collisional
heating rate $\Delta M_{\rm coll}v_{\rm rel}^{2}/P_{\rm QPE}$ deposited inside
the star. The flare timing pattern observed in some of the known QPE systems
indicate mild orbital eccentricities $e\approx 0.1$, suggesting significant
tidal heating, possibly inflating the star’s radius to about twice its
unperturbed, main-sequence size. Repeated ablation of the star’s outer layers
may also indirectly increase the effective collision area. The bulk of the
stripped material, $\Delta M_{\star}$, escapes the star at velocities
comparable to $v_{\rm esc,\star}\approx\sqrt{GM_{\star}/R_{\star}}$. During
the time between two consecutive collisions, the ablated material propagates a
distance of order $v_{\rm esc,\star}P_{\rm QPE}\approx R_{\star}$, where we
assumed $P_{\rm QPE}\gtrsim P_{\rm QPE,min}$. In an average sense, the star is
thus shrouded by a substantial layer of ablated material, which itself can
intersect the disk. Given that typically $\Delta M_{\star}\gg\Delta M_{\rm
coll}$ and that the extent of the bulk of stripped material only doubles in
radius, the effective surface density of the impacting shroud exceeds the
disk’s $\Sigma$. Hence, the ablated material is not expected to be
substantially decelerated and incorporated into the disk in just a single
collision, but rather over multiple orbital periods. This suggests that a
complex flow composed of liberated disk material, ablated stellar debris, and
possibly tidally stripped stellar material all evolve under the combined
gravity of the SMBH and the star itself, governing the timing and emission
properties of the resulting flares. In this case, the effective stellar radius
may not be limited by tidal gravity to the Roche lobe size, $R_{\rm RL}\approx
0.5\,r_{0}(M_{\star}/M_{\bullet})^{1/3}\approx
1.1\,R_{\odot}\,\mathcal{P}_{\rm QPE,4}^{2/3}\mathcal{M}_{\star}^{1/3}$.

As discussed in Sec. 4.1, evolution of the EMRI’s orbital semi-major axis or
eccentricity is expected to be very slow (e.g., $\tau_{\rm mig}\gtrsim 10^{3}$
yr). Even in the presence of rapid nodal precession, the timing of disk
collisions and hence the pattern of the twice-per-orbit periodic eruptions
should be relatively steady over timescales of years if the accretion disk
remains flat with its midplane aligned with the black hole spin axis (e.g.,
Fig. 1). However, if such alignment cannot occur or is incomplete, then the
disk (or just the collision annulus) will also precess about the black hole
spin axis, on a characteristic timescale comparable or shorter than $\tau_{\rm
prec}$ (Eq. (3)). This would give rise to a more complex temporal spacing of
the eruptions (e.g., Franchini et al. 2023), which may be required to explain
those QPE sources which exhibit less regular timing properties. Indeed,
temporary disk-spin alignment or warping could occur as substantial spin-
misaligned angular momentum is added to the disk by stellar mass-ablation at
the beginning of a limit-cycle (when $\Sigma$ is minimal), on a timescale
shorter than the precession period (see also Xiang-Gruess et al. 2016). We
speculate this could explain why the repetition pattern of the eruptions from
GSN 069, after their reappearance following the second quiescent-outburst,
temporarily differed (it remains to be seen how well the repetition pattern
will settled back into that observed prior to the outburst; G. Minutti,
private communication).

A key property that enables QPEs to be detected is the higher temperatures of
the eruptions $k_{\rm B}T_{\rm QPE}\sim 100-200$ eV (Miniutti et al., 2019;
Giustini et al., 2020; Arcodia et al., 2021; Chakraborty et al., 2021; Arcodia
et al., 2022; Miniutti et al., 2023b; Webbe & Young, 2023) relative to those
of the quiescent disk ($k_{\rm B}T\lesssim 50$ eV). LM23 found that the value
of $k_{\rm B}T_{\rm QPE}$ radiated by the debris cloud created from a disk-
star collision is sensitive to the disk’s midplane density, insofar the rate
of photon production within the hot shocked disk material scales as the
density squared. The fact that our thermally-unstable models reveal that
$\Sigma$ can vary by up to several orders of magnitude over its evolution
(Figs. 3, 4) therefore also implies variability in the eruptions
temperatures$-$and hence in the detectability of QPE activity$-$over
timescales of years to decades. This may lead to windows of time where QPE
eruptions temporarily “disappear”, if local thermal-instabilities indeed
manifest. Such behavior is consistent with the cessation of QPE eruptions in
GSN069 around the time of the second quiescent outburst (Miniutti et al.,
2023b). However, even for the most extreme limit cycles present in the
$\beta=0$ case, the disk spends over half of its time within a factor of a few
of $\Sigma_{\rm max}$ (right panel of Fig. 3).

###  4.4 Accretion Activity in Galactic Nuclei and QPE Demographics

From Eq. (12), we see that fall-back accretion after a TDE remains above a
given Eddington ratio $\dot{M}_{\rm TDE}/\dot{M}_{\rm Edd}$ for a timescale:

$t_{\rm Edd}\simeq t_{\rm TDE}\left(\frac{M_{\rm\star,TDE}}{3t_{\rm
TDE}\dot{M}_{\rm TDE}}\right)^{3/5}\approx 6.8\,{\rm yr}\,\left(\frac{t_{\rm
TDE}}{30\,{\rm
d}}\right)^{2/5}\left(\frac{M_{\rm\star,TDE}}{M_{\odot}}\right)^{3/5}\left(\frac{\dot{M}_{\rm
TDE}}{0.1\dot{M}_{\rm Edd}}\right)^{-3/5}.$ (52)

By contrast, we have shown that a disk fed by a stellar-EMRI can maintain an
accretion rate $\dot{M}_{\rm acc}\sim 0.1\dot{M}_{\rm Edd}$ (Eqs. (31),(36))
for a timescale $\tau_{\rm dest}\sim 10^{2}-10^{3}$ yr (Eqs. (32),(A8),(39)),
i.e. a factor of $\sim 10-10^{2}$ times longer than $t_{\rm Edd}$. Thus, even
if stellar EMRIs occur $\sim 100$ times less frequently than TDEs in galactic
nuclei (e.g., Linial & Sari 2023), their contribution to accretion activity in
otherwise inactive galactic nuclei (i.e., those without large-scale AGN) can
be comparable to that of TDEs.333The tidal disruption of giant stars has been
argued to contribute $\sim 10\%$ of the total TDE rate (MacLeod et al., 2012)
and to last for tens of hundreds of years, comparable to the disk lifetimes
estimated in this work. However, the susceptibility of their accretion streams
to disruption by dense gas in the nucleus (Bonnerot et al., 2016) remains an
uncertainty. However, unlike a TDE, such a source of ionizing UV/X-ray
radiation lasting hundreds or thousands of years, will create a more
spatially-extended narrow-line region in the nuclei of galaxies hosting
interacting star/disk systems. This may be compatible with the estimated long
duration of the accretion activity in QPE sources (e.g., Wevers et al. 2023;
Patra et al. 2023).

The per galaxy abundance of QPE sources has been recently constrained by
Arcodia et al. (2024b) to be of order $\sim 10^{-5}\,\rm gal^{-1}$, which
corresponds to a volumetric rate of $0.6\times 10^{-7}(\tau_{\rm dest}/10\,\rm
yr)^{-1}\;Mpc^{-3}\,yr^{-1}$. Thus, the relatively long QPE source lifetimes
found by our calculations $\tau_{\rm dest}\sim 10^{2}-10^{3}$ yr imply that a
lower formation rate of QPE sources is needed to explain their observed
occurrence rate in galactic nuclei than if the collision target is a
comparatively short-lived TDE disk. This is relevant if rare conditions, such
as a stellar EMRI with a mass/radius larger than the Sun (and correspondingly
larger cross sections for interaction with the disk), or low inclination
angles $\chi\ll 1$, are needed to produce detectable eruptions. However, such
long-lived self-sustained disks obviously cannot be produced the disk-
colliding body is a compact object, such as an intermediate-mass black holes
(IMBH), instead of a nondegenerate star. While IMBH-disk collisions can
generate QPE-like eruptions with similar properties to star-disk collisions
(LM23), they were already disfavored as QPE sources because of their
comparative rarity to stars and constraints based on the implied SMBH growth
rate via IMBH cannibalism (LM23). The inability of black holes to create self-
sustained disks with which to interact, would only further reduce their
appearance in the QPE population.

Although the sample size is still small, QPE sources appear to exhibit a
preference for occurring in low-mass galaxies (e.g., Wevers et al. 2022)
hosting correspondingly low-mass black holes, compared even to TDE host
galaxies. This observational finding does not have a clear dynamical
explanation in stellar-EMRI QPE models insofar that the TDE and EMRI rate
should scale similarly with SMBH mass (however, see Lu & Quataert 2023; King
2023; LM23). However, insofar that the presence of a gaseous disk is a
necessary condition for QPE activity, our disk evolution models may offer new
insights into QPE demographics. Firstly, the predicted EMRI destruction time
scales inversely with the SMBH mass (e.g., $\tau_{\rm dest}\propto
M_{\bullet,6}^{-2/3}$ for $\beta=0$; Eq. (32)) such that QPE sources resulting
from disk interaction will be longer-lived and hence more abundant at any
given time in galactic nuclei which contain lower-mass SMBH. However, this
effect may be in part compensated in a flux-limited sample if the eruptions
are indeed more luminous from stars orbiting higher SMBH masses (see Eq.
(51)). A second reason that QPE sources might preferentially occur around low-
mass SMBH relates to the requirement of a short precession time to maintain a
large disk-star misalignment angle and hence a long-lived stellar-collision-
fed disk (Fig. 1). This condition is more challenging to satisfy for larger
black hole masses because $\tau_{\rm prec}\propto M_{\bullet}$.

##  5 Summary/Conclusions

Motivated by questions about the origins and lifetimes of the black hole
accretion flows in QPE systems, we have developed a model for the coupled
evolution of a stellar EMRI on an inclined quasi-circular orbit and the
gaseous disk (modeled as a single annulus centered near their interaction
radius). The model accounts for mass stripped from the star and thermal energy
injected into the disk as the result of periodic disk-star collisions.

Our main conclusions can be summarized as follows:

  * •

While a transient gaseous disk (e.g., from an independent TDE) is necessary to
initiate the process, stellar stripping soon becomes the dominant source of
mass feeding the disk and the subsequent evolution “forgets” about the initial
state. Qualitative features of the disk’s evolution$-$particularly the
presence or absence of limit-cycles$-$depend on the uncertain efficiency with
which the mass stripped from the star heats the disk, as encapsulated within
our model by the dimensionless parameter $\beta\simeq\dot{E}_{\rm
abl}/|\dot{M}_{\star}|v_{\rm k}^{2}$. Future radiation hydrodynamic
simulations of disk-star collisions are necessary to estimate the value of
$\beta$ as a function of the system parameters such as the disk/star
inclination angle.

  * •

For inefficient collision/ablation heating ($\beta<3/16$) under radiation-
dominated conditions, the disk is susceptible to the standard thermal
instability. The instability manifests as limit-cycle behavior, in which the
disk starts at low surface density $\Sigma\sim\Sigma_{\rm min}$ in the gas-
pressure dominated radiatively-cooled “low” state, where it builds up mass
gradually from stellar stripping. Eventually, once the surface density reaches
a critical value $\Sigma\sim\Sigma_{\rm max}$, the disk transitions on the
thermal time to an advection-dominated “high” state (Figs. 2,3). In the
geometrically-thick high state, the disk rapidly “flushes” until the high
branches vanishes at $\Sigma<\Sigma_{\rm min}$ and the cycle starts again
after the disk quickly cools back to the low branch. The total duration of the
cycle is of order years to decades (Eq. (29), (A6)).

  * •

As $\beta$ increases, collisional/ablation heating becomes relevant in
offsetting radiative cooling in the low state and the boundaries of the limit
cycles shrink (Fig. 4; Appendix A). For sufficiently large $\beta$ above a
critical threshold $\beta_{\rm stable}=3/16$ (Appendix B), collision/ablation
heating is sufficient to completely stabilize the disk’s thermal evolution,
even under radiation-dominated conditions, resulting in a unique surface
density, scale-height and mass accretion rate (Sec. 3.2).

  * •

Regardless of the value of $\beta\lesssim 1$, the time-averaged mass-accretion
rates are characteristically $\langle\dot{M}_{\rm acc}\rangle\sim
0.03-0.3\dot{M}_{\rm Edd}$ for QPE-like systems (Eqs. (31), (36)). These
rates, and the characteristic timescales over which they are predicted to vary
due to changes in the disk properties at the collision radius ($\tau_{\rm
cyc}$ or $\tau_{\nu}$), broadly agree with the observed quiescent emission of
QPE sources. In particular, it is natural to expect both rising and decaying
quiescent X-ray light curves over timescales of years to decades, e.g. without
the need to invoke multiple full- or partial-TDEs in the same galactic
nucleus.

  * •

In disk-star collision models for QPEs, time-dependence in the disk properties
at the collision radius should naturally translate into secular timescale
changes in the eruption properties, such as their luminosity, duration, and
temperature (Fig. 3). For thermally-unstable disks, modulation of the eruption
temperatures (relative to the comparatively stable quiescent disk temperature)
due to large fluctuations in the disk midplane density (which controls photon
production in the collision ejecta), could cause QPE sources to appear or
disappear from detectability on timescales of years to decades. A reduction in
the QPE temperature due to the higher disk density at the peak of the cycle
may explain the temporary disappearance of detectable eruptions from GSN069
during its quiescent 2019 outburst (Miniutti et al., 2019, 2023b), before
their recent reappearance (Miniutti et al., 2023a).

  * •

Table 2 and Fig. 6 summarize the key timescales of the problem. For the
observed range of QPE periods $P_{\rm QPE}\sim 2-20$ hr, the star is
completely ablated by disk collisions on a timescale of $\tau_{\rm dest}\sim
10^{2}-10^{3}$ yr, depending weakly on the free parameters of the problem
(Eqs. (A8),(39)). However, the stellar lifetime may depend more sensitively on
our simplified $\alpha-$viscosity disk model, which may over-estimate the
midplane density (and hence under-estimate $\tau_{\rm dest}$) compared to
radiation GRMHD simulations (e.g., Jiang et al. 2019). The destruction time
also depends on our assumed prescription for stellar mass ablation (Eq. (6)),
which though supported by numerical simulations performed in the regime $h\sim
R_{\star}$, requires further elucidation through targeted studies (e.g., Yao &
Quataert, in prep), including both the $h\ll R_{\star}$ and $h\gg R_{\star}$
limits. Across most of the parameter space, the star is destroyed before its
orbit will experience significant migration due to gas drag or gravitational
interaction with the disk.

  * •

Despite their comparatively rare formation rate (e.g., relative to TDEs), the
longer durations stellar collision-fed disks spend accreting at a substantial
fraction of the Eddington accretion rate may have a number of consequences.
They should produce radially-extended narrow-line emission regions, consistent
with those observed from QPE hosts (Wevers et al., 2022; Patra et al., 2023),
in $\sim 0.1\%$ of low-mass galaxies (for an assumed EMRI rate of $10^{-6}$
yr-1 per galaxy). Longer QPE source lifetimes alleviate constraints on the
fraction of stellar EMRIs that need to produce detectable eruptions in order
to reproduce the observed QPE source population.

  * •

Our model assumes that the accretion disk remains on averaged aligned with the
black hole spin axis, which is justified provided that the stellar orbit
precesses sufficiently rapidly compared to the timescale over which the disk
surface density evolves (Fig. 1). While this condition is trivially satisfied
for thermally-stable disks with long viscous timescales
$\tau_{\nu}\gg\tau_{\rm prec}$, for the limit-cycle case the more stringent
condition $\tau_{\rm prec}\lesssim\tau_{\rm cyc}/10$ requires a combination of
a low-mass and/or rapidly-spinning SMBH. Even temporary changes in the star-
disk inclination angle resulting from the angular momentum of the deposited
stripped mass (e.g., Xiang-Gruess et al. 2016) could result in short-lived
changes to the QPE eruption recurrence time pattern, similar that observed in
GSN 069 following their recent reappearance.

We thank Riccardo Arcodia, Margherita Giustini, Giovanni Miniutti, and Eliot
Quataert for comments on the manuscript. IL acknowledges support from a
Rothschild Fellowship and The Gruber Foundation. BDM was supported, in part,
by the National Science Foundation (grant No. AST-2009255). The Flatiron
Institute is supported by the Simons Foundation.

##  Appendix A Analytic Estimates of Thermally Unstable Disk Evolution

As in the $\beta=0$ case (Sec. 3.1.1), the results for finite
collision/ablation heating ($0\ll\beta<\beta_{\rm stable}$) can be understood
analytically. Balancing now radiative cooling with both accretion and
collision/ablation heating ($\dot{E}_{\rm abl}+\dot{E}_{\rm
acc}\approx\dot{E}_{\rm rad}$) gives the following expression for the disk
surface density:

$\Sigma_{\dot{E}_{\rm rad}=\dot{E}_{\rm abl}+\dot{E}_{\rm
acc}}\approx\tilde{A}T^{8/3}\left(1-\tilde{B}T^{4}\right)^{1/3}\,,$ (A1)

where

$\tilde{A}=\left(\frac{\pi}{18\beta\eta}\frac{ca^{2}}{\kappa}\left(\frac{R_{\star}}{r_{0}}\right)^{-4}\left(\frac{M_{\bullet}}{M_{\star}}\right)^{-1}\frac{r_{0}^{2}}{v_{\rm
k}^{5}}\right)^{1/3};\tilde{B}=\frac{3r_{0}\alpha a\kappa}{8v_{\rm k}c}.$ (A2)

Here we have assumed $p_{\rm rad}\gg p_{\rm gas}$ and
$\kappa\approx\kappa_{\rm es}$.

As the disk surface density and temperature $T\propto\Sigma^{3/8}$ rise during
each cycle, thermal stability is lost once the accretion heating term (first
term on the right hand side of Eq. (A1)) becomes comparable to
collision/ablation heating. This occurs at a maximum density

$\Sigma_{\rm
max}=\tilde{A}\left(\frac{4}{27\tilde{B}^{2}}\right)^{1/3}=\frac{4}{9}\left(\frac{2\pi}{3}\right)^{1/3}\frac{1}{\kappa_{\rm
es}}\left(\frac{c}{v_{\rm
k}}\right)\left(\frac{1}{\beta\eta\alpha^{2}}\left(\frac{r_{0}}{R_{\star}}\right)^{4}\left(\frac{M_{\star}}{M_{\bullet}}\right)\right)^{1/3}\approx\\\
1.9\times 10^{4}\,{\rm
g\,cm^{-2}}\;(\beta/0.1)^{-1/3}\eta_{-3}^{-1/3}\alpha_{-1}^{-2/3}\mathcal{M}_{\star}^{1/3}\mathcal{R}_{\star}^{-4/3}M_{\bullet,6}^{-2/9}\mathcal{P}_{\rm
QPE,4}^{11/9}\,,$ (A3)

as is attained at temperature $T=(2/(3\tilde{B}))^{1/4}$. The time spent at
around $\Sigma\simeq\Sigma_{\rm max}/2$ each phase, is then approximately

$\tau_{\Sigma}\approx\frac{\pi r_{0}^{2}\Sigma_{\rm
max}}{|\dot{M}_{\star}|}\approx\frac{1}{2}P_{\rm
QPE}\,\eta^{-1}\left(\frac{r_{0}}{R_{\star}}\right)^{4}\left(\frac{M_{\star}}{M_{\bullet}}\right)\left(\frac{h}{r}\right)\,,$
(A4)

with

$\frac{h}{r}=\frac{aT^{4}r_{0}}{3\Sigma v_{\rm
k}^{2}}=\frac{2}{9}\frac{ar_{0}}{\tilde{B}\Sigma v_{\rm
k}^{2}}=2\left(\frac{4}{9\pi}\right)^{1/3}\left(\frac{\alpha}{\beta\eta}\left(\frac{r_{0}}{R_{\star}}\right)^{4}\left(\frac{M_{\star}}{M_{\bullet}}\right)\right)^{-1/3}\approx\\\
8.8\times
10^{-3}\;(\beta/0.1)^{1/3}\eta_{-3}^{1/3}\alpha_{-1}^{-1/3}\mathcal{R}_{\star}^{4/3}\mathcal{M}_{\star}^{-1/3}\mathcal{P}_{\rm
QPE,4}^{-8/9}M_{\bullet,6}^{-1/9}\,.$ (A5)

The time spent at any surface density, $\Sigma<\Sigma_{\rm max}$,
$\tau_{\Sigma}\propto h/r\propto\Sigma^{1/2}$ is again dominated by large
$\Sigma$. The total cycle duration can thus be estimated as:

$\tau_{\rm cyc}\approx\tau_{\rm\Sigma_{\rm max}}\approx P_{\rm
QPE}\,\eta^{-2/3}\,\left(\frac{r_{0}}{R_{\star}}\right)^{8/3}\left(\frac{M_{\star}}{M_{\bullet}}\right)^{2/3}\left(\frac{\beta}{\alpha}\right)^{1/3}\left(\frac{4}{9\pi}\right)^{1/3}\approx\\\
3.4\,{\rm
yr}\;\eta_{-3}^{-2/3}\alpha_{-1}^{-1/3}(\beta/0.1)^{1/3}\mathcal{M}_{\star}^{2/3}\mathcal{R}_{\star}^{-8/3}\mathcal{P}_{\rm
QPE,4}^{25/9}M_{\bullet,6}^{2/9}\,.$ (A6)

Similar to Eqs. (31), (32), the mass-averaged accretion rate and stellar
destruction time are given by:

$\langle\dot{M}_{\rm acc}\rangle\approx\frac{\pi r_{0}^{2}\Sigma_{\rm
max}}{\tau_{\rm cyc}}\approx 0.07\,\dot{M}_{\rm
Edd}\;\eta_{-3}^{1/3}\alpha_{-1}^{-1/3}(\beta/0.1)^{-2/3}\mathcal{R}_{\star}^{4/3}\mathcal{M}_{\star}^{-1/3}\mathcal{P}_{\rm
QPE,4}^{-2/9}M_{\bullet,6}^{-7/9}$ (A7) $\tau_{\rm
dest}\simeq\frac{9}{16}\frac{M_{\star,0}}{\langle\dot{M}_{\rm
acc}\rangle}\approx 310\,{\rm
yr}\,\eta_{-3}^{-1/3}\alpha_{-1}^{1/3}(\beta/0.1)^{2/3}\mathcal{R}_{\star}^{-4/3}\mathcal{M}_{\star}^{4/3}\mathcal{P}_{\rm
QPE,4}^{2/9}M_{\bullet,6}^{-2/9},$ (A8)

where the factor $\approx 9/16$ accounts for the dependence $\tau_{\rm
dest}\propto\mathcal{M}_{\star}^{16/9}$ as the star loses mass adiabatically
($\mathcal{R}_{\star}\propto\mathcal{M}_{\star}^{-1/3}$; Eq. (8)). Eq. (A8)
depends weakly on most of the parameters, accounting also for the fact that
lower-main sequence stars roughly obey
$\mathcal{R}_{\star}\propto\mathcal{M}_{\star}$.

##  Appendix B Thermal Stability with Collision/Ablation Heating

Neglecting radial advection, the disk cools through radiation, at the rate
(Eq. (17))

$\dot{E}_{-}=\dot{E}_{\rm rad}\approx\frac{4\pi}{3}r_{0}^{2}\sigma_{\rm
SB}\frac{T^{4}}{\kappa\Sigma}=\frac{1}{4}L_{\rm
Edd}\left(\frac{h}{r}\right)\,,$ (B1)

where in the second equality we have assumed $p_{\rm rad}\gg p_{\rm gas}$ and
$\kappa\approx\kappa_{\rm es}$. The total heating rate, including viscous
accretion and collision/ablation heating, can be written (Eq. (17)):

$\displaystyle\dot{E}_{+}=\dot{E}_{\rm acc}+\dot{E}_{\rm
abl}=\left(\frac{3}{8}|\dot{M}_{\rm acc}|+\beta|\dot{M}_{\star}|\right)v_{\rm
k}^{2}\,.$ (B2)

We have neglected heating from the disk material directly shocked by the star
because $|\dot{M}_{\star}|\gg\dot{M}_{\rm coll}$ in our fiducial models (Figs.
2).

Holding $\Sigma$ constant (since it varies on the longer, viscous timescale),
the criterion for thermal stability can be written (e.g., Piran 1978):

$\left(\frac{\partial\log{\dot{E}_{+}}}{\partial
h}\right)<\left(\frac{\partial\log{\dot{E}_{-}}}{\partial h}\right)\,,$ (B3)

or equivalently,

$\frac{2\frac{3}{8}\dot{M}_{\rm
acc}-\beta\dot{M}_{\star}}{\frac{3}{8}\dot{M}_{\rm
acc}+\beta\dot{M}_{\star}}<1\,.$ (B4)

In steady-state, $\dot{M}_{\rm acc}\approx\dot{M_{\star}}$, and thus stability
requires $\beta>\beta_{\rm stable}=3/16$. This can be understood because (at
fixed $\Sigma$), while accretion heating is unstable $\dot{E}_{\rm acc}\propto
h^{2}$ (e.g., Lightman & Eardley 1974), collision/ablation heating
$\dot{E}_{\star}\propto h^{-1}$ is stabilizing and hence can dominate for
sufficiently large $\beta$.

## References

  * Abramowicz et al. (1988) Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646, doi: 10.1086/166683
  * Arcodia et al. (2024a) Arcodia, R., Liu, Z., Merloni, A., et al. 2024a, arXiv e-prints, arXiv:2401.17275, doi: 10.48550/arXiv.2401.17275
  * Arcodia et al. (2021) Arcodia, R., Merloni, A., Nandra, K., et al. 2021, Nature, 592, 704, doi: 10.1038/s41586-021-03394-6
  * Arcodia et al. (2022) Arcodia, R., Miniutti, G., Ponti, G., et al. 2022, A&A, 662, A49, doi: 10.1051/0004-6361/202243259
  * Arcodia et al. (2024b) Arcodia, R., Merloni, A., Buchner, J., et al. 2024b, arXiv e-prints, arXiv:2403.17059.  https://arxiv.org/abs/2403.17059
  * Armitage et al. (1996) Armitage, P. J., Zurek, W. H., & Davies, M. B. 1996, ApJ, 470, 237, doi: 10.1086/177864
  * Bahcall & Wolf (1976) Bahcall, J. N., & Wolf, R. A. 1976, ApJ, 209, 214, doi: 10.1086/154711
  * Bardeen & Petterson (1975) Bardeen, J. M., & Petterson, J. A. 1975, ApJ, 195, L65, doi: 10.1086/181711
  * Bonnerot et al. (2016) Bonnerot, C., Rossi, E. M., & Lodato, G. 2016, MNRAS, 458, 3324, doi: 10.1093/mnras/stw486
  * Chakraborty et al. (2021) Chakraborty, J., Kara, E., Masterson, M., et al. 2021, ApJ, 921, L40, doi: 10.3847/2041-8213/ac313b
  * Chakraborty et al. (2024) Chakraborty, J., Arcodia, R., Kara, E., et al. 2024, arXiv e-prints, arXiv:2402.08722, doi: 10.48550/arXiv.2402.08722
  * Chen et al. (1995) Chen, X., Abramowicz, M. A., Lasota, J.-P., Narayan, R., & Yi, I. 1995, ApJ, 443, L61, doi: 10.1086/187836
  * Dai & Blandford (2013) Dai, L., & Blandford, R. 2013, MNRAS, 434, 2948, doi: 10.1093/mnras/stt1209
  * Davis et al. (2010) Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52, doi: 10.1088/0004-637X/713/1/52
  * Dexter & Begelman (2019) Dexter, J., & Begelman, M. C. 2019, MNRAS, 483, L17, doi: 10.1093/mnrasl/sly213
  * Di Matteo et al. (2002) Di Matteo, T., Perna, R., & Narayan, R. 2002, ApJ, 579, 706, doi: 10.1086/342832
  * Done et al. (2007) Done, C., Gierliński, M., & Kubota, A. 2007, A&A Rev., 15, 1, doi: 10.1007/s00159-007-0006-1
  * Franchini et al. (2023) Franchini, A., Bonetti, M., Lupi, A., et al. 2023, A&A, 675, A100, doi: 10.1051/0004-6361/202346565
  * Gallego-Cano et al. (2018) Gallego-Cano, E., Schödel, R., Dong, H., et al. 2018, A&A, 609, A26, doi: 10.1051/0004-6361/201730451
  * Generozov & Perets (2023) Generozov, A., & Perets, H. B. 2023, MNRAS, 522, 1763, doi: 10.1093/mnras/stad1016
  * Gezari et al. (2006) Gezari, S., Martin, D. C., Milliard, B., et al. 2006, ApJ, 653, L25, doi: 10.1086/509918
  * Giustini et al. (2020) Giustini, M., Miniutti, G., & Saxton, R. D. 2020, A&A, 636, L2, doi: 10.1051/0004-6361/202037610
  * Hills (1975) Hills, J. G. 1975, Nature, 254, 295, doi: 10.1038/254295a0
  * Hirose et al. (2009) Hirose, S., Blaes, O., & Krolik, J. H. 2009, ApJ, 704, 781, doi: 10.1088/0004-637X/704/1/781
  * Hopkins et al. (2016) Hopkins, P. F., Torrey, P., Faucher-Giguère, C.-A., Quataert, E., & Murray, N. 2016, MNRAS, 458, 816, doi: 10.1093/mnras/stw289
  * Ingram et al. (2021) Ingram, A., Motta, S. E., Aigrain, S., & Karastergiou, A. 2021, MNRAS, 503, 1703, doi: 10.1093/mnras/stab609
  * Ivanov et al. (1998) Ivanov, P. B., Igumenshchev, I. V., & Novikov, I. D. 1998, ApJ, 507, 131, doi: 10.1086/306324
  * Jacquemin-Ide et al. (2023) Jacquemin-Ide, J., Rincon, F., Tchekhovskoy, A., & Liska, M. 2023, arXiv e-prints, arXiv:2311.00034, doi: 10.48550/arXiv.2311.00034
  * Jiang et al. (2019) Jiang, Y.-F., Blaes, O., Stone, J. M., & Davis, S. W. 2019, ApJ, 885, 144, doi: 10.3847/1538-4357/ab4a00
  * 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
  * Kaur et al. (2023) Kaur, K., Stone, N. C., & Gilbaum, S. 2023, MNRAS, 524, 1269, doi: 10.1093/mnras/stad1894
  * King (2020) King, A. 2020, MNRAS, 493, L120, doi: 10.1093/mnrasl/slaa020
  * King (2022) —. 2022, MNRAS, 515, 4344, doi: 10.1093/mnras/stac1641
  * King (2023) —. 2023, MNRAS, 520, L63, doi: 10.1093/mnrasl/slad006
  * King et al. (2007) King, A. R., Pringle, J. E., & Livio, M. 2007, MNRAS, 376, 1740, doi: 10.1111/j.1365-2966.2007.11556.x
  * Kormendy & Richstone (1995) Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581, doi: 10.1146/annurev.aa.33.090195.003053
  * Krolik & Linial (2022) Krolik, J. H., & Linial, I. 2022, ApJ, 941, 24, doi: 10.3847/1538-4357/ac9eb6
  * Lightman & Eardley (1974) Lightman, A. P., & Eardley, D. M. 1974, ApJ, 187, L1, doi: 10.1086/181377
  * Linial & Metzger (2023) Linial, I., & Metzger, B. D. 2023, arXiv e-prints, arXiv:2303.16231, doi: 10.48550/arXiv.2303.16231
  * Linial & Metzger (2024) —. 2024, ApJ, 963, L1, doi: 10.3847/2041-8213/ad2464
  * Linial & Quataert (2024) Linial, I., & Quataert, E. 2024, MNRAS, 527, 4317, doi: 10.1093/mnras/stad3470
  * Linial & Sari (2017) Linial, I., & Sari, R. 2017, MNRAS, 469, 2441, doi: 10.1093/mnras/stx1041
  * Linial & Sari (2023) —. 2023, ApJ, 945, 86, doi: 10.3847/1538-4357/acbd3d
  * Liu et al. (2015) Liu, Z.-W., Tauris, T. M., Röpke, F. K., et al. 2015, A&A, 584, A11, doi: 10.1051/0004-6361/201526757
  * Lu & Quataert (2023) Lu, W., & Quataert, E. 2023, MNRAS, 524, 6247, doi: 10.1093/mnras/stad2203
  * MacLeod et al. (2012) MacLeod, M., Guillochon, J., & Ramirez-Ruiz, E. 2012, ApJ, 757, 134, doi: 10.1088/0004-637X/757/2/134
  * MacLeod & Lin (2020) MacLeod, M., & Lin, D. N. C. 2020, ApJ, 889, 94, doi: 10.3847/1538-4357/ab64db
  * MacLeod et al. (2013) MacLeod, M., Ramirez-Ruiz, E., Grady, S., & Guillochon, J. 2013, ApJ, 777, 133, doi: 10.1088/0004-637X/777/2/133
  * Metzger et al. (2012) Metzger, B. D., Rafikov, R. R., & Bochkarev, K. V. 2012, MNRAS, 423, 505, doi: 10.1111/j.1365-2966.2012.20895.x
  * Metzger & Stone (2017) Metzger, B. D., & Stone, N. C. 2017, ApJ, 844, 75, doi: 10.3847/1538-4357/aa7a16
  * Metzger et al. (2022) Metzger, B. D., Stone, N. C., & Gilbaum, S. 2022, ApJ, 926, 101, doi: 10.3847/1538-4357/ac3ee1
  * Milosavljević et al. (2006) Milosavljević, M., Merritt, D., & Ho, L. C. 2006, ApJ, 652, 120, doi: 10.1086/508134
  * Miniutti et al. (2023a) Miniutti, G., Giustini, M., Arcodia, R., et al. 2023a, A&A, 674, L1, doi: 10.1051/0004-6361/202346653
  * Miniutti et al. (2023b) —. 2023b, A&A, 670, A93, doi: 10.1051/0004-6361/202244512
  * Miniutti et al. (2019) Miniutti, G., Saxton, R. D., Giustini, M., et al. 2019, Nature, 573, 381, doi: 10.1038/s41586-019-1556-x
  * Nelson & Papaloizou (2000) Nelson, R. P., & Papaloizou, J. C. B. 2000, MNRAS, 315, 570, doi: 10.1046/j.1365-8711.2000.03478.x
  * Pan et al. (2022) Pan, X., Li, S.-L., Cao, X., Miniutti, G., & Gu, M. 2022, ApJ, 928, L18, doi: 10.3847/2041-8213/ac5faf
  * Pasham et al. (2024) Pasham, D. R., Coughlin, E. R., Zajacek, M., et al. 2024, arXiv e-prints, arXiv:2402.09690.  https://arxiv.org/abs/2402.09690
  * Patra et al. (2023) Patra, K. C., Lu, W., Ma, Y., et al. 2023, arXiv e-prints, arXiv:2310.05574, doi: 10.48550/arXiv.2310.05574
  * Payne et al. (2021) Payne, A. V., Shappee, B. J., Hinkle, J. T., et al. 2021, ApJ, 910, 125, doi: 10.3847/1538-4357/abe38d
  * Phinney (1989) Phinney, E. S. 1989, in The Center of the Galaxy, ed. M. Morris, Vol. 136, 543
  * Piran (1978) Piran, T. 1978, ApJ, 221, 652, doi: 10.1086/156069
  * Pringle (1981) Pringle, J. E. 1981, ARA&A, 19, 137, doi: 10.1146/annurev.aa.19.090181.001033
  * Pringle & Rees (1972) Pringle, J. E., & Rees, M. J. 1972, A&A, 21, 1
  * Quintin et al. (2023) Quintin, E., Webb, N. A., Guillot, S., et al. 2023, A&A, 675, A152, doi: 10.1051/0004-6361/202346440
  * Rees (1988) Rees, M. J. 1988, Nature, 333, 523, doi: 10.1038/333523a0
  * Salpeter (1964) Salpeter, E. E. 1964, ApJ, 140, 796, doi: 10.1086/147973
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shen & Matzner (2014) Shen, R.-F., & Matzner, C. D. 2014, ApJ, 784, 87, doi: 10.1088/0004-637X/784/2/87
  * Sheng et al. (2021) Sheng, Z., Wang, T., Ferland, G., et al. 2021, ApJ, 920, L25, doi: 10.3847/2041-8213/ac2251
  * Shu et al. (2018) Shu, X. W., Wang, S. S., Dou, L. M., et al. 2018, ApJ, 857, L16, doi: 10.3847/2041-8213/aaba17
  * Soltan (1982) Soltan, A. 1982, MNRAS, 200, 115, doi: 10.1093/mnras/200.1.115
  * Sorathia et al. (2013) Sorathia, K. A., Krolik, J. H., & Hawley, J. F. 2013, ApJ, 777, 21, doi: 10.1088/0004-637X/777/1/21
  * Stone et al. (2013) Stone, N., Sari, R., & Loeb, A. 2013, MNRAS, 435, 1809, doi: 10.1093/mnras/stt1270
  * Stone et al. (2017) Stone, N. C., Küpper, A. H. W., & Ostriker, J. P. 2017, MNRAS, 467, 4180, doi: 10.1093/mnras/stx097
  * Stone & Metzger (2016) Stone, N. C., & Metzger, B. D. 2016, MNRAS, 455, 859, doi: 10.1093/mnras/stv2281
  * Suková et al. (2021) Suková, P., Zajaček, M., Witzany, V., & Karas, V. 2021, ApJ, 917, 43, doi: 10.3847/1538-4357/ac05c6
  * Tagawa & Haiman (2023) Tagawa, H., & Haiman, Z. 2023, MNRAS, 526, 69, doi: 10.1093/mnras/stad2616
  * van Velzen et al. (2019) van Velzen, S., Stone, N. C., Metzger, B. D., et al. 2019, ApJ, 878, 82, doi: 10.3847/1538-4357/ab1844
  * Webbe & Young (2023) Webbe, R., & Young, A. J. 2023, MNRAS, 518, 3428, doi: 10.1093/mnras/stac3318
  * Wevers et al. (2022) Wevers, T., Pasham, D. R., Jalan, P., Rakshit, S., & Arcodia, R. 2022, A&A, 659, L2, doi: 10.1051/0004-6361/202243143
  * Wevers et al. (2023) Wevers, T., Coughlin, E. R., Pasham, D. R., et al. 2023, ApJ, 942, L33, doi: 10.3847/2041-8213/ac9f36
  * Xian et al. (2021) Xian, J., Zhang, F., Dou, L., He, J., & Shu, X. 2021, ApJ, 921, L32, doi: 10.3847/2041-8213/ac31aa
  * Xiang-Gruess et al. (2016) Xiang-Gruess, M., Ivanov, P. B., & Papaloizou, J. C. B. 2016, MNRAS, 463, 2242, doi: 10.1093/mnras/stw2130
  * Zalamea et al. (2010) Zalamea, I., Menou, K., & Beloborodov, A. M. 2010, MNRAS, 409, L25, doi: 10.1111/j.1745-3933.2010.00930.x
  * Zanazzi & Lai (2019) Zanazzi, J. J., & Lai, D. 2019, MNRAS, 487, 4965, doi: 10.1093/mnras/stz1610
  * Zhao et al. (2022) Zhao, Z. Y., Wang, Y. Y., Zou, Y. C., Wang, F. Y., & Dai, Z. G. 2022, A&A, 661, A55, doi: 10.1051/0004-6361/202142519
  * Zhou et al. (2024) Zhou, C., Huang, L., Guo, K., Li, Y.-P., & Pan, Z. 2024, arXiv e-prints, arXiv:2401.11190, doi: 10.48550/arXiv.2401.11190