# Measuring eccentricity and gas-induced perturbation from gravitational waves
of LISA massive black hole binaries

Mudit Garg,1 Andrea Derdzinski,1,2,3 Shubhanshu Tiwari,4 Jonathan Gair,5 and
Lucio Mayer1
1Department of Astrophysics, University of Zurich, Winterthurerstrasse 190,
CH-8057 Zürich, Switzerland
2Department of Life and Physical Sciences, Fisk University, 1000 17th Avenue
N., Nashville, TN 37208, USA
3Department of Physics & Astronomy, Vanderbilt University, 2301 Vanderbilt
Place, Nashville, TN 37235, USA
4Physik-Institut, Universität Zürich, Winterthurerstrasse 190, 8057 Zürich,
Switzerland
5Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
Am Mühlenberg 1, Potsdam 14476, Germany E-mail<EMAIL_ADDRESS>
(Received / Accepted)

###### Abstract

We assess the possibility of detecting both eccentricity and gas effects
(migration and accretion) in the gravitational wave (GW) signal from LISA
massive black hole binaries (MBHBs) at redshift $z=1$. Gas induces a phase
correction to the GW signal with an effective amplitude ($C_{\rm g}$) and a
semi-major axis dependence (assumed to follow a power-law with slope $n_{\rm
g}$). We use a complete model of the LISA response, and employ a gas-corrected
post-Newtonian in-spiral-only waveform model TaylorF2Ecc. By using the Fisher
formalism and Bayesian inference, we explore LISA’s ability to constrain
$C_{\rm g}$ together with the initial eccentricity $e_{0}$, the total
redshifted mass $M_{z}$, the primary-to-secondary mass ratio $q$, the
dimensionless spins $\chi_{1,2}$ of both component BHs, and the time of
coalescence $t_{c}$. We find that simultaneously constraining $C_{\rm g}$ and
$e_{0}$ leads to worse constraints on both parameters with respect to when
considered individually. Assuming a standard thin viscous accretion disc, for
$M_{z}=10^{6}~{}{\rm M}_{\sun}$, $q=8$, $\chi_{1,2}=0.9$, and $t_{c}=4$ years,
we can confidently measure (with a relative error of $<50$ per cent) an
Eddington ratio as small as ${\rm f}_{\rm Edd}\sim 0.1$ for a circular binary
while for an eccentric system only ${\rm f}_{\rm Edd}\gtrsim 1$ can be
inferred. The minimum measurable eccentricity is $e_{0}\gtrsim 10^{-2.75}$ in
vacuum and $e_{0}\gtrsim 10^{-2}$ in the presence of a circumbinary disc. A
weak environmental perturbation (${\rm f}_{\rm Edd}\lesssim 1$) to a circular
binary can be mimicked by an orbital eccentricity during in-spiral, implying
that an electromagnetic counterpart would be required to confirm the presence
of an accretion disc.

###### keywords:

methods: data analysis – methods: statistical – black hole physics –
gravitational waves – accretion, accretion discs.

††pubyear: 2024††pagerange: Measuring eccentricity and gas-induced
perturbation from gravitational waves of LISA massive black hole binaries–12

##  1 Introduction

The prospect of the observation of gravitational waves (GWs) in the mHz band
in the 2030s looks promising following the adoption by ESA of the Laser
Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017; Barack et al.
2019) and with other projects, such as TianQin (Wang et al., 2019) and Taiji
(Gong et al., 2021), being developed. One of the primary expected
extragalactic sources for LISA are massive black hole binaries (MBHBs) with
primary-to-secondary mass ratios $q\lesssim 10$ and total masses between
$10^{4}~{}{\rm M}_{\sun}$ and $10^{8}~{}{\rm M}_{\sun}$, which LISA will be
able to observe up to redshift $z\sim 20$ (Amaro-Seoane et al., 2017). Another
expected source are intermediate/extreme mass ratio inspirals (I/EMRIs; Babak
et al. 2017; Amaro-Seoane 2018) with $q\gtrsim 10^{2}$, which can be observed
up to $z\lesssim 2$. MBHBs, with their high signal-to-noise ratios (SNRs;
Amaro-Seoane et al. 2017), provide exciting opportunities to not only measure
source properties with high accuracy but also place constraints on the
properties of their environments.

The main formation channel for MBHBs is via galaxy mergers (Begelman et al.,
1980). To shrink these binaries from a large scale to the coalescence phase
within a Hubble time requires an environmental perturbation that could come
from either gas or stars (see, e.g. Amaro-Seoane et al., 2023). In this paper,
we will be primarily concerned with the dynamical effects of gas, as they can
non-negligibly perturb both the semi-major axis and eccentricity of MBHBs in
the LISA regime more strongly than stellar interactions, given the tight
separations. Therefore, when we refer to an environment we will always mean a
gas accretion disc. MBHs are often observed to be accompanied by an accretion
disc at the center of active galactic nuclei (AGN) galaxies, especially beyond
$z\gtrsim 1$ and up to $z\lesssim 7$ (Padovani et al., 2017). Therefore, as
galaxy mergers trigger gas inflow and AGN activity (Mayer, 2013), MBHBs can be
driven to coalescence by a surrounding gas reservoir. For the near-equal mass
MBHBs considered here, we expect the accretion disc to take the form of a
circumbinary disc (CBD; D’Orazio et al. 2016).

GWs can be an important tool to not only measure the source properties but
also probe imprints of the environment in which the binary is evolving. In the
coalescence phase, gas mainly affects the binary via migration torques and
mass accretion. The detectability of the imprint on the emitted waveform of
these effects (see, e.g. Barausse et al., 2014; Garg et al., 2022; Caputo et
al., 2020) is dependent on the details of the gas inflow. While most of the
works on environmental measurements from GWs have focused on I/EMRIs (Levin,
2007; Barausse et al., 2014; Derdzinski et al., 2019, 2021; Speri et al.,
2023; Cole et al., 2023), recent works have indicated that measurements of gas
effects on more equal-mass MBHBs are also possible (Garg et al., 2022;
Dittmann et al., 2023; Tiede et al., 2023). However, these studies for MBHBs
have focussed on the detectable accumulated dephasing in the GW waveform
caused by gas, modelled using only Newtonian-order terms. This makes it
optimistic as well as impossible to confidently pin down gas as a sole cause
for this dephasing in the absence of an electromagnetic (EM) counterpart,
since either higher-order post-Newtonian (PN) terms, eccentricity, or other
environmental effects can also produce similar dephasings (see, e.g. Zwick et
al., 2023). Therefore, it becomes crucial to identify the region of parameter
space in which we could confidently identify a gas accretion disc as the
environment of a MBHB, utilizing only the observed emitted GWs.

While shrinking the MBHB from a large scale, gas can also excite eccentricity
that can be measurable up to $\sim 10^{-2.75}$ in the LISA band one year
before the merger (Garg et al., 2024), despite partial circularization due to
GWs (Peters & Mathews, 1963; Peters, 1964). Depending on whether the CBD is
prograde or retrograde and extremely or moderately thin, and whether the
binary is equal-mass or unequal-mass, we can expect different eccentricities
in the LISA band. Therefore, measurement of eccentricities can provide
evidence towards certain disc configurations even if gas effects themselves
become negligible near coalescence. However, if there are measurable gas
effects, then there can be an interplay between them and the eccentricity when
performing parameter estimation.

This work considers eccentric binaries of two aligned spinning MBHs embedded
in a CBD. We aim to consider eccentricity as well as gas parameters during
parameter estimation for either a one-year or four-year observation window. To
be close to realistic data analysis methodologies, we use high-order post-
Newtonian eccentric waveforms with aligned spin corrections to the circular
part, we model LISA’s motion and use the time delay interferometry response
model, which will be needed to cancel the laser noise, and we consider both
analytical and numerical techniques to assess the achievable constraints on
the parameters of interest.

The paper is structured as follows. In Section 2, we explain our methodology
to include gas-induced corrections in the GW phase for eccentric MBHBs.
Section 3 studies modeling of environmental effects from CBD simulations to
get the leading-order dephasing from different gas effects. We summarize our
parameters of interest and waveform model in Section 4. In Section 5, we
analytically compute errors on different parameters using the Fisher matrix
formalism. We summarize our results from Bayesian inference in Section 6. In
Section 7, we study if a wrong template can mimic an injected signal. We
discuss our findings in Section 8 and summarize the key takeaways of this work
in Section 9.

##  2 The stationary phase approximation

Let us consider two spinning BHs in vacuum at redshift $z$, with a redshifted
total binary mass $M_{z}$ and a primary-to-secondary mass ratio $q\geq 1$,
revolving around each other in an eccentric orbit with their dimensionless
spins111The dimensionless spin of a BH of mass $m$ and spin angular momentum
$J$ is $\chi\equiv cJ/Gm^{2}$, where $G$ is the gravitational constant and $c$
is the speed of light. $\chi_{1,2}$ aligned to the angular momentum of the
binary. This is equivalent to the motion of a BH of reduced mass $\eta M_{z}$,
where $\eta\equiv q/(1+q)^{2}$ is the symmetric mass ratio, in an elliptical
orbit around a black hole of mass $M_{z}$, fixed at the focus.222See Tables 1
and 7 for definitions of commonly used variables and terms. This orbit has a
LISA frame semimajor axis $a$, eccentricity $e$, and orbital angular frequency
$\Omega$.

Due to the quadrupole nature of GW emission, the GW emission from a small
eccentricity ($\lesssim 0.1$) and near-equal mass system is dominated by the
$n=2$ eccentric harmonic, which is twice the orbital frequency:

$f=\frac{1}{\pi}\left(\frac{GM_{z}}{a^{3}}\right)^{\frac{1}{2}},$ (1)

In the small-eccentricity limit, we can also approximate the orbital angular
frequency as $\Omega=\pi f$ at all PN orders.

At the Newtonian-order,333Newtonian-order terms are denoted by the superscript
(0). the orbital averaged GW-driven semi-major axis decay rate is (Peters &
Mathews, 1963; Peters, 1964)

$\displaystyle\dot{a}^{(0)}_{\rm GW}$
$\displaystyle=-\frac{64}{5}\frac{G^{3}}{c^{5}}\frac{\eta
M_{z}^{3}}{a^{3}}\text{F}(e),$ (2)

where

$\displaystyle\text{F}(e)$
$\displaystyle=\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right)(1-e^{2})^{-\frac{7}{2}}.$
(3)

The stationary phase approximation (SPA), which holds for slowly varying phase
and amplitude over an orbital period (Cutler & Flanagan, 1994), is valid for
the inspiral part of the GW signal. The SPA phase can be expressed as

$\displaystyle\psi(f)$ $\displaystyle=2\pi ft_{c}-\phi_{c}+2\pi
f\frac{da}{\dot{a}}-2\pi\int f\frac{da}{\dot{a}},$ (4)

where $t_{c}$ and $\phi_{c}$ are the time and phase of coalescence,
respectively.

For simplicity, we only consider circular Newtonian-order amplitude with the
quadrupole mode (i.e., $(2,2)$ mode) for all cases. The phase is more
sensitive than the amplitude to minor corrections arising from small
eccentricity ($e\lesssim 0.1$) (Moore et al., 2016), and the same should be
true for weak environmental effects.

