Probing Accretion Physics with Gravitational Waves

Gravitational-wave observations of extreme mass ratio inspirals (EMRIs) offer the opportunity to probe the environments of active galactic nuclei (AGN) through the torques that accretion disks induce on the binary. Within a Bayesian framework, we study how well such environmental effects can be measured using gravitational wave observations from the Laser Interferometer Space Antenna (LISA). We focus on the torque induced by planetary-type migration on quasicircular inspirals, and use different prescriptions for geometrically thin and radiatively efficient disks. We find that LISA could detect migration for a wide range of disk viscosities and accretion rates, for both $\alpha$ and $\beta$ disk prescriptions. For a typical EMRI with masses $50M_\odot+10^6M_\odot$, we find that LISA could distinguish between migration in $\alpha$ and $\beta$ disks and measure the torque amplitude with $\sim 20\%$ relative precision. Provided an accurate torque model, we also show how to turn gravitational-wave measurements of the torque into constraints on the disk properties. Furthermore, we show that, if an electromagnetic counterpart is identified, the multimessenger observations of the AGN EMRI system will yield direct measurements of the disk viscosity. Finally, we investigate the impact of neglecting environmental effects in the analysis of the gravitational-wave signal, finding 3$\sigma$ biases in the primary mass and spin, and showing that ignoring such effects can lead to false detection of a deviation from general relativity. This work demonstrates the scientific potential of gravitational observations as probes of accretion-disk physics, accessible so far through electromagnetic observations only.

Gravitational-wave observations of extreme mass ratio inspirals (EMRIs) offer the opportunity to probe the environments of active galactic nuclei (AGN) through the torques that accretion disks induce on the binary. Within a Bayesian framework, we study how well such environmental effects can be measured using gravitational wave observations from the Laser Interferometer Space Antenna (LISA). We focus on the torque induced by planetary-type migration on quasicircular inspirals, and use different prescriptions for geometrically thin and radiatively efficient disks. We find that LISA could detect migration for a wide range of disk viscosities and accretion rates, for both α and β disk prescriptions. For a typical EMRI with masses 50M⊙ + 10 6 M⊙, we find that LISA could distinguish between migration in α and β disks and measure the torque amplitude with ∼ 20% relative precision. Provided an accurate torque model, we also show how to turn gravitational-wave measurements of the torque into constraints on the disk properties. Furthermore, we show that, if an electromagnetic counterpart is identified, the multimessenger observations of the AGN EMRI system will yield direct measurements of the disk viscosity. Finally, we investigate the impact of neglecting environmental effects in the analysis of the gravitational-wave signal, finding 3σ biases in the primary mass and spin, and showing that ignoring such effects can lead to false detection of a deviation from general relativity. This work demonstrates the scientific potential of gravitational observations as probes of accretion-disk physics, accessible so far through electromagnetic observations only.

I. INTRODUCTION
Extreme mass ratio inspirals (EMRIs) are a primary target of the future space-borne gravitational-wave (GW) detector LISA [1]. A typical EMRI involves a stellar-mass compact object with mass O(10)M ⊙ orbiting a massive black hole (MBH) with mass O(10 6 )M ⊙ . While emitting in the LISA band, the compact object will complete as many as hundreds of thousands of orbits around the MBH, making EMRIs ideal sources to test the nature of black holes (BHs), general relativity (GR), and the astrophysics of galactic nuclei [2][3][4][5][6].
Mergers involving MBHs commonly occur within gasrich environments. Interactions with the gas can help binaries -including MBH binaries and EMRIs -form and become more compact. Based solely on observations of actively accreting MBHs, it is expected that 1%-10% of the EMRIs observed by LISA will reside in the accretion disk of active galactic nuclei (AGNs) [7]. However, the formation of EMRIs involves a wide range of dynamical processes and timescales, so that predictions for LISA in the literature still vary by orders of magnitude, even in the absence of gas [4][5][6]8]. Recently, Refs. [9,10] argued that there could be even more EMRIs taking place in AGNs than previously estimated, and that LISA could detect 10 − 10 4 EMRIs per year from dense accretion * lorenzo.speri@aei.mpg.de disks, compared to only 1 − 10 2 per year from relatively gas-free environments (see also Ref. [11]).
The presence of an accretion disk can considerably modify the orbital trajectory of an EMRI emitting in the LISA band. The modification comes in the form of torques, originating either from hydrodynamical effects such as "accretion winds" [12,13], or through purely gravitational torques [14] such as dynamical friction [13,[15][16][17], which can take the form of planetary-style migration if the density wakes produce torques through shocks [12,18,19]. Through these processes, the disk exchanges energy and angular momentum with the system.
It is expected that the parameters characterizing EMRI signals will be measured with extreme precision by LISA, thanks to the large number of orbits that can be observed in the LISA band. For this reason, it is reasonable to believe that accretion-disk effects which are estimated to be detectable [12,13,[15][16][17][19][20][21] can also be used to measure accretion-disk properties, such as the accretion rate and disk viscosity. This would open up the possibility of studying the properties of accretion disks through purely gravitational means.
In this work, we perform the first quantitative study of the measurability of accretion-disk effects using EMRI observations. In particular, we model migration in the disk with two distinct analytical prescriptions -the α [22] and β [23] -to account for expected uncertainties in the underlying disk physics. Our analysis uses stateof-the-art waveforms [24][25][26][27] and is performed within a realistic framework for EMRI parameter estimation. For the first time, we show which accretion-disk quantities can be effectively constrained via gravitational waves. To achieve this, we perform a full Bayesian inference study over the parameter space of a circular-equatorial EMRI system. This setup is the most general in the scenario of interest: while generic EMRI trajectories can be eccentric and inclined, EMRIs formed in accretion disks are likely to have circularized and aligned with (or be born in) the disk [9,11,[28][29][30][31]. We show that, with an agnostic torque model employed for parameter estimation, LISA could detect and characterize migration in the EMRI GW signal.
We find that migration could be detected in both α and β disks. Our results validate earlier studies [12,19] within a fully Bayesian setting, and extend the discussion to the measurability of the disk parameters. Compared to Refs. [12,19], we find that migration can be detected and also characterized when assuming the α disk model with realistic disk parameters. Furthermore, for our reference EMRI system, we find that the GW observation alone can distinguish between disk prescriptions and constrain regions of the disk parameter space. Assuming a torque model, we show that the GW measurement can be combined with electromagnetic observations to infer the viscosity of the host disk. Finally, we investigate the degree to which we expect key EMRI parameters to be biased if one ignores torques in EMRI parameter inference, finding that the primary's mass and spin may be biased. We also demonstrate that not accounting for such torques can lead to the false detection of a GR deviation.
While our work focuses on migration in thin accretion disks, our torque parametrization could be used to study other effects from the environment or modifications of GR, provided that the induced loss of angular momentum or energy can be written as a simple power law of the orbital separation. This paper is organized as follows. In Sec. II we present the models for the accretion-disk structure and the migration torque used throughout our work, including a very general phenomenological torque model suited for GW parameter estimation. In Sec. III we describe the vacuum waveform model and how we modify it to account for an environmental torque. We also summarize our methods for parameter estimation and the properties of our reference EMRI. Finally, in Sec. IV we present the results: the projected constraints on the amplitude of environmental torques (Sec. IV A), a detailed study of an EMRI signal with a detectable disk torque (including a discussion of multimessenger implications, Sec. IV B), and a study of the bias induced on the EMRI parameters and beyond GR deviations when environmental effects are present in the signal but ignored in the GW analysis (Sec. IV C). We discuss future prospects in Sec. V. We work in units in which G = c = 1, unless explicitly stated.

II. ACCRETION-DISK EFFECTS
We begin by describing the accretion-disk prescriptions, torque model and assumptions used in this work. We denote the primary and secondary masses with M 1 and M 2 , respectively. We assume that the MBH is accreting at a ratio f Edd ≡Ṁ 1 /Ṁ Edd,1 = ϵṀ 1 /L Edd,1 of the Eddington rateṀ Edd,i = 2.536 × 10 −8 M i (ϵ/0.1) yr −1 .
Here ϵ denotes the efficiency of conversion of mass energy into luminosity in the disk. The accretion efficiency and the Eddington ratio enter all disk quantities in the combination ϵ −1 f Edd .
The astrophysics of accretion disks is notoriously uncertain. Sophisticated numerical simulations are often required to produce realistic disks [32], and to describe the rich phenomenology of the torques that disks generate [20]. However, a fully numerical approach is intractable for this work's scope, which requires models that are fast to generate. Thus, we adopt analytical models for the disk and its torques. We expect future analyses of real data to be conducted in conjunction with numerical simulations [20,33].
We employ radiatively efficient, Newtonian, stationary, thin accretion-disk models, considering both the α and β disk prescriptions for the standard Shakura-Sunyaev solutions for the inner (radiation-pressure dominated) region of the disk [22,34]. For α disks the viscous stress is proportional to the total pressure t rϕ = α(p gas +p rad ) [22], while for β disks t rϕ = αp gas [23]. There is a long-standing question regarding the stability and realism of these analytic solutions [32,[35][36][37][38]. In particular, β disks have raised concerns for being only loosely motivated by the absence of thermal instabilities, which appear in analytical solutions of α disks. Despite these instabilities, the α disk model is still considered a good approximation for realistic, turbulent accretion flows in the radiation-dominated regime (although still far from reproducing all the features of radiation magnetohydrodynamic simulations, see Ref. [32]). In both cases α is the viscosity, which parametrizes the complex (and uncertain) magnetohydrodynamic features of accretion disks. Viscosity in AGN disks is typically believed to be around α = 0.1, and as low as α = 0.01 [32,39]. However, smaller values of α have not been excluded yet.
The Shakura-Sunyaev disk with α and β viscosity prescriptions predicts the following density profiles and scale height [22,23]  A derivation of these quantities from the properties of the disk can be found in our Appendix A, and with more details in [40] and [12]. The corresponding disk densities are obtained as ρ = Σ/2H, and are inversely proportional to the viscosity and the accretion rate. The disk models that we employ are valid for geometrically thin disks only [41], where H ≪ r. This condition is satisfied for the entire EMRI evolution and all disk parameters explored in this work. Observations of the inner regions of AGN accretion disks around 10 6 M ⊙ black holes have found electron densities log n e [cm −3 ] ≃ 16 − 18.6 [42]. In our disk models, these correspond to Σ = n e m p × 2H ≃ 2 × 10 5 kg/m 2 − 3 × 10 7 kg/m 2 . The presence of an accretion disk induces environmental torques that modify the trajectory of the compact object orbiting the MBH. For EMRI sources in the LISA band, such torques are suppressed with respect to the effect of GW emission, but are potentially observable [12,16,17,19,20,33]. Because of this suppression, we can treat environmental torques as a perturbative effect, with the total torque acting on the secondary given by the standard, gravitational torque (see Sec. III), plus the environmental one:L =L GW +L disk .
Environmental torques come in a variety of forms, with different dependencies on the disk parameters and the location of the EMRI secondary in the disk. Moreover, analytic models for the torques (and the underlying disks, as discussed above) still suffer from systematic uncertainties. In order to search for environmental effects in the GW signal, we need a flexible, agnostic model. Motivated by torque models in the literature [12,16,19], we propose the following parametrization: where we scale the torque byL (0) GW = − 32 5 M 2 /M 1 (r/M 1 ) −7/2 , the leading order circularorbit GW torque, and where r is the distance between the compact object and the MBH. In this model, the torque is simply characterized by an amplitude A and a radial slope n r , which corresponds to a post-Newtonian (PN) order −n r PN. We further justify the torque model (2.2) in the following.
The approximation Eq. (2.2) breaks down within the inner edge of the disk, where accretion-disk torques should wane due to the low gas density. Modeling this phase is however not necessary: our waveform model also ignores the plunge of the EMRI system, because its duration is only a small fraction of the full inspiral [43][44][45].
Environmental effects in AGN disks may include the disk gravitational potential, the formation of density waves (dynamical friction or planetlike migration, depending on the scales involved), winds, and mass accretion on the primary or the secondary [12,13,[15][16][17]19]. In this work, we focus on migration because, in preliminary explorations of the dephasing of the EMRI inspiral induced by the aforementioned effects, we find migration to be the dominant one, in agreement with previous works [12,16,17,19].
In analogy with planet-planetary disk interactions [46,47], a compact object spiralling through a disk can produce spiral density waves in the surrounding gas, which can, in turn, excite resonances in the disk. Some waves move inward of the secondary's orbit, and some outward. The latter extract angular momentum from the orbit, which shrinks in size. The secondary is said to have undergone "type I" migration in this case. The EMRI secondary can also undergo "type II" migration, if a gap opens in the disk at the radius of the secondary's orbit.
In type-I migration, it is sufficient to approximate the disk response with linear perturbation theory [48]. We use the analytic formula for type-I migration for a Newtonian orbiter in a three-dimensional, isothermal disk from Ref. [49], where η = M 2 /M 1 is the mass ratio and γ = 3/2 (α disk) or γ = −3/5 (β disk). This takes precisely the form Eq. (2.2), with amplitude and powers given in Table I. This torque is negative (i.e., the migration is inward) and decreasing along the inspiral, as can be seen substituting the expressions for the surface density and disk height Eq. (2.1). Being proportional to the disk density and inversely proportional to the disk height Eq. (2.1), the migration torque Eq. (2.3) scales with inverse powers of the disk parameters α and f Edd . This formula was derived in the context of planetary disks, where it successfully reproduces numerical studies [49][50][51] (see also the review [52]). Unfortunately, numerical simulations have not yet reproduced the realistic conditions of an EMRI migrating in a nuclear accretion disk, such as relativistic inspirals and turbulent, radiating, nonisothermal disks. However, hydrodynamical simulation for higher mass-ratio systems (intermediate mass ratio inspirals) on a post-Newtonian trajectory [20,33] found that the analytic formula of Ref. [49] successfully reproduces the torque amplitude (within an order of magnitude) and its dependence on the radius.
The type-I torque described above does not take into account the effect of gas on horseshoe or corotation orbits near the inspiraling body. This is sometimes referred to as "type-III" migration, and has been shown to lead to rapid runaway migration in the context of protoplanetary dynamics [53]. This corotation torque should only be significant when the mass in the disk near the EMRI orbit is comparable or greater than the secondary mass. In our system, the local disk mass is orders of magnitude smaller than the secondary mass while the system emits in the LISA band [12]. If our EMRI enters the type-III regime, our type-I model is likely be a conservative choice for the amplitude of the torque.
Because of the turbulent nature of AGN disks and the rapid evolution of gas around the inspiraling body, migration can oscillate over time and average to both negative and positive values [20]. Here we use the positive, constant overall torque amplitude predicted in Ref. [49]. Oscillations around the average migration torque were found to be suppressed for smaller mass ratios [20] and could be modeled separately; see e.g., Ref. [21].
Type-II migration takes place when the orbiter's tidal torques remove gas from the orbit faster than viscous forces refill it, creating a lower-density annular gap in the disk. The conditions for gap opening in the two disk models we consider are given in Ref. [12], Eqs. (42)- (46). While migration torques are of type-I in most of the parameter space relevant for EMRI GW observations, the type-II regime may also occur for some disk-parameter and EMRI radii. In this regime, previous work on EMRIs [12,16] used the quasistationary approximation of Syer and Clarke [54]. However, contrary to the standard paradigm, recent simulations have shown that gas can flow through the gap and that long timescales might be required for the torque to reach the quasistationary estimate [55][56][57][58]. To avoid overestimating the torque and introducing additional model dependence in our analysis when entering the type-II regime, we use the type-I formula whenever migration is active, as suggested by recent numerical simulations of migration in intermediate mass-ratio inspirals [20,33].
Migration can also be quenched at certain radii in the disk (see condition (113) in [12]), a transition which could not be described by Eq. (2.2). We verify that the EMRI does not cross the quenching radius in the LISA band across the parameter space explored in this work.
To summarize, in the following we use the type-I migration torque of Table I to simulate the impact of the disk on the compact object, and we recover it in terms of the agnostic parametrization of Eq. (2.2). This formulation is crucial for three reasons. Firstly, it reduces the parameter space dimension and allows us to perform parameter estimation efficiently. It compresses the information encoded in f Edd , α, ϵ, C, n r , n α , n f Edd , n M1 into the amplitude A and slope n r , which characterize the size of the effect and the disk type, respectively. Secondly, this parametrization allows us to extract the astrophysical parameters a posteriori. Once the disk type is determined through the slope n r , we can infer the astrophysical parameters from the measured amplitude A → (f Edd , α, ϵ) and Table I. Furthermore, if a more accurate model for the mapping A → (f Edd , α, ϵ) is available, we can just update the mapping without repeating the slow parameter estimation procedure. Lastly, the constraints obtained with this parametrization can also be used to test general relativity. In fact, for some fixed value of n r (or −n r PN order) the amplitude A can represent the size of a GR deviation through the parametrized post-Einsteinian expansion [59] commonly used to test general relativity with gravitational waves [60].

III. EMRI WAVEFORMS
The GW emission from EMRI systems will only provide measurements of their parameters with unprecedented precision [61][62][63][64] if their waveforms are accurately modeled. EMRI waveform models rely on perturbative solutions in which the Einstein equations are expanded about the limit of small mass ratio η = M 2 /M 1 ≪ 1 [65]. In this limit, the orbital evolution of the compact object in vacuum is governed by the Kerr geodesic equations with a forcing term, called the gravitational self-force [66,67]. The self-force takes into account the finite size and mass of the body and its backreaction on the background Kerr spacetime. As a result, the orbit of the compact object slowly decays into the MBH due to the emission of GWs. The presence of an accretion disk can further enhance the decay of the compact object's orbit.
In this work we employ the FastEMRIWaveform (few) waveform model package [24][25][26][27] and modify the trajectory evolution to take into account the presence of environmental effects. We model the evolution of Kerr circular-equatorial orbits at adiabatic order [68] by interpolating the Teukolsky fluxesL GW using the Teukolsky package in the Black Hole Perturbation Toolkit [69]. We model the disk-induced effects by writing the rate of change of angular momentum aṡ L =L GW +L disk , with the environmental contribution from Eq. (2.2). This implementation provides a fast and accurate adiabatic trajectory that can be fed into a waveform-generation formalism through the augmented analytic kludges [70].
We assume that the disk rotation axis is aligned with the spin of the primary BH. While this might not be the case for the outer regions of a disk, we expect the inner regions to be aligned by the Bardeen-Petterson effect [71][72][73]. We also restrict our analysis to circular and planar orbits (i.e., in the plane of the disk), as disk-induced density waves are expected to damp the EMRI inclination and eccentricity long before it enters the LISA band [9,[28][29][30][31]. Compact bodies are also likely to form preferentially within the disk, giving rise to planar EMRIs [11]. We conservatively assume prograde orbits, as retrograde orbits can suffer from even larger disk-induced torques [13,74,75]. Compact objects formed in the disk are expected to be on prograde orbits [76,77], and prograde EMRIs can be seen to much greater distances, so we expect these to dominate detected LISA events [64]. However, our implementation is flexible enough to model generic orbits once prescriptions for disk-induced effects become available for this scenario.
A quasicircular equatorial EMRI waveform is described by the masses of the two bodies M 1 , M 2 ; the dimensionless spin parameter a of the primary; the initial phase and radius, Φ 0 and r 0 ; the azimuthal and polar angles, θ K and ϕ K , describing the orientation of the spinangular momentum vector of the MBH; the polar and azimuthal sky location angles, θ S and ϕ S , given in the Solar System barycenter frame; and the luminosity distance d L . The presence of accretion effects introduces two additional parameters, A and n r . We will refer to (θ S , ϕ S , θ K , ϕ K , d L ) as extrinsic parameters, and to (M 1 , M 2 , a, r 0 , Φ 0 ) as intrinsic parameters.
The speed of generation of few waveforms allows us to perform Bayesian analyses with Markov chain Monte Carlo (MCMC) methods. In this work, we use our version of parallel-tempered MCMC, based on both emcee [78] and ptemcee [79] with five adaptive temperatures and 16 walkers per temperature. The sampler was run until chains were longer than 50τ , whereτ is the average autocorrelation time, determined across chains [80]. This ensures the convergence of the samples. We sample using a standard Gaussian likelihood Here we assume stationary Gaussian noise, we define the data stream s, GW signal h(λ) with parameters λ, and the inner product weighted by the LISA power spectral density (PSD), S n (f ) [81,82] and tilde means Fourier transform. In all the studies reported, we inject noise-free data streams since this yields results averaged over noise realizations. A noise injection would not affect the shape of the posteriors, but only shift the recovered parameters away from the injected ones. Since our conclusions are based on the measurement precision obtained by the width of the posteriors, this choice does not affect our results. We assume uniform priors in all parameters, restricted to a narrow range around the true parameters. All posteriors shown in this work have a support that is much tighter than the prior ranges.
Because of the high computational cost of the analysis we focus on a reference EMRI with masses M 1 = 10 6 M ⊙ and M 2 = 50M ⊙ , and primary spin a = 0.9. Black holes of 50M ⊙ or more have already been observed in binaries by LIGO and Virgo [83] Collaborations, although with a lower rate compared to lighter BHs. Our choice of a relatively heavy secondary is motivated by the fact that black holes are expected to grow via accretion when originating in an AGN disk [11,84] (for possible formation scenarios see Ref. [85,86]). Depending on the accretion rate and the time available for growth, a much larger secondary mass could be reached, in principle. However, we decided conservatively to set M 2 = 50M ⊙ to be within the range M 2 ∈ [10, 60]M ⊙ reported in the LISA Science Requirements Document [87], and because a larger secondary mass would only improve our results. Since the SNR grows linearly with M 2 , a larger secondary mass would provide tighter constraints on environmental effects. Furthermore, due to the negative-PN nature of the disk torques, for a fixed four years inspiral the compact object could explore farther regions from the central MBH, yielding tighter constraints on environmental effects.
We set the EMRI initial radius at the beginning of the LISA observation to r 0 = 15.482 M 1 (= 7.45 × 10 −7 pc), such that the compact object plunges into the MBH after 4 years. We assume that LISA will observe the signal in its entirety, consistently with the nominal science mission duration [87] and neglecting gaps in the detector's data. We fix the other parameters to randomly chosen values of (Φ 0 , θ S , ϕ S , θ K , ϕ K ) = (3.0, 0.54, 5.36, 1.73, 3.2) and a luminosity distance d L = 1.456 Gpc [88], chosen to give a signal-to-noise ratio SNR= ⟨h|h⟩=50. A detailed investigation of how different configurations affect the detectability of environmental effects with EMRIs is beyond the scope of this work. We stress that this choice of parameters is typical for observable EMRIs and well within the range provided in the LISA Science Requirements Document [87], where a "Golden" EMRI system has SNR> 50, spin a > 0.9, and a prograde orbit.

IV. EMRI PARAMETER ESTIMATION WITH ENVIRONMENTAL TORQUES
There are three ways in which environmental effects could be relevant to GW observations of EMRIs: 1. If environmental effects are absent or too weak to affect the waveform, EMRI observations could be used to set an upper limit on the torque amplitudes. The same analysis can be used to forecast the detectability of a given effect, and it only requires knowledge of the torque's radial dependence (Sec. IV A); 2. If environmental effects are strong enough and can be modeled with a simple power of the radius, their presence can be detected (Sec. IV B). Provided reliable torque models, we can use such detections to constrain some of the properties of the disk, especially in coordination with electromagnetic observations; 3. If we analyze GW data ignoring environmental effects, and these are relatively strong, EMRI parameter estimation will be biased (Sec. IV C).
We discuss each of these scenarios in the following sections.  Table I for the dependence on accretion rate f Edd and viscosity α). We consider a torque detectable if its amplitude is above the 95% percentile line (dashed black line). Assuming an accretion-disk model, this implies a bound on combinations of α and f Edd (for ϵ = 0.1). We also tick the amplitudes at which the phase difference between a disk-embedded (A ̸ = 0) and vacuum system (A = 0) reaches ∆Φ = 1, 100 at the end of the inspiral. The EMRI is observed by LISA 4 years before plunge with SNR= 50, and has central black hole spin a = 0.9 and masses M1 = 10 6 M⊙, M2 = 50M⊙.
Assuming the disk models and migration torque from Sec. II, we ask, when are environmental effects large enough to be detectable?
The dephasing ∆Φ between vacuum and matterinfluenced waveforms is a commonly used metric for detectability, with ∆Φ ≳ 1 a reference threshold. However, a significant dephasing is a necessary, but not sufficient condition for detectability. Using this threshold as a proxy can lead to overestimating detectability, implying that it must be used as a qualitative, but not quantitative metric. To better assess the detectability of environmental effects, we perform a Bayesian parameter estimation over the EMRI parameters and the torque amplitude for a vacuum signal. This approach is similar to tests of GR with parametrized post-Einsteinian expansions of the phase [89]. We sample over the intrinsic EMRI parameters as well as the torque amplitude, fixing the slope n r and the extrinsic parameters for computational efficiency. Sampling over the extrinsic parameters does not affect our results, because these are not strongly correlated with the amplitude A (see Sec. IV B).
The resulting posterior for A has a variance that carries information about the detectability of the effect (slope) selected. In this work, we assume that any effect with an amplitude falling inside the 95% (symmetric) bound is 2σ-consistent with noise, under our assumed (flat) prior. Therefore any effect with A < A 95% is considered indistinguishable from an EMRI in vacuum, and viceversa. In reality, the indistinguishability of the two hypotheses should be assessed by calculating the Bayes factors between the disk and vacuum (possibly eccentric and inclined) hypotheses when injecting different amplitudes. This is a computationally expensive analysis which is beyond the scope of this work.
We show the results obtained for our reference EMRI in Fig. 1 for the two slopes predicted for the migration torque by the two disk models (Table I). We find that the symmetric 95% bounds for β and α disks are A 95% = {3.5, 1.8} × 10 −6 . The constraint on A becomes tighter as the slope increases, consistent with the fact that larger slopes correspond to more negative PN orders probed by the long inspiral.
In Fig. 1, we also indicate at what amplitude the torque induces a dephasing of 1 and 10 2 rad. This shows torques may not be detected even for dephasings much larger than a radian, confirming that the previously adopted requirement ∆Φ ≳ 1 overestimates the detectability of environmental effects.
A more accurate estimate of detectability consists in requiring A > A 95% . For migration in β disks, we find that this implies For typical parameters α = f Edd = ϵ = 0.1 the effect is detectable, as found in [19]. The dependence on the disk accretion rate and viscosity implies that lower values would lead to more observable effects. Referring to the same typical parameters, Ref. [12] found that migration in the more realistic α disks could not be observed. In this case, we find the condition This implies that, while α = f Edd = ϵ = 0.1 are still excluded, the effect could be detected in disks with lower accretion rates and/or viscosity. For instance, for values α ≲ 0.05 the effect is measurable for f Edd ≲ 0.01, ϵ = 0.1. These values are within what is expected from global simulations [32] and x-ray observations of AGNs [90]. The conditions (4.2) and (4.1) are based on the "typical" parameters for EMRI in accretion disks used in our signal injection. These detection thresholds would degrade, for instance, if LISA observed a shorter portion of the inspiral.
While here we present the bounds for the radial slopes predicted for the migration torque in two disk models, this analysis can be easily generalized to other effects at different slopes, or PN orders. In Fig. 6 of Appendix B, we present the symmetric 95% bound on the torque amplitude as a function of the slope for −4 ≤ n r ≤ 10. . Migration in β disks predicts a slope nr = 59/10 (green dash-dotted line), which is excluded at more than 95% (black dashed lines). The expected degeneracy between the amplitude and slope parameters is shown as a blue dotted line. The EMRI configuration is the same as in Fig. 1 and the complete corner plot can be found in Appendix C in Fig. 7. The contours show the 1σ,2σ,3σ Gaussian credible regions.
We find that for n r ≳ 3 the upper limit follows the relation log 10 A 95% = −4.63 − 0.14 n r . This can be readily used as an approximate bound for other effects.

B. Detection of environmental torques
We now consider the more optimistic scenario in which disk effects are above the detection threshold. Can we measure and characterize environmental torques? Provided a reliable model, can we infer the parameters of the disk hosting the EMRI?
We again investigate this scenario by performing a Bayesian analysis for our reference EMRI. We limit our study to the more realistic α disks, injecting a signal with slope n r = 8 and amplitude A = 1.92 × 10 −5 . Choosing α = 0.03 as in Ref. [33], this amplitude corresponds to a value f Edd = 0.005 (with ϵ = 0.1), and to a surface density Σ α ≈ 3.6 × 10 5 kg/m 2 at r = 10 M 1 , consistent with observations [42]. However, our agnostic procedure means the constraint applies to all combinations of α, f Edd and ϵ resulting in the same amplitude through Eq. (2.2).
The full posterior can be found in Appendix C in Fig. 7. As expected for typical EMRI observations, the intrinsic parameters are measured with ∼ 10 −5 relative precision. The sky localization error is ∆Ω = 1.8 deg 2 [91,92] and the relative luminosity distance error is 6% [93]. The comoving volume error for this source is ≈ 5 × 10 −5 Gpc 3 , which means that this source would be promising for follow-up electromagnetic campaigns.
The marginalized posteriors of the environmental parameters A and n r are shown in Fig. 2. The posterior of A is inconsistent with A = 0 at more than 3σ, as expected from our previous discussion. Moreover, the measurement of n r can be used as a discriminator between disk (or torque) models. In our example, we can distinguish that the injection is due to migration in α and not β disks, since the latter predicts a slope n r = 59/10 (dash-dotted green line) that is excluded by the posterior at more than 95%.
We find a strong degeneracy between the amplitude A and slope n r , with correlation coefficient −0.97. This means that the waveform template does not change significantly if we vary the amplitude A and slope n r at the same time along such a correlation. The reason for the degeneracy is that the environmental torque has the largest impact on the waveform at the beginning of the inspiral, namely when the radius roughly coincides with the initial separation r ≈ r 0 . As a result, the torque is well approximated by a constantL disk (r 0 ) ≈ A r nr 0 . This gives rise to the observed correlation, as shown in Fig 2 (blue dotted line). This degeneracy is characteristic of this model and, therefore, could be used as a model consistency check.
The correlation coefficients of the amplitude and slope with the intrinsic parameters are of the order ∼ 0.6 − 0.7, whereas with the sky localization and distance are of order ∼ 0.2 and 0.08, respectively. As expected, extrinsic parameters are not strongly correlated with the accretion parameters.
Our analysis shows that a sufficiently strong disk torque can be detected with a model-independent template, as long as the torque can be described as a power law of the radius. Provided that we have a trusted model for the torque slope and amplitude, the latter carries further information about the effect, which can be used in conjunction with measurements of the mass to extract the values of α and f Edd consistent with the observation. This is represented in the top left panel of Fig. 3  We also show the undetectable region according to the criterion in Eq. (4.2) (black-shaded region). We find that purely GW observations can single out a narrow region of the parameter space (α, f Edd ), or equivalently (Σ α , f Edd ). Note that this parameter space region is determined through the aforementioned relation. Therefore, if there was an analytic prescription defining the validity region of Eq. (2.2), it would be possible to further constrain the parameter space of (α, f Edd ). Extremely low (large) values of α (Σ α ) should not be misinterpreted as actual possible constraints, but only as a consequence of the degenerate measurement of α, f Edd (Σ α , f Edd ) consistent with the amplitude torque A.
In order to break the degeneracy and fully characterize the disk, we need electromagnetic observations. If the host AGN for this EMRI is identified in a follow-up campaign, electromagnetic observations across the spectrum could be used to determine the bolometric luminosity of the central engine. The bolometric luminosity, together with the GW measurement of the central BH mass and spin, could be used to determine the underlying efficiency and accretion rate of the AGN. Finally, this multimessenger determination of the accretion rate f Edd could be used to extract the disk viscosity α from the joint posterior provided by the GW analysis. We give a concrete example of this procedure for f Edd = 0.01 and α = 0.00375 consistent with amplitude A = 1.92 × 10 −5 and a surface density Σ α ≈ 1.44 × 10 6 kg/m 2 at r = 10 M 1 , well within observational limits [42]. We use the numerical fits provided in [94] to relate the (observable) bolometric luminosity to the intrinsic accretion rate and radiative efficiency [95]. In the lower panels of Fig. 3, we show the constraints obtained on the viscosity (or on the disk surface density) when the AGN accretion rate is inferred with error σ f Edd [96]. This example showcases the potential for multimessenger observations of EMRIs in accretion disks.
Electromagnetic detection and host association will be somewhat challenging for our reference EMRI. In our error volume, there could be up to O(10 − 100) black holes with a mass measurement from EM observations consistent with M 1 within 0.5 dex (as estimated from the black hole mass function [97]), too many for unique identification of the source. However, other properties of the source (e.g., the inclination of the MBH spin) could be used to simplify the identification problem. Accounting for band-dependent bolometric corrections [98], we can also estimate whether electromagnetic missions contemporary to LISA will be able to detect the source. In the X-rays, the future Athena wide field imager, with field of view ∆Ω = 0.4 deg 2 [99], will require for this source a total integration time of ∼ 2.8 days. The near-infrared or optical emission of this system is also within the sensitivity limit of the near-infrared camera instrument on the James Webb Space Telescope [100]. As for the radio emission, the low-mass MBHs (≲ 10 7 M ⊙ ) typical of EMRI systems detectable by LISA are not expected to host powerful radio jets. However, radio emissions may be observed from, e.g., synchroton emission or even relativistic protojets [101].

C. Biased parameter estimation from environmental effects
Finally, we consider the case in which the EMRI GW signal is analyzed ignoring environmental effects and ask: how would ignoring environmental torques affect the inference of the EMRI parameters and tests of general relativity?
We investigate this aspect by analyzing the same GW signal from a migrating EMRI injected in the previous section. This time, we analyze the data using two waveform templates that do not allow for environmental torques ("vacuum template" and "GR deviation template") [102].
Firstly, we search for a signal using the vacuum template. We perform several runs where the MCMC walkers are allowed to explore a parameter space with priors extending up to 5% around the true value [103]. When using an incorrect template, we are not guaranteed to find a maximum likelihood point when trying to match the full signal present in four years of LISA data. In fact, the migration torque we consider here is strong enough that we cannot find any match , i.e. any maximum likelihood point.
In the last part of the inspiral the orbital decay due to GW emission is stronger than the disk-torque dissipation. Therefore, we expect to match the signal with a vacuum template when considering only the portion of the data closer to the plunge. We refine our search by considering the last 3.5, 3, 2.5, and, 2 years of data. We find a maximum likelihood only when we analyze the last two years of data. This maximum likelihood L ≈ 0.07, is approximately 14 times smaller than the one obtained with the correct template ("migration template"). Figure 4 shows the posterior for primary mass and spin recovered by the analysis of the last 2 years of data with the vacuum template. For reference, we also show the posterior distribution using the template matching the injection ("migration template"). We find that the vacuum-template posteriors are significantly biased, i.e. they are shifted 3σ away from the true values. In and spin a for the analysis of the same GW signal affected by environmental effects injected in Fig. 2 and analyzed with three different waveform models (templates). The contours show the 1,2,3-sigma Gaussian credible regions and the black dashed lines show the 95% credible intervals. The migration template (dark red) takes into account the presence of environmental effects and is able to match the full four-year inspiral, whereas the vacuum (black) and the GR deviation (blue) templates can match only the last two years of inspiral and recovers a biased value of primary mass and spin. The vacuum template is 3-sigma biased, whereas the 1σ contours from the GR deviation posteriors are consistent with the injected parameters. particular, unaccounted (inward) migration leads to an overestimation of the mass and spin of the primary, as it increases the rate of inspiral.
The absolute size of these biases is small and would not adversely impact any conclusions about the astrophysics of the sources. However, if a similar bias occurred on a parameter that characterizes a deviation from GR, it could shift the inferred value of that parameter away from zero and possibly lead to false detection of a GR deviation.
We verify this by searching for a deviation from GR in the two years of data with a standard parametrized-PN model. In particular, we allow for a deviation from GR coming from a time-varying gravitational constant [104][105][106][107]. This deviation will manifest itself in the waveform as a −4PN term and, therefore, can be accounted for in our agnostic model of Eq. for the primary mass and spin shown Fig. 4 is slightly shifted from the injected parameters and broader than the two posteriors obtained with the migration and vacuum templates. However, this broadening makes it consistent at 1-sigma with the injected parameter values. The posterior for the amplitude of the GR deviation shown in Fig. 5 is centered around A GR deviation = 3.35 +1.6 −1.6 × 10 −5 and it is 1.87σ inconsistent with GR, where the amplitude is expected to be zero. This demonstrates the degeneracy between disk effects and modifications of gravity in EMRI signals, and should motivate further studies into how to test GR with systems potentially affected by the environment.
We expect to perform exquisitely sensitive tests of GR with EMRI observations [16,64,105,108,109], but our analysis suggests it will be important to allow for additional environmental perturbations when carrying out these tests. Additionally, if environmental effects are ignored, the residuals between the template and the signal might affect parameter estimation of other sources [110].

V. DISCUSSION
In this paper, for the first time, we quantitatively study how to measure accretion-disk-induced torques with GW observations from EMRIs. The analysis we carry out assumes the binary is affected by migration in radiatively efficient and geometrically thin disks [12,19], and it is based upon a realistic waveform-generation formalism for EMRI parameter estimation within a fully Bayesian framework.
We investigate three different scenarios in which torques from accretion disks could affect EMRI parameter inference. In the first scenario, we measure how well accretion-disk effects are bounded to zero if they are absent, and we find that the migration torque amplitude can be constrained with a precision of ∼ 5 × 10 −6 . This allows us to infer when environmental effects are strong enough to be detectable. Interpreting the constraints as coming from migration torques, we confirm previous estimates in Ref. [19] that migration is observable assuming β disk prescriptions, but we also point out that the same is true for a wide range of accretion-disk parameters with α disks. Our work is the first realistic assessment of the detectability of an environmental effect in an EMRI, identifying the region of parameter space that GW observations can realistically probe.
In the second scenario, we analyze the GW signal of a typical migrating EMRI in an α disk. We find that we can distinguish between different disk prescriptions and constrain the amplitude of the environmental effects with ∼ 20% relative error. Using the (agnostic) measurements of torque amplitude and mass, we infer 2D marginalized posteriors for the disk viscosity and accretion rate. Moreover, assuming a multimessenger measurement of the bolometric luminosity with 10% precision, we show that the viscosity can be measured with 50% precision.
In the third and final scenario, we investigate the size of biases in EMRI parameter estimates caused by ignoring a strong migration torque. Our proof-ofprinciple analysis shows that the size of the bias that one should expect from reasonable migration torques will not significantly affect the inference of the astrophysics of galactic nuclei.
However, we demonstrated that these biases adversely affect tests of general-relativity with EMRIs [16], and that unmodeled environmental effects can lead to a false "detection" of a deviation from GR. This emphasizes the importance of including environmental effects when performing tests of GR. Should a population of EMRIs be detected, we also expect beyond-GR effects to be universal across the population, unlike environmental effects. Bayesian model selection pipelines should therefore be able to tell the two models apart.
Our work highlights the science potential of EMRIs embedded in accretion disks and the need for accurate torque models. In this work, we use prescriptions for EMRI migration that are designed for planetary (type-I) migration in a 3D isothermal disk [49]. This model has several limitations: for instance, it does not account for the fact that migration torques can change significantly close to the inner edge of the disk [111], where migration can halt altogether [112]. In general, EMRI migration differs from planetary migration in that the secondary object inspirals rapidly due to GW emission [33]. Targeted numerical simulations, although limited by the wide range of scales and timescales involved, will therefore be crucial to accurately model migration torques for GW observations. The first such simulations (in 2D) have shown promising results [20,33]. Similarly to what happened in planetary science, we could see in the next decade a progression from 2D to 3D simulations and the inclusion of more and more physics (radiation, magnetic fields, temperature and entropy gradients, thinner disks, etc.).
Another important takeaway point from this work concerns EMRI search and parameter inference strategies. Phenomenological models capable of capturing a host of environmental effects are likely to be needed in future analyses. Our work discusses a possible way of doing this. Our analysis should be considered as a proof-of-principle study of the impact of environmental effects on inference on EMRIs. The code used in this work will be available in the near future as an extension of the few package [26].
While this work relies on a reference EMRI configuration, we expect the detectability of environmental effects to improve when the small compact object explores farther regions around the MBH. This is due to the negative PN order of the effects considered in this work, which become dominant over gravitational emission at low frequencies. For fixed inspiral length and primary mass, a larger secondary mass (due to accretion [11,84] and/or hierarchical mergers [30,113] in the disk) would not only have higher SNR, but also be observable at larger radii. Therefore, intermediate massive black hole binary systems might be the best sources for detecting environmental effects. However, accurate waveforms for such systems are not available yet, and the environmental effects might require different modeling from the one presented here [20].
In our work, we have assumed circular equatorial orbits. A realistic, agnostic analysis of an EMRI signal with LISA would allow for non-zero eccentricity and inclination, even when looking for accretion-disk torques (or at least perform Bayesian comparisons between the two hypotheses). There is also a possibility that eccentricity may evolve as a result of environmental torques, as seen in binaries with a circumbinary disk [114]. Since EMRI observations are able to measure eccentricity with 10 −5 relative precision, the study of such a scenario would require accurate eccentric vacuum trajectories, and a reliable model for the environment-driven evolution. We also did not consider the possibility of searching for the migrating signal with an eccentric or inclined vacuum waveform. In this case, we expect to be able to tell if the EMRI was truly eccentric by studying the harmonic content of the signal. We plan to investigate the effects of eccentricity and inclination in future work, when accurate generic orbits around Kerr BHs become available.
While our study shows that accretion-disk properties can be resolved with EMRIs if observed for 4 or possibly more years (up to 10), it remains to be seen if the same holds true when considering different source classes detectable by LISA [115], or when taking into account competing torques, such as from dark-matter spikes [116][117][118][119][120][121][122][123], hierarchical triples [124], or modifications to GR [16,108,125,126]. A detailed study of the distinguishability between these different effects will be the subject of future work. In our analysis we derived the constraints LISA could put on the amplitude of two disk-induced effects, which predict different torque powers n r . Other beyond-vacuum effects might also manifest with a specific power law-like dependence on the orbital separation. Here, we explore how the constraints change as a function of n r for our reference EMRI source. We show the results in Fig. 6 in terms of the symmetric 95% bounds on the amplitude A 95% . We find that for n r ≳ 3 the bound can be fitted with a straight line in log-space as follows, log 10 A 95% = −4.63 − 0.14 n r .
(B1) Similar results are found in Fig. 8 and 9 of Ref. [132] and in Fig. 2 of Ref. [105], where the bounds are set on a different amplitude parameter. In future work, we plan to investigate how to map our parametrization to the parametrized post-Einsteinian expansions [89].