Interacting Stellar EMRIs as Sources of Quasi-periodic Eruptions in Galactic Nuclei

The first EMRI is assumed to be of a star (or brown dwarf or planet) of mass M₁ = m₁ M_⊙ and radius R₁ = r₁ R_⊙. If M₁ is overflowing its Roche lobe onto the SMBH of mass M_• = 10⁶ M_•,6 M_⊙, then its semimajor axis is given by

$$\begin{eqnarray}\begin{array}{rcl}{a}_{1} & = & {R}_{\mathrm{RL}}\simeq 2.17{R}_{1}{\left(\displaystyle \frac{{M}_{\bullet }}{{M}_{1}}\right)}^{1/3}\\ & \simeq & 1.0\,\mathrm{au}\,{M}_{\bullet ,6}^{1/3}\displaystyle \frac{{r}_{1}}{{m}_{1}^{1/3}}\approx 1.0\mathrm{au}\,\displaystyle \frac{{M}_{\bullet ,6}^{1/3}}{{\tilde{\rho }}_{1}^{1/3}},\end{array}\end{eqnarray} \tag{ 1 }$$

where ${\tilde{\rho }}_{1}\equiv {\rho }_{1}/{\rho }_{\odot }$ is the mean density ρ₁ of M₁ normalized to the solar value ρ_⊙. The semimajor axis must also exceed that of the innermost stable circular orbit (ISCO),

$$\begin{eqnarray}&&{R}_{\mathrm{ISCO}}\approx 0.06\,\mathrm{au}({R}_{\mathrm{ISCO}}/6{R}_{{\rm{g}}}){M}_{\bullet ,6},\end{eqnarray} \tag{ 2 }$$

where the gravitational radius R_g ≡ GM_•/c² and R_ISCO/R_g varies from 1 to 9 as the dimensionless SMBH spin a_• varies from +1 to −1. We see that R_RL > R_ISCO for SMBH masses in the range estimated from QPE host galaxies (M_• ≲ 10⁶ M_⊙), for all physically allowed values of M₁, R₁ corresponding to brown dwarfs or nondegenerate stars.

The first EMRI M₁ undergoes RLOF evolution on the timescale set by gravitational wave radiation (e.g., MS17),

$$\begin{eqnarray}&&{\tau }_{\mathrm{GW}}\approx 1.3\times {10}^{6}\,\mathrm{yr}\,\chi {M}_{\bullet ,6}^{-2/3}{m}_{1}^{-1}{\tilde{\rho }}_{1}^{-4/3}.\end{eqnarray} \tag{ 3 }$$

The dimensionless factor χ = 1 in the case of free inspiral and χ = 3/(3p − 1) ≈ 2–4 if M₁ is undergoing RLOF, where ${R}_{\star }\propto {M}_{\star }^{p}$ and the given range of χ corresponds to p ≈ 0.6–0.8 for a range of stellar masses and thermal states (e.g., Linial & Sari 2017; MS17). Assuming stable mass transfer, this results in a mass-accretion rate $\langle \dot{M}\rangle \sim {M}_{1}/{\tau }_{\mathrm{GW}}$ and corresponding accretion luminosity,

$$\begin{eqnarray}&&\langle L\rangle \approx 0.1\langle \dot{M}\rangle {c}^{2}\approx 4\times {10}^{39}\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\chi }^{-1}{M}_{\bullet ,6}^{2/3}{m}_{1}^{2}{\tilde{\rho }}_{1}^{4/3}.\end{eqnarray} \tag{ 4 }$$

For typical values ${m}_{1}\sim {\tilde{\rho }}_{1}\sim 1$, this is several orders of magnitude too small to explain time-averaged QPE luminosities, demonstrating why single EMRI models are challenged. Single EMRI models that invoke denser stars ${\tilde{\rho }}_{1}\gg 1$, like white dwarfs (Zalamea et al. 2010; King 2020) or helium cores (Zhao et al. 2021), can produce higher 〈L〉, but these run into their own challenges with respect to rates (Section 5.2).

The second EMRI is a star of mass M₂ = m₂ M_⊙ and radius R₂ = r₂ R_⊙ on an orbit of semimajor axis a₂. We assume that both EMRIs have nearly circularized their orbits due to energy loss via gravitational wave (GW) emission. A strong interaction between two consecutive EMRIs will only occur if their orbits approach one another because the rate of gravitational-wave-driven orbital decay of M₂ is faster than that of M₁. We thus require M₂ ≳ M₁ for an interaction. If both M₁ and M₂ are filling their Roche radii, then they must possess roughly equal mean densities due to their common semimajor axes near the point of strongest interaction, i.e., ${M}_{1}/{R}_{1}^{3}\approx {M}_{2}/{R}_{2}^{3}$, and hence we also require R₂ ≳ R₁.

Once the orbits of the two EMRIs approach within a separation Δa ≡ a₂ − a₁ of several stellar radii, strong tidal interactions occur between them. At this point, the EMRIs share a roughly common semimajor axis a₁ ≃ a₂ = a and orbital period,

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{orb}} & \simeq & 2\pi {\left(\displaystyle \frac{{a}^{3}}{{{GM}}_{\bullet }}\right)}^{1/2}\approx 8.8\,\mathrm{hr}\,{\tilde{\rho }}^{-1/2}{\left(\displaystyle \frac{a}{{R}_{\mathrm{RL}}}\right)}^{3/2},\end{array}\end{eqnarray} \tag{ 5 }$$

where $\tilde{\rho }$ is the mean density of either star with respective Roche radius R_RL (Equation (1)).

To simplify the analysis below, due to the slower evolution of M₁, we approximate its orbit as being fixed during its interaction with M₂. The number of orbits required for M₂ to migrate inward radially, via gravitational wave emission, by a distance δ a ≪ a₂ is given by

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{GW}} & \sim & \displaystyle \frac{{\tau }_{\mathrm{GW}}}{{T}_{\mathrm{orb}}}\left(\displaystyle \frac{\delta a}{{a}_{2}}\right)\\ & \approx & 6\times {10}^{6}\displaystyle \frac{\chi }{{M}_{\bullet ,6}}\displaystyle \frac{{r}_{2}}{{m}_{2}{\tilde{\rho }}_{2}^{1/2}}{\left(\displaystyle \frac{a}{{R}_{\mathrm{RL}}}\right)}^{3/2},\end{array}\end{eqnarray} \tag{ 6 }$$

where for τ_GW we use Equation (3), replacing M₁ with M₂.

Each close flyby will result in the removal of mass from one or both stars and an accretion-powered flare, such that the QPE recurrence time T_QPE is roughly the time between flybys, T_fly (however, see Section 3.5). The specific mechanism of the mass removal is described below. There are two cases to consider, depending on whether both EMRIs are orbiting in the same direction ("co-orbiting" case) or in opposite directions ("counter-orbiting" case).

In the counter-orbiting case, close passages occur twice per orbital period (Equation (5)),

$$\begin{eqnarray}&&{T}_{\mathrm{fly}}\approx \displaystyle \frac{{T}_{\mathrm{orb}}}{2}\approx 4.4\,\mathrm{hr}\,{\tilde{\rho }}^{-1/2},\mathrm{Counter}-\mathrm{orbiting}.\end{eqnarray} \tag{ 7 }$$

where $\tilde{\rho }$ is the mean density of the star or stars undergoing RLOF (in our fiducial scenario, at least M₂).

Figure 2 shows the mean density of stars in different evolutionary stages (WDs, brown dwarfs, and stars at different phases of the main sequence) as a function of their mass, compared to the minimum density compatible with the observed values of T_QPE = T_fly according to Equation (7) assuming a = R_RL. For example, eRO-QPE1(eRO-QPE2) require ρ/ρ_⊙ ≳ 0.06(2.7) to match the observed eruption periods T_QPE ≈ 18.5(2.7) hr. eRO-QPE2 is consistent with a brown dwarf/planet or ZAMS stars of mass 3 × 10⁻³ ≲ M₁ ≲ 0.7M_⊙. The longer period of eRO-QPE1 is not compatible with a ZAMS star undergoing RLOF, but is compatible with an evolved star of mass ≳2M_⊙.

Zoom In Zoom Out Reset image size

The ZAMS/HAMS/TAMS stellar properties shown in Figure 2 assume thermal equilibrium, which is not a good approximation when the stars are losing mass at a high rate (e.g., Linial & Sari 2017; see discussion at the end of Appendix A). Insofar that thermal timescale mass-loss will cause a star to inflate, the lines in Figure 2 represent an upper limit on the inferred density of the RLOFing star (lower limit on T_QPE) at a given stellar mass.