Since both eccentricity and environmental interactions affect the phase
evolution of the source, we describe the cumulative phase of an event by its
individual contributions: $\psi_{\rm TF2}$ is the phase a circular inspiral
will accumulate in vacuum, solely due to emission of GWs; $\Delta\psi_{\rm
TF2Ecc}$ represents the phase correction to an event’s waveform when orbital
eccentricity alters its GW emission; and $\Delta\psi_{\rm gas}$ represents the
phase correction that is a consequence of environmental interaction, which
further speeds up or slows down the inspiral. For the latter, we adopt the
interaction with a gas disk as our fiducial environmental effect. As we
discuss later on, $\Delta\psi_{\rm TF2Ecc}$ and $\Delta\psi_{\rm gas}$ have
negligible cross-terms. We consider terms up to $3.5$PN order (Buonanno et
al., 2009) for $\psi_{\rm TF2}$ with aligned spin corrections also up to
$3.5$PN order (Arun et al., 2009; Mishra et al., 2016), $3$PN and
$\mathcal{O}(e^{2})$ order for $\Delta\psi_{\rm TF2Ecc}$ (Moore et al., 2016),
and the leading-order correction from gas in $\Delta\psi_{\rm gas}$. The total
SPA phase can be expressed as

$\psi=2\pi ft_{c}-\phi_{c}+\psi_{\rm TF2}+\Delta\psi_{\rm
TF2Ecc}+\Delta\psi_{\rm gas}.$ (5)

The overall semi-major axis evolution rate can be written down assuming no
cross-term between GWs and gas effects444Thus far, most hydrodynamical
simulations show that the gas torques and accretion rates are not strongly
affected by GW-inspiral (Tang et al., 2018; Derdzinski et al., 2019, 2021).
However, there may be exceptions within the parameter space, which is yet to
be fully explored.:

$\displaystyle\dot{a}$ $\displaystyle=\dot{a}_{\rm GW}+\dot{a}_{\rm gas}.$ (6)

Given expected traditional disk model properties, we assume $\dot{a}_{\rm
gas}\ll\dot{a}_{\rm GW}$ in the near-coalescence phase. Therefore, Eq. (4) can
be expanded to separate the different contributions (all the integration
constants are absorbed into $t_{c}$ and $\phi_{c}$):

$\displaystyle\psi_{\rm TF2}+\Delta\psi_{\rm TF2Ecc}$ $\displaystyle=2\pi
f\int\frac{da}{\dot{a}_{\rm GW}}-2\pi\int f\frac{da}{\dot{a}_{\rm GW}},$ (7)
$\displaystyle\Delta\psi_{\rm gas}$ $\displaystyle=-2\pi f\int
da\frac{\dot{a}_{\rm gas}}{\dot{a}^{2}_{\rm GW}}+2\pi\int
da~{}f\frac{\dot{a}_{\rm gas}}{\dot{a}^{2}_{\rm GW}}.$ (8)

The GR phasing contributions, $\psi_{\rm TF2}$ and $\Delta\psi_{\rm TF2Ecc}$
are well known in the literature. The gas contribution, $\Delta\psi_{\rm
gas}$, appears here in its current form for the first time, although different
forms of the same expressions do exist in the literature (see, e.g. Yunes et
al., 2011).

In the next section, we consider the evolution of an MBHB in the presence of
an accretion disc to get an estimate of $\Delta\psi_{\rm gas}$.

##  3 Modeling effects from a circumbinary disc

The torque exerted by a near-equal mass MBHB typically carves out a central
cavity in the inner accretion disc that leads to the formation of a
circumbinary disc (CBD; D’Orazio et al. 2016). However, streams of gas still
flow into this cavity, feeding mini-disks that accrete on to the binary as
well as adding an additional torque component apart from an outer CBD (Farris
et al., 2014). Non-axisymmetric features in this gas configuration exert a
gravitational torque on the binary, which can lead to inspiral or outspiral of
the MBHB depending on various disc-binary parameters. Typically a binary
shrinks before GWs take over to drive it to merger, if we have a retrograde
system (Tiede & D’Orazio, 2024), or a prograde disc with an unequal mass
binary (Duffell et al., 2020), a sufficiently thin disc (Tiede et al., 2020),
or a moderately eccentric system (D’Orazio & Duffell, 2021; Siwek et al.,
2023). In 2D CBD studies, this gas torque is usually expressed as a function
of the accretion rate ($\propto\dot{M}a^{2}\Omega$) and in the circular
limit555All quantities in the circular limit have bar on the top. as
$\bar{\Gamma}_{\rm CBD}=\xi\dot{M}a^{2}\Omega$, where the fudge factor $\xi$
depends upon the disc parameters and binary mass ratio (Duffell et al., 2020;
Garg et al., 2022). This expression is consistent with the viscous torque
estimate (Lin & Papaloizou, 1986). Furthermore, $\xi$ tends to be positive
(expand binary) for an equal-mass binary (D’Orazio & Duffell, 2021) and
negative (shrink binary) for $q\gtrsim 10$ (Cuadra et al., 2009; Moody et al.,
2019; Muñoz et al., 2019; Muñoz et al., 2020; Duffell et al., 2020; Dittmann
et al., 2023; Tiede et al., 2020; Tiede & D’Orazio, 2024). However, depending
on the thermodynamic assumptions, we can also have a negative torque for a
circular near-equal mass system (Bourne et al., 2023). Therefore, in this
work, we will survey $\xi$ values that cover all realistic possibilities.

We consider the impact of both migration and mass accretion on the MBHB
evolution in the following sections.

###  3.1 Migration

Studying the influence of gas-induced migration on the GW inspiral of the MBHB
near merger has a lot of challenges. Gas effects on a tight near-equal mass
circular MBHB have been simulated extensively, but only in the regime where GW
inspiral is not important. Therefore, extrapolating results from these studies
to near coalescence ($\dot{a}_{\rm gas}\ll\dot{a}_{\rm GW}$) could potentially
lead to errors. A few studies on circular extreme-mass ratio BHBs embedded in
a gas disc by Tang et al. (2018); Derdzinski et al. (2019, 2021) find that gas
effects do not change due to GW emission near the merger. However, they only
consider Newtonian-order terms (Peters, 1964) without including higher-order
relativistic corrections. The inclusion of eccentricity could further
exacerbate these problems. Hence, to encompass modeling uncertainties
regarding gas migration effects on the embedded eccentric MBHB in the LISA
band, we can write down a generic power law 666We found any cross terms
between gas and eccentricity to be heavily suppressed in the phase, see
Appendix B.

$\displaystyle\dot{\bar{a}}_{\text{mig}}=$
$\displaystyle{\mathcal{A}}\left(\frac{a}{GM_{z}/c^{2}}\right)^{n_{\rm
g}}\dot{\bar{a}}^{(0)}_{\rm GW},$ (9)

where dimensionless $\mathcal{A}$ and $n_{\rm g}$ are assumed to be constants
for the duration of the binary coalescence time, $t_{c}$. Effectively, gas
corrections to the SPA phase due to migration will enter at the $-n_{\rm g}$
PN order in this parametrization. Given that gas effects become increasingly
negligible compared to GW emission towards the merger, we can safely assume
$n_{\rm g}>0$. Moreover, we will only consider Newtonian-order hydrodynamical
simulations777We note that relativistic corrections to the gas and binary
motion can result in changes to the gas dynamics and resulting torque
(Berentzen et al., 2009; Liu, 2021), but here we focus on detecting gas with a
more generic parameterization. to study $\dot{a}_{\text{mig}}$.

For $\dot{a}_{\rm gas}=\dot{\bar{a}}_{\rm mig}$ in Eq. (8), we get the
leading-order dephasing from migration by assuming a CBD torque fiducial model
(see Appendix C):

$\displaystyle\Delta\psi_{\rm mig}$ $\displaystyle=\psi^{(0)}_{\rm TF2}C_{\rm
mig}v^{-8},$ (10) $\displaystyle C_{\rm mig}$ $\displaystyle=1.66\times
10^{-16}\frac{\xi}{1.0}\frac{{\rm f}_{\rm
Edd}}{1.0}\frac{0.1}{\epsilon}\left(\frac{\eta}{0.25}\right)^{-2}\frac{M_{z}}{10^{5}{\rm
M}_{\sun}},$

where $v\equiv(GM_{z}\pi f/c^{3})^{\frac{1}{3}}$ is the characteristic
velocity of the binary, $\psi^{(0)}_{\rm TF2}\equiv(3/128\eta)v^{-5}$ is the
leading-order phase term, ${\rm f}_{\rm Edd}$ is the Eddington ratio,888We can
express the accretion rate as $\dot{M}={\rm f}_{\rm Edd}\dot{M}_{\text{Edd}}$,
where $\dot{M}_{\text{Edd}}=2.22\times 10^{-8}M(0.1/\epsilon){\rm yr}^{-1}$.
and $\epsilon$ is the radiative efficiency, which we always set to $0.1$.

###  3.2 Mass accretion

The increase in mass of either BH due to mass accretion can also affect the
binary evolution. Heavier BHs will coalesce faster due to stronger GW emission
as per Eq. (2). Preferential accretion onto one of the BHs will change the
center of mass equilibrium, and accretion of gas with linear or angular
momentum will change a BH’s momentum and spin, respectively. Each of these
effects produces an additional effective torque, which can be comparable to or
much weaker than the gravitational component discussed in Section 3.1. We
neglect these components in this study, given that we only consider accretion
onto the total binary and not onto its individual components, nor its small-
scale gas configuration. We note however that the inclusion of these effects
may justify an increase in the torque amplitude $C_{\rm mig}$ in Eq. (10).

The phase contribution due to mass accretion is (see Appendix C):

$\displaystyle\Delta\psi_{\rm acc}$ $\displaystyle=\psi^{(0)}_{\rm TF2}C_{\rm
acc}v^{-8},$ $\displaystyle C_{\rm acc}$ $\displaystyle=-7.82\times
10^{-17}\frac{{\rm f}_{\rm
Edd}}{1.0}\frac{0.1}{\epsilon}\left(\frac{\eta}{0.25}\right)^{-1}\frac{M_{z}}{10^{5}{\rm
M}_{\sun}}.$

###  3.3 Effective dephasing due to gas

Both phasing terms from migration and accretion have the same $v^{-8}$
frequency dependence with respect to the vacuum GW phase, which implies that
gas corrections to the SPA phase enter at the $-4$PN order. The amplitudes
($C_{\rm mig}$ and $C_{\rm acc}$) of both effects have similar dependencies on
the binary-disc parameters. However, migration dephasing becomes stronger than
accretion for unequal mass ratios.

The overall phasing contribution from gas can be expressed as

$\Delta\psi_{\rm gas}=10^{-16}C_{\rm g}\psi^{(0)}_{\rm TF2}v^{-2n_{\rm g}},$
(11)

where $C_{\rm g}$ and $n_{\rm g}$ depend on the underlying gas model. For our
CBD torque fiducial model, we have:

$\displaystyle C_{\rm g}$ $\displaystyle=\frac{\xi}{1.0}\frac{{\rm f}_{\rm
Edd}}{1.0}\frac{0.1}{\epsilon}\left(\frac{\eta}{0.25}\right)^{-2}\frac{M_{z}}{10^{5}{\rm
M}_{\sun}},$ (12) $\displaystyle n_{\rm g}$ $\displaystyle=4.$ (13)

where we model $C_{\rm g}$ based on the migration dephasing in Eq. (10) and
any uncertainties about the dependence on the binary-disc parameters are
folded into $\xi$, which can either be positive or negative.

##  4 Parameter space, LISA response, and time delay interferometry

