© The Authors 2023

1. Introduction

X-ray quasi-periodic eruptions (QPEs) are an extreme and rather puzzling X-ray variability phenomenon that has recently been detected from the nuclei of nearby galaxies. They are fast bursts in the soft X-ray band, repeating every few hours, superimposed to an otherwise stable quiescent X-ray level that is consistent with emission from a radiatively efficient accretion flow around relatively low-mass massive black holes (MBHs). During these bursts, the X-ray count rate increases by up to two orders of magnitude, depending on the considered energy band. QPEs have thermal-like X-ray spectra with temperature evolving from kT ≃ 50 − 80 eV to ≃100 − 250 eV (and back) in about one to a few hours, with a typical duty cycle of 10 − 30%, depending on the specific source. The peak X-ray luminosity is ≃10^{42 − 43} erg s⁻¹. Following their discovery in GSN 069 (Miniutti et al. 2019), QPEs have been observed so far from the nuclei of other four galaxies: RX J1301.9+2747 (Giustini et al. 2020), eRO-QPE1, and eRO-QPE2 (Arcodia et al. 2021), as well as XMMSL1 J024916.6−04124 (Chakraborty et al. 2021). Only 1.5 QPEs were detected in the latter source, but all observed properties indicate it is nevertheless a robust QPE candidate.

We may classify QPE host galaxies as relatively low-mass post-starburst galaxies (Arcodia et al. 2021; Wevers et al. 2022), a population that is very similar to that of the hosts of the preferred tidal disruption events (TDEs, French et al. 2020), with which they also share a preference towards low-mass MBHs of 10^{5 − 6.7} M_⊙ (Wevers et al. 2022). Two out of the five QPE sources known to date have been associated with X-ray detected TDEs (Miniutti et al. 2023a; Chakraborty et al. 2021), and a new QPE candidate is associated with the X-ray decay following an optically detected TDE (Quintin et al. 2023), further strengthening the QPE-TDE connection. Despite very little X-ray obscuration, none of the QPE galaxies show the typical optical or UV broad emission lines associated with unobscured, type 1 active galactic nuclei (AGNs), although the observed narrow-line ratios suggest the presence of an ionising AGN-like continuum (Wevers et al. 2022). The lack of broad emission lines may be due to an AGN that has switched off sometime in the past leaving relic narrow emission lines or be perhaps associated with TDE-like accretion discs that are likely too compact to support a mature AGN-like broad line region.

Several theoretical models have been proposed to explain QPEs, based on modified disc instabilities (Raj & Nixon 2021; Pan et al. 2022; Kaur et al. 2022), gravitational lensing (Ingram et al. 2021), mass transfer from one or more bodies in different configurations (King 2020, 2022; Chen et al. 2022; Metzger 2022; Zhao et al. 2022; Krolik & Linial 2022; Lu & Quataert 2022; Wang et al. 2022; Linial & Sari 2023), and collisions between an orbiting secondary object and the accretion accretion flow that is formed around the primary MBH (Suková et al. 2021; Xian et al. 2021; Linial & Metzger 2023; Tagawa & Haiman 2023)¹. Although these models can explain some of the features observed in the light curves of the confirmed QPE sources (mainly GSN 069) there is no comprehensive model able to reproduce the variety of behaviours observed in the QPEs of different sources. Here, our efforts are aimed at building a comprehensive and flexible model that can reproduce these variabilities and provide a possible emission mechanism that is compatible with the observations.

We propose a new model in which QPEs are produced by extreme mass ratio inspirals (EMRIs). Those are binaries in which a lighter companion (usually a stellar mass black hole, BH) adiabatically inspirals onto a MBH under the effect of gravitational wave (GW) back-reaction. In the specific scenario we propose, the observed X-ray flares emerge as the companion crosses an accretion disc that is rigidly precessing around the central MBH. This is somewhat similar to the model proposed by Xian et al. (2021), but it also includes the relativistic Lense-Thirring (Lense & Thirring 1918) precession of the accretion disc. We include both apsidal and nodal precession of the EMRI companion up to 3.5PN order beyond Newtonian dynamics. Our proposed model has the flexibility to reproduce the wide variety of QPE light curves including luminosity variations and timing irregularities, owing to the various orbital and precession frequencies involved in the aforementioned processes and their interplay. We emphasise that our aim, at this stage, is not to apply our model to a specific source for parameter inference, but to show that the different amplitudes and timing properties of the QPE phenomenon can be qualitatively reproduced.

The paper is organised as follows. In Sect. 2, we describe our model in detail and we present the resulting light curves in Sect. 3. We then draw our conclusions and discuss future improvements of the model in Sect. 4.

2. Methods

Using a semi-analytic approach, we modelled a system composed of a MBH with mass, M₁, and spin parameter, χ, surrounded by a misaligned and rigidly precessing accretion disc, whose properties are compatible with a disc generated by a stellar TDE (see Sect. 2.2). The MBH is further orbited by a stellar-mass BH (the EMRI companion, hereafter) with a mass, M₂, that intersects the accretion disc (up to three times per orbit; see Sect. 2.3). The mass ratio of the system q = M₂/M₁ is assumed to be q ≪ 1 throughout the paper – a feature that allows us to primarily characterise the binary system as an EMRI (see Fig. 1 and Table 1).

Fig. 1. Schematic view of the model where the central big black dot represents the primary MBH, the small black dot represents the EMRI companion, the grey area represents the tilted disc and the EMRI-disc crossings are shown with yellow stars. The vertical arrow shows the MBH spin vector. The misalignment between the EMRI orbit and the disc is enhanced in this cartoon with the purpose of showing the impacts more clearly. However, we use smaller misalignments to produce the QPEs light curves. The blue and yellow arrows show the EMRI and gas velocity vectors, respectively, while the yellow dashed line only mimics an orbit inside the disc.

2.1. EMRI post-Newtonian trajectory

We integrated the equations of motion of the EMRI system considering post-Newtonian (PN) dynamics and retaining corrections up to 3.5 PN order. Integer terms such as 1, 2, or 3 PN represent conservative terms, while 2.5 and 3.5 PN orders are dissipative terms accounting for GW emission. Finally, the 1.5 (conservative) PN term takes into account the effect of the MBH spin, specifically: the leading order spin-orbit coupling. We consider the equations of motion derived using a Lagrangian formalism and referred to the relative position and velocities of the binary components. Practically, we evolve forward in time the relative separation, r, and velocity, v, of the EMRI components according to:

$$\begin{array}{c}\hfill \frac{d^2\mathit{r}}{\mathrm{d}{t}^2}=-\frac{\mathit{GM}}{r^2}((1+\mathcal{A})\mathit{n}+\mathcal{B}\mathit{v})+{\mathit{C}}_{1.5}+\mathcal{O}\left(\frac{1}{c^8}\right),\end{array}$$(1)