Next, consider the co-orbiting case. Here, the inner star M₁ must "chase" the outer one M₂ due to its slightly shorter orbital period ΔT_orb ≪ T_orb, as results from their small semimajor axis difference, Δa ≪ a₁, a₂. The greater number of orbits required for a close passage in this case, N_fly ≈ T_orb/∣ΔT_orb∣ ≈(2/3)(a/Δa) ≫ 1, results in a larger time interval between collisions,

$$\begin{eqnarray}&&{T}_{\mathrm{fly}}\approx {N}_{\mathrm{fly}}{T}_{\mathrm{orb}}\approx 10.8\,{\rm{d}}\,\displaystyle \frac{{M}_{\bullet ,6}^{1/3}}{{m}_{2}^{1/3}{\tilde{\rho }}^{1/2}}{\left(\displaystyle \frac{{\rm{\Delta }}a}{5{R}_{2}}\right)}^{-1},\mathrm{Co}-\mathrm{orbiting},\end{eqnarray} \tag{ 8 }$$

where we have assumed M₂ is overflowing its Roche lobe and Δa is normalized to a characteristic value ∼ 5R₂ necessary for a strong interaction (Equation (15) below). The large value of T_fly in the co-orbiting case is challenging to reconcile with the short observed QPE periods T_QPE ∼ hours, unless the colliding stars are WDs (r₂ ∼ 0.01; ρ ≳ 10⁵ g cm⁻³). For this reason, we favor the counter-orbiting case. However, our results to follow would apply equally to the co-orbiting case, and the latter may be relevant for longer-period AGN variability (Section 5.1).

During the interval between close encounters, T_fly, the orbital separation Δa decreases due to the gravitational wave inspiral of M₂ by an amount Δa_GW, obtained by setting N_GW =N_fly = 1/2 and N_GW = N_fly ≃ (2/3)(a/Δa) in the counter-orbiting and co-orbiting cases, respectively:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{\Delta }}{a}_{\mathrm{GW}}}{{R}_{2}} & \approx & 8\times {10}^{-8}\displaystyle \frac{{r}_{2}^{2}}{\chi }{M}_{\bullet ,6}{\tilde{\rho }}^{3/2},\,\mathrm{Counter}-\mathrm{orbiting},\\ & \approx & 5\times {10}^{-6}\displaystyle \frac{{r}_{2}}{\chi }{M}_{\bullet ,6}^{4/3}{\tilde{\rho }}^{7/6}{\left(\displaystyle \frac{{\rm{\Delta }}a}{5{R}_{2}}\right)}^{-1},\,\mathrm{Co}-\mathrm{orbiting},\end{array}\end{eqnarray} \tag{ 9 }$$

where we have used Equation (6) with δ a = Δa_GW and have assumed both stars are overflowing their Roche lobes. For stellar parameters {ρ ∼ 0.1−10 g cm⁻³, r₂ ∼ 0.1−1} and {ρ ≳ 10⁵ g cm⁻³, r₂ ∼ 0.01} necessary to match the QPE timescales in the counter-orbiting and co-orbiting cases, respectively, we have Δa_GW/R₂ ≪ 1. The two stars will thus be subject to many strong flybys prior to any direct contact between their surfaces.

A close passage between M₁ and M₂ can generate mass loss from one or both stars exceeding their rate of steady mass transfer onto the SMBH. Mass loss can in principle arise either from a direct physical collision between the stars ("hydrodynamical" mass loss), or as the result of tidal forces impacting the rate of mass-transfer onto the SMBH ("tidal" mass-loss). In both cases, the more compact lower-mass star M₁ will preferentially remove mass from the more dilute outer layers of M₂. For this reason and others related to the geometry of the Roche surface (see below), the bulk of this discussion focuses on mass loss from M₂. Furthermore, we focus on tidal instead of hydrodynamical mass loss because: (1) as we show below, it becomes significant once Δa shrinks to a few stellar radii; and (2) many such close flybys occur before the first physical collision (Equation (9)). The latter point contrasts with the noncoplanar case, for which many more orbits separate the close encounters and physical collisions are more relevant (MS17).

As we show in Appendix A, the gravitational influence of a close passage from M₁ is to briefly shrink the Hill radius r_H of M₂, according to:

$$\begin{eqnarray}&&\displaystyle \frac{{r}_{{\rm{H}}}}{{r}_{{\rm{H}},0}}\equiv 1-\epsilon ;\,\,\epsilon \simeq \displaystyle \frac{{M}_{1}}{3{M}_{2}}{\left(\displaystyle \frac{{\rm{\Delta }}a}{{r}_{{\rm{H}},0}}\right)}^{-2},\end{eqnarray} \tag{ 10 }$$

where ${r}_{{\rm{H}},0}\simeq {({M}_{2}/3{M}_{\bullet })}^{1/3}$ is the usual (unperturbed) Hill radius (Equation (1)).

Insofar that r_H,0 ≃ R₂ if M₂ is filling its Roche lobe and losing mass through the inner Lagrange point L₁, then the close passage of M₁ causes the Roche surface of R₂ to penetrate below its photosphere by an additional factor Δr ≃ R₂. To the extent that Δr exceeds the atmosphere scale height H ∼ (10⁻⁴ − 10⁻³)R₂ near the photosphere of M₂, this increases its mass-loss rate through L₁ by a large factor for the brief time interval τ_fly ∼ (Δa/a)T_QPE the two stars spend close to each other.

In Appendix A, we estimate the mass-loss Δm_fly from M₂ per close passage, following the formalism of Ginzburg & Quataert (2021). We find (Equation (A11))

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{m}_{\mathrm{fly}} & \sim & 5\times {10}^{-8}{M}_{\odot }{\tilde{\kappa }}^{-1}\displaystyle \frac{{m}_{1}^{4.5}}{{r}_{2}^{3.5}}\\ & & \times \displaystyle \frac{{m}_{2}^{1/3}}{{M}_{\bullet ,6}^{1/3}}\left(\displaystyle \frac{{T}_{\mathrm{QPE}}}{10\,\mathrm{hr}}\right){\left(\displaystyle \frac{{T}_{\mathrm{eff}}}{{10}^{4}{\rm{K}}}\right)}^{-4}{\left(\displaystyle \frac{{\rm{\Delta }}a}{5{R}_{2}}\right)}^{-8},\end{array}\end{eqnarray} \tag{ 11 }$$

where $\tilde{\kappa }$ is the photosphere opacity (normalized to the electron scattering opacity) and T_eff ∼ 10⁴ K the surface temperature, where we have assumed an n = 3 polytrope for the outer envelope structure. These surface properties are expected due to the strong influence of irradiation of the star by the luminous SMBH accretion flow, which usually overwhelms its internal nuclear luminosity (Appendix A).

The total mass loss from M₂ during the time the two stars spend separated by any distance ∼ Δa is given by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{\Delta }}m}{{M}_{2}} & \sim & {N}_{\mathrm{GW}}({\rm{\Delta }}a)\displaystyle \frac{{\rm{\Delta }}{m}_{\mathrm{fly}}}{{M}_{2}}\\ & \sim & 1.2{M}_{\bullet ,6}^{-4/3}\displaystyle \frac{\chi }{\tilde{\kappa }}\displaystyle \frac{{m}_{1}^{4.5}}{{m}_{2}^{2}{r}_{2}^{1.5}}\left(\displaystyle \frac{{T}_{\mathrm{QPE}}}{10\,\mathrm{hr}}\right)\\ & & \times {\left(\displaystyle \frac{{T}_{\mathrm{eff}}}{{10}^{4}{\rm{K}}}\right)}^{-4}{\left(\displaystyle \frac{{\rm{\Delta }}a}{5{R}_{2}}\right)}^{-7},\end{array}\end{eqnarray} \tag{ 12 }$$

where we have used Equation (6) for N_GW with a = R_RL and δ a = Δa.

We thus see that M₂ will be completely destroyed (Δm ≳ M₂) once gravitational wave radiation reduces the orbital separation Δa below a critical value

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{\Delta }}{a}_{\mathrm{dest}}}{{R}_{2}} & \approx & 5{M}_{\bullet ,6}^{-4/21}\displaystyle \frac{{\chi }^{1/7}}{{\tilde{\kappa }}^{1/7}}\displaystyle \frac{{m}_{1}^{9/14}}{{r}_{2}^{3/28}{m}_{2}^{2/7}}\\ & & \times {\left(\displaystyle \frac{{T}_{\mathrm{QPE}}}{10\,\mathrm{hr}}\right)}^{1/7}{\left(\displaystyle \frac{{T}_{\mathrm{eff}}}{{10}^{4}\,{\rm{K}}}\right)}^{-4/7},\end{array}\end{eqnarray} \tag{ 13 }$$

which we note is a weak function of the relevant parameters. The destruction of M₂ will occur gradually, over a timescale

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{dest}} & \sim & {N}_{\mathrm{GW}}({\rm{\Delta }}{a}_{\mathrm{dest}}){T}_{\mathrm{orb}}\\ & \approx & 3\times {10}^{4}\,\mathrm{yr}\displaystyle \frac{\chi }{{M}_{\bullet ,6}}\displaystyle \frac{{r}_{2}^{7/2}}{{m}_{2}^{11/6}}\left(\displaystyle \frac{{\rm{\Delta }}{a}_{\mathrm{dest}}}{5{R}_{2}}\right).\end{array}\end{eqnarray} \tag{ 14 }$$

Most of the total mass of M₂ accreted by the SMBH will therefore occur when the per-flyby mass loss Δm_fly is near the critical value

$$\begin{eqnarray}&&\begin{array}{l}\,{\rm{\Delta }}{m}_{\mathrm{dest}}\,\sim \,{M}_{2}\displaystyle \frac{{T}_{\mathrm{QPE}}}{{t}_{\mathrm{dest}}}\\ \sim 4\times {10}^{-8}{M}_{\odot }\displaystyle \frac{{M}_{\bullet ,6}}{\chi }\displaystyle \frac{{m}_{2}^{17/6}}{{r}_{2}^{7/2}}\left(\displaystyle \frac{{T}_{\mathrm{QPE}}}{10\,\mathrm{hr}}\right){\left(\displaystyle \frac{{\rm{\Delta }}{a}_{\mathrm{dest}}}{5{R}_{2}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 15 }$$

We thus find that Δm_dest and t_dest are constrained to lie within a couple orders of magnitude of ∼10⁻⁸ M_⊙ and ∼10⁴ yr, respectively (for the allowed ranges of m₂ and M_•).

The above expressions refer to mass loss from the inner L₁ point of M₂ due to tidal interactions with M₁. In the case when both stars are undergoing RLOF, M₁ can also experience enhanced mass loss through its the outer L₂ Lagrange point, due to the gravitational force of M₂. However, because of the significant radial separation between the unperturbed L₁ and L₂ points ${\rm{\Delta }}{R}_{{\rm{L}}1,{\rm{L}}2}/{R}_{1}\approx {(2/3)({M}_{1}/{M}_{\bullet })}^{1/3}\sim {10}^{-2}$ (Linial & Sari 2017) relative to the photosphere scale height of M₂ (to which the mass-loss rate is extremely sensitive), mass loss from M₁ during the flyby will generally be smaller than that from M₂.

The estimates presented so far assume perfectly circular orbits. While our scenario invokes quasi-circular EMRI orbits, some residual eccentricity e ≪ 1 may be present (Equation (28)). This residual eccentricity can modulate the QPE peak luminosity, L_peak, in an observable way. Since L_peak ∝ Δm_fly ∝ Δa⁻⁸ ∝ (1 − e)⁻⁸ (Equation (11)), even a residual eccentricity of e ∼ 10⁻² (e ∼ 10⁻³) suffices to change L_peak by a factor of 10 (by a factor of 2), possibly contributing to the large observed variation in QPE amplitudes within a single source. See Appendix A for more details.

To zeroth order, the QPE period equals the interval between flybys, i.e., T_QPE = T_fly (Equations (7) and (8)). The gaseous disk generated by the stripped mass will accrete onto the SMBH, powering X-ray emission, nominally on the viscous time, τ_visc, at the circularization radius, r_circ. Associating the viscous time with the QPE flare duration,

$$\begin{eqnarray}\begin{array}{rcl} & & {\tau }_{\mathrm{QPE}}\sim {\tau }_{\mathrm{visc}}\sim {\left.\displaystyle \frac{{r}^{2}}{\nu }\right|}_{{r}_{\mathrm{circ}}}\approx \displaystyle \frac{1}{\alpha }\displaystyle \frac{{T}_{\mathrm{orb}}}{2\pi }{\left(\displaystyle \frac{h}{r}\right)}^{-2}\\ & \approx & \displaystyle \frac{3.2}{{\alpha }_{0.1}}{T}_{\mathrm{QPE}}{\left(\displaystyle \frac{h}{r}\right)}^{-2}{\left(\displaystyle \frac{{r}_{\mathrm{circ}}}{{R}_{\mathrm{RL}}}\right)}^{3/2},\mathrm{Counter}-\mathrm{orbiting}\\ & \approx & \displaystyle \frac{0.054}{{\alpha }_{0.1}}{T}_{\mathrm{QPE}}{\left(\displaystyle \frac{h}{r}\right)}^{-2}\left(\displaystyle \frac{{\rm{\Delta }}a}{5{R}_{2}}\right)\displaystyle \frac{{m}_{1}^{1/3}}{{M}_{\bullet ,6}^{1/3}}\\ & & \times {\left(\displaystyle \frac{{r}_{\mathrm{circ}}}{{R}_{\mathrm{RL}}}\right)}^{3/2},\mathrm{Co}-\mathrm{orbiting},\end{array}\end{eqnarray} \tag{ 16 }$$

where ν = α c_s h is the kinematic viscosity, h the vertical aspect ratio, c_s = hΩ_K the sound speed, ${{\rm{\Omega }}}_{K}={({{GM}}_{\bullet }/{r}^{3})}^{1/2}$, and α = 0.1α_0.1 the viscosity parameter.

The duty cycle τ_QPE/T_QPE ∼ 0.1−0.4 inferred from observations of QPEs (e.g., Table 1) is difficult to satisfy in the counter-orbiting case based on Equation (16) if r_circ ∼ R_RL and h/r ≪ 1. However, note that: (1) the accretion rate from an initially thin ring of material typically peaks at ∼1/10 of t_visc as measured at the ring radius (e.g., Pringle 1981); (2) disk material formed from the collision will be hot and may find itself in a slim-disk-like state (e.g., Abramowicz et al. 1988) with h/r ∼ 1, if the accretion luminosity is indeed approaching the Eddington value, L_Edd ∼ 10⁴⁴ M_•,6 erg s⁻¹, as may be achieved depending on the SMBH mass; (3) if both stars lose significant mass from the interaction, then due to the opposing specific angular momenta of the counter-orbiting stellar orbits, the disk that forms from the mixture of debris will circularize at radii r_circ < R_RL; (4) systems with T_QPE > τ_QPE would not exhibit strong X-ray periodicity and hence would be observationally selected against in QPE searches.

In this section we develop a simplified but illustrative toy model for QPE light-curve evolution. Following standard procedures (e.g., Pringle 1981; Metzger et al. 2008), we numerically solve the diffusion equation for the time-dependent evolution of disk surface density Σ(r) assuming an initially narrow, δ-function distribution of mass at the radius R₀ ≈ a (a "spreading ring solution"). We assume a kinematic viscosity law of the form ν ∝ r^1/2, appropriate for a disk with aspect ratio h/r ∼ 1 (as is required to reproduce observed τ_QPE; see above). We calculate the emission in various X-ray energy bands by assuming blackbody emission at the local equilibrium temperature (obtained by balancing viscous heating with radiative cooling; we neglect color corrections due to electron scattering) and integrating over disk radii, starting from the assumed inner boundary (taken to be the ISCO radius of a nonspinning black hole). We treat the initial viscous timescale at R₀, t_visc,0, and its dependence on the initial disk mass Δm as a free parameter of the problem. Standard α − disk models predict t_visc,0 ∝ 1/ν₀ ∝ Δm² and a viscosity law that deviates from our assumed ν ∝ r^1/2; however, the same disk models are known to be thermally and viscously unstable (e.g., Lightman & Eardley 1974) in the small radial range relevant for us. Given the lack of consensus in the literature on viscous stability of radiation-dominated disks, we adopt the idealized ν ∝ r^1/2 parameterization for simplicity.

Panel (a) of Figure 3 shows three different light curves assuming a 10⁶ M_⊙ SMBH and an initial viscous time t_{visc, 0} = 0.1 days. The three colors correspond to three different initial masses in the spreading ring (labeled in the figure). All light curves are for photons with h ν = 300 eV, i.e., soft X-rays emitted from the Wien tail of the disk's multicolor blackbody spectrum. Consequently, the light curves transition to an exponential decline in individual X-ray bands, even though the bolometric luminosity (and mass accretion rate) falls off as a power law. Panel (b) of Figure 3 presents the achromatic behavior fundamental to any quasi-thermal spreading disk model: for observations on the Wien tail, high-energy observing bands show steeper decays than (comparatively) low-energy bands. This behavior qualitatively reproduces the achromatic evolution seen in Figure 2 of Miniutti et al. (2019).

Zoom In Zoom Out Reset image size

We have seen that achieving a match between theoretical and observed values of T_QPE is only possible if both stars share a common orbital plane. However, even if this is true at one moment in time, it may not be true later, due to the effect of nodal precession from the SMBH spin.

If we assume that both stars are misaligned from the SMBH equatorial plane by an angle I, then the maximum distance between the two orbits is ${d}_{\max }\approx a\sin I\sin ({{\rm{\Omega }}}_{1}-{{\rm{\Omega }}}_{2}$), where Ω₁ and Ω₂ are the nodal angles of each orbit with respect to a reference direction in the SMBH equatorial plane, and the approximate equality here reflects the assumption that Ω₁ − Ω₂ ≪ 1 (i.e., the orbits are nearly coplanar). At leading post-Newtonian order, nodal precession is driven by Lense-Thirring frame dragging, with the nodal shift per orbit for a circular orbit given by

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\Omega }}=4\pi {\chi }_{\bullet }{\left(\displaystyle \frac{a}{{R}_{{\rm{g}}}}\right)}^{-3/2},\end{eqnarray} \tag{ 17 }$$

where 0 ≤ χ_• ≤ 1 is the dimensionless spin magnitude of the SMBH (Merritt et al. 2010). Differential nodal precession will cause initially coplanar orbits to precess into a 3D configuration, so long as their semimajor axes a₁ and a₂ = a₁ + Δa differ slightly. After a time t, two initially co-aligned orbits will achieve a nodal separation

$$\begin{eqnarray}\begin{array}{rrl}\,{{\rm{\Omega }}}_{1}-{{\rm{\Omega }}}_{2} & = & \left(\displaystyle \frac{{\rm{\Delta }}{{\rm{\Omega }}}_{1}}{{T}_{\mathrm{orb},1}}-\displaystyle \frac{{\rm{\Delta }}{{\rm{\Omega }}}_{2}}{{T}_{\mathrm{orb},2}}\right)t\\ \, & \mathop{\approx }\limits_{{\rm{\Delta }}a\ll 1}\,6{\chi }_{\bullet }\displaystyle \frac{{G}^{2}{M}_{\bullet }^{2}t}{{c}^{3}{a}^{3}}\displaystyle \frac{{\rm{\Delta }}a}{a},\qquad \end{array}\end{eqnarray} \tag{ 18 }$$

where T_orb,1 and T_orb,2 are the orbital periods of M₁ and M₂, respectively.

The assumption of coplanarity will break down (and QPEs will turn off, for a time) once ${d}_{\max }\gtrsim {\rm{\Delta }}{a}_{\mathrm{dest}}\sim 5{R}_{2}$ (Equation (13)). In the small precession limit, ${{\rm{\Omega }}}_{1}-{{\rm{\Omega }}}_{2}\approx \sin ({{\rm{\Omega }}}_{1}-{{\rm{\Omega }}}_{2})$, initially coplanar orbits will cease producing QPEs after a time

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{prec}} & \approx & \displaystyle \frac{1}{6}\displaystyle \frac{1}{{\chi }_{\bullet }\sin I}\displaystyle \frac{{c}^{3}{a}^{3}}{{G}^{2}{M}_{\bullet }^{2}}\\ & \approx & 100\,\mathrm{days}\,{M}_{\bullet ,6}^{-2}\left(\displaystyle \frac{0.1}{{\chi }_{\bullet }I}\right){\left(\displaystyle \frac{a}{\mathrm{au}}\right)}^{3}\\ & \approx & 100\,\mathrm{days}\,{M}_{\bullet ,6}^{-1}{\tilde{\rho }}^{-1}\left(\displaystyle \frac{0.1}{{\chi }_{\bullet }I}\right),\end{array}\end{eqnarray} \tag{ 19 }$$

where in the final line we have taken a = R_RL (Equation (1)).

Precession can thus lead to a long-term modulation in the QPE activity on the timescale T_prec. For modest values of the SMBH spin and/or inclination, Equation (19) shows that T_prec is approaching the active timescale of known QPE systems—for example, eRO-QPE1 has been seen to persist for at least T_active ≳ 400 days (R. Arcodia 2021, private communication), although the low inferred stellar density in this case $\tilde{\rho }\lesssim 0.1$ (Figure 2) acts to increase T_prec. On the other hand, observations of GSN069 in 1990 (Shu et al. 2018; Miniutti et al. 2019) and eRO-QPE2 in 2014 (Arcodia et al. 2021) rule out QPE emission at the level of the present-day quiescent flux, consistent with a scenario in which precession recently brought these systems into alignment.

If future observations demonstrate that some QPEs do not turn off on the timescale ∼ T_prec, then they must arise not merely from coplanar EMRI pairs, but from pairs that lie within the SMBH equatorial plane, at least to within an angle ∼ Δa_dest/a ∼ 5R₂/a ∼ 10⁻². This has implications for the required EMRI rate in different formation channels (Section 3.7).

Finally, we note that a baseline quiescent level of SMBH X-ray accretion activity would be expected, even at times when the EMRI orbits are not aligned to enable strong periodic tidal stripping and QPE emission; this is due to the elevated mass-transfer rate of M₂, which results from it being overinflated as a result of the most recent period of tidally enhanced mass loss (Appendix A).

The maximum mass stripped per EMRI flyby is given by Δm_dest (Equation (15)). The predicted range of Δm_dest ∼10⁻⁹ − 10⁻⁷ M_⊙ is broadly consistent with the radiated X-ray energy of QPE flares (Arcodia et al. 2021), and is relatively insensitive to the free parameters (e.g., m₁, m₂, r₂, T_eff, χ). The predicted positive correlation between Δm_dest ∝ T_QPE is consistent with the flare luminosity of eRO-QPE1 (T_QPE ≈ 18.5 hr) being an order of magnitude higher than that of eRO-QPE2 (T_QPE ≈ 2.4 hr). Indeed, this trend of increasing peak luminosity with QPE period is shared by all four known QPEs (Arcodia et al. 2021).

The maximum duration of SMBH activity from EMRI tidal stripping is set by the timescale for strong encounters to destroy one or both stars. The destruction of M₂ occurs on the timescale t_dest ∼ 10³ − 10⁵ yr (Equation (14)), consistent with the upper limit on the AGN activity age T_active ≲ 10³ − 10⁴ yr in eRO-QPE1 and eRO-QPE2 based on the lack of narrow-line emission from the nuclei of their host galaxies (Arcodia et al. 2021; Section 3.6). Although the destruction of M₂ can take place over tens of thousands of years or longer, we note that flyby-powered QPEs may not be visible throughout this entire interval, due to precession of the EMRI orbital planes by the SMBH spin (on a timescale of T_prec ∼ months−years; Section 3.3).

Another effect that could potentially reduce the lifetime of the interacting EMRI system is ablation of the stars due to interaction with the gaseous accretion disk. Ablation will be particularly strong in the counter-orbiting case in which M₁ orbits in the opposite direction of the disk seeded by mass loss from M₂. In Appendix B, we estimate the ablation timescale of the star, t_abl (Equation (B3)). For typical gaseous disk properties (e.g., h/r ≳ 0.1, L_X ∼ 10⁴² erg s⁻¹), we find that t_abl ≳ t_dest ∼ 10³ − 10⁵ yr (Equation (14)) and hence gas ablation is unlikely to destroy the stars faster than their own self-interaction. Nevertheless, t_abl may be comparable to the lifetime of a pre-existing AGN or the radial migration time of the two EMRIs to their interaction radius (Section 4). In the case of a pre-existing AGN, destruction of the stellar EMRI could in principle also occur due to interaction with a relativistic jet from the SMBH (e.g., Zajaček et al. 2020) when the EMRI orbit is misaligned with the plane of the AGN disk and crosses the jet axis.

The QPE source GSN 069 (Miniutti et al. 2019) exhibited an ≈8% increase in its period T_QPE over several months of observations, while RXJ 1301.9+2747 exhibited a ≈50% change in the peak-to-peak interval between two consecutive bursts (Giustini et al. 2020). While the interacting EMRI scenario so far described would predict slow changes in T_QPE, due to evolution of the stellar orbits (from gravitational wave emission or angular momentum transfer with gaseous material), the timescale for significant changes ∼ t_dest is much longer than these observed changes.

However, one must consider possible variation in the delay between the release of gas by the (strictly) periodic tidal flyby, and the subsequent accretion onto the SMBH (see Figure 3). The interval between flares is actually the sum of the flyby period and the viscous timescale of the gaseous disk, i.e.,

$$\begin{eqnarray}&&{T}_{\mathrm{QPE}}\simeq {T}_{\mathrm{fly}}+{\tau }_{\mathrm{QPE}}.\end{eqnarray} \tag{ 20 }$$

The viscous timescale depends on the scale height of the gaseous disk τ_QPE ∝ (h/r)⁻² (Equation (16)), which in turn depends sensitively on the accretion rate. Larger accretion rates tend to lead to thicker disks (larger h/r) and hence shorter τ_QPE. Such a scenario would nominally predict shorter T_QPE following larger-amplitude flares, consistent with the observed trend of increasing T_QPE and decreasing L_QPE in GSN 069 (Miniutti et al. 2019). On the other hand, due to uncertainty in the applicability of the α − formalism in the context of various instabilities that may afflict radiation-dominated accretion flows, it is challenging to make more definitive predictions for these correlations.

Several mechanisms could give rise to stochastic or secular evolution in the amount of mass loss per flyby. If either EMRI were to possess a small eccentricity, modulations in the stellar separation between flybys will lead to large variations in the tidally stripped mass loss (Appendix A). This point is illustrated using our toy model for QPE light curves (Section 3.2) in Figure 3, panel (c). This panel makes a specific assumption for the relationship between accretion time and mass loss (t_visc,0 ∝ Δm²), and finds that ∼5% aperiodicities (as measured from peak to peak) can be obtained for a factor ≈1.5 variation in Δm, as is expected for residual eccentricities e ∼ 10⁻³ (Appendix A). Larger aperiodicities could be produced by larger residual eccentricities.

Tidal forces from the companion could also excite periodic oscillations in the mass-losing star, rendering the amount of mass loss sensitive to the phase (amplitude) of the oscillation at the time of the flyby. The accretion timescale τ_QPE (and hence T_QPE) is also sensitive to the circularization radius of the gaseous debris, which depends on its (complex) interaction with the pre-existing gaseous disk.

While the nuclei of the eROSITA QPE hosts appear to be inactive (Arcodia et al. 2021), the first two QPEs occurred in galaxies with active nuclei possessing narrow emission-line regions (Miniutti et al. 2019; Giustini et al. 2020). It is thus of interest to ask whether these galaxies are "intrinsically" active due to a pre-existing AGN disk, or whether the long-lived phase of accretion due to the interacting EMRIs studied here could power their activity.

Accretion of the total stellar mass M₂ ∼ M_⊙ over the active duration t_dest ∼ 10³ − 10⁴ yr (Equation (14)) will release a total energy E ∼ 0.1M₂ c² ≈ 2 × 10⁵³ m₂ erg in UV/X-ray radiation, sufficient to ionize M_ion ∼ (E/_Ryd)m_p ∼ 7 × 10⁶ m₂ M_⊙ of hydrogen, where _Ryd ≃ 13.6 eV is the Rydberg energy and m_p the proton mass. This is broadly consistent with the inferred masses of Seyfert 2 narrow-line regions, while the radial extent of the predicted transient narrow-line region ∼ ct_dest ∼ 0.3 − 3 kpc is also typical (e.g., Vaona et al. 2012). It thus appears possible that a system of interacting EMRIs, if caught sufficiently late in their evolution, could generate its own transient narrow-line region. However, the presence of a pre-existing gaseous AGN disk can help facilitate the migration of circular EMRIs into galactic nuclei (Section 4.2), and hence the preferential occurrence of QPEs in intrinsic AGN environments might also be expected.

Here, we provide rough estimates for the rate of circular EMRI formation needed to explain the observed QPE population. Motivated by the blind nature of the eROSITA survey, we focus on the QPE discoveries in otherwise quiescent galactic nuclei. ⁵ The first EMRI M₁ undergoes RLOF evolution over a timescale τ_GW ∼ 1−10 Myr (Equation (3)) for M₁ ∼0.3–1. This is longer than the interval between consecutive EMRIs if the latter occur at a per-galaxy rate ${\dot{N}}_{\mathrm{EMRI}}\,\gtrsim 1/{\tau }_{\mathrm{GW}}\sim {10}^{-7}-{10}^{-6}$ yr⁻¹.

The co-moving volume within the redshift z = 0.0505 of the most distant eROSITA source (eRO-QPE1) is ${ \mathcal V }\approx 0.04$ Gpc³. Using the local density of Milky Way (MW)–like galaxies of ${{ \mathcal R }}_{\star }\sim 6\times {10}^{6}$ Gpc⁻³ as a proxy for potential QPE hosts, and assuming a QPE active lifetime τ_dest (Equation (14)), the number of QPEs in the survey can be estimated as

$$\begin{eqnarray}\begin{array}{rl} & {N}_{\mathrm{QPE}}\,\sim \,{\dot{N}}_{\mathrm{EMRI}}{f}_{\mathrm{cop}}{ \mathcal V }{\tau }_{\mathrm{dest}}\\ & \,\sim \,2.4\left(\displaystyle \frac{{\dot{N}}_{\mathrm{EMRI}}}{{10}^{-7}\,{\mathrm{yr}}^{-1}}\right)\left(\displaystyle \frac{{f}_{\mathrm{cop}}}{{10}^{-2}}\right)\left(\displaystyle \frac{{\tau }_{\mathrm{dest}}}{{10}^{4}\,\mathrm{yr}}\right),\end{array}\end{eqnarray} \tag{ 21 }$$

where f_cop is the fraction of the EMRIs that are coplanar to within the range of mutual inclination $i\lesssim {i}_{\max \,}\sim 5{R}_{2}/a\sim {10}^{-2}$ that permit interactions of the type required to generate strong periodic mass loss.

The relevant value of f_cop depends on the EMRI formation channel. While ${f}_{\mathrm{cop}}\sim {i}_{\max }\sim {10}^{-2}$ for EMRIs that arrive with an isotropic distribution of inclination angles, we could expect f_cop ∼ 1 for EMRIs that arrive by migrating through a gaseous AGN disk (Section 4.2).

However, this simple rate estimate is complicated by general relativistic nodal precession of EMRI orbits inclined with respect to the SMBH spin (Section 3.3). While precession will bring even initially misaligned EMRIs into temporary alignment (on the precession timescale T_prec of months to years; Equation (19)), the limited duty cycle of the alignment compensates by increasing the required rate by a factor $\sim {i}_{\max }^{-1}\sim {10}^{2}$, so one is back to f_cop ∼ 10⁻² as in the isotropic case. On the other hand, if the observed QPE population exhibits no evidence for precession (e.g., as a "turn-off" of the QPE signal on a timescale ∼ T_prec), then the required "double coincidence," namely that the EMRI orbits must be aligned both with each other as well as with the SMBH spin, acts to increase the required rate by a larger factor ${i}_{\max }^{-2}\sim {10}^{4}$, and hence one has an effective value of f_cop ∼ 10⁻⁴ entering Equation (21).

In summary, the number of QPEs N_QPE = 2 detected thus far by eROSITA (Arcodia et al. 2021) requires a rate ${\dot{N}}_{\mathrm{EMRI}}\gtrsim {10}^{-7}$ yr⁻¹ (for f_cop = 10⁻²; precessing case) and ≳10⁻⁵ yr⁻¹ (for f_cop = 10⁻⁴; nonprecessing case). The next section explores different circular EMRI channels and to what extent they can generate these rates.

Our calculations thus far have neglected changes in the orbital separation between M₁ and M₂ that arise if the angular momentum of the stripped debris during flybys is transferred back into one or the other orbit. In the most extreme case, in which 100% of the angular momentum is transferred back to M₂, the increase in orbital separation between M₁ and M₂ due to a loss of mass Δm_fly is given by (MS17)

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{a}_{m}}{{R}_{2}}\approx 2\displaystyle \frac{a}{{R}_{2}}\left(\displaystyle \frac{{\rm{\Delta }}{m}_{\mathrm{fly}}}{{M}_{2}}\right)\approx 4\times {10}^{-4}\displaystyle \frac{{M}_{\bullet ,6}^{1/3}}{{m}_{2}^{1/3}}\left(\displaystyle \frac{{\rm{\Delta }}{m}_{\mathrm{fly}}}{{10}^{-6}{M}_{2}}\right).\end{eqnarray} \tag{ 22 }$$

For ejecta masses Δm_fly ∼ 10⁻⁶ M_⊙ (Equation (15)), Δa_m can greatly exceed the compensating decrease in the stellar orbital separation between flybys driven by gravitational wave emission, Δa_GW/R₂ ∼ 10⁻⁷ (Equation (9)).

If Δa_m ≫ Δa_GW, then the mass-loss rate per flyby will be regulated to a value smaller than what we have calculated by neglecting angular momentum added back to the orbit. Indeed, in the limit of fully conservative mass transfer, the time-averaged accretion luminosity is regulated to a value 〈L〉 ∼ 10³⁹ − 10⁴¹ erg s⁻¹ (Equation (4)), which, as already mentioned, is too small to explain observed QPE luminosities, 〈L_QPE〉 ≳ 10⁴² erg s⁻¹.

However, there are several reasons why one would expect mass transfer onto the SMBH to be highly nonconservative in this environment. For instance, the vertical scale height of the gaseous disk h/r ≳ 0.1 (as is necessary to explain the short durations of QPE flares; Section 3.1) greatly exceeds M₂'s Hill radius ∼ R₂/a ∼ 0.01 (Equation (1)). Geometrically thick, trans- or super-Eddington accretion disks are thought to launch powerful outflowing disk winds; if these winds are magnetized, they can carry away significant amounts of angular momentum, rendering mass transfer highly nonconservative.

Furthermore, even if the mass transfer is fully conservative and the accretion luminosity is fixed to the value 〈L〉 ≪ 〈L_QPE〉 (Equation (4)), limited periods of much higher mass transfer are allowed as long as they are compensated by much longer periods at lower $\dot{M}$. This very situation occurs as a result of the EMRI orbital precession (Section 3.3), which results in an active duty cycle ∼5ΔR₂/a ∼ 0.01 for flyby-induced flares when the orbits are aligned in a common plane. Increasing 〈L〉 (Equation (4)) by a factor ∼100 during the comparatively brief coplanar activity window results in accretion luminosities∼10⁴²–10⁴³ erg s⁻¹, consistent with the time-averaged QPE luminosities. Thus, although the question of whether the mass transfer is fully conservative or nonconservative may change detailed predictions of time evolution in our model, interacting EMRI systems (when active) can produce mass-transfer rates in broad agreement with QPE observations.

4.1.1. General Constraints

Stars can not be formed in situ on scales of ∼ R_RL ∼ 10² R_g;they must be delivered from larger radii. If this occurs through high-eccentricity migration from some initial semimajor axis a₀, then GW emission is the only way the orbit can circularize. Internal tidal dissipation would blow up the star long before circularization, as GM_•/R_RL ≫ GM_•/a₀.

The time required to substantially circularize an orbit of initial eccentricity e₀ ≈ 1 via GW emission is (Peters 1964)

$$\begin{eqnarray}&&{T}_{\mathrm{GW}}^{\mathrm{ecc}}\approx \displaystyle \frac{24\sqrt{2}}{85}\displaystyle \frac{{a}_{0}^{1/2}{q}_{0}^{7/2}{c}^{5}}{{G}^{3}{M}_{\bullet }{M}_{1}},\end{eqnarray} \tag{ 23 }$$

where as before we consider a star with initial mass M₁ and radius R₁, and an initial pericenter q₀ = a₀(1 − e₀). If orbital perturbations do not affect the star's angular momentum on timescales $\gtrsim {T}_{\mathrm{GW}}^{\mathrm{ecc}}$, then the star can circularize and eventually become an RLOF EMRI. However, rapid perturbations to the stellar orbital angular momentum will either move it to larger pericenter q (aborting the circularization process) or smaller q (resulting in the disruption of the star if $q\lt {R}_{{\rm{t}}}\,\equiv {R}_{1}{({M}_{\bullet }/{M}_{1})}^{1/3}$, the parabolic tidal disruption radius). The most generic source of angular momentum perturbations is two-body nonresonant relaxation, ⁶ which operates on a star of eccentricity e on the timescale

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{AM}} & \approx & {T}_{{\rm{r}}}(1-{e}^{2})\approx \displaystyle \frac{0.3}{\mathrm{ln}{\rm{\Lambda }}}\displaystyle \frac{{\sigma }^{3}(r)(1-{e}^{2})}{{G}^{2}n(r)\langle {m}^{2}\rangle }\\ & \approx & \displaystyle \frac{0.75}{(3-\gamma ){\left(1+\gamma \right)}^{3/2}}\displaystyle \frac{{M}_{\bullet }^{1/2}\langle m\rangle {r}_{\mathrm{infl}}^{3-\gamma }}{{G}^{1/2}\langle {m}^{2}\rangle \mathrm{ln}{\rm{\Lambda }}}{a}^{\gamma -5/2}q.\end{array}\end{eqnarray} \tag{ 24 }$$

Here, we have used the standard energy relaxation time T_r for an assumed power-law stellar density profile n(r) ∝ r^−γ , which has a 1D velocity dispersion $\sigma (r)=\sqrt{{{GM}}_{\bullet }/r}/(1+\gamma )$. We have defined: the influence radius r_infl inside of which the enclosed stellar mass equals M_•; the first (〈m〉) and second (〈m²〉) moments of the stellar present-day mass function (PDMF); and the Coulomb logarithm $\mathrm{ln}{\rm{\Lambda }}\approx \mathrm{ln}(0.4{M}_{\bullet }/\langle m\rangle )$.

Requiring ${T}_{\mathrm{AM}}\gt {T}_{\mathrm{GW}}^{\mathrm{ecc}}$ for a stellar EMRI can be translated into the condition that q₀ < q_GW, where

$$\begin{eqnarray}&&{q}_{\mathrm{GW}}\approx 3.2{R}_{{\rm{g}}}{\xi }^{-2/5}{\zeta }^{2/5}{\left(\displaystyle \frac{{a}_{0}}{{r}_{\mathrm{infl}}}\right)}^{-(6-2\gamma )/5},\end{eqnarray} \tag{ 25 }$$

and we have defined $\xi \equiv {\left(1+\gamma \right)}^{3/2}(3-\gamma )\,\mathrm{ln}\,{\rm{\Lambda }}$ and ζ ≡ M₁〈m〉/〈m²〉. We note that ξ ∼ 10 always, and ζ ∼ 1 usually, although ζ ≪ 1 for very low-mass target stars and/or PDMFs rich in stellar mass black holes.

For GW circularization to be possible, we also require that q_GW > R_t, which places an upper limit on the initial semimajor axis: a₀ < a_GW, where

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{\mathrm{GW}}}{{r}_{\mathrm{infl}}}\approx {\left(0.0028\right)}^{\tfrac{1}{6-2\gamma }}{\left(\displaystyle \frac{\zeta }{\xi }\right)}^{\tfrac{1}{3-\gamma }}{\left(\displaystyle \frac{{R}_{{\rm{g}}}}{{R}_{1}}\right)}^{\tfrac{5}{6-2\gamma }}{\left(\displaystyle \frac{{M}_{\bullet }}{{M}_{1}}\right)}^{\tfrac{-5}{18-6\gamma }}.\end{eqnarray} \tag{ 26 }$$

Assuming the constraints a₀ < a_GW and q₀ < q_GW indeed hold, we can now compute the residual eccentricity, ${e}_{\mathrm{res}}$, left over at the beginning of RLOF. We do this by making use of the Peters (1964) constant of motion ${c}_{0}=a{(1-{e}^{2}){e}^{-12/19}(1+121{e}^{2}/304)}^{-870/2299}$, assuming that the initial 1 − e₀ ≪ 1, and the final residual (i.e., beginning of RLOF) eccentricity ${e}_{\mathrm{res}}\ll 1$. This yields

$$\begin{eqnarray}\begin{array}{rcl}{e}_{\mathrm{res}} & \approx & 0.22{\left(\displaystyle \frac{\xi }{\zeta }\right)}^{19/30}{\left(\displaystyle \frac{{q}_{0}}{{q}_{\mathrm{GW}}}\right)}^{-19/12}{\left(\displaystyle \frac{{R}_{1}}{{R}_{{\rm{g}}}}\right)}^{19/12}\\ & & \times {\left(\displaystyle \frac{{M}_{\bullet }}{{M}_{1}}\right)}^{19/36}{\left(\displaystyle \frac{{a}_{0}}{{r}_{\mathrm{infl}}}\right)}^{19(3-\gamma )/30}.\end{array}\end{eqnarray} \tag{ 27 }$$

Setting a₀ = a_GW, we thus obtain a simple expression for the maximum value of the residual eccentricity:

$$\begin{eqnarray}&&{e}_{\mathrm{res}}\leqslant {e}_{\mathrm{res}}^{\max }\approx 0.034{\left(\displaystyle \frac{{q}_{0}}{{q}_{\mathrm{GW}}}\right)}^{-19/12}.\end{eqnarray} \tag{ 28 }$$

Interestingly, this value is large enough to produce substantial variation of the QPE amplitude between eruptions (Equation (A12)).

4.1.2. Single-star Scattering

4.1.3. Hills Mechanism

The challenge of producing stellar EMRIs from two-body relaxation alone has motivated past work to consider the Hills mechanism (Hills 1988). If a binary star with initial internal semimajor axis A_bin and total mass M_bin approaches the SMBH on a highly radial orbit, it will be tidally detached if its external pericenter q_bin is smaller than the binary detachment radius ${R}_{\mathrm{bin}}\approx {A}_{\mathrm{bin}}{\left({M}_{\bullet }/{M}_{\mathrm{bin}}\right)}^{1/3}$. One binary component will be ejected as a hypervelocity star (see Koposov et al. 2020 for a recent example from our own Galactic center), while the other star will become bound to the SMBH with pericenter q₀ ∼ q_bin and a semimajor axis a₀ ≳ a_Hills, where

$$\begin{eqnarray}&&{a}_{\mathrm{Hills}}=\displaystyle \frac{{A}_{\mathrm{bin}}}{2}{\left(\displaystyle \frac{{M}_{\bullet }}{{M}_{\mathrm{bin}}}\right)}^{2/3}.\end{eqnarray} \tag{ 29 }$$

In practice, a₀ ∼ a_Hills usually, although in a minority of cases a₀ ≫ a_Hills will occur.

The Hills mechanism is an appealing way to produce tightly bound stars with a₀ < a_GW, because the binaries can approach the SMBH from arbitrarily far away: the post-detachment a₀ of the bound star is determined primarily by A_bin. Indeed, if we require that a₀ < a_GW, we find that

$$\begin{eqnarray}&&\displaystyle \frac{{A}_{\mathrm{bin}}}{{r}_{\mathrm{infl}}}\lt 2{\left(0.0028\right)}^{\tfrac{1}{6-2\gamma }}{\left(\displaystyle \frac{\zeta }{\xi }\right)}^{\tfrac{1}{3-\gamma }}{\left(\displaystyle \frac{{R}_{{\rm{g}}}}{{R}_{1}}\right)}^{\tfrac{5}{6-2\gamma }}{\left(\displaystyle \frac{{M}_{\bullet }}{{M}_{1}}\right)}^{-\tfrac{11-2\gamma }{18-6\gamma }},\end{eqnarray} \tag{ 30 }$$

where we have approximated M_bin ≈ M₁. For Bahcall–Wolf cusps (γ = 7/4 ) and main-sequence stars, this requirement (for ξ = 10; ζ = 1 ),

$$\begin{eqnarray}&&{A}_{\mathrm{bin}}\lesssim 0.6{R}_{\odot }{M}_{\bullet ,6}{r}_{1}^{2}{m}_{1}\left(\displaystyle \frac{{r}_{\mathrm{infl}}}{0.1\,\mathrm{pc}}\right),\end{eqnarray} \tag{ 31 }$$

becomes quite restrictive and limits us to considering extremely tight binaries with A_bin ≲ R_⊙.

The rate of such "Hills EMRIs" is quite uncertain. The total rate of binary separations in simple, spherically symmetric models for an MW-like galactic nucleus is∼10⁻⁵ yr⁻¹(f_B/0.1), where f_B is the binary fraction (Yu & Tremaine 2003). This rate can increase by one to two orders of magnitude if the regions outside the SMBH influence radius have a strongly triaxial geometry (Merritt & Poon 2004) or contain massive perturbers such as giant molecular clouds (Perets et al. 2009) or nuclear spiral arms (Hamers & Perets 2017). However, only ∼1% of post-detachment stars will successfully evolve into a quasi-circular EMRI, as it is much more common for the post-detachment a > a_GW (Amaro-Seoane et al. 2012).

At the order-of-magnitude level, we thus expect the rate of quasi-circular EMRIs sourced from the total nuclear stellar population in MW-like galaxies to be ∼10⁻⁷ − 10⁻⁵ yr⁻¹. This is 1−3 orders of magnitude higher than the single-star scattering rate and broadly consistent with that needed to explain the observed eROSITA QPEs (Section 3.7) for typical source lifetimes τ_dest ∼ 10⁴ yrs.

An alternative route to producing the observed eROSITA QPE sample is to rely on secular dynamics in nuclear stellar disks. The S stars in the center of the MW may be the bound byproducts of the Hills mechanism, operating on binaries originating in a subpc disk of stars (e.g., Madigan et al. 2014; Generozov & Madigan 2020). Secular torques produced by global eccentricity features of the disk can quickly excite binaries to radial orbits that are vulnerable to the Hills mechanism. The production of ≈100 S stars in the last ≈10⁷ yr indicates a time-averaged binary detachment rate of ∼10⁻⁵ yr⁻¹ in MW-like galaxies, similar to the minimum rates calculated above. One appeal of the disk scenario, however, is that a significant fraction of the bound binary components may orbit the SMBH in roughly the same orbital plane. In particular, secular eccentricity excitations of disk members are accompanied by inclination excitation, which Wernke & Madigan (2019) find results in ∼10%–20% of the orbits at disruption being inclined in a narrow range of angles centered around ±180° relative to the stellar disk.

Another source for producing circular EMRIs is via inward migration of stars embedded within a gaseous AGN disk (Levin 2007). Stars may form in situ in AGN disks (Sirko & Goodman 2003), or alternatively may be captured by gas drag (Syer et al. 1991). Regardless of their origin, once embedded within the disk (typically at radii ≫1 au), they will migrate inward, primarily due to "Type I" torques (Goldreich & Tremaine 1980).

While Type I migration is generally inward, the gaseous torque can flip sign in small regions of AGN disks, in a way that depends on the accretion rate, viscosity, and radial location within the disk (Paardekooper et al. 2010). In particular, radial zones can develop where the net migration torque is positive, causing outward migration. The interface between an inner region exerting positive torque and an exterior region exerting a negative torque is a "migration trap" where many stellar mass objects can accumulate (Bellovary et al. 2016). We estimate the location of these traps in low-mass AGN using a Shakura & Sunyaev (1973) disk model with realistic opacities (Rogers et al. 1996) and computing Type I torques following Paardekooper et al. (2010). As shown in Figure 4, migration traps naturally develop on scales ∼10²⁻³ R_g for SMBH masses M_• ∼ 3 × 10⁵ − 3 × 10⁶ M_⊙ characteristic of QPE host galaxies.

Zoom In Zoom Out Reset image size

In active galaxies, migration traps on ∼au (i.e., ∼10² R_g) scales offer a natural mechanism for producing the coplanar EMRI pairs needed to generate QPEs. The simplest version of this mechanism involves a star (or compact object) parked in the migration trap, following a Type I inspiral that is prograde with respect to the gas. The subsequent retrograde inspiral ⁷ of an outer star eventually triggers its RLOF, with the QPEs beginning once the two stars approach sufficiently closely. One uncertainty in this formation channel concerns the modifications to RLOF for stars embedded in AGN gas. For stars that are on prograde orbits, these modifications are likely modest, as the subsonic AGN gas will simply add an extra pressure force to the outer boundary condition of the star. For stars on retrograde orbits, however, the supersonic ram pressure of the AGN gas may play a role in ablating the stellar atmosphere (Appendix B).

Even in the absence of retrograde orbiters, migration traps can accumulate large chains of prograde stars trapped in mean motion resonances (Secunda et al. 2019). After the AGN episode ends, this chain of stars could begin to undergo GW migration inward, setting off a sequence of co-orbiting EMRIs (Section 5.1). This could result in a preference for QPEs or other periodic nuclear sources in post-starburst galaxies that have undergone major mergers and associated AGN activity in the relatively recent past (indeed, RX J1301.9+2747 is hosted by a post-starburst galaxy; Giustini et al. 2020).

Counter-orbiting EMRIs are favored as the origin of the recently discovered X-ray QPEs, due to the need to generate intervals between close flybys of less than a day (see Equation (7)). However, co-orbiting collisions should also occur, which for the same stellar parameters (ρ ∼ 0.1 − 10 g cm⁻³) would predict longer QPE periods of T_QPE ∼ days to months and commensurately larger average ejecta masses Δm_dest ∝ T_QPE (Equation (15)).

A candidate for such a long-period QPE is the periodically flaring AGN, ASASSN-14ko, which exhibits outbursts at regular intervals of around 114 days (Payne et al. 2021). During the rise of its ouburst in May 2020, ASASSN-14ko exhibited UV-bright, thermal spectral energy distribution similar to tidal disruption events. However, the X-ray flux decreased by a factor of ≈4 at the beginning of the outburst before returning to its quiescent flux after ∼8 days. The large inferred black hole mass M_• ∼ 10⁸ M_⊙ for ASASSN-14ko would require two stars with ρ ≲ 0.1ρ_⊙ for RLOF to occur outside the ISCO radius. For the same parameters, Equation (8) predicts a QPE period of T_QPE ∼ 100 d, consistent with the activity period in ASASSN-14ko. The EMRI destruction time for such a massive SMBH can be relatively relatively short, ${t}_{\mathrm{destr}}\sim 100\,\mathrm{yr}$ (Equation (14)), and so appreciable evolution of the system could be observable.

Zalamea et al. (2010) consider a scenario in which a WD on an eccentric orbit undergoes periodic RLOF (or equivalently, partial tidal disruption) onto an SMBH, feeding gas onto the SMBH and powering a quasi-periodic string of flares. A variant of this scenario was proposed in King (2020), and both could, in principle, appear as QPEs (Arcodia et al. 2021). This eccentric single-star partial disruption scenario requires a high-eccentricity orbit for a typical WD density ρ ∼ 3 × 10⁵ ρ_⊙, as the orbital period at this object's Roche radius is only ∼1 minute (Equation (5)). Reproducing observed QPE periods of ≈2–19 hours would thus require orbits with eccentricity e ≳ 0.9–0.99.

One challenge to this scenario is the difficulty of putting a single star on a significantly eccentric orbit with a semimajor axis of ∼1 au. King (2020) invokes the tidal disruption of the envelope of a red giant star, which leaves behind a degenerate core on an orbit with appropriate parameters. This version of the single-star scenario has two major challenges:

• 1.  
  The magnitude of the "core kick" is far too low to put a surviving core on an orbit with period ∼ T_QPE. The origin of core kicks lies in the deviation tensor (the third-order expansion of the gravitational potential), which encodes the asymmetries of the tidal field (Brassart & Luminet 2008; Cheng & Evans 2013). As this is simply the next-order expansion of the gravitational potential, beyond the second-order tidal expansion, we can estimate the magnitude of the specific orbital energy perturbation to the surviving core on dimensional grounds (in analogy to the reasoning of Stone et al. 2013) as

  $$\begin{eqnarray}&&\delta \epsilon \sim \displaystyle \frac{{{GM}}_{\bullet }}{{R}_{{\rm{t}}}}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{{\rm{t}}}}\right)}^{2}\displaystyle \frac{{\rm{\Delta }}M}{{M}_{\mathrm{core}}}\sim \displaystyle \frac{{{GM}}_{\star }}{{R}_{\star }}\displaystyle \frac{{\rm{\Delta }}M}{{M}_{\mathrm{core}}}.\end{eqnarray} \tag{ 32 }$$

  Here, we have considered the partial disruption of a star with initial mass M_⋆ and radius R_⋆, which leaves behind a surviving core of mass ${M}_{\mathrm{core}}={M}_{\star }-{\rm{\Delta }}M$ following a close encounter near or inside the tidal radius ${R}_{{\rm{t}}}={R}_{\star }{({M}_{\bullet }/{M}_{\star })}^{1/3}$. This order-of-magnitude estimate is in good agreement with numerical hydrodynamic simulations of core kicks from partial disruptions (Manukian et al. 2013; Gafton et al. 2015). For the disruption of a red giant (M_⋆ ∼ 1M_⊙; R_⋆ ∼ 200R_⊙) with initial specific orbital energy _RG, it predicts a final specific energy for the bound core of _RG ± δ , where δ ∼ 10⁻⁸ c². Since δ ≪ GM_•/R_RL ∼ 0.01c², there is no way for the leftover core to have an orbital period comparable to T_QPE.
• 2.  
  The sign of the core kick is likely positive-definite, i.e., δ > 0, such that surviving cores always become more loosely bound (rather than more tightly bound). In the partial disruption simulations of Faber et al. (2005), Manukian et al. (2013), Gafton et al. (2015), significant mass loss in a partial disruption is always associated with a positive energy kick to the surviving core.

The production mechanisms and rates of eccentric WD EMRIs are not considered extensively by Zalamea et al. (2010), although they suggest two-body scattering of a single WD or the Hills mechanism separating a binary with at least one WD component. Both of these possibilities are disfavored on rates grounds. For single-star scattering, we have already seen (Section 4.1.2) that rates are negligibly small, and they become even smaller when considering (relatively uncommon) WDs.

This leaves the Hills mechanism as the favored way to produce stars with $a={(2\pi )}^{-2/3}{({{GM}}_{\bullet })}^{1/3}{T}_{\mathrm{QPE}}^{2/3}$. Recalling that the post-separation semimajor axis of the bound star a₀ ≳ a_Hills, in order for a single star to be born into an orbit with a period equal to T_QPE, the semimajor axis of the original binary must obey

$$\begin{eqnarray}&&{A}_{\mathrm{bin}}\lesssim 0.1{R}_{\odot }{\left(\displaystyle \frac{{T}_{\mathrm{QPE}}}{10\,\mathrm{hr}}\right)}^{2/3}{\left(\displaystyle \frac{{M}_{\bullet }}{{10}^{5}{M}_{\odot }}\right)}^{-1/3}{\left(\displaystyle \frac{{M}_{\mathrm{bin}}}{1.2{M}_{\odot }}\right)}^{2/3}.\end{eqnarray} \tag{ 33 }$$

This disfavors the tidal detachment of a main-sequence or helium star binary (Zhao et al. 2021) as capable of generating QPEs as short as 10 hr. A binary composed of two WDs could satisfy Equation (33), but any such tight binary has its own problem: a short lifetime. The GW inspiral time of an equal-mass WD binary is only

$$\begin{eqnarray}&&{T}_{\mathrm{GW}}\approx 3\times {10}^{4}\,\mathrm{yr}\,{\left(\displaystyle \frac{{A}_{\mathrm{bin}}}{0.1{R}_{\odot }}\right)}^{4}{\left(\displaystyle \frac{{M}_{\mathrm{bin}}}{1.2{M}_{\odot }}\right)}^{-3},\end{eqnarray} \tag{ 34 }$$

rendering such systems exceedingly rare.

Quantitatively, the rate of Hills separation of double WD binaries can be written as ${\dot{N}}_{\mathrm{Hills}}^{\mathrm{WD}}\sim {\dot{N}}_{\mathrm{Hills}}{f}_{\mathrm{WD}}{f}_{\mathrm{hard}}$, where ${\dot{N}}_{\mathrm{Hills}}\sim {10}^{-5}-{10}^{-3}\,{\mathrm{yr}}^{-1}$ in an MW-type galaxy (Section 4.1.3). Here, f_WD is the fraction of all tidally detached binaries comprised of two WDs, and is likely ≲ 0.002 (the total WD number fraction for a Salpeter IMF and an old stellar population). The factor f_hard is the fraction of all double WD binaries with A_bin less than the critical value given by Equation (33); lifetime arguments imply that f_hard < T_GW/T_H ∼4 × 10⁻⁶, where T_H is the Hubble time. Taken together, the total rate at which the Hills mechanism deposits single WDs onto sufficiently short-period orbits to explain the observed QPEs is ${\dot{N}}_{\mathrm{Hills}}^{\mathrm{WD}}\lesssim {10}^{-10}\,{\mathrm{yr}}^{-1}$ per MW-type galaxy, orders of magnitude below what is required by observations.

The rates problem is further exacerbated by the short predicted WD lifetime once Roche overflow starts. Mass transfer onto the SMBH is likely to be unstable due to the inverted mass–radius relation of WDs and the nonconservative nature of mass transfer in highly eccentric binaries. In the fiducial example given by (Zalamea et al. 2010, their Figure 2), the WD only loses a fraction ∼10⁻⁶ − 10⁻⁸ of its mass (as required to explain QPE amplitudes) for a few hundred orbits. This short lifetime ∼1 month is in tension with archival X-ray detections of RX J1301.9+2747 and GSN 0691 going back decades (Miniutti et al. 2019; Giustini et al. 2020). If the QPE lifetime is τ_dest ≲ 1 yr, then the required WD EMRI rate to explain the eROSITA QPE sample is ${\dot{N}}_{\mathrm{EMRI}}\sim {10}^{-5}$ gal yr⁻¹ (Equation (21)), five orders of magnitude larger than we have estimated above for the Hills mechanism. Finally, the small disk radius predicted in WD RLOF scenarios is at odds with the ultraviolet dimming observed concurrently with the X-ray QPEs from XMMSL1 J024916.6-041244, which supports the gaseous disk production occurring on a radial scale ∼10³ R_g (Chakraborty et al. 2021).