To study the evolution of a MBHB embedded in an accretion disc, we mainly
consider the binary-disc parameters summarized in Table 1, which are defined
in the LISA frame. There are eight intrinsic parameters (first eight rows of
Table 1) and six extrinsic parameters (last six rows of Table 1). We further
divide the intrinsic parameters into five intrinsic-merger
$\\{M_{z},q,\chi_{1},\chi_{2},t_{c}\\}$ and three intrinsic-inspiral
$\\{e_{0},C_{\rm g},n_{\rm g}\\}$ parameters, due to their relative importance
in different phases of the GW source evolution.

$\boldsymbol{\theta}$ |  Definition | Units
---|---|---
$M_{z}$ |  Total redshifted mass | ${\rm M}_{\sun}$
$q$ |  Mass ratio | Dimensionless
$\chi_{1,2}$ |  Spin parameters of both BHs | Dimensionless
$t_{c}$ |  Time of coalescence | year
$e_{0}$ |  Initial eccentricity | Dimensionless
$C_{\rm g}$ |  Environmental amplitude | Dimensionless
$n_{\rm g}$ |  Environmental semi-major axis power-law relative to GWs | Dimensionless
$D_{\rm L}$ |  Luminosity distance | Mpc
$\phi_{c}$ |  Phase at coalescence | Radian
$\imath$ |  Inclination | Radian
$\lambda$ |  Ecliptic latitude | Radian
$\beta$ |  Ecliptic longitude | Radian
$\psi$ |  Initial polarization angle | Radian
Table 1: Parameters of interest in the LISA frame.

We add the gas-induced phasing term in Eq. (11) to the TaylorF2Ecc (Moore et
al., 2016) phase and then modify lisabeta (Marsat et al., 2021) to incorporate
this waveform model. lisabeta takes into account LISA’s motion and computes
the time delay interferometry (TDI) response of the detector. We use the LISA
sensitivity curve including galactic confusion noise from Marsat et al.
(2021).

Based on realistic astrophysical expectations for the different parameters in
Table 1, we generate waveforms until the innermost stable circular orbit over
these parameter grids:

$\displaystyle M_{z}\in$
$\displaystyle\\{10^{4.5},10^{5},10^{5.5},10^{6}\\}~{}{\rm M}_{\sun},$
$\displaystyle q\in$ $\displaystyle\\{8,1.2\\}\iff\eta\in\\{0.1,0.25\\},$
$\displaystyle\chi_{1,2}\in$ $\displaystyle\\{0.9\\},$ $\displaystyle
e_{0}\in$
$\displaystyle\\{0,10^{-3},10^{-2.75},10^{-2.5},10^{-2.25},10^{-2},10^{-1.75},10^{-1.5},$
$\displaystyle 10^{-1.25},0.1\\},$ $\displaystyle t_{c}\in$
$\displaystyle\\{1~{}{\rm year},4~{}{\rm years}\\},$ $\displaystyle C_{\rm
g}\in$
$\displaystyle\\{-10^{3},-10^{2},-10^{1},-10^{0},-0.1,0,0.1,10^{0},10^{1},10^{2},10^{3}\\},$
$\displaystyle n_{\rm g}\in$ $\displaystyle\\{4\\}.$ (14)

The extrinsic fiducial parameters are $z=1$, which corresponds to $D_{\rm
L}=6791.3$ Mpc for the best-fit Planck Collaboration et al. (2020) cosmology,
and all angles are set to $0.5$ radians. The systems we consider here spend at
least 4 years in the LISA band before merging. We limit $e_{0}\geq 10^{-3}$,
since eccentricities below that are unmeasurable for MBHBs (Garg et al.,
2024). Considering that BHs in gas are expected to be highly spinning (see,
e.g. Reynolds, 2021), we chose $\chi_{1,2}=0.9$ as our fiducial case. However,
we find that our constraints on intrinsic-inspiral parameters are almost
independent of the exact spin magnitude. The SNRs of the sources in the grid
specified in (4) range from $\sim 150$ to $\sim 2500$ at $z=1$ with tiny
differences for the two times of coalescence.999See Fig. 2 of Garg et al.
(2024) for SNRs measured in a one year observation for our systems of interest
at different redshifts. In the next section, we compute relative errors on the
measurement of parameters.

##  5 Fisher formalism

###  5.1 Fisher matrix measurement of parameters for a $M_{z}=10^{5}~{}{\rm
M}_{\sun}$ and $q=8.0$ MBHB with $t_{c}=4$ years

Using the Fisher formalism (Vallisneri, 2008), we compute the expected error
($\sigma^{\rm Fisher}_{\theta}$) on each parameter $\theta$ and set a
threshold on the relative error $\delta\theta[\%]\equiv 100*\sigma^{\rm
Fisher}_{\theta}/|\theta|<50$ to define when a parameter is considered to be
measured. We use $50$ per cent as a threshold since we are interested in
determining whether the parameter differs from zero, i.e., whether an effect
is present, rather than wanting an accurate measurement of a large effect. We
first consider a MBHB with $M_{z}=10^{5}~{}{\rm M}_{\sun}$ and $q=8.0$ and
with a coalescence time of $t_{c}=4$ years, which has an SNR of $377.85$ and
generates $19159$ GW cycles in the LISA band. We study both circular and
eccentric (with $e_{0}=0.01$ or $e_{0}=0.1$) cases. In this section, we always
keep the intrinsic-merger parameters $\\{M_{z},q,\chi_{1},\chi_{2},t_{c}\\}$
free and we show relative errors in vacuum in Table 2. Unsurprisingly,
including eccentricity as a free parameter increases the uncertainties on all
the intrinsic-merger parameters. Furthermore, given that we expect our
fiducial CBD system to have the environmental power-law $n_{\rm g}=4$ and
$n_{\rm g}$ has a non-Gaussian posterior, based on Fig. 8, we always keep it
fixed to its fiducial value.

$\boldsymbol{\delta\theta[\%]}$ | $\boldsymbol{\rm Circular}$ | $\boldsymbol{e_{0}=0.01}$ | $\boldsymbol{e_{0}=0.1}$
---|---|---|---
$\delta M_{z}[\%]/10^{-2}$ | $3.41$ | $3.97$ | $3.97$
$\delta q[\%]/10^{-2}$ | $7.30$ | $8.51$ | $8.52$
$\delta\chi_{1}[\%]/10^{-2}$ | $6.52$ | $6.58$ | $6.56$
$\delta\chi_{2}[\%]/10^{-1}$ | $9.04$ | $9.56$ | $9.51$
$\delta t_{c}[\%]/10^{-7}$ | $2.67$ | $3.08$ | $3.05$
$\delta e_{0}[\%]/10^{0}$ | $-$ | $1.41$ | $0.01$
Table 2: The relative uncertainties, computed using the Fisher formalism, on
the intrinsic-merger parameters and eccentricity for a vacuum system with
$M_{z}=10^{5}~{}{\rm M}_{\sun}$, $q=8.0$, $\chi_{1,2}=0.9$, and $t_{c}=4$
years between three eccentricities: $e_{0}=0,~{}e_{0}=0.01$ and $e_{0}=0.1$.

####  5.1.1 Circular MBHB

We compute relative errors in the presence of an environment for a circular
system. In Fig. 1, we show results for $\delta C_{\rm g}[\%]$ as a function of
$C_{\rm g}$. We find that $\delta C_{\rm g}[\%]\approx 5.34\times
10^{2}/|C_{\rm g}|$, almost independently of the sign of the environmental
amplitude. The absolute error $\sigma^{\rm Fisher}_{C_{\rm g}}$ is independent
of the value of $C_{\rm g}$, as expected due to the phase correction from the
gas being linearly dependent on $C_{\rm g}$. Moreover, we can well-measure
$|C_{\rm g}|\gtrsim 10$.

Figure 1: The relative error on the measurement of the environmental effective
amplitude $\delta C_{\rm g}[\%]$ as a function of $C_{\rm g}$, denoted by a
red ‘+’ symbol, assuming zero eccentricity and fixed environmental power-law
$n_{\rm g}$. The dashed line represents the circular case. The $C_{\rm g}=0$
case is computed using $C_{\rm g}=10^{-3}$ for numerical reasons. We denote
the $50$ per cent well-measured threshold by a solid gray line.

The relative errors on the intrinsic-merger parameters with the inclusion of
gas are higher than the values in vacuum given in Table 2, but are independent
of the magnitude of $C_{\rm g}$. Therefore, the change in uncertainties is due
to having an extra free parameter in the model.

####  5.1.2 Eccentric MBHB

In the presence of an environment, we compute relative uncertainties for
injected signals with eccentricities $e_{0}=0.01$ and $e_{0}=0.1$. In Fig. 2,
we show results for $\delta C_{\rm g}[\%]$ and $\delta e_{0}[\%]$ as a
function of $C_{\rm g}$. We find that $\delta C_{\rm g}[\%]\approx 6.35\times
10^{3}/|C_{\rm g}|$, irrespective of the eccentricity value, as expected, but
higher than the circular case shown in Fig. 1. The eccentricity uncertainty,
$\delta e_{0}[\%]$, is nearly independent of the environmental amplitude, but
reaches higher values than for the vacuum case shown in Table. 2. Both
increases can be attributed to having an extra free intrinsic-inspiral
parameter. $|C_{\rm g}|\gtrsim 10^{2}$ is well-measured for our fiducial
system.

Figure 2: Same as in Fig. 1, but now including relative errors $\delta e_{0}$
for systems with initial eccentricities $e_{0}=0.01$ and $e_{0}=0.1$, denoted
by a blue ‘$\Ydown$’ symbol. We indicate the $e_{0}=0.01$ and $e_{0}=0.1$
cases by dot-dashed and solid lines, respectively. By including eccentricity,
the relative errors on $C_{\rm g}$ increases by almost an order of magnitude
with respect to Fig. 1.

The relative errors on intrinsic-merger parameters are again almost
independent of the magnitude of $C_{\rm g}$ but are higher than the values for
eccentric vacuum systems shown in Table 2, because of the extra free parameter
in the model.

####  5.1.3 Summary

The major conclusions that can be drawn from this Section are:

  * •

All relative errors are nearly independent of the sign of the environmental
amplitude $C_{\rm g}$.

  * •

The relative errors on the intrinsic-merger parameters,
$\\{M_{z},q,\chi_{1},\chi_{2},t_{c}\\}$, and eccentricity, $e_{0}$, are
independent of the value of $C_{\rm g}$, i.e., they remain almost the same for
any value of $C_{\rm g}$. Moreover, they only increase slightly with respect
to their values in vacuum when the environmental effect is included in the
model. This is mainly because of the inclusion of an extra free parameter.

  * •

The relative errors on $C_{\rm g}$ are higher in the eccentric case than in
the circular case. However, the change is the same for both $e_{0}=0.01$ and
$e_{0}=0.1$. Again, this is because of the inclusion of an extra degree of
freedom in the eccentric environmental model with respect to the circular
environmental one.

  * •

For all the scenarios explored here, amplitudes less than $C_{\rm g}]\lesssim
10$ are not well measured.

In the next section, we will extend the Fisher results to cover the full
parameter space defined in (4). Based on the results in this section, we will
only consider $C_{\rm g}\geq 1$ and choose $C_{\rm g}=100$ and $e_{0}=0.1$ as
our fiducial intrinsic-inspiral parameters.

###  5.2 Fisher matrix measurement of parameters for systems of interest

We would now like to explore the parameter space of Eq. (4). For this, we will
make matrix plots for the environmental amplitude $C_{\rm g}$ and the initial
orbital eccentricity $e_{0}$ by always considering a fixed environmental
power-law $n_{\rm g}=4$. We are again keeping the intrinsic-merger parameters
$\\{M_{z},q,\chi_{1},\chi_{2},t_{c}\\}$ free.

####  5.2.1 $\delta C_{\rm g}[\%]$

We want to study the uncertainty in the environmental amplitude $C_{\rm g}$ as
a function of $C_{\rm g}$ for both circular and eccentric systems. Based on
the findings of Section 5.1, we only need to consider positive values of
$C_{\rm g}\in\\{10^{3},10^{2},10^{1},10^{0}\\}$ and choose one eccentricity
(here $e_{0}=0.1$). In Fig. 3, we show $\delta C_{\rm g}[\%]$ for various
different choices of the intrinsic-merger parameters and strength of the
environmental effect. Since $M_{z}$ defines if the inspiral part of the signal
is in the low or high sensitivity region of the LISA frequency band, it
affects the uncertainty on $C_{\rm g}$. Moreover, $q$ and $t_{c}$ set the
number of GW cycle and so they also affect the measurement of the
environmental amplitude. Based on the parameter dependence in Eq. (12), for
$M_{z}=10^{6}{\rm M}_{\sun}$, $q=8.0$ or $\eta=0.1$, and $\xi=10$, we can
confidently measure ${\rm f}_{\rm Edd}\gtrsim 0.1$ for circular systems and
${\rm f}_{\rm Edd}\gtrsim 1$ for eccentric MBHBs. These constraints on ${\rm
f}_{\rm Edd}$ depends upon setting $\xi=10$, which is consistent with an
extrapolation of the simulation results of Dittmann & Ryan (2022), who found
increasingly strong torques on binaries the more thin the disc is in their
Fig. 3.

Figure 3: The relative Fisher errors on the environmental amplitude, $\delta
C_{\rm g}[\%]$, as a function of $|C_{\rm g}|$ and the total mass $M_{z}$. The
environmental power-law $n_{\rm g}$ is fixed to $4$ in all panels. In all
eight panels, we vary $M_{z}$ from $10^{4.5}$ to $10^{6}~{}{\rm M}_{\sun}$ and
$C_{\rm g}$ from $10^{3}$ to $1$. For all the left panels, we set $t_{c}=1$
year and assume $t_{c}=4$ years for all the right panels. We have either
$q=1.2$ (in the first and third rows) or $q=8.0$ (in the second and fourth
rows). In the top four panels, we assume circular orbits and in the bottom
four panels we set $e_{0}=0.1$, and allow the initial eccentricity $e_{0}$ to
also be a free parameter. We draw a solid blue line to identify the region of
well-measured parameters which have relative errors below $50$ per cent.
Moreover, we suppress all errors beyond $100$ per cent. We can at best well
measure $C_{\rm g}\gtrsim 10$ for circular binaries with $q=8$ and $t_{c}=4$
years and are not able to constrain $C_{\rm g}$ at all for eccentric systems
with $q=1.2$ and $t_{c}=1$ year.

####  5.2.2 $\delta e_{0}[\%]$

Now we want to compute the relative errors on the initial eccentricity,
$e_{0}$, as a function of $e_{0}$, with or without the environmental effect.
In Fig. 4, we show $\delta e_{0}[\%]$ in vacuum and in the presence of an
environment with $C_{\rm g}=100$. The measurement on $e_{0}$ strongly depends
upon $M_{z}$ for a similar reason to the dependence of $\delta C_{\rm g}[\%]$
upon $M_{z}$. It also varies weakly with $q$ and $t_{c}$, since limited
information is contained in the GW cycles observed from the very early
inspiral.

Figure 4: The relative error on the initial eccentricity, $\delta e_{0}[\%]$,
as a function of $e_{0}$ and $M_{z}$. In all eight panels we vary
eccentricities from $10^{-3}$ to $0.1$ and masses from $10^{4.5}$ to
$10^{6}~{}{\rm M}_{\sun}$. In the top four panels, we assume vacuum and in the
bottom four panels we have $C_{\rm g}=100$ and fixed $n_{\rm g}=4$. All other
parameter choices are the same as in Fig. 3. In vacuum, we are able to measure
at minimum $e_{0}\gtrsim 10^{-2.75}$. While in the presence of an environment,
at best we can measure $e_{0}\gtrsim 10^{-2}$.

For $t_{c}=1$ year we find a higher minimum measurable eccentricity in vacuum
than what is reported in Garg et al. (2024) because of the larger number of
free intrinsic-merger parameters $\\{\chi_{1,2},t_{c}\\}$ in the analysis.

In the next section, we will carry out Bayesian inference to verify the
results of this section and demonstrate the validity of our Fisher analysis.

##  6 Bayesian inference

We use the exact same Bayesian setup outlined in Section 4 of Garg et al.
(2024) to find posterior distributions for a few systems of interest. Salient
features include a zero noise realization and the Fisher
initialization101010starting walkers inside a multivariate Gaussian with its
mean given by the injected values and covariance by the Fisher matrix. to
accelerate the parallel tempering Markov-chain Monte Carlo (PTMCMC) runs using
the ptmcmc111111https://github.com/JohnGBaker/ptmcmc sampler. We set uniform
priors for all parameters. The resulting posteriors can be used to cross-
verify the Fisher matrix results reported in Section 5.

First, we show posteriors for our fiducial system in Fig. 5. We can infer that
all the parameters are well-recovered and the posteriors almost overlap with
the Fisher matrix predictions. We conclude that our Fisher results in Section
5 are robust.121212The posterior for $t_{c}=1$ year is similarly Gaussian as
in Fig. 5, although with higher covariances as expected.

Figure 5: Posterior distributions (solid black) for a gas-embedded eccentric
MBHB with parameters $M_{z}=10^{5}~{}{\rm M}_{\sun}$, $q=8.0$,
$\chi_{1,2}=0.9$, and $t_{c}=4$ years with effective gas-amplitude $C_{\rm
g}=100$ and fixed power-law $n_{\rm g}=4$. The blue lines mark the injected
values, the middle dashed line shows the median of the distribution, and the
two extreme vertical dashed lines indicate the $90$ per cent symmetric
credible interval. The two-dimensional contours of the posteriors indicate
$68$, $95$, and $99$ per cent credible intervals. We also show the Fisher
matrix predictions (dashed red) for comparison.

Next, we compare posteriors on the environmental amplitude for a circular
system to those for two eccentric systems in Fig. 6. The $C_{\rm g}$
posteriors in all cases peak around the injected value of $100$. Also, as
expected from Fig. 2, the shape of the posteriors for both non-zero
eccentricities are almost the same. In the circular case, we have no support
for $C_{\rm g}=0$. The broadening of the posterior is dominated by the extra
degree of freedom in the model, rather than by the the presence of
eccentricity in the signal itself.

Figure 6: The environmental amplitude posteriors for three eccentricities:
$e_{0}=0$ (solid black), $e_{0}=0.01$ (dashed red), and $e_{0}=0.1$ (dot-
dashed green). In all cases, we use the same template to inject the signal and
recover it.

Last, we compare $e_{0}$ and $C_{\rm g}$ posteriors between fixed intrinsic-
merger parameters, free intrinsic-merger parameters, and free intrinsic-merger
and extrinsic parameters in Fig. 7.131313See Fig. E2 of Garg et al. (2024) for
posteriors of the extrinsic parameters. We find that all posteriors peak
around the injected values with smallest uncertainties for the fixed
intrinsic-merger variables and largest uncertainties for the free extrinsic
variables. These results can be attributed to having fewer or more free
parameters. For fixed intrinsic-merger parameters, we have no support for
$C_{\rm g}=0$. Also, the errors on the intrinsic parameters are almost
independent of the inclusion or not of extrinsic parameters in the model.

Figure 7: $C_{\rm g}$ and $e_{0}$ posteriors’ comparison for three cases:
fixed intrinsic-merger parameters (solid black), free intrinsic-merger
parameters (our fiducial case; dashed red) or varying all parameters (dot-
dashed green) in Table 1 except $n_{\rm g}$.

In the next section, we compute biases and Bayes factors to compare different
templates when fitting a given data.

##  7 Fitting a wrong template

In this section, we examine the consequences of fitting a wrong template to an
injected signal. We compute Bayes factors of the correct template with respect
to the wrong template. We always set $\\{M_{z},q,\chi_{1},\chi_{2},t_{c}\\}$
to their fiducial values.

We evaluate Bayes factors by taking the ratio between the evidence ($Z$) of
fitting a true template and a false template to a given signal:

$\displaystyle\mathcal{B}=\frac{Z_{\rm true}}{Z_{\rm false}}.$ (15)

To reduce statistical variance during PTMCMC, we take the average of Bayes
factors from two sets of independent runs, while we report the recovered
parameter from the first set of runs.

###  7.1 Injecting only either eccentricity or gas perturbation

We fit an eccentric template to a gas-perturbed circular injected signal to
recover $e_{0}$, and a gas-perturbed circular template (fixed $n_{\rm g}$) to
an injected eccentric signal to recover $C_{\rm g}$.

In Table 3, we show the recovered values of $e_{0}$ and Bayes factors when
analysing an injected environmental circular signal for positive and negative
values of $C_{\rm g}$ with an eccentric template. The stronger the negative
$C_{\rm g}$, the higher the recovered eccentricity from its peaked Gaussian
posterior. There are also larger biases on the intrinsic-merger parameters
(see Appendix E). Since eccentricity can not explain an effect that slows down
the inspiral caused by positive $C_{\rm g}$ (an outward torque), it is
restricting $e_{0}$ to small values $\lesssim 10^{-3.5}$ with again large
biases and one-sided Gaussian posteriors for eccentricity. For weak
environmental effects, $C_{\rm g}\in\\{-10^{2},-10^{1},10^{1}\\}$, the Bayes
factors are inconclusive. However, for strong environmental effects we can
confidently accept the true template. Therefore, an EM counterpart will be
crucial to be sure about the presence of a weak environmental effect.

$\boldsymbol{C_{g,\rm inj}}$ | $\boldsymbol{\log_{10}e_{0}}$ | $\boldsymbol{\ln\mathcal{B}}$ | Strength
---|---|---|---
$-10^{3}$ | $-1.7$ | $\pagecolor{green}34.6$ | Decisive
$-10^{2}$ | $-2.1$ | $\pagecolor{lime}1.4$ | Inconclusive
$-10^{1}$ | $-2.6$ | $\pagecolor{pink}-2.7$ | Inconclusive
$10^{1}$ | $-3.9$ | $\pagecolor{lime}2.3$ | Inconclusive
$10^{2}$ | $-4.3$ | $\pagecolor{green}54.5$ | Decisive
$10^{3}$ | $-4.2$ | $\pagecolor{green}212.0$ | Decisive
Table 3: The injected environmental amplitude $C_{0,\rm inj}$, recovered
$e_{0}$, Bayes factor $\ln\mathcal{B}$ and the strength ov evidence for the
true model when fitting an eccentric vacuum model to a circular environmental
signal. Here $\ln\mathcal{B}\gtrsim 5$ (green) represents definitive evidence
in favor of the true model, $1\lesssim\ln\mathcal{B}\lesssim 3$ (lime) means
the true model is weakly preferred, and $-3\lesssim\ln\mathcal{B}\lesssim-1$
(pink) implies the false model is weakly preferred (Taylor, 2021). We term
$-3\lesssim\ln\mathcal{B}\lesssim 3$ as inconclusive and
$\ln\mathcal{B}\gtrsim 5$ as decisive. We recover negligible eccentricities
for positive migration ($C_{\rm g}>0$) since any significant $e_{0}$ cannot
explain slower inspiral.

In Table 4, we show the recovered value of $C_{\rm g}$ and Bayes factors, when
fitting a gas perturbed circular template to an injected eccentric signal for
various small eccentricities. We get narrow-peaked Gaussian posteriors for
$C_{\rm g}$. Even for $e_{0,\rm inj}=10^{-2.5}$, we find Bayes factors in
favor of the true model, which becomes stronger for higher eccentricity.

$\boldsymbol{\log_{10}e_{0,\rm inj}}$ | $\boldsymbol{C_{\rm g}}$ | $\boldsymbol{\ln\mathcal{B}}$ | Strength
---|---|---|---
$-2.5$ | $-19.2$ | $\pagecolor{lime}2.5$ | Inconclusive
$-2.0$ | $-192.7$ | $\pagecolor{green}6.0$ | Decisive
$-1.5$ | $-2307.4$ | $\pagecolor{green}58.7$ | Decisive
Table 4: Results obtained when analysing an injected eccentric vacuum signal
with a circular environmental template. The columns show the injected
eccentricity, $e_{0,\rm inj}$, recovered $C_{\rm g}$, Bayes factor,
$\ln\mathcal{B}$, and the strength of evidence in favour of the true model.

###  7.2 Injecting both eccentricity and gas perturbation

In this section, we fit either only a vacuum eccentric or an environmental
circular template to a gas-perturbed eccentric injected signal. Therefore,
recovering $e_{0}$ in the former case and $C_{\rm g}$ in the latter case,
respectively.

In Table 5, we show the recovered value of $e_{0}$ when fitting a vacuum
eccentric signal to an injected eccentric environmental signal. The posteriors
for $e_{0}$ are narrow-peaked Gaussians. We find that for a strong
environmental perturbation ($|C_{\rm g}|=10^{3}$) PTMCMC heavily favors the
true model, while the results are inconclusive for ($|C_{\rm g}|=10^{2}$).
This again emphasizes the need to have a complementary EM signal to have
definite proof of an environment.

$\boldsymbol{\log_{10}e_{0,\rm inj}}$ | $\boldsymbol{C_{g,\rm inj}}$ | $\boldsymbol{\log_{10}e_{0}}$ | $\boldsymbol{\ln\mathcal{B}}$ | Strength
---|---|---|---|---
$-2$ | $-10^{3}$ | $-1.6$ | $\pagecolor{green}33.3$ | Decisive
$-2$ | $-10^{2}$ | $-1.9$ | $\pagecolor{pink}-1.1$ | Inconclusive
$-2$ | $10^{2}$ | $-2.2$ | $\pagecolor{pink}-1.0$ | Inconclusive
$-2$ | $10^{3}$ | $-2.4$ | $\pagecolor{green}151.8$ | Decisive
Table 5: The injected $e_{0,\rm inj}$ and $C_{0,\rm inj}$, recovered $e_{0}$,
Bayes factor, and the strength of evidence in favour of the true model, when
fitting an injected eccentric environmental signal with a vacuum eccentric
template.

In Table 6, we show the recovered value of $C_{\rm g}$ when studying an
injected eccentric environmental signal with a gas-perturbed circular template
and find that while the true model is preferred over the false model, the
Bayes factors are not decisive.

$\boldsymbol{\log_{10}e_{0,\rm inj}}$ | $\boldsymbol{C_{g,\rm inj}}$ | $\boldsymbol{C_{\rm g}}$ | $\boldsymbol{\ln\mathcal{B}}$ | Strength
---|---|---|---|---
$-2$ | $-10^{3}$ | $-1194.5$ | $\pagecolor{lime}2.2$ | Inconclusive
$-2$ | $-10^{2}$ | $-293.9$ | $\pagecolor{lime}1.8$ | Inconclusive
$-2$ | $10^{2}$ | $-92.2$ | $\pagecolor{lime}1.8$ | Inconclusive
$-2$ | $10^{3}$ | $810.4$ | $\pagecolor{lime}1.7$ | Inconclusive
Table 6: Results obtained when analysing an injected eccentric environmental
signal with a circular environmental template. The columns show the injected
$e_{0,\rm inj}$ and $C_{g,\rm inj}$, recovered $C_{\rm g}$, Bayes factors
$\ln\mathcal{B}$, and the strength of evidence in favour of the true model.

##  8 Discussion

Most previous studies have assumed that near-merger MBHBs will be in a
circular binary since the expectation is that the orbital eccentricity would
be radiated away due to GW emission (Peters & Mathews, 1963; Peters, 1964)
before the source reaches the LISA band ($e_{\rm LISA}\approx 0$). However,
this has been challenged by recent suites of high-resolution hydrodynamical
simulations that embed an eccentric MBHB in a thin accretion disc (Zrake et
al., 2021; D’Orazio & Duffell, 2021; Siwek et al., 2023; Tiede & D’Orazio,
2024). These studies extrapolate to the LISA band using Peters (1964) orbital
evolution equations to find $e_{\rm LISA}\sim 10^{-3}$ for the case of
prograde discs, while $e_{\rm LISA}\sim 10^{-2}$ for retrograde CBDs given the
stronger boosting of eccentricity. Similarly, the interaction with a stellar
potential can lead to $e_{\rm LISA}\sim 10^{-3.5}$ (Gualandris et al., 2022).
In recent work Garg et al. (2024), we found that a gas-induced eccentricity as
low as $10^{-2.75}$ could be observable for MBHBs in a vacuum. As a
consequence, neglecting eccentricity could lead to potential biases in the
inference of the intrinsic-merger parameters,
$\\{M_{z},q,\chi_{1},\chi_{2},t_{c}\\}$. These studies also emphasize that a
gaseous environment can significantly alter the inspiral GW waveform of MBHBs
not only via observable gas-induced dephasing (Garg et al., 2022; Dittmann et
al., 2023), but also by exciting small but measurable eccentricities.

Usually gas perturbations for MBHBs are measured in terms of a cumulative
phase shift in the GW phase using Newtonian waveforms (Garg et al., 2022;
Dittmann et al., 2023). While this is a good first step to estimate if gas
could leave an observable imprint, it makes it impossible to relate this phase
shift unequivocally to the presence of a CBD disc. This is because either
higher PN corrections, a small eccentricity, or other environmental influences
such as a dark matter spike or third body interaction can mimic this effect
(Zwick et al., 2023). To be confident that this phase shift is most likely
from gas, we need to measure both its effective amplitude ($C_{\rm g}$) and
its power-law slope ($n_{\rm g}$). Since gas, via either migration or
accretion, induces a phase correction at the $-4$PN order, we not only need
high SNR but also numerous cycles in the inspiral phase to detect it. This is
why EMRIs are traditionally preferred for this kind of study (e.g., by Speri
et al. 2023), as a $q\gtrsim 10^{4}$ circular system spends $\sim 10^{5}$ GW
cycles in band during a four-year LISA observation window, which allows for
constraints on both $C_{\rm g}$ and $n_{\rm g}$ even for a moderate Eddington
ratio of ${\rm f}_{\rm Edd}=0.1$. What MBHBs lack in the number of cycles,
they make up for in their high SNRs. For this reason, we find that we would be
able to confidently measure (relative error $<50$ per cent) $C_{\rm g}$ for a
MBHB with $M_{z}=10^{6}~{}{\rm M}_{\sun}$ and $q=8$, if ${\rm f}_{\rm Edd}\sim
0.1$ for a circular system or ${\rm f}_{\rm Edd}\sim 1$ for an eccentric
system, assuming a strong gravitational torque ($\xi\sim 10$) from an
extremely thin accretion disc.

Estimating the strength of a gas torque ($\xi$) on the MBHB near-merger is
challenging, since most simulations study the system in the regime where GWs
are not dominant. Current simulations for non-inspiraling MBHBs predict
$\xi\lesssim 1$ for a moderately thin disc. There have been studies of the
effects of gas on the GW phase of an EMRI in the LISA band (Derdzinski et al.,
2019, 2021; Nouri & Janiuk, 2023), albeit using Newtonian-order
approximations, but nothing yet for measuring gas torques on coalescing MBHBs
aside from Dittmann et al. (2023).141414There are interesting studies on the
relativistic accretion flows onto merging MBHBs which show complex accretion
flows, however these are typically not evolved for long enough to measure
accurate torques (see e.g. Avara et al. 2023; Ennoggi et al. 2023; Gutiérrez
et al. 2022). Many of the simulations with MBHBs that focus on the long-term
evolution of the binary neglect magnetic fields and radiative transfer,
further increasing the uncertainty in the estimated torques on realistic
systems. Therefore, it is not clear whether $\xi$ becomes stronger, weaker, or
remains the same as a binary approaches merger. For this reason, in this study
we have remained agnostic about the torque strength, and considered a wide-
range of $\xi$ which allows for a $|C_{\rm g}|$ as high as $10^{3}$, in the
case that future studies find stronger gas torques in the LISA regime.

The effective gas-amplitude $C_{\rm g}$ in Eq. (12) depends on various binary-
disc parameters including the accretion rate onto the binary. Our
parameterization implies that a super-Eddington accretion rate may be
responsible for a high value of $C_{\rm g}$, _as well as_ a nonlinear
dependence of the torque strength $\xi$ in the regime of extremely thin
accretion disks which remains poorly explored. In other words, a degeneracy
exists between the disc parameters (in addition to $f_{\rm Edd}$ and radiative
efficiency) and the resulting torque ($\xi$). This can be broken by future
simulations of these systems that explore more representative parameters for
luminous AGN systems. Furthermore, inspiraling MBHBs are naturally expected to
produce bright EM counterparts, which can provide valuable constraints on
$f_{\rm Edd}$ and binary eccentricity. The presence of gas allows for a
characteristic X-ray emission during the inspiral (see, e.g. Haiman, 2017; Dal
Canton et al., 2019; Mangiagli et al., 2022; Cocchiararo et al., 2024).
Detection of such counterparts will not only confirm the presence of an
accretion disk, but also provide a narrower prior on the disc parameters from
an independent channel. This would help to break the degeneracies between disc
parameters such as ${\rm f}_{\rm Edd}$, $\xi$, and the radiative efficiency.
The combination of this knowledge with measurements of environmental
parameters from the GW signal will provide the strongest constraints on
accretion disc structure and binary-disc interaction.

The concurrent measurement of gas-amplitude and initial eccentricity leads to
worse constraints on both parameters. While it reduces the minimum measured
$C_{\rm g}$ from $\sim 10^{2}$ for a circular orbit to $\sim 10^{3}$ for an
eccentric system, it raises the eccentricity threshold from $e_{0}\sim
10^{-2.75}$ in vacuum to $\sim 10^{-2}$ in the presence of an accretion disc.
Since $e_{\rm LISA}$ depends more on the disc parameters than the binary
parameters for a near-equal mass MBHB (Zrake et al., 2021), we expect that
higher ${\rm f}_{\rm Edd}$ or $\xi$ will lead to a higher eccentricity in the
LISA band. Therefore, a concurrent measurement allows the possibility of
measuring both inspiral effects to place, in-turn, stronger limits on the disc
properties.

We also find that a weak environmental perturbation ($|C_{\rm g}|\lesssim 10$)
around a circular binary can be confused with an eccentric MBHB in a vacuum.
This is not true for the vice versa case, as a vacuum eccentric signal can not
be confused with a circular environmental system for small $e_{0}\gtrsim
10^{-2.5}$. This further emphasizes the importance of finding an EM
counterpart to conclusively determine that an accretion disc is present.
Otherwise, the parameters estimated from the inspiral part of the phase will
include biases which can not be completely mitigated by information from the
merger-ringdown part of the GW signal. A wrong analysis of the inspiral phase
could thus lead to a false detection of eccentricity, and it could even
possibly mimic a deviation from GR (Gair et al., 2013). This can have far-
reaching consequences, such as raising doubts on the validity of GR, or
leading to a biased distribution of MBHB parameters, which would be used to
disentangle SMBH formation and growth channels.

There are a number of caveats in the current work. The TaylorF2Ecc model does
not include spin-eccentricity cross-terms in the phase, which should be
negligible for the small eccentricities we study here. Moreover, we only
consider the Newtonian GW amplitude without eccentric and gas-induced
corrections. However, the inclusion of higher PN orders and eccentric-
environmental cross-terms should only help to improve the measurements of the
parameters. While the accuracy of TaylorF2Ecc reduces towards our cutoff at
the innermost stable circular orbit, we expect the results will not change if
an earlier cutoff is used (Garg et al., 2024). From the astrophysical
perspective, the phase correction due to the gas perturbation given by Eq.
(11) is a relatively simple model. It ignores any gas torque fluctuations
during the binary orbit (Zwick et al., 2022) that, in the very early inspiral
phase, produce a secular phase shift that mimics our model and can deviate
from our fiducial CBD torque (=$\xi\dot{M}\Omega r^{2}$) even if torque
fluctuations themselves are negligible after orbital average. Moreover, once
the orbital averaged value of $\dot{a}_{\rm gas}$ does not well-approximate
the semi-major axis decay rate towards the merger, these fluctuations can not
be ignored. However, this time-domain effect is cumbersome to include
consistently in the frequency-domain waveforms considered here, and we leave
this to future work. Also, including higher PN terms in the gas phase would
require even higher resolution hydrodynamical simulations so that we could go
beyond the Newtonian equations of motion, whose rich features may help us find
additional terms. Finally, we have always kept the radiative efficiency
$\epsilon$ as $0.1$ even for our highly spinning MBHs. Even if the dependence
of gas torque on $\epsilon$ remains the same near merger, it can vary for the
same spin magnitude between prograde or retrograde. We have kept it fixed, and
assume that any uncertainty is folded into $\xi$.

##  9 Conclusion

In this paper, we considered GWs from eccentric MBHBs embedded in a CBD to
estimate if both eccentricity and gas imprints (migration and accretion) could
be concurrently measured from the emitted gravitational waves observed by
LISA. We study systems of interest at redshift $z=1$ with highly spinning
($\chi_{1,2}=0.9$) individual BH masses, $M_{z}$, between
$10^{4}$–$10^{6}~{}{\rm M}_{\sun}$ and a primary-to-secondary mass ratio,
$q\in[1.2,8]$, such that the MBHBs spend at least four years in the LISA band
before merging. We considered both one-year and four-year times of coalescence
($t_{c}$) to study the measurability of the parameters for systems with
initial eccentricity, $e_{0}$, from $10^{-3}$–$10^{-1}$ and effective gas
amplitude, $C_{\rm g}$, between $-10^{3}$ and $10^{3}$. We assumed that the
power-law ($n_{\rm g}$) scaling of the environmental effect with semi-major
axis is $n_{\rm g}=4$, as expected from both migration and accretion (see Sec.
3). Due to the high expected SNRs, $\sim 150$–$2500$, we found that LISA
observations should be able to place constraints on the intrinsic-inspiral
variables, $e_{0}$ and $C_{\rm g}$, as well as on the intrinsic-merger binary
parameters, $\\{M_{z},q,\chi_{1,2},t_{c}\\}$. To account for LISA’s motion
around the Sun and model the time delay interferometry response, we used the
lisabeta software, and included the dephasing due to gas described by Eq. (11)
into the TaylorF2Ecc waveform model. We surveyed the parameter space
analytically using the Fisher formalism and then studied a few cases using
Bayesian inference. Finally, we assessed whether a weak environmental imprint
could be confused with a small eccentric signal using GW data alone. We
itemize our main findings below.

  * •

Since the gas correction to the GW phase depends linearly on $C_{\rm g}$ in
Eq. (11), the absolute error on $C_{\rm g}$ is independent of its magnitude
(see Fig. 1 and Section 5.1).

  * •

The cross terms between eccentricity and gas are negligible (see Appendix C).
Therefore, when constraining $C_{\rm g}$ and $e_{0}$ simultaneously, the
relative uncertainties on $e_{0}$ and $C_{\rm g}$ are independent of their
exact magnitude (see Fig. 2).

  * •

Using the Fisher formalism, we found that the relative errors on $e_{0}$ and
$C_{\rm g}$ are almost independent of the exact spin magnitude but have the
strongest dependence on $M_{z}$ out of the intrinsic-merger parameters (see
Figs 3 and 4). The constraints on $C_{\rm g}$ and $e_{0}$ depend upon $q$ and
$t_{c}$ (which set the number of GW cycles) strongly and weakly, respectively.
This is due to the rapid versus slow evolution in the very early inspiral of
$C_{\rm g}$ and $e_{0}$, respectively.

  * •

$C_{\rm g}\gtrsim 100$ is constrained to $<50$ per cent relative error for
circular binaries and $C_{\rm g}\gtrsim 1000$ for eccentric systems. This
translates to confidently measuring ${\rm f}_{\rm Edd}$ to around $0.1$ and
$1$ respectively in a four-year observation window for a $M_{z}=10^{6}~{}{\rm
M}_{\sun}$ and $q=8.0$ MBHB embedded in an extremely thin disc (see Fig. 3 and
the connection between $C_{\rm g}$ and accretion properties in Eq. 12).

  * •

We should be able to measure eccentricities as low as $10^{-2.75}$ in vacuum
and as low as $10^{-2}$ in the presence of an accretion disc (see Fig. 4).

  * •

Bayesian inference verifies the results of the Fisher results, showing
posteriors that overlapped with the prediction from the Fisher matrix (see
Fig. 5 and that peaked at the same value of $C_{\rm g}$ with or without
eccentricity (see Fig. 6).

  * •

Sampling extrinsic parameters does not affect the recovery of $C_{\rm g}$ and
$e_{0}$ (see Fig. 7).

  * •

An eccentric vacuum template can mimic a circular environmental signal for a
weak injected gas amplitude, $|C_{g,\rm inj}|\lesssim 10$. However, vice versa
does not hold for small $e_{0,\rm inj}\gtrsim 10^{-2.5}$ (see Tables 3 and 4).

  * •

An injected eccentric environmental signal could be confused with a vacuum
eccentric signal for $|C_{g,\rm inj}|\lesssim 100$ and $e_{0,\rm inj}=0.01$.
Similarly, the same signal could be mimicked by a circular gas-perturbed
template even for $|C_{g,\rm inj}|\sim 10^{3}$ (see Tables 5 and 6).

  * •

A stronger environmental perturbation or a higher eccentricity leads to
proportionally large biases on the intrinsic-merger parameters if fitted with
a wrong model (see Appendix E).

## Data availability statement

The data underlying this article will be shared on reasonable request to the
authors.

## Acknowledgements

MG, AD, and LM acknowledge support from the Swiss National Science Foundation
(SNSF) under the grant 200020_192092. ST is supported by the SNSF Ambizione
Grant Number: PZ00P2-202204. We acknowledge John G. Baker and Sylvain Marsat
for providing us access to their lisabeta software. The authors also
acknowledge use of the Mathematica software (Wolfram Research Inc., 2021),
NumPy (Harris et al., 2020), and inspiration drawn from the GWFAST package
(Iacovelli et al., 2022) regarding the python implementation of TaylorF2Ecc.

## References

  * Amaro-Seoane (2018) Amaro-Seoane P., 2018, Phys. Rev. D, 98, 063018
  * Amaro-Seoane et al. (2017) Amaro-Seoane P., et al., 2017, arXiv e-prints, p. arXiv:1702.00786
  * Amaro-Seoane et al. (2023) Amaro-Seoane P., et al., 2023, Living Reviews in Relativity, 26, 2
  * Arun et al. (2009) Arun K. G., Buonanno A., Faye G., Ochsner E., 2009, Phys. Rev. D, 79, 104023
  * Avara et al. (2023) Avara M. J., Krolik J. H., Campanelli M., Noble S. C., Bowen D., Ryu T., 2023, arXiv e-prints, p. arXiv:2305.18538
  * Babak et al. (2017) Babak S., et al., 2017, Phys. Rev. D, 95, 103012
  * Barack et al. (2019) Barack L., et al., 2019, Classical and Quantum Gravity, 36, 143001
  * Barausse et al. (2014) Barausse E., Cardoso V., Pani P., 2014, Phys. Rev. D, 89, 104059
  * Begelman et al. (1980) Begelman M. C., Blandford R. D., Rees M. J., 1980, Nature, 287, 307
  * Berentzen et al. (2009) Berentzen I., Preto M., Berczik P., Merritt D., Spurzem R., 2009, ApJ, 695, 455
  * Bourne et al. (2023) Bourne M. A., Fiacconi D., Sijacki D., Piotrowska J. M., Koudmani S., 2023, arXiv e-prints, p. arXiv:2311.17144
  * Buonanno et al. (2009) Buonanno A., Iyer B. R., Ochsner E., Pan Y., Sathyaprakash B. S., 2009, Phys. Rev. D, 80, 084043
  * Caputo et al. (2020) Caputo A., Sberna L., Toubiana A., Babak S., Barausse E., Marsat S., Pani P., 2020, ApJ, 892, 90
  * Cocchiararo et al. (2024) Cocchiararo F., Franchini A., Lupi A., Sesana A., 2024, arXiv e-prints, p. arXiv:2402.05175
  * Cole et al. (2023) Cole P. S., Bertone G., Coogan A., Gaggero D., Karydas T., Kavanagh B. J., Spieksma T. F. M., Tomaselli G. M., 2023, Nature Astronomy, 7, 943
  * Cuadra et al. (2009) Cuadra J., Armitage P. J., Alexander R. D., Begelman M. C., 2009, MNRAS, 393, 1423
  * Cutler & Flanagan (1994) Cutler C., Flanagan É. E., 1994, Phys. Rev. D, 49, 2658
  * D’Orazio & Duffell (2021) D’Orazio D. J., Duffell P. C., 2021, ApJ, 914, L21
  * D’Orazio et al. (2016) D’Orazio D. J., Haiman Z., Duffell P., MacFadyen A. I., Farris B. D., 2016, Mon. Not. Roy. Astron. Soc., 459, 2379
  * Dal Canton et al. (2019) Dal Canton T., Mangiagli A., Noble S. C., Schnittman J., Ptak A., Klein A., Sesana A., Camp J., 2019, ApJ, 886, 146
  * Derdzinski et al. (2019) Derdzinski A. M., D’Orazio D., Duffell P., Haiman Z., MacFadyen A., 2019, MNRAS, 486, 2754
  * Derdzinski et al. (2021) Derdzinski A., D’Orazio D., Duffell P., Haiman Z., MacFadyen A., 2021, MNRAS, 501, 3540
  * Dittmann & Ryan (2022) Dittmann A. J., Ryan G., 2022, MNRAS, 513, 6158
  * Dittmann et al. (2023) Dittmann A. J., Ryan G., Miller M. C., 2023, ApJ, 949, L30
  * Duffell et al. (2020) Duffell P. C., D’Orazio D., Derdzinski A., Haiman Z., MacFadyen A., Rosen A. L., Zrake J., 2020, ApJ, 901, 25
  * Ennoggi et al. (2023) Ennoggi L., Campanelli M., Zlochower Y., Lopez Armengol F., Mewes V., Ji L., Noble S., 2023, in APS April Meeting Abstracts. p. N09.009
  * Farris et al. (2014) Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2014, ApJ, 783, 134
  * Gair et al. (2013) Gair J. R., Vallisneri M., Larson S. L., Baker J. G., 2013, Living Reviews in Relativity, 16, 7
  * Garg et al. (2022) Garg M., Derdzinski A., Zwick L., Capelo P. R., Mayer L., 2022, MNRAS, 517, 1339
  * Garg et al. (2024) Garg M., Tiwari S., Derdzinski A., Baker J. G., Marsat S., Mayer L., 2024, MNRAS, 528, 4176
  * Gong et al. (2021) Gong X., Xu S., Gui S., Huang S., Lau Y.-K., 2021, in , Handbook of Gravitational Wave Astronomy. Springer Singapore, p. 24, doi:10.1007/978-981-15-4702-7_24-1
  * Gualandris et al. (2022) Gualandris A., Khan F. M., Bortolas E., Bonetti M., Sesana A., Berczik P., Holley-Bockelmann K., 2022, MNRAS, 511, 4753
  * Gutiérrez et al. (2022) Gutiérrez E. M., Combi L., Noble S. C., Campanelli M., Krolik J. H., López Armengol F., García F., 2022, ApJ, 928, 137
  * Haiman (2017) Haiman Z., 2017, Phys. Rev. D, 96, 023004
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Iacovelli et al. (2022) Iacovelli F., Mancarella M., Foffa S., Maggiore M., 2022, ApJS, 263, 2
  * Levin (2007) Levin Y., 2007, MNRAS, 374, 515
  * Lin & Papaloizou (1986) Lin D. N. C., Papaloizou J., 1986, ApJ, 309, 846
  * Liu (2021) Liu W., 2021, MNRAS, 504, 1473
  * Mangiagli et al. (2022) Mangiagli A., Caprini C., Volonteri M., Marsat S., Vergani S., Tamanini N., Inchauspé H., 2022, Phys. Rev. D, 106, 103017
  * Marsat et al. (2021) Marsat S., Baker J. G., Canton T. D., 2021, Phys. Rev. D, 103, 083011
  * Mayer (2013) Mayer L., 2013, Classical and Quantum Gravity, 30, 244008
  * Mishra et al. (2016) Mishra C. K., Kela A., Arun K. G., Faye G., 2016, Phys. Rev. D, 93, 084054
  * Moody et al. (2019) Moody M. S. L., Shi J.-M., Stone J. M., 2019, ApJ, 875, 66
  * Moore et al. (2016) Moore B., Favata M., Arun K. G., Mishra C. K., 2016, Phys. Rev. D, 93, 124061
  * Muñoz et al. (2019) Muñoz D. J., Miranda R., Lai D., 2019, ApJ, 871, 84
  * Muñoz et al. (2020) Muñoz D. J., Lai D., Kratter K., Miranda R., 2020, ApJ, 889, 114
  * Nouri & Janiuk (2023) Nouri F. H., Janiuk A., 2023, arXiv e-prints, p. arXiv:2309.06028
  * Padovani et al. (2017) Padovani P., et al., 2017, A&ARv, 25, 2
  * Peters (1964) Peters P. C., 1964, PhD thesis, California Institute of Technology
  * Peters & Mathews (1963) Peters P. C., Mathews J., 1963, Physical Review, 131, 435
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Reynolds (2021) Reynolds C. S., 2021, ARA&A, 59, 117
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
  * Siwek et al. (2023) Siwek M., Weinberger R., Hernquist L., 2023, MNRAS, 522, 2707
  * Speri et al. (2023) Speri L., Antonelli A., Sberna L., Babak S., Barausse E., Gair J. R., Katz M. L., 2023, Physical Review X, 13, 021035
  * Tang et al. (2018) Tang Y., Haiman Z., MacFadyen A., 2018, MNRAS, 476, 2249
  * Taylor (2021) Taylor S. R., 2021, arXiv e-prints, p. arXiv:2105.13270
  * Tiede & D’Orazio (2024) Tiede C., D’Orazio D. J., 2024, MNRAS, 527, 6021
  * Tiede et al. (2020) Tiede C., Zrake J., MacFadyen A., Haiman Z., 2020, ApJ, 900, 43
  * Tiede et al. (2023) Tiede C., D’Orazio D. J., Zwick L., Duffell P. C., 2023, arXiv e-prints, p. arXiv:2312.01805
  * Vallisneri (2008) Vallisneri M., 2008, Phys. Rev. D, 77, 042001
  * Wang et al. (2019) Wang H.-T., et al., 2019, Phys. Rev. D, 100, 043003
  * Wolfram Research Inc. (2021) Wolfram Research Inc. 2021, Mathematica, Version 13.0.0, https://www.wolfram.com/mathematica
  * Yunes et al. (2011) Yunes N., Kocsis B., Loeb A., Haiman Z., 2011, Phys. Rev. Lett., 107, 171103
  * Zrake et al. (2021) Zrake J., Tiede C., MacFadyen A., Haiman Z., 2021, ApJ, 909, L13
  * Zwick et al. (2022) Zwick L., Derdzinski A., Garg M., Capelo P. R., Mayer L., 2022, MNRAS, 511, 6143
  * Zwick et al. (2023) Zwick L., Capelo P. R., Mayer L., 2023, MNRAS, 521, 4645

##  Appendix A Definition of different terms

In Table 7, we summarize various quantities that are used repeatedly in the
main text.

Terms |  Definition
---|---
$\eta$ |  Symmetric mass ratio $q/(1+q)^{2}$
$v$ |  Characteristic velocity $(GM_{z}\pi f/c^{3})^{\frac{1}{3}}$
$\psi^{(0)}_{\rm TF2}$ |  leading-order circular phase contribution $(3/128\eta)v^{-5}$
$\Delta\psi_{\rm gas}$ |  leading-order gas-induced phasing correction $10^{-16}C_{\rm g}\psi^{(0)}_{\rm TF2}v^{-2n_{\rm g}}$
Intrinsic-merger parameters |  $\\{M_{z},q,\chi_{1},\chi_{2},t_{c}\\}$
Intrinsic-inspiral parameters |  $\\{e_{0},C_{\rm g},n_{\rm g}\\}$
Extrinsic parameters |  $\\{D_{\rm L},\imath,\phi_{c},\lambda,\beta,\psi\\}$
Fiducial parameters in the LISA frame |  $M_{z}=10^{5}~{}{\rm M}_{\sun},q=8,\chi_{1,2}=0.9,t_{c}=4~{}{\rm years},e_{0}=0.1,C_{\rm g}=100,n_{\rm g}=4,D_{\rm L}=6791.3~{}{\rm Mpc},\\{\imath,\phi_{c},\lambda,\beta,\psi\\}=0.5~{}{\rm radians}$
Table 7: Definition of different variables and terms in the main text.

##  Appendix B Gas-eccentricity cross terms?

To consider cross-terms between gas and eccentricity, we can modify
$\dot{a}_{\text{mig}}$ in Eq. (9) and add terms up to $\mathcal{O}(e^{2})$:

$\displaystyle\dot{a}_{\text{mig}}=$
$\displaystyle{\mathcal{A}}(1+\mathcal{A}_{1}e+\mathcal{A}_{2}e^{2})\left(\frac{a}{GM_{z}/c^{2}}\right)^{n_{\rm
g}}\dot{\bar{a}}^{(0)}_{\rm GW},$ (16)

where the dimensionless parameters
$\\{{\mathcal{A}},\mathcal{A}_{1},\mathcal{A}_{2},n_{\rm g}\\}$ are assumed to
be constants for the observation window and this parametrization should be
valid in the low-eccentricity limit.

There are only a few high-resolution hydrodynamical studies which have
considered both gas and eccentricity in the case of near-equal mass binary
systems embedded in a CBD. Muñoz et al. (2019); Muñoz et al. (2020); Zrake et
al. (2021); D’Orazio & Duffell (2021); Siwek et al. (2023) study prograde
orbits, and Tiede & D’Orazio (2024) focus on retrograde orbits. A recent study
by Siwek et al. (2023) has focused on the eccentricity evolution of unequal-
mass binaries. They all assume that the binary is accreting at the Eddington
rate (i.e., ${\rm f}_{\rm Edd}=1$) and with a radiative efficiency
$\epsilon=0.1$, and Shakura & Sunyaev (1973) viscosity coefficient
$\alpha=0.1$. However, the works listed here do not focus on the GW-dominated
regime (i.e., $\dot{a}_{\rm GW}>\dot{a}_{\rm gas}$), which adds further
uncertainty about how valid their results are in the LISA band. Still, we
include their results here as a starting point for understanding binary semi-
major axis and eccentricity coupling in gas. The values of
$\\{\xi,\mathcal{A}_{1},\mathcal{A}_{2}\\}$ inferred from simulations are
given in Table 8.

From | Motion | $\boldsymbol{q}$ | $\boldsymbol{\xi}$ | $\boldsymbol{\mathcal{A}_{1}}$ | $\boldsymbol{\mathcal{A}_{2}}$
---|---|---|---|---|---
DD21 | pro | $1.0$ | $0.31$ | $-13.61$ | $87.07$
SWH23 | pro | $1.0$ | $0.22$ | $-22.07$ | $66.76$
SWH23 | pro | $2.0$ | $0.16$ | $-24.13$ | $93.01$
SWH23 | pro | $10.0$ | $-0.05$ | $-68.09$ | $385.55$
TD23 | retro | $1.0$ | $-1.25$ | $0$ | $0$
Table 8: Coefficients for $\dot{a}_{\text{gas}}$ in Eq. (19), from D’Orazio &
Duffell (2021) (DD21), Siwek et al. 2023 (SWH23), and Tiede & D’Orazio (2024)
(TD23) for ${\rm f}_{\rm Edd}=1.0$, $\epsilon=0.1$, and $\alpha=0.1$. These
approximations are valid for $e\lesssim 0.15$.

For the gas-induced migration torque described by the parameterized form in
Eq. (16), the leading-order phase contribution can be computed assuming
$\dot{a}_{\rm GW}=\dot{a}^{(0)}_{\rm GW}$ in Eq. (7):

$\displaystyle\Delta\psi_{\rm mig}$ $\displaystyle=-\psi^{(0)}_{\rm
TF2}20{\mathcal{A}}v^{-2n_{\rm g}}\bigg{[}\frac{1}{(n_{\rm g}+4)(2n_{\rm
g}+5)}+$ $\displaystyle\frac{\mathcal{A}_{1}e_{0}}{(n_{\rm
g}+4+\frac{19}{12})(2n_{\rm
g}+5+\frac{19}{6})}\left(\frac{v_{0}}{v}\right)^{\frac{19}{6}}+$
$\displaystyle\left(\mathcal{A}_{2}-\frac{157}{12}\right)\frac{e_{0}^{2}}{(n_{\rm
g}+4+\frac{19}{6})(2n_{\rm
g}+5+\frac{19}{3})}\left(\frac{v_{0}}{v}\right)^{\frac{19}{3}}\bigg{]}.$ (17)

Note that even if there is no eccentricity-dependent term in $\dot{a}_{\rm
mig}$ (i.e. $\mathcal{A}_{1}=\mathcal{A}_{2}=0$), $\Delta\psi_{\rm mig}$ still
has an $e_{0}^{2}$ dependence due to $\dot{a}_{\rm GW}^{2}$ in the denominator
of Eq. (8).

Both $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ are suppressed by powers of
$e_{0}(v_{0}/v)^{19/6}$ and multiplied by ${\mathcal{A}}$, which is itself
extremely small ($\sim\mathcal{O}(10^{-16})$) as per Eq. (19). Therefore,
measurements of $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ either requires them
to be extremely large, or we need to observe the signal at a much lower
frequency than what LISA can measure. Current simulations suggest that
$\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ are at most $\mathcal{O}(10^{2})$ as
per Table 8.151515In this work we focus only on the orbital-averaged value of
the torque and neglect fluctuations that can have a higher magnitude (Zwick et
al., 2022). Hence, we drop the cross-terms in Eq. (10) with
$e_{0}(v_{0}/v)^{19/6}$ to recover the circular limit in Eq. (10).

##  Appendix C $\Delta\psi_{\rm gas}$ for a CBD model

###  C.1 Migration

Assuming a thin CBD torque fiducial model, we have in the circular limit:

$\displaystyle\bar{\Gamma}_{\rm CBD}$
$\displaystyle=\xi\dot{M}a^{2}\Omega=\frac{d}{dt}(\eta
M_{z}a^{2}\Omega)\implies\dot{\bar{a}}_{\rm
mig}=\xi\frac{\dot{M}}{M_{z}}\frac{2}{\eta}a,$

where in the first line we have ignored the mass accretion term when taking
the time-derivative (see Section 3.2). Therefore, we have

$\displaystyle\dot{a}_{\text{mig}}$
$\displaystyle={\mathcal{A}}\left(\frac{a}{GM_{z}/c^{2}}\right)^{n_{\rm
g}}\dot{\bar{a}}_{\rm GW},$ (19) $\displaystyle{\mathcal{A}}$
$\displaystyle=-8.64\times 10^{-16}\frac{\xi}{1.0}\frac{{\rm f}_{\rm
Edd}}{1.0}\frac{0.1}{\epsilon}\left(\frac{\eta}{0.25}\right)^{-2}\frac{M_{z}}{10^{5}{\rm
M}_{\sun}},$ $\displaystyle n_{\rm g}$ $\displaystyle=4.$

Plugging $\mathcal{A}$ and $n_{\rm g}$ into Eq. (B) in the circular limit
gives us $\Delta\psi_{\rm mig}$ in Eq. (10).

###  C.2 Mass accretion

Assuming both companion BHs are accreting at the same constant ${\rm f}_{\rm
Edd}$ during the observation window, we can write down the mass evolution as

$\displaystyle\dot{M}$ $\displaystyle={\rm f}_{\rm Edd}{\dot{M}}_{\rm
Edd}=2.22\times 10^{-8}\frac{{\rm f}_{\rm
Edd}}{\tau}\frac{0.1}{\epsilon}M_{z},$ $\displaystyle\implies M_{z}$
$\displaystyle=M_{z,0}\exp{\left(2.22\times 10^{-8}{\rm f}_{\rm
Edd}\frac{0.1}{\epsilon}t\right)},$ (20)

where $M_{0}$ is the initial mass.

Under the SPA, a phase shift due to mass accretion in the circular limit can
be expressed as (Caputo et al., 2020):

$\displaystyle\Delta\psi_{\rm acc}=\psi^{(0)}_{\rm
TF2}\frac{25}{256}\frac{{\rm f}_{\rm
Edd}}{\tau}\frac{0.1}{\epsilon}\frac{GM_{z}\eta^{-1}}{c^{3}}\left(\frac{1}{3}v_{0}^{-8}-\frac{15}{26}v^{-8}\right),$
(21)

where we have absorbed terms independent of $f$ into $\phi_{c}$ and terms
proportional to $f$ in $t_{c}$. We have replaced $M_{z,0}$ with $M_{z}$ after
computations, since the LISA observation duration we consider is at maximum
$4~{}{\rm years}$, implying $2.22\times 10^{-8}{\rm f}_{\rm Edd}t\ll 1$ in Eq.
(C.2). $\Delta\psi_{\rm acc}$ can be re-expressed as

$\displaystyle\Delta\psi_{\rm acc}$ $\displaystyle=\psi^{(0)}_{\rm
TF2}\mathcal{\hat{C}}_{\rm acc}+\psi^{(0)}_{\rm TF2}C_{\rm acc}v^{-8},$ (22)
$\displaystyle\mathcal{\hat{C}}_{\rm acc}$ $\displaystyle=4.52\times
10^{-17}\frac{{\rm f}_{\rm
Edd}}{1.0}\frac{0.1}{\epsilon}\left(\frac{\eta}{0.25}\right)^{-1}\frac{M_{z}}{10^{5}{\rm
M}_{\sun}}\left(\frac{a_{0}}{GM/c^{2}}\right)^{4},$ $\displaystyle C_{\rm
acc}$ $\displaystyle=-7.82\times 10^{-17}\frac{{\rm f}_{\rm
Edd}}{1.0}\frac{0.1}{\epsilon}\left(\frac{\eta}{0.25}\right)^{-1}\frac{M_{z}}{10^{5}{\rm
M}_{\sun}},$

where $a_{0}$ is the initial separation, which cannot be more than
$\sim\mathcal{O}(100GM_{z}/c^{2})$ to have the MBHB merge within the LISA
observation window. Therefore, $C_{I}$ is at maximum $\mathcal{O}(10^{-6})$
even for optimistic accretion rates (${\rm f}_{\rm Edd}$ $\lesssim 100$) for
near equal-mass MBHBs. Hence, the first term in Eq. (22) can be dropped due to
being negligible with respect to the leading-order SPA phase contribution
$\psi^{(0)}_{\rm TF2}$ and we get the accretion dephasing in Eq. (22).

##  Appendix D Posteriors for free environmental power-law

In Fig. 8, we show posteriors for our fiducial system when the environmental
power-law exponent, $n_{\rm g}$, is also a free parameter.

Figure 8: Same as in Fig. 5 but for free $n_{\rm g}$.

##  Appendix E Bias due to ignoring gas-perturbation or eccentricity

Our goal in this section is to quantify the potential shift in the measured
intrinsic-merger parameters (bias) in case that either eccentricity or
environment is neglected during analysis of a detected signal. To compute the
bias (denoted by $\Delta$) induced on merger parameters due to fitting a wrong
template, we take a difference of maximum likelihood values (denoted by
$\hat{}$ on top) when recovering an injected signal between a wrong template
(denoted by false) and a right template (denoted by true) and divide it by the
standard deviation of the given parameter when recovering with the true
template.

$\displaystyle\Delta\theta$ $\displaystyle=\frac{\hat{\theta}_{\rm
false}-\hat{\theta}_{\rm true}}{\sigma^{\theta}_{\rm true}}.$ (23)

This way of computation should also minimize statistical uncertainties in the
two models during PTMCMC. Moreover, we take average of biases from two sets of
independent runs to reduce statistical variance even further. We show results
in the following tables for the same systems of interest that were used in
Tables 3, 4, 5, and 6, respectively. The stronger that either gas-amplitude or
initial eccentricity are, the larger the bias induced in the intrinsic-merger
parameters.

$\boldsymbol{C_{g,\rm inj}}$ | $\boldsymbol{\Delta M_{z}}$ | $\boldsymbol{\Delta q}$ | $\boldsymbol{\Delta\chi_{1}}$ | $\boldsymbol{\Delta\chi_{2}}$ | $\boldsymbol{\Delta t_{c}}$
---|---|---|---|---|---
$-10^{3}$ | $-6.1$ | $-6.2$ | $-0.3$ | $2.0$ | $-4.2$
$-10^{2}$ | $-1.2$ | $-1.2$ | $-0.2$ | $0.5$ | $-0.9$
$-10^{1}$ | $-0.1$ | $-0.1$ | $0.0$ | $0.0$ | $-0.0$
$10^{1}$ | $-0.9$ | $-0.9$ | $-0.3$ | $0.5$ | $0.7$
$10^{2}$ | $-5.0$ | $-5.0$ | $-1.9$ | $3.0$ | $-4.3$
$10^{3}$ | $-11.3$ | $-11.4$ | $-2.3$ | $5.3$ | $-8.9$
Table 9: The injected environmental amplitude $C_{g,\rm inj}$ and biases on intrinsic-merger parameters due to fitting an vacuum eccentric template to a circular environmental signal. $\boldsymbol{\log_{10}e_{0,\rm inj}}$ | $\boldsymbol{\Delta M_{z}}$ | $\boldsymbol{\Delta q}$ | $\boldsymbol{\Delta\chi_{1}}$ | $\boldsymbol{\Delta\chi_{2}}$ | $\boldsymbol{\Delta t_{c}}$
---|---|---|---|---|---
$-2.5$ | $0.4$ | $0.4$ | $0.2$ | $-0.3$ | $0.4$
$-2$ | $2.1$ | $2.1$ | $0.2$ | $-0.8$ | $1.5$
$-1.5$ | $9.1$ | $9.2$ | $0.2$ | $-2.9$ | $6.2$
Table 10: The injected eccentricity $e_{0,\rm inj}$ and biases on intrinsic-merger parameters due to fitting a circular environmental template to a vacuum eccentric signal. $\boldsymbol{\log_{10}e_{0,\rm inj}}$ | $\boldsymbol{C_{g,\rm inj}}$ | $\boldsymbol{\Delta M_{z}}$ | $\boldsymbol{\Delta q}$ | $\boldsymbol{\Delta\chi_{1}}$ | $\boldsymbol{\Delta\chi_{2}}$ | $\boldsymbol{\Delta t_{c}}$
---|---|---|---|---|---|---
$-2$ | $-10^{3}$ | $5.8$ | $5.8$ | $0.4$ | $-2.0$ | $-3.9$
$-2$ | $-10^{2}$ | $0.9$ | $0.9$ | $-0.1$ | $-0.2$ | $0.6$
$-2$ | $10^{2}$ | $-1.1$ | $-1.1$ | $-0.0$ | $0.4$ | $-0.8$
$-2$ | $10^{3}$ | $6.5$ | $6.5$ | $0.8$ | $-2.6$ | $4.6$
Table 11: The injected eccentricity $e_{0,\rm inj}$ and environmental amplitude $C_{g,\rm inj}$, and biases on intrinsic-merger parameters due to fitting a vacuum eccentric template to an eccentric environmental signal. $\boldsymbol{\log_{10}e_{0,\rm inj}}$ | $\boldsymbol{C_{g,\rm inj}}$ | $\boldsymbol{\Delta M_{z}}$ | $\boldsymbol{\Delta q}$ | $\boldsymbol{\Delta\chi_{1}}$ | $\boldsymbol{\Delta\chi_{2}}$ | $\boldsymbol{\Delta t_{c}}$
---|---|---|---|---|---|---
$-2$ | $-10^{3}$ | $1.8$ | $1.8$ | $0.4$ | $-0.9$ | $1.5$
$-2$ | $-10^{2}$ | $1.7$ | $1.7$ | $0.1$ | $-0.7$ | $1.2$
$-2$ | $10^{2}$ | $1.3$ | $1.4$ | $0.1$ | $-0.6$ | $1.0$
$-2$ | $10^{3}$ | $1.9$ | $1.9$ | $0.5$ | $-1.1$ | $1.6$
Table 12: The injected eccentricity $e_{0,\rm inj}$ and environmental
amplitude $C_{g,\rm inj}$, and biases on intrinsic-merger parameters due to
fitting a circular environmental template to an eccentric environmental
signal.