where n = r/|r|, while M is the binary total mass. The coefficients 𝒜 and ℬ can be found in Blanchet (2014) and encode the interaction due to the non-spin terms, while the coefficient C_1.5 introduces the leading order spin-orbit interaction and can be retrieved from Bohé et al. (2013). We employed a C++ implementation of the Bulirsch-Stoer algorithm to integrate the above differential equations. More details can be found in Bonetti et al. (2016), except that here we do not include the interaction with a third body or the hardening of the binary due to the presence of a stellar cluster, with the latter assumption being justified at such small separations. The main advantage of this approach, compared to the one adopted by Xian et al. (2021), is that we are not limited by the test particle approximation and therefore in the choice of the EMRI companion mass.

2.2. Disc model

Based on the growing observational evidence that QPEs arise after a TDE-like precursor (Chakraborty et al. 2021; Miniutti et al. 2023a), we considered a disc formed by a tidal disruption of a star by the primary MBH and modelled it following Franchini et al. (2016). The main difference between the current model and Franchini et al. (2016) concerns the accretion rate onto the central MBH, that in the current investigation is set to reproduce the observed nuclear activity when QPEs are not occurring (the “quiescence” luminosity, hereafter).

More in detail, the disc extends all the way to the innermost stable circular orbit (ISCO), which is a function of the primary BH spin χ (Bardeen et al. 1972), and we assumed the disc to be prograde with respect to the spin of the MBH, i.e. the z component of its angular momentum is positive.

We assumed the mass of the disc to be M_d = f M_⊙ with f ≤ 4 and to be distributed with a power law surface density profile Σ = Σ₀(R/R_g)^−p out to a radius R_out = 200 − 400 R_g, depending on the primary mass, where R_g = GM₁/c² is the MBH gravitational radius. The normalization of the density profile can be therefore computed in a straightforward way as:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_0=\frac{M_{\mathrm{d}}(2-p)}{2\pi R_{\mathrm{g}}^2}({\left(\frac{R_{\mathrm{out}}}{R_{\mathrm{g}}}\right)}^{2-p}-{\left(\frac{R_{\mathrm{ISCO}}}{R_{\mathrm{g}}}\right)}^{2-p}),\end{array}$$(2)

with p as the power law index.

We define the disc aspect ratio as (Strubbe & Quataert 2009; Franchini et al. 2016)

$$\begin{array}{c}\hfill \frac{H}{R}=\frac{3}{2}{(2\pi )}^{1/2}{\eta }^{-1}K{(R)}^{-1}{\dot{m}}_1{\left(\frac{R}{R_{\mathrm{g}}}\right)}^{-1},\end{array}$$(3)

where

$$\begin{array}{c}\hfill K(R)=\frac{1}{2}+{[\frac{1}{4}+6{\left(\frac{{\dot{m}}_1}{\eta }\right)}^2{\left(\frac{R}{R_{\mathrm{g}}}\right)}^{-2}]}^{1/2}.\end{array}$$(4)

The K(R) coefficient takes into account the modified structure for a slim disc. Here, η = 0.1 is the accretion efficiency and ṁ₁ = Ṁ₁/Ṁ_1,Edd, namely, the mass accretion rate of the central MBH normalised to its Eddington limit.

We assumed that the angular momentum of the disc is misaligned with respect to the MBH spin and that the warp induced by the Lense-Thirring (LT) effect (Lense & Thirring 1918) propagates as a bending wave allowing the disc to rigidly precess around the MBH (Franchini et al. 2016). We modelled the disc viscosity according to the Shakura & Sunyaev α-prescription (Shakura & Sunyaev 1973). In the assumed bending waves regime, α is lower than H/R, so we can reasonably assume that the alignment timescale is long enough to neglect any disc-spin alignment effects (Stone & Loeb 2012; Franchini et al. 2016).

Under these assumptions², the disc is predicted to precess rigidly with a frequency Ω_p that is the angular-momentum-weighted average of the Lense-Thirring precession Ω_LT(R) frequency over the extent of the disc:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{p}}=\frac{{\int }_{R_{\mathrm{ISCO}}}^{R_{\mathrm{out}}}{\mathrm{\Omega }}_{\mathrm{LT}}(R)L(R)2\pi R\mathrm{d}R}{{\int }_{R_{\mathrm{ISCO}}}^{R_{\mathrm{out}}}L(R)2\pi R\mathrm{d}R},\end{array}$$(5)

with the Lense-Thirring frequency given by:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{LT}}(R)=\frac{c^3}{2G{M}_1}\frac{4\chi {\left(\frac{R}{R_{\mathrm{g}}}\right)}^{-3/2}-3{\chi }^2{\left(\frac{R}{R_{\mathrm{g}}}\right)}^{-2}}{{\left(\frac{R}{R_{\mathrm{g}}}\right)}^{3/2}+\chi },\end{array}$$(6)

and with L(R) = ΣΩ(R)R² being the angular momentum modulus at a radius, R. We assumed that the gas orbits the MBH on circular Keplerian orbital thus, the gas orbital frequency and its velocity modulus are simply given by Ω(R) = (GM₁/R³)^1/2 and |v_gas(R)| = (GM₁/R)^1/2, respectively.

We stress that we have assumed a disc model broadly consistent with the tidal disruption of a star by the primary MBH. However, the disc properties are poorly constrained parameters and we can therefore vary them to reproduce different observed timing behaviours in QPE sources.

We assumed that the disc is precessing rigidly around the primary MBH, as there has been (to date) no precise criterion for a disc to precess rigidly other than order of magnitude estimates that compare the sound crossing time to the viscous time, with the latter being largely unconstrained. If the annulus of the disc that corresponds to the crossing does precess with its own frequency, then this will probably cause a different modulation in the light curves. We leave the investigation of a differentially precessing disc annuli to a future work.

We further note that precession is not crucial when comparing our model with single-epoch QPE light curves, as it induces a modulation on timescales that are longer than the typical observation duration (∼1.5 days, except in the case of eRO-QPE1). Future, multi-epoch observations of QPE sources are needed to detect unambiguously the putative modulations induced by disc precession.

2.3. Disc crossing

The EMRI companion crosses the disc at different radii between one and three times per orbit around the MBH since the encounter can occur at pericentre and along the precessing disc line of nodes (similar to the model of non periodic outburst in Be/X-ray binaries proposed by Martin & Franchini 2021).

The apsidal and nodal precession of the EMRI companion further complicate the time series of disc crossings making the interval between subsequent episodes change during time, as similarly observed during the different campaigns in the light curve of RX J1301.9+2747 (Giustini et al. 2020, and in prep.).

In our model, the EMRI evolution is computed from the numerical integration of Eq. (1), thus the characterisation of the disc crossings has to be done numerically as well. We stress that in our model, the disc is represented by a rigidly rotating plane that does not influence the EMRI evolution. In order to completely characterise the crossings and infer the numerical values of the physical quantities involved, we proceed as follows. As a first step, we identify the times and positions at which the EMRI companion crosses the disc plane. Specifically, those are found along the EMRI trajectory when the coordinates r₂ = (x₂, y₂, z₂) of M₂ satisfy the relation w_xx₂ + w_yy₂ + w_zz₂ = 0. Here, w = (w_x, w_y, w_z) represents the vector normal to the disc plane, which changes with time according to³

$$\begin{array}{cc}\hfill w_x& =sin({\mathrm{\Omega }}_{\mathrm{p}}t)sin(i_{\mathrm{disc}}),\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill w_y& =cos({\mathrm{\Omega }}_{\mathrm{p}}t)sin(i_{\mathrm{disc}}),\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill w_z& =cos(i_{\mathrm{disc}}).\hfill \end{array}$$(9)

Given the crossing location, we can determine the velocity vector of the gas (v_gas) at each passage. We calculate the components of r₂ in the disc plane and then the polar angle ϕ_d = arctan(y_2, d/x_2, d) in such a frame. Since the gas is assumed to be in circular motion, we know that the gas velocity has to be perpendicular to the projection of r₂ in the disc plane r_2, d, thus the gas velocity in the reference plane of the disc is given by:

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathrm{gas},\mathrm{d}}={\left(\frac{G{M}_1}{R_{\mathrm{cross}}}\right)}^{1/2}(-sin({\varphi }_{\mathrm{d}}),cos({\varphi }_{\mathrm{d}}),0),\end{array}$$(10)

where R_cross is the radius at which the crossing occurs. Finally, the velocity vector above is rotated back in the initial frame and the relative velocity between the gas and M₂ is evaluated from the vector difference, v_rel = v_gas − v₂. We assumed the relative velocity does not change during the crossing as this occurs on a short timescale. For instance, for GSN 069, this is 2H/v₂ ∼ 200 s.

We note that the nodal LT precession of the disc does in principle make it possible for the precession velocity to eventually point in the opposite direction with respect to the EMRI velocity. However, in all of our cases, the precession of the disc is slower than the orbital period and the precession velocity always points in the same direction as the orbital one. In other words, we never have a full precession of the disc within an EMRI orbital period.

The crossing times determined by the above procedure are referred to the reference frame of the binary. In order to relate the time at which the crossing occurs to the arrival time of emitted photons at the observer, we included the relevant delays, specifically: the Roemer, Shapiro, and Einstein delays. These delays take into account the fact that, once photons are released by the EMRI companion crossing the disc, they have to propagate towards a distant observer, thus they are influenced by the light-travel time within the binary system (Roemer delay Δ_R(τ)) and the non-zero curvature generated by the primary companion which bends the travel path (Shapiro Δ_S(τ) and Einstein delays Δ_E(τ)). Specifically, we express the arrival time to the distant observer as (see e.g. Poisson & Will 2014, Sect. 10.3.6) as:

$$\begin{array}{c}\hfill t_a=\tau +{\mathrm{\Delta }}_{\mathrm{R}}(\tau )+{\mathrm{\Delta }}_{\mathrm{S}}(\tau )+{\mathrm{\Delta }}_{\mathrm{E}}(\tau ),\end{array}$$(11)

where the three delays are given by the following expression:

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{R}}(\tau )& =|{\mathit{r}}_{\mathrm{obs}}-{\mathit{r}}_2|/c,\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{S}}(\tau )& =\frac{2G{M}_1}{c^3}ln\left(\frac{|{\mathit{r}}_{\mathrm{obs}}-{\mathit{r}}_1|+({\mathit{r}}_{\mathrm{obs}}-{\mathit{r}}_1)·\mathit{k}}{r+\mathit{r}·\mathit{k}}\right),\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{E}}(\tau )& =\frac{M_2+2M_1}{M}\sqrt{\frac{a^3}{\mathit{GM}}}\frac{G{M}_1}{a{c}^2}sin(u).\hfill \end{array}$$(14)

In the above equations, M = M₁ + M₂ is the mass of the binary, r_obs denotes the observer position, r₁ and r₂ represent the position vectors of the primary and secondary BH with respect to the centre of mass of the system, while the unit vector k is given by:

$$\begin{array}{c}\hfill \mathit{k}=\frac{{\mathit{r}}_{\mathrm{obs}}-{\mathit{r}}_2}{|{\mathit{r}}_{\mathrm{obs}}-{\mathit{r}}_2|}·\end{array}$$(15)

Here a, e, and u respectively represent the orbit semi-major axis, eccentricity and eccentric anomaly, respectively.

2.4. Emission mechanism

Since we are interested in modelling the light curve that results from the disc crossings, we need to investigate possible emission mechanisms that can reproduce the amplitude, duration, and spectral properties of QPEs.

The emission mechanisms that have thus far been proposed in the literature to explain QPEs, within the disc-crossing model, consider either accretion onto the secondary object or the emission generated by the shock as the companion passes through the accretion flow (Suková et al. 2021; Xian et al. 2021). We can safely rule out the first hypothesis since, in order to produce the observed luminosities of ∼10⁴² − 10⁴³ erg s⁻¹, the accretion onto the EMRI companion would greatly exceed its Eddington limit. As for the shock mechanism, this has been investigated and modelled extensively in Nayakshin et al. (2004) where the authors considered a star that passes through an accretion disc surrounding a MBH, generating a shock due to the compression of the disc material in front of and aside of it. However, some of the assumptions in their work are tailored to the specific case of Sgr A^* and are not well suited to the case of the radiatively efficient, high-mass accretion rate discs that are relevant for QPE sources.

Here, we instead consider the emission from a gas cloud that is pulled out of the disc as a consequence of the crossing of the EMRI companion, broadly following the approach outlined in Pihajoki (2016). We assume the gas cloud to have an initial radius, R_in, that is of the order of the influence radius of the secondary object that crosses the disc, namely:

$$\begin{array}{c}\hfill R_{\mathrm{in}}\sim R_{\inf }=\frac{G{M}_2}{c_{\mathrm{s}}^2+v_{\mathrm{rel}}^2},\end{array}$$(16)

where c_s is the disc sound speed, while v_rel is the modulus of the relative velocity (see previous section for details). We consider the cloud to be optically thick (we assume its density to be equal to the disc density at the crossing radius) and with a post-shock temperature of ∼10⁶ K, which corresponds to the observed QPE emission of ∼100 − 120 eV. The exact value of the cloud initial (post shock) temperature, T₂, depends on the source considered and it is determined by using the Rankine-Hougoniot condition for shocks in a radiative pressure dominated gas disc:

$$\begin{array}{c}\hfill T_2={[1+\frac{8}{7}({\mathcal{M}}_{\mathrm{e},1}^2-1)]}^{1/4}{T}_1,\end{array}$$(17)

where T₁ is the pre-shock disc gas temperature and ℳ_e, 1 is the effective pre-shock Mach number, given by:

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{e},1}={\left(3\gamma \frac{{\mathcal{M}}_1^2}{4R_{\mathrm{p}}}\right)}^{1/2},\end{array}$$(18)

where ℳ₁ = v_rel/c_s is the pre-shock Mach number, R_p is the ratio between the radiation and gas pressure before the shock occurs, and γ = 4/3 is the adiabatic index. The pre-shock temperature T₁ is determined by the properties of the disc and is obtained from the relation linking the total pressure to the gas sound speed, namely⁴:

$$\begin{array}{c}\hfill P=P_{\mathrm{rad}}+P_{\mathrm{gas}}=\frac{\rho c_{\mathrm{s}}^2}{\gamma },\end{array}$$(19)

with the radiation and gas pressure respectively given by

$$\begin{array}{cc}\hfill P_{\mathrm{rad}}& =\frac{4}{3}\frac{{\sigma }_{\mathrm{sb}}{T}_1^4}{c},\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill P_{\mathrm{gas}}& =\frac{\rho k_{\mathrm{B}}{T}_1}{m_{\mathrm{p}}\mu },\hfill \end{array}$$(21)

with σ_sb being the Stefan-Boltzmann constant, k_B the Boltzmann constant, μ = 1 the mean molecular weight, and m_p the proton mass.

We note that our model implies that T₂/T₁ is a decreasing function of mass accretion rate ṁ₁, asymptotically reaching T₂/T₁ at high ṁ₁. As QPEs only carry a small fraction of the disc bolometric luminosity (Miniutti et al. 2023a), they become undetectable once T₂/T₁ ≃ 1. Our result appears therefore to be consistent with observations of GSN 069 where QPEs are not detected above a disc luminosity threshold that corresponds to an Eddington ratio of ≃0.4 ± 0.2 (Miniutti et al. 2023b).

We assumed the temperature of the cloud to increase from T₁ to T₂ instantaneously, giving us a fast rise in the QPE profile. This is an assumption we made for simplicity, as we are mostly interested in reproducing the timing and peak luminosity properties rather than the exact individual QPE profile (see Sect. 3). This cloud of gas might take some time to exit from the disc and the gas would also have to thermalise before emitting at T₂. This assumption can likely be relaxed by considering the effects of photon diffusion and disc crossing timescales⁵ that induce a rise on a slightly shorter timescale compared to the decay.

We assumed the gas cloud to emit black body radiation while it adiabatically expands after leaving the dense and pressure confining accretion disc environment. The actual properties of the hot atmosphere around the accretion disc are not observationally constrained and not derivable from first principles in our model. However, QPEs can only be detected if their peak temperature, T₂, is significantly higher than the pre-shock disc emission, T₁. Observationally, QPEs are characterised by T₂ ≃ 2 − 3 × T₁. We therefore tailored the expansion of the ejected cloud so that the initial temperature drops by a factor of 2 − 3 during the typical QPE duration assuming adiabatic expansion. The radius of the cloud changes over time as:

$$\begin{array}{c}\hfill R(t)=R_{\mathrm{in}}+\frac{2R_{\mathrm{in}}}{\mathrm{\Delta }{t}_{\mathrm{QPE}}}t,\end{array}$$(22)

where Δt_QPE is the duration of the QPE starting from the time of the impact, while the temperature at which the black body spectrum is emitted changes with time as T_exp = T₂(R_in/R(t))⁶.

Assuming the cloud initial radius to be R_in ∼ 10¹¹ cm and for it to expand by a factor of 2−3 in ∼1 h, we infer an expansion velocity v_exp ∼ 400 km s⁻¹. This is consistent with the estimate v_exp ∼ v_rel(R_in/H) (Arnett 1980; Pihajoki 2016), with v_rel being the relative velocity between the EMRI companion and the disc, which we infer in the case of GSN 069 from our model. We note that the assumed expansion factor is consistent with the observed data of GSN 069 (e.g. lower panel of Fig. 18 in Miniutti et al. 2023a).

The emitted luminosity in the soft-X ray band 0.2 − 2 keV is therefore:

$$\begin{array}{c}\hfill L_{\mathrm{X}}=4\pi R{(t)}^2{\displaystyle {\int }_{0.2\mathrm{keV}}^{2\mathrm{keV}}\frac{2h{\nu }^3}{c^2}\frac{\mathrm{d}\nu }{e^{h\nu /k_{\mathrm{B}}{T}_{\exp }}-1}·}\end{array}$$(23)

As discussed in Sect. 2.5, we specifically consider the case in which the secondary is a BH of mass, M₂. We assume the initial cloud radius R_in to be either determined by the secondary influence radius, R_inf, or by the disc half-height at the crossing radius, H_cross, depending on which one is the smallest. For our choice of parameters we have R_inf ≲ H_cross for all the sources we consider.

We then construct the light curve by assuming that the emission starts at the time of the impact. We then assign a time length that corresponds to the observed QPE duration (therefore, different for each source and tailored to the specific observed burst duration for that source) and an amplitude that follows the decay in the luminosity according to Eq. (23).

As an example, Fig. 2 shows the spectrum of the radiation emitted in the band 10 eV −2 keV by a 10⁶ K cloud that detaches from the disc as a consequence of the EMRI companion passage through it. The cloud initial size is R_in = R_inf and it is assumed to expand by a factor of 3 in radius over the 1 h duration of the observed QPE signal. The exact value depends on R_inf, which in turn depends on the relative velocity between the EMRI companion and the disc – and, thus, ultimately on the source considered.

Fig. 2. Spectrum emitted by a cloud of initial size Rinf and initial temperature T = 106 K that expands by a factor of 3 in size over one hour.

2.5. The secondary object

The QPE timing properties and the initial temperature of the expanding cloud producing them, T₂, are set by the disc properties and by the EMRI dynamics. Thus, these parameters are independent of the nature and mass of the secondary, as long as q ≪ 1. On the other hand, the QPE (peak) luminosity depends on the initial size of the expanding cloud, R_in, as shown in Eqs. (22) and (23). Thus, only impacts that form clouds with R_in ≃ 10¹¹ cm are able to produce L_X ≃ 10⁴² − 10⁴³ erg s⁻¹ at QPE peak, as observed.

For compact companions such as BHs, neutron stars, and white dwarfs, the initial size of the expanding cloud is R_in ∼ R_inf and it thus depends on companion mass and relative velocity at impact, v_rel. It is easy to show that BHs with masses ≪100 M_⊙ never reach the required R_in even for very low v_rel. The same applies to neutron stars and white dwarfs that can therefore be safely ruled out with respect to the framework of the simple emission mechanism we have considered. Slightly less massive BHs, that is, around and above 30−40 M_⊙ are still viable candidates given that the luminosity of the flares is also significantly influenced by the disc properties (e.g. denser clouds produce brighter flares in our model). On the other hand, the lifetime of the EMRI system, regulated by the GW emission, becomes too short for significantly more massive BHs to be consistent with the recurrence rate of the QPEs. Thus, BHs with mass ∼100 M_⊙⁷ can potentially produce clouds with sufficiently large R_in provided that v_rel is low enough; that is to say, if the secondary orbit is prograde with respect to the accretion disc and the misalignment between the orbit and disc planes is not too large. This is precisely the situation we assume in our model.

In general, we cannot completely rule out that the companion responsible for the QPEs is a star, as has been considered in the very recent works by Linial & Metzger (2023) and Tagawa & Haiman (2023). A typical main sequence star has indeed a radius comparable to the required R_in and, thus, it could (in principle) produce QPEs of a high enough luminosity with the same mechanism outlined in Sect. 2.4. However, assuming R_in to be on the order of the stellar radius would make the initial cloud size independent of v_rel, making it difficult to produce the modulation of the luminosity observed in almost every QPEs source. The initial cloud temperature for stars-disc collisions is expected to be significantly lower than the typical QPE one, but Linial & Metzger (2023) have shown that Comptonisation in the highly ionised medium can actually produce soft X-ray emission. On the other hand, such a scenario poses two additional constraints to consider. First, the pericentre distance of the secondary, constrained from the observed recurrence of QPEs has to be larger than the tidal radius R_t ≈ (M₁/M_⊙)^1/3 R_⊙, namely, the radius at which the star is torn apart by the tidal field of the primary MBH. In principle, if QPEs are only produced by impacts, the star also needs to avoid overfilling its Roche lobe at pericentre – otherwise, any mass transfer will likely contribute to the QPE phenomenon as well. Second, we have to compare the energy of the impact with the star’s binding energy to determine whether the star can actually survive the crossing.

Considering the orbital parameters we obtained (as described in Sect. 3 below and which are basically independent on the nature and mass of the secondary for q ≪ 1), we note that the secondary pericentre distance is well within the tidal radius of a main sequence star for RX J1301.9+2747, and only marginally outside for GSN 069. We further note that in three (GSN 069, eRO-QPE2, and RX J1301.9) out of the four QPE sources considered here, a solar-type star that orbits the primary MBH with the inferred separations would overfill its Roche lobe, therefore transferring mass to the primary. We caution that these estimates depend on the primary MBH masses, which are not well constrained, and on the complex physics driving TDEs and mass transfer; however, we point out that these considerations must be taken into account when discussing the possible nature of the EMRI companion. We have also computed the impact energy for different sources and compared it with a solar mass star binding energy. We find the impact energy to be (at best) comparable to the binding energy. This would suggest that the star is unlikely to survive many passages through the disc without significantly modifying its structure. Linial & Metzger (2023) discarded BH companions based on the influence radius size argument. However, we have noted that the estimate in their Eq. (A.1) considers the orbital velocity instead of the relative velocity with respect to the disc. This is valid only under the assumption that the companion velocity is much larger than the disc velocity and this does indeed lead to an influence radius smaller by about two orders of magnitude, as compared to the case of prograde orbits with modest misalignments that we study here.

In summary, QPEs induced by impacts between the secondary EMRI component and the accretion disc around the primary MBH can (in principle) be produced either by massive stellar BHs (which we know to exist in virtue of LIGO/Virgo findings; The LIGO Scientific Collaboration 2021) in prograde orbits (i.e. with low v_rel), or by stars. We focus here on the former scenario for the reasons given above and since this avoids complications due to the detailed treatment of the star structure (and its evolution due to the impacts) and the possible onset of mass transfer events. As already mentioned, the timing of QPEs is essentially independent from the nature and mass of the companion, so that the orbital parameters derived in Sect. 3 can be considered to be general, while the individual QPE luminosity does actually depend on the secondary nature and on the envisaged detailed emission mechanism.

Finally, another relevant aspect to consider when discussing the nature of the secondary object that triggers the QPE is that when assuming a typical evolved stellar mass function (e.g. Kroupa 2001), the number density of stellar BHs is ∼10⁻³ that of normal stars. However, stellar BHs are heavier than normal stars and segregate closer to the MBH owing to dynamical friction. This mass segregation is generally completed in a fraction of the relaxation time and results in the well-known Bahcall & Wolf (1977) cusp: stellar black holes distribute into a steep cusp with density ∝r^−1.75, while lighter stars follow a milder power-law cusp ∝r^−1.35 (Bortolas 2022). The cusp extends within a radius close to the influence radius of the central MBH. It has been shown that typical nuclear star clusters hosting the relatively low-mass MBHs involved in QPEs feature short relaxation times (shorter than a Gyr, Broggi et al. 2022); thus, they are very likely to be segregated and dynamically relaxed. More specifically, we can safely assume mass segregation and the Bahcall & Wolf (1977) cusp is in place only below ∼1 pc (which is close to the influence radius of the primary MBH). Thus, at a 1 pc scale, we can assume the number density of stellar BHs to be ∼10⁻³ that of stars and extrapolate their number density at smaller scales, considering that stellar BHs and stars have distributions with different profiles (r^−1.75 vs. r^−1.35). It is easy to see that – at radii of the order of the inferred semi-major axis at which QPEs are occurring (i.e. ∼100 R_g ∼ 10⁻⁶ − 10⁻⁵ pc as presented in the next section) – the number density of stellar black holes is comparable to that of stars. This consideration, together with the fact that stars may be tidally destroyed at such small separations or consumed by the repeated interactions with the disc, seems to support the scenario in which a stellar BH is generating the QPEs.

Another interesting feature of QPEs is that, as mentioned in the introduction, their host galaxies seem to overlap with the typical hosts of tidal disruption events. These galaxies feature a rate of TDEs that is ∼10⁻⁴ per year or even greater (e.g. French et al. 2020; Hammerstein et al. 2021). As mentioned above, QPEs may involve the remnant disc of a tidal disruption event. Assuming the nuclear stellar cluster of the host galaxy is mass-segregated and dynamically relaxed, the event rates of EMRIs can be shown to be ∼10⁻² − 10⁻³ times the event rates of tidal disruption events (Bortolas et al., in prep.). Yet EMRIs may take ∼10⁵ years to complete their inspiral, as GWs are initially relatively inefficient at inducing their orbital decay, while tidal disruptions happen promptly and remain observable over timescales of years at most. A quantitative assessment of this statement will be presented in an upcoming study.

3. Results

We used the model developed in Sect. 2 to produce the synthetic light curves of four observed QPE sources: GSN 069, eRO-QPE2, eRO-QPE1, and RX J1301.9+2747. The value of the primary mass M₁ is crucial to correctly match the QPEs timing properties, as it sets the scale of the system, and we constrained it through the M–σ relation (Wevers et al. 2022). The derived M₁ is therefore associated with significant uncertainties. The semi-major axis of the EMRI system was then constrained by the observed QPE recurrence period for each source. We further assumed the secondary to be a M₂ = 100 M_⊙ Schwarzschild black hole for all sources. We note that varying values of M₂ changes the radius of influence R_inf and, therefore, the relative velocity would need to be adjusted to match the QPE peak luminosities. Changing the relative velocity is likely to lead to slightly different orbits and, consequently, timing properties. The orbital eccentricity, the MBH spin, and the disc parameters, as well as the misalignment angles of the EMRI orbit and disc with respect to the plane perpendicular to the MBH spin, do not have stringent observational constraints. Therefore, they must be tuned to reproduce the observed light curves.

Our aim is to qualitatively reproduce the recurrence times and relative X-ray luminosities of QPEs in individual sources and to account for the diversity of QPE properties in the different objects. In principle, our model can be used to infer the system parameters accurately from the observed (often multi-epoch) available light curves and we plan on building a Bayesian inference algorithm to set a meaningful posterior on the model parameters. This will be the subject of a follow-up paper. Here, we do not focus on the individual QPE profile, preferring to defer a detailed study of the involved micro-physics to a future work. For each source, we also provide an estimate of the timescale over which GWs drive the EMRI to merge.

3.1. Synthetic light curves

3.1.1. GSN 069

GSN 069 is the most studied source that has shown QPEs in different epochs. Here, we modelled five peaks observed during an XMM-Newton observation performed on 31 May 2019, described as XMM5 in Miniutti et al. (2023a). The quiescence bolometric luminosity is 2 − 5 × 10⁴³ erg s⁻¹ for this source, implying ṁ₁ ∼ 0.1 for a 10⁶ M_⊙ primary black hole. We assumed the MBH to be slowly spinning with spin χ = 0.1. The secondary BH orbits the primary on a slightly eccentric and inclined orbit (i_EMRI = 10° with respect to the x − y plane), with a semi-major axis of a = 160 R_g, giving an orbital period of 18 h, eccentricity e = 0.1, yielding a GW timescale of T_gw ∼ 2 × 10⁴ yr. We assumed that the disc contains 4 M_⊙, is misaligned by i_disc = 5° with respect to the x − y plane and precesses rigidly with a period of ∼125 days. Based on the detected quiescence level, we assumed ṁ₁ = 0.1. As shown in Fig. 3, our model nicely reproduces not only the observed alternating long-short recurrence times, but also the associated alternating strong-weak QPEs (Miniutti et al. 2023a). Our model also predicts a peak temperature of ∼100 eV, in agreement with the observed values, namely, 90−95 eV (Miniutti et al. 2023a).

Fig. 3. Comparison between observed and predicted light curve. Upper panel: 0.2 − 2 keV quiescence-subtracted X-ray luminosity light curve from the XMM-Newton observation XMM5 of GSN 069 (Miniutti et al. 2023a). Lower panel: synthetic light curve obtained with the parameters given in Sect. 3.1.1.

3.1.2. eRO-QPE2

For this less massive source, the bolometric luminosity of the quiescence is ∼5 − 6 × 10⁴¹ erg s⁻¹, corresponding to ṁ₁ ∼ 0.01 for a M₁ = 10⁵ M_⊙ primary MBH (Arcodia et al. 2021). We assumed the primary to have a spin of χ = 0.5. The secondary BH orbits the primary on an almost circular inclined orbit with a semi-major axis of a = 320 R_g, giving an orbital period of 5 h, eccentricity e = 0.05 (yielding T_gw ∼ 3000 yr) and inclination of i_EMRI = 15° with respect to the x − y plane.

We note that the disc mass in this case is likely too low to be compatible with the tidal disruption of a solar mass star. However, we cannot rule out the possibility that part of the disc mass has already been accreted onto the MBH.

We assumed that the disc contains 0.01 M_⊙ and that is misaligned by i_disc = 10° with respect to the x − y plane, while precessing rigidly on a period of ∼5.7 days. Based on the detected quiescence level, we assumed ṁ₁ = 0.04. As shown in Fig. 4, the observed alternating long-short recurrence times and strong-weak QPE amplitudes are remarkably well reproduced, as is the case for GSN 069. The peak temperatures we obtained, namely ∼170 eV, are slightly below the ∼215 eV inferred from observations, although these are, in fact, subject to large uncertainties (Arcodia et al. 2021).

Fig. 4. Comparison between observed and predicted light curve. Upper panel: 0.2 − 2 keV quiescence-subtracted X-ray luminosity light curve from one of the XMM-Newton observations of eRO-QPE2. Lower panel: synthetic light curve obtained with the parameters listed in Sect. 3.1.2.

3.1.3. eRO-QPE1

The quiescence level was not clearly detected in this source, but it is estimated to be below 1.6 × 10⁴¹ erg s⁻¹, broadly corresponding to ṁ₁ ∼ 0.05 − 0.001 for a primary MBH with M₁ = 10^5.8 M_⊙ (Arcodia et al. 2021). For this source, we assumed the MBH to have a spin of χ = 0.65, while the secondary orbit is taken to be slightly eccentric and inclined, with a semi-major axis of a = 355 R_g, giving an orbital period of 40 h, eccentricity e = 0.05 (with T_gw ∼ 2 × 10⁵ yr) and inclination i_EMRI = 20° with respect to the x − y plane. We assumed that the disc contains 4 M_⊙, is misaligned by i_disc = 5° with respect to the x − y plane and precesses rigidly on a ∼7.5 day timescale. Based on the detected quiescence level, we assumed ṁ₁ = 0.05.

We show the light curve in Fig. 5. We note that there is a hint of a superposed modulation on a ∼150 h period in the data of this QPE source. Having reproduced a similar modulation, we argue that it is linked to the disc precession frequency (∼6 days). We inferred a peak temperature of ∼190 eV, broadly consistent with the observed value of ∼260 eV and again subject to large uncertainties (Arcodia et al. 2021). Here, we note that since our model predicts the QPE occurrence based on the interplay of three different precession frequencies, it could (in principle) reproduce more complex behaviours, such as in the case shown in Arcodia et al. (2022). In particular, those authors found one single QPE peak to be reasonably well fitted by multiple flare profiles (see left panel of their Fig. 2). A more eccentric EMRI could reproduce peaks of different amplitudes much closer in time, according to our model. However, the proper modelisation of the eRO-QPE1 light curve presented in Arcodia et al. (2022) warrants a separate, more thorough investigation.

Fig. 5. Comparison between observed and predicted light curve. Upper panel: 0.2 − 2 keV X-ray luminosity light curve from a NICER ∼250 h-long monitoring of eRO-QPE1. The quiescent level is undetected by NICER. Lower panel: synthetic light curve obtained with the parameters listed in Sect. 3.1.3.

3.1.4. RX J1301.9+2747

RX J1301.9+2747 has a quiescence level similar to that of GSN 069, namely, ∼3 − 6 × 10⁴³ erg s⁻¹, corresponding to ṁ₁ ∼ 0.1 for a primary MBH with M₁ = 10^6.3 M_⊙ (Giustini et al. 2020). We assumed the primary MBH spin to be χ = 0.5. The secondary orbits the primary on an eccentric and inclined orbit with a semi-major axis of a = 50 R_g, giving an orbital period of ∼6 h, eccentricity of e = 0.4 (implying the smallest GW timescale of the four QPEs, i.e. T_gw ∼ 400 yr) and inclination of i_EMRI = 6° with respect to the x − y plane. We assumed that the disc is characterised by 4M_⊙, misaligned by i_disc = 5° with respect to the x − y plane, and that is precesses rigidly with a period of ∼19.8 days. Based on the detected quiescence level, we assumed ṁ₁ = 0.5. The peak temperature our model predicts is in the range of ∼180 − 450 eV. This is up to a factor of 3−4 higher than the observed value, namely, 110−130 eV (Giustini et al. 2020).

Figure 6 shows that the synthetic QPE amplitudes are compatible with the observed ones within a factor of ∼2 and the complex timing properties of the light curve are reasonably reproduced. However, we note that the observed QPE relative amplitude of the peak pairs is not exactly reproduced by our model. These differences, together with those in the peak temperatures, will be better addressed by both a multi-epoch fitting and a further refinement of the emission mechanism modelisation.

Fig. 6. Comparison between observed and predicted light curve. Upper panel: 0.2 − 2 keV quiescence-subtracted X-ray luminosity from one of the XMM-Newton observations of RX J1301.9+2747. Lower panel: synthetic light curve obtained with the parameters listed in Sect. 3.1.4.

3.2. Effect of the model free parameters on the light curves

There are essentially five main parameters that we vary in order to reproduce the synthetic light curves. Here, we describe the effects of each on the timing and burst properties. We show the main effects in Fig. 7 for the case of GSN 069.

Fig. 7. Effect of changing the model parameters on the synthetic light curves for the specific case of GSN 069. We remind that our best-inferred values for the eccentricity of this source is e = 0.1, MBH spin χ = 0.1, and inclination iEMRI = 10° (shown by grey dashed lines). In each panel, we only vary one parameter at a time. Top left panel: MBH spin changed to χ = 0.9. Top right panel: EMRI eccentricity increased to e = 0.6. Bottom left panel: retrograde EMRI with iEMRI = 160° (note the different y-scale). Bottom right panel: disc precession is switched off.

The mass of the primary MBH sets the semi-major axis of the EMRI orbit and determines the size of the accretion disc, therefore (together with the MBH spin) the rigid precession period. The size of the disc is also determined by the primary MBH spin value through the ISCO. In particular the ISCO is smaller, namely, closer to the primary, for higher spin values; therefore, the disc rigidly precesses on shorter timescales. This is because the Lense-Thirring precession frequency at each radius is higher and the average precession frequency increases as the Lense-Thirring precession is stronger at smaller radii. Changing the MBH spin then results in a modulation of the peak amplitude on a timescale that is shorter for higher spins, as shown in the top left panel of Fig. 7. This modulation can also be seen in Fig. 5 and corresponds to the disc precession period.

Another parameter that we varied between different sources is the disc mass. This parameter affects the absolute amplitude of the peaks without influencing either the timing or the relative amplitude between subsequent peaks. We do not explicitly show its effect on the light curve here.

The eccentricity of the EMRI changes substantially the timing properties of the light curve and the amplitude pattern between subsequent peaks. Specifically, when the eccentricity is large the crossings near pericentre are closer in time, generating a light curve characterised by closely spaced double peaks separated by larger time intervals. Furthermore, since the eccentric motion causes the binary to change velocity along its orbit, at each crossing, the velocity of the secondary object is quite different from that of the gas in the disc, assumed to rotate in circular motion at all radii. This means that in general eccentric EMRIs are characterised by higher relative velocities, which correspond to lower influence radii and, consequently, lower luminosities. This is shown in the top right panel of Fig. 7.

Changing the relative inclination between the EMRI companion and the disc modifies not only the peak luminosities, but also the relative amplitude of subsequent peaks. In particular, if the EMRI companion is orbiting in a retrograde orbit with respect to the disc, the luminosities remain much lower than the observed values as the relative velocity between the EMRI companion and the disc at the crossing is significantly higher (by roughly an order of magnitude) than in the prograde case. This is shown in the bottom left panel of Fig. 7.

Finally, we have taken the disc not to be precessing in the case of GSN 069. We find that if the disc position is fixed in time, the peak luminosities are lower by ∼30%, while the relative amplitude of subsequent peaks decreases by about a factor of 2. We show this in the bottom right panel of Fig. 7. The main reason for this change is that the crossing radii are different compared to the precessing case and this essentially translates into different relative velocities.

4. Discussion and conclusions

We have built a complete model to reproduce the behaviour of the sources that produce QPEs confirmed to date: GSN 069, eRO-QPE1, eRO-QPE2, and RX J1301.9+2747. We assumed the QPEs to be triggered when an EMRI companion (mass ratio q ∼ 10⁻⁴ − 10⁻³) crosses the tilted and rigidly precessing accretion disc that surrounds the primary MBH. Since QPEs appear to be connected with candidate TDEs (either partial or total), we modelled the accretion disc following Franchini et al. (2016). Therefore, the disc contains a mass on the order of 1 M_⊙ and it is geometrically thick, allowing it to precess as a solid body around the primary. In our model, the emission comes from a cloud of gas that is pulled out of the disc by the passage of the companion and cools down through adiabatic expansion, emitting ∼10⁴² − 10⁴³ erg s⁻¹ in the soft X-ray band.

We assumed the disc to extend up to ∼6 R_t, where R_t is the tidal disruption radius of the primary MBH. This might be an overestimation of the disc’s outer radius, as the circularisation radius for the debris of a tidally disrupted star is expected to be ∼2 R_t. However, to date, there have been no complete numerical simulations to validate this assumption and the disc is expected to spread viscously. Viscous spreading of the disc may actually explain why QPEs are only observed with a delay ≳4.5 yr with respect to the initial TDE-like X-ray flare in GSN 069 (Miniutti et al. 2023a), as QPE emission only starts when the outer disc reaches the pre-existing EMRI companion orbit.

We have specifically addressed the case of a BH companion, but we cannot exclude that QPEs are produced by star-disc impacts instead (Linial & Metzger 2023; Tagawa & Haiman 2023), although with the caveats and limitations discussed in Sect. 2.5. Once a secondary BH is considered, observable QPEs are only produced in a relatively narrow range of BH masses on the order of ∼100 M_⊙ and only for prograde secondary orbits with relatively small inclination between the orbital and disc planes. Significantly lower mass BHs, retrograde secondary orbits, or a large misalignment between the secondary and disc planes all produce X-ray luminosities that are undetectable against the quiescence emission as the influence radius is R_inf ≪ 10¹¹ cm. This is either because of the smaller BH size or because of higher v_rel (or both). Moreover, the initial temperature of the cloud would be significantly different and the black body emission peak would not be compatible with the observed ∼100 eV. On the other hand, the EMRI lifetime for significantly higher BH masses could be too short to be reasonably accessible to observations. Therefore, the observed QPE sources might well represent the high-luminosity tail of the QPE population with most systems being actually below the current detection threshold or too short-lived to be detected. In this case, the overall rate of EMRIs could be significantly higher than the one that can be estimated from the observed QPE rate, possibly representing a crucial mechanism for the mass growth of MBHs in the mass range 10⁵ − 10⁷ M_⊙. Furthermore, the fact that even a small number of QPE sources can give a lower limit on the EMRI rates constitutes a relevant prediction for LISA (Amaro-Seoane et al. 2023).

We set Δt_QPE to match the observed duration of the QPEs of each individual source. We do not link the duration to the physical properties of the impact. Therefore our model can (in principle) allow for a wider range of Δt_QPE. The constraints on the duration interval would depend on the physical properties of the expanding cloud of gas as well as on the properties of the gas surrounding the accretion disc. We do, however, note that the duration of the QPE is broadly consistent with the expected expansion velocity that we would get from our theoretical model – at least for the cases of GSN 069, eRO-QPE1, and eRO-QPE2. We also note that we fixed the expansion factor to 3 based on the observations of GSN 069 (Miniutti et al. 2023a), as the properties of the gas outside the accretion disc (which determine the timescale for the adiabatic expansion) are unconstrained. However, QPEs become undetectable once the contrast between the cloud temperature and that of the accretion disc emission is ∼1 which implies an expansion factor of 2 − 3 in the other QPE sources as well. Thus, we refrain from discussing the details of the physical mechanism that drives the expansion of the cloud producing QPEs, deferring a physical treatment to a follow-up paper in which both the QPE rise properties (that we neglect here), decay timescales, and physical properties will be studied in detail.

We find that only mutual inclinations (i.e. between the EMRI and the disc) smaller than ∼15° and prograde orbits are capable of reproducing the observed luminosities. This would imply, assuming isotropy, that roughly 2% of the EMRI-disc systems fall into this region of the parameter space. This is clearly a small fraction of the putative number of possible systems, but considering the scarcity of the sources where QPEs have been identified so far, we believe that this constraint can still be reconciled with current observations. However, we can reasonably expect higher inclination systems to progressively align the disc outer parts with the EMRI companion (Miller & Krolik 2013), therefore transiting through the region of inclinations that would give rise to detectable flares and alleviating a possible tension. Finally, since the number of expected EMRI is not very well constrained, we can reasonably expect that for each EMRI detected as QPE there would likely be a population of EMRIs that do not produce QPEs that would be strong enough to be detected.

Another relevant aspect that should be discussed concerns the orientation of the systems considered with respect to the line of sight. Given the primary masses and semi-major axes of the modelled sources, we do not expect timing delays (see Sect. 2.3) to be signficant enough to critically change the picture for certain line of sight directions over others. The same can be said for the lensing generated by the MBHs, which (due to the lower masses and relatively wide orbits) seems to produce quite small deflection angles (∼0.01 − 0.1 radians in all considered cases) in general. Microlensing effects by the orbiting EMRI companions are expected to be even smaller as lensing is only efficient for a very limited range of possible geometrical configurations. We cannot, however, exclude the occasional magnification of one QPE that can occur when the emitting cloud, the MBH, and the observer are close to be perfectly aligned.

The proposed model has five unconstrained parameters: MBH spin value, orbital eccentricity, disc mass, EMRI, and disc inclination. All five may (in principle) be tuned to reproduce the different observed behaviours. For simplicity, we have fixed the values of M₁ and M₂ (and therefore a), but all the eight parameters can (in principle) be determined accurately using multi-epoch observations with several QPEs. We aim to build a Bayesian inference algorithm based on our model in order to place meaningful posterior on these parameters from the observed QPE light curve. This will be the subject of a follow up paper.

Although the emission process we consider is simple and subject to further refinements and modifications, we have shown that our model can reproduce the diversity of QPE recurrence times and behaviour in the four best-studied QPE sources. Moreover, the model is consistent with the observed spectral properties of QPEs and with their evolution during QPE decay and is also able to account naturally for the varying QPE amplitudes and recurrence times in individual erupters.

Acknowledgments

We thank the anonymous referee for the useful comments and suggestions. A.F., E.B. and A.S. acknowledge financial support provided under the European Union’s H2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691). M.B. acknowledges support provided by MUR under grant “PNRR – Missione 4 Istruzione e Ricerca – Componente 2 Dalla Ricerca all’Impresa – Investimento 1.2 Finanziamento di progetti presentati da giovani ricercatori ID:SOE_0163” and by University of Milano-Bicocca under grant “2022-NAZ-0482/B”. R.A. received support for this work by NASA through the NASA Einstein Fellowship grant No. HF2-51499 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. E.B. acknowledges support from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158). M.G. is supported by the “Programa de Atracción de Talento” of the Comunidad de Madrid, grant number 2018-T1/TIC-11733. G.M. was partly supported by grant n. PID2020-115325GB-C31 funded by MCIN/AEI/10.13039/501100011033.

References

All Tables

All Figures

Fig. 1. Schematic view of the model where the central big black dot represents the primary MBH, the small black dot represents the EMRI companion, the grey area represents the tilted disc and the EMRI-disc crossings are shown with yellow stars. The vertical arrow shows the MBH spin vector. The misalignment between the EMRI orbit and the disc is enhanced in this cartoon with the purpose of showing the impacts more clearly. However, we use smaller misalignments to produce the QPEs light curves. The blue and yellow arrows show the EMRI and gas velocity vectors, respectively, while the yellow dashed line only mimics an orbit inside the disc.

In the text

	Fig. 2. Spectrum emitted by a cloud of initial size Rinf and initial temperature T = 106 K that expands by a factor of 3 in size over one hour.
In the text

	Fig. 3. Comparison between observed and predicted light curve. Upper panel: 0.2 − 2 keV quiescence-subtracted X-ray luminosity light curve from the XMM-Newton observation XMM5 of GSN 069 (Miniutti et al. 2023a). Lower panel: synthetic light curve obtained with the parameters given in Sect. 3.1.1.
In the text

	Fig. 4. Comparison between observed and predicted light curve. Upper panel: 0.2 − 2 keV quiescence-subtracted X-ray luminosity light curve from one of the XMM-Newton observations of eRO-QPE2. Lower panel: synthetic light curve obtained with the parameters listed in Sect. 3.1.2.
In the text

	Fig. 5. Comparison between observed and predicted light curve. Upper panel: 0.2 − 2 keV X-ray luminosity light curve from a NICER ∼250 h-long monitoring of eRO-QPE1. The quiescent level is undetected by NICER. Lower panel: synthetic light curve obtained with the parameters listed in Sect. 3.1.3.
In the text

	Fig. 6. Comparison between observed and predicted light curve. Upper panel: 0.2 − 2 keV quiescence-subtracted X-ray luminosity from one of the XMM-Newton observations of RX J1301.9+2747. Lower panel: synthetic light curve obtained with the parameters listed in Sect. 3.1.4.
In the text

Fig. 7. Effect of changing the model parameters on the synthetic light curves for the specific case of GSN 069. We remind that our best-inferred values for the eccentricity of this source is e = 0.1, MBH spin χ = 0.1, and inclination iEMRI = 10° (shown by grey dashed lines). In each panel, we only vary one parameter at a time. Top left panel: MBH spin changed to χ = 0.9. Top right panel: EMRI eccentricity increased to e = 0.6. Bottom left panel: retrograde EMRI with iEMRI = 160° (note the different y-scale). Bottom right panel: disc precession is switched off.

In the text