The effect of interstellar medium on LVK’s black holes

Gravitational radiation alone is not efficient in hardening the orbit of a wide binary black hole (BBH). By employing a toy model for the interstellar medium (ISM) surrounding BBHs, here we discuss the effect of this baryonic medium on BBH dynamics. Depending on the BBH’s mass, we show that a binary surrounded by an isotropic cold neutral medium (i.e., an asymptotic temperature 𝑇 ∞ ≈ 100 K) with a time-averaged particle density of ⟨ 𝑛 𝐻 ⟩ = O( 1 ) cm − 3 can play a significant role in hardening the binary orbit over a O( 10 9 ) yr time scale. Additionally, this causes the black hole’s mass to grow at a rate ∝ 𝑚 2 . We thus discuss the impact of the ISM on the LIGO-Virgo-KAGRA (LVK) observables and quantify the properties of the ISM under which the latter could act as an additional important pathway for driving a subset of LVK’s BBH mergers.


INTRODUCTION
Emission of gravitational radiation alone is not efficient in hardening the orbit of a binary black hole (BBH).The merger time of a binary can be calculated entirely in the weak field limit, and for a circular BBH (with equal mass components) to merge within a Hubble time, the binary separation  has to satisfy (cf.Peters 1964 or Eq.40 later) where  0 ,  is the Hubble's constant and the mass of the component black hole (in  ⊙ ), respectively.For a typically stellar mass black hole with  ≤ 50 ⊙ , Eq. 1 requires  ≲ 69 ⊙ which is much smaller than the radial extent of most stars during the later phase of their evolution.However, the LIGO-Virgo-KAGRA (LVK) network of ground-based gravitational wave detectors has already observed close to eighty BBH mergers (Abbott et al. 2021(Abbott et al. , 2023)).
Although multiple 1D binary stellar evolution studies are able to generate tight BBHs from progenitor binary stars with relatively large initial separation (e.g., Dominik et al. 2012;Bavera et al. 2021;Zevin et al. 2021;Broekgaarden et al. 2022;van Son et al. 2022;Briel et al. 2023), this involves either limiting the mass accretion onto the primary black hole -during the phase of stable stellar mass transfer -to a sub-Eddington rate (e.g., van Son et al. 2020;Zevin & Bavera 2022) or fine-tuning certain free parameters during the common envelope evolution (CEE) phase of stellar evolution which lacks a full understanding (e.g., Gallegos-Garcia et al. 2021;Klencki et al. 2021;Marchant et al. 2021).In the case of stable mass transfer between binaries, limiting accretion to a sub-Eddington value causes the mass transfer to (at times) become non-conservative.This results in mass and angular momentum loss from the binary ★ sohanghodla9@gmail.comand thus a gradual hardening of the orbit.However, it is not clear if such an assumption is justified, with studies suggesting that super-Eddington accretion on black holes could be readily achievable (e.g., Sądowski & Narayan 2016;Ghodla & Eldridge 2023).On the other hand, CEE is a highly nonlinear process, and to attain a better understanding, 3D numerical simulation might be necessary.Meanwhile, conducting comprehensive multidimensional magnetohydrodynamic simulations of the CEE phase across a diverse range of spatiotemporal parameter spaces remains computationally too expensive (e.g., Ricker & Taam 2008, 2012;Chamandy et al. 2018;Moreno et al. 2022;Lau et al. 2022).As such, presently, one has to resort to 1D studies, which, at times, may make ad hoc choices for the associated free parameters to produce BBHs with favorable orbital separation (e.g., Breivik et al. 2020;Riley et al. 2022).
It is also disputed whether CEE could readily occur in close enough stellar binaries.In particular, in the detailed 1D treatment of binary interaction, mass transfer can remain stable (e.g., see Gallegos-Garcia et al. 2021;Klencki et al. 2021;Marchant et al. 2021 for analyses using Mesa models and Briel et al. 2023 using Bpass), which reduces the scope of CEE.As the latter plays a key role in hardening the orbit of the progenitors of BBHs, one can expect a reduced BBH merger rate estimate.Thus the question arises: what other potential mechanism could help to harden a binary orbit?Three-body interactions (e.g., Ziosi et al. 2014;Dorozsmai et al. 2024) or the efficient inspiral of compact remnants in the disks of active galactic nuclei (e.g., Leigh et al. 2018;Kaaz et al. 2023) may offer such alternatives.Separately, gravitational wave mergers may also occur via the dynamical capture of black holes in dense astrophysical environments (e.g., Kulkarni et al. 1993;Liu & Bromm 2020;Gamba et al. 2023), which may not require an efficient hardening mechanism.
Since BBHs reside in a sea of interstellar medium (ISM), here we explore the effect of this surrounding baryonic medium on the orbital dynamics and merger properties of BBHs.Presently, most population studies assume that BBHs (born from twin stars) live in isolation (i.e., vacuum).While this assumption is justified for a short evolutionary timescale, for longer timescales of O (10 9 ) yr (typical for LVK BBHs), this may strictly not be true.This is because, over such a large timescale, the minuscule impact of ISM interaction has a sufficient time to accumulate.As such, this work aims to quantify the characteristics of the baryonic medium surrounding the BBHs, under which it could play a important role in driving the LVK BBH mergers.
The ISM typically is an extremely tenuous fluid, containing nonrelativistic (charged and neutral) matter, relativistic charged particles known as cosmic rays, and magnetic fields (e.g., see Ferrière 2001; Draine 2011 for a review).A large fraction of the ISM in disk galaxies is in the form of neutral atomic hydrogen (HI) that can be broken into two parts: the cold neutral medium (CNM) with an average number density of   = 20-50 cm −3 and temperature  ≈ 100 K and the warm natural medium (WNM) with   = 0.6 cm −3 and  = 5000 K (e.g.Draine 2011).Under optimal circumstances, these two phases can coexist in pressure equilibrium, resulting in a twophase medium (Field et al. 1969;Wolfire et al. 1995, 2003, also see McKee & Ostriker 1977) containing relatively concentrated regions of CNM embedded in a diffused WNM.
Given that the progenitors of stellar black holes form in dense pockets of the H 2 regions, which themselves are the result of atomicto-molecular gas conversion, the surroundings of these progenitor stars would presumably be rich in neutral hydrogen.However, because of the intense radiation emitted by the black hole forming OB-type star, much of the ISM surrounding the BBHs immediately after formation (including the circumstellar medium created by the progenitor's stellar winds and mass ejection during mass transfer) should be composed of warm neutral and warm ionized hydrogen (the latter also known as a warm ionized medium, WIM).The abundance of the latter type will depend on the amount of ionizing photons emitted by the progenitor stars and the strength of the supernova explosion (if any).
There is also the possibility that the supernova explosion expels a substantial portion of the ISM to large distances, thus creating a scarcity of interstellar matter around the subsequently formed black holes.However, the absence of core-collapse supernova observations from stars with mass ≳ 18 ⊙ suggests that most black holes might form from a failed supernova explosion (Smartt 2015).When a successful explosion does occur, a rough estimate (using energy conservation) shows that matter displaced to a radial distance  can fall back on the BBH in time where  is the total mass of the BBH.For  = 10 ⊙ and  = 1 pc, this gives  ≈ 3 Myr, which is negligible compared to the lifetime of many potential LVK BBHs, i.e., O (10 9 ) yr.Furthermore, assuming no significant reheating of the surrounding ISM occurs post BBH formation, the time taken by the WNM to radiatively cool down can be estimated as where  = 3  /2 is its thermal energy,  being the gas temperature and C its mean cooling rate.The latter is defined as (e.g., Draine 2011) where Λ is the cooling function that depends on the gas temperature and metallicity.For a medium constituting of HI gas, we assume (Wolfire et al. 2003) (5) For typical values for the number density and temperature of the WNM (n  = 0.6 cm −3 ,  = 5000 K), this yields  cool ≈ 2 Myr (although this could be much higher in strongly metal deficient environments where cooling would be driven by Lyman- transition which does not depend on the metal abundance, e.g., chapter 30 in Draine 2011).A similar time frame also exists for the cooling of WIM (e.g., see Draine 2011).Since we are interested in BBHs that live for O (10 9 ) yr, as such in this paper, we assume that, for the most part, the BBH surroundings always have some presence of a generic cold ISM1 .We further assume that these BBHs have a negligible velocity w.r.t.their surrounding medium given their presumably common source of origin (the case for a high relative velocity is discussed later in Section 6).The validity of these assumptions is a major source of uncertainty in the present work, and therefore, the formalism employed should only be viewed as a toy model of the ISM around BBHs.Under these assumptions for the ISM, we show that the properties of the LVK's BBHs (including their mass and merger rates) residing in such a medium could be significantly affected.On the other hand, if the assumptions of our toy model are strongly violated, one may safely assume the ISM to have a negligible impact on the BBH dynamics.We also show that the accretion of WNM is inefficient and typically does not affect the BBH dynamics.Hence, the accretion of CNM will remain the focus of this paper.In all this, we ignore the effect of dark matter interaction with the BBHs and leave this analysis for a future study.
The remainder of this text is organized as follows.In Section 2, we discuss the effect of ISM accretion on the BBH mass and orbital period.In Section 3, we quantify the rate of the binary's orbital decay due to gravitational radiation emission, dynamical friction, and ISM accretion on the BBH.In Section 4, we then calculate the resulting BBH merger time incorporating all these means of orbital decay.In Section 5, we conduct a population study that allows us to quantify the effect of ISM on a generic cosmological population of BBHs and compare the result to the LVK observations.We end with a brief discussion in Section 6 followed by a conclusion in Section 7.

Effect on a black hole
Consider a Schwarzschild black hole immersed in an isotropic surrounding of ISM.Such a fluid can be approximated with an ideal gas equation of state with pressure  and density  obeying the polytropic relation Due to its gravitational influence the black hole would slowly accrete mass from the environment.The resulting mass accretion rate can be modeled using the Bondi-Hoyle accretion formula as (Bondi & Hoyle 1944) where  ∞ is the density,  ∞ is the speed of the sound, and  ∞ is the velocity of the black hole w.r.t. the ISM.The subscript "∞" indicates that the parameters are locally measured by an observer who is asymptotically far away from the black hole.This for us is the distance where the black hole's influence on the ISM becomes negligible.Additionally, the numerical factor (Bondi 1952) ; (1) =  3/2 /4 ≈ 1.12 (8) and where we have set  = 1 in the final expression for  ∞ .

Mass evolution of the black hole
For a black hole comoving w.r.t. to its surrounding ISM, this implies  ∞ = 0 with the resulting expression in Eq. 7 known as Bondi accretion rate (Bondi 1952).Solving Eq. 7 for  ∞ = 0 thus gives us the mass evolution of the black hole as where   is the mass of the black hole at some initial time   and ∞ (e.g., see Fig. 1 for an illustration).Eq. 10 suggests that the mass  of the black hole immersed in an ISM environment would approach infinity in a finite time Δ :=  −   .The latter can be calculated by setting the denominator in Eq. 10 to zero, yielding Interestingly, even for a modest (time-averaged) value of ⟨  ⟩ = 1 cm −3 and  ∞ = 100 K, this implies that the mass of a black hole with   ≳ 8  ⊙ would diverge within a Hubble time (e.g., Fig. 1).On the other hand, if the surrounding constitutes a medium with a higher temperature, the black hole mass growth is strongly restricted, e.g.see Fig. B1.In particular, for  ∞ ≳ 1000 K, one can calculate that a black hole would experience negligible mass growth for the quoted values of ⟨  ⟩.Thus, in this work, we focus only on a low-temperature environment.

An upper bound on ⟨𝑛 𝐻 ⟩
Even when the ISM has a low enough temperature, the supply of matter around the black hole would be limited, which would then set an upper limit on the mass of the black hole.Requiring that the latter remain finite allows us to set an upper bound on the value of ⟨  ⟩.For example, if we are interested in time-averaging over a Hubble time, then setting Δ = 1/ 0 gives where   (in the RHS expression) is assumed to be in units of  ⊙ .As expected, ⟨  ⟩ is directly proportional to  0 ,  ∞ and inversely proportional to   .

Mass evolution of the BBH
Now, let us consider a binary system composed of two Schwarzschild black holes embedded in an ISM.As earlier, we assume the binary as a whole to be coming w.r.t. to its surrounding ISM.However, the binary would also rotate around its center of mass, thus giving rise to a non-zero  ∞ for the individual black holes.Nevertheless, the gas is still asymptotically at rest w.r.t. the center of mass of the BBH.
Let us introduce the length scale known as the transonic or Bondi radius of the BBH configuration with total mass .This corresponds to the radial distance from the center of mass of the BBH at which the accretion flow first becomes supersonic.Taking  as the semi-major axis of the binary, if  ≪   (which is always the case for LVK black holes 2 ), then accretion at a large distance from the BBH can still be treated as Bondi accretion on a central compact object with mass .On the other hand, the accretion dynamics near the binary could be complex.However, numerical simulations of binaries embedded in a gaseous environment with  ≪   suggest that the accretion rate on the binary (as a whole) is enhanced compared to accretion on two masses with the same total mass  (Farris et al. 2010;Comerford et al. 2019;Kaaz et al. 2019).In this work, we adopt a conservative approach, assuming each black hole in the binary accretes independently at the rate it would accrete in isolation as determined by setting  ∞ = 0 in Eq. 7.

Impact on the BBH's orbital dynamics
The orbital angular momentum of an isolated binary can be expressed as where ,  and  =  1  2 /( 1 +  2 ) are the eccentricity, semimajor axis, and the reduced mass of the binary, respectively ( 1 ,  2 being the component mass and  =  1 +  2 ).If mass accretion remains isotropic -implying that the matter contains zero angular momentum at the radial asymptote, the inflowing material does not contribute to .Therefore,  should remain conserved (barring the loss in gravitational radiation).For simplicity, in this work, we set  = 0 which reduces  to the separation between the black holes.This is a conservative estimate as a more eccentric orbit decays faster by the emission of gravitational radiation, leading to a more rapid inspiral.
To understand the orbital evolution of such a mass-accreting BBH, it would be sufficient to study how  evolves with time.Using Eq. 10, we can calculate the time-evolution of the expression  2  as where in deriving the inequality, we have assumed  1 ≥  2 and hence  1 ≥  2 .Thus, from Eq. 14 and 16, an estimate for the evolution of  for an ISM accreting binary can be written as where   is the initial binary separation and the factor in the denominator is what causes  to become time dependent.For the rest of this paper, we assume the equality in the above equation to hold.This implies a less rapid orbital decay.We note that accretion of material can also occur from the circumstellar medium around the BBH, which is the leftover matter from the earlier episode of stellar evolution.In this case, the material might contain angular momentum, unlike the above-mentioned assumption.We discuss its impact on the BBH's orbital dynamics in Section 6.

Feedback on the ISM
The previous discussion highlights that the ISM temperature (at least for  >   ) must remain low for it to significantly impact the evolution of BBHs (e.g., see Fig. 1 and B1).However, there exists a potential scenario where the energy released during the accretion process could heat up the surrounding ISM, thereby increasing its temperature.

If the infalling matter does not form an accretion disk
In the region with  ≤   , the isotropically infalling matter moves at a supersonic rate, thus preventing it from radiating energy through viscous dissipation.Therefore, this should not increase the temperature of the gas.Meanwhile, the binary components with nonzero  ∞ can stir up their surrounding gas due to their gravitational influence, creating local density variations around the orbit.This can induce turbulence, thus causing the gas to radiate and heat the surrounding ISM.However, as mentioned previously, numerical studies of mass accretion in gravitationally bound systems (Farris et al. 2010;Comerford et al. 2019;Kaaz et al. 2019) suggest that this does not significantly affect the mass accretion rate.This may be attributed to the fact that in Eq. 7, the mass accretion rate on the BBH is determined by the gas properties near  ≳   , which is much larger than the separation of the black holes.Once the gas falls past   , it is descending at a supersonic velocity (which further increases with infall).As a result, the gas should become more or less bound to the binary.

If the infalling matter forms an accretion disk
It is possible that the infalling gas borrows sufficient angular momentum from the BBH's orbit such that a fraction of it is assimilated into an accretion disk around the component masses.The characteristics of such an accretion disk will depend on the magnitude of the mass accretion rate on the individual black holes (e.g., Frank et al. 2002;Yuan & Narayan 2014).For Bondi accretion on a black hole with mass  (which holds at least at larger scales), the accretion rate scales as Normalizing the above w.r.t. the Eddington accretion rate yields the fraction where  is the accretion efficiency and  Edd is the Eddington accretion rate and takes the form Above, in calculating the last term, we assumed that the opacity  of the matter is purely driven by the accretion of ionized hydrogen.
Let us further assume that only a fraction  of the accreted mass is assimilated into the accretion disk (with the rest getting accreted directly).Then, for the parameter values of   = 1 cm −3 and  ∞ = 100 K, Eq. 19 suggest that the mass assimilation rate in the disk w.r.t. the Eddington rate should should evolve as For illustrative purposes if we consider  ∈ [0.1, 1],  ∈ [3 ⊙ , 40 ⊙ ] and  = 1/16 3 , then one finds that for such a choice of   and  ∞ ,  ∈ (6.3 × 10 −5 , 0.01).Such values of  span two solutions of the accretion flow, namely the radiatively inefficient accretion disks (e.g., Narayan & Yi 1994;Blandford & Begelman 1999;Yuan & Narayan 2014) and the Shakura & Sunyaev (1973) accretion disks.In the case of radiatively inefficient accretion flows, the temperature of matter is almost virial.The inflowing plasma develops a two-temperature state, where the ions contain most of the energy, but it is the electrons that are the more efficient cooling source.The ions and electrons do not efficiently thermalize; thus, the disk develops a two-temperature plasma.As the disk cannot cool, nearly all of the energy generated through viscous dissipation will be advected by the black hole (Narayan & Yi 1994;Yuan & Narayan 2014) or a large fraction of the infalling matter will be lost in winds (Blandford & Begelman 1999;Yuan & Narayan 2014) (will part of the lost matter not significantly contributing towards  either).Such an accretion flow may form a corona at the inner edge of the accretion disks and thus glow in hard X-rays.However, the absorption cross-section of (hard) X-rays for low-metallicity matter is relatively small (Wilms et al. 2000).As such, they might not significantly heat the ISM surrounding the black holes, implying that a steady accretion rate could be maintained.Additionally, if the disk develops winds, the gas will also remove with it its acquired angular momentum, thus helping in hardening the binary without increasing the mass of the black hole.
In the scenario where the accretion flow becomes radiatively efficient (that is,  ≳ 0.01, e.g., Yuan & Narayan 2014), it will begin to radiate black-body-like radiation.For such a case, ISM heating could be significant and further mass will be prevented from efficiently assimilating in the disk.In such a scenario, the accretion process may still proceed but with a duty cycle.Initially, mass is assimilated in the disk, which eventually reaches the black hole on a vicious timescale.This mass then heats up the surrounding ISM, thus halting further mass accretion.The ISM then gradually cools over a few million years, allowing for an accretion disk to form again, thus continuing the cycle.As this may reduce the mass accretion rate, thus the BBH's evolution may not be affected significantly.We note that large  values are more likely for massive black holes.
The efficiency at which the infalling matter becomes part of the accretion disk rather than getting directly accreted remains a major caveat in the current work.In the following, we assume that  or  is small, which then allows us to ignore the effect of ISM heating in the vicinity of the black holes.The validity of the results shown in the following should be interpreted within the domain of this assumption.

Decay due to gravitational radiation with time-varying masses
The orbital decay of a BBH due to the emission of gravitational radiation in the weak field regime is governed by the time variation of the quadrupole moment tensor where  is the position of the point mass.The rate of radiative loss of binary's orbital energy can then be written as (Peters 1964) where the subscript "GW" stands for gravitational wave.For a circular orbit, this reduces to We note that due to time-varying masses, the first term on the RHS of the above equation is at least larger than the one given in Peters (1964) by a factor of  1  4 2 (cf.Eq. 16) which is implicit in the mass term .This results in a relatively larger loss of energy in gravitational radiation, which further increases monotonically with time.Additionally, the last term also appears due to the time-varying mass of the black holes in Eq. 22.As this term only has a sub-leading effect, it is thus ignored.If we now write the total orbital energy of the binary as Here, we consider a particle density of ⟨  ⟩ = 1 cm −3 .The two line styles for the blue curve (i.e., dynamical friction) apparently overlap due to the low resolution on the y-axis.
then Eq. 24 and 25 together give us the orbital decay rate due to gravitational wave emission as with  now being time-dependent.

Decay due to accretion of (asymptotically) zero angular momentum matter
As the accreted material adds mass but not angular momentum to the BBH, the orbit of the binary would gradually decay.While it is possible that the motion of the BBH could potentially stir up the ISM (thus injecting angular momentum into the surrounding media), the source of this angular momentum ultimately originates from the BBH's orbit.Eventually, the surrounding material will be accreted by the BBH.It is thus fair to assume that, outside of the angular momentum loss in gravitational radiation, the orbital angular momentum of the mass-accreting binary remains conserved.Consequently, we use Eq. 17 to calculate the orbital decay rate due to isotropic mass accretion as where the subscript "AM" implies that the decay occurs due to the conservation of the BBH's angular momentum.

Decay due to dynamical friction
As discussed in Section 2.2.1, even a BBH comoving with the ISM will still have components with  ∞ ≠ 0. Thus due to their interaction with the ISM, the black holes in the binary would consequently experience a gravitational drag force, also known as dynamical friction (Chandrasekhar 1943).Under the assumption that the surrounding medium constitutes a cold ISM with  ∞ = 100 K, all the black holes of interest would be moving at a supersonic speed w.r.t. the background medium.
The force due to dynamical friction on a mass moving from a uniform background medium can be approximated as (Ostriker 1999) where I = ln 1 − 1/M 2 /2 + ln ( max / min ).Here, M =  ∞ / ∞ is the Mach number of the problem with M ≫ 1 for supersonic motion.
Additionally,  max ,  min correspond to the minimum and maximum impact parameters of the interaction.We choose  min equal to the accretion radius of the black hole (also known as the Bondi-Hoyle-Lyttleton radius), i.e.,  min = 2/ 2 ∞ , while  max remains a free parameter in the present work and we set  max =  ∞ Δ, Δ being the time elapsed since the black hole was born.This allows us to calculate the energy dissipation rate as Using the above in conjunction with Eq. 25 then yields In the case of a black hole in a binary, the immediate surroundings of the black hole will differ from that assumed in deriving Eq. 28 (e.g., see Section 2.3).In particular, for such black holes, a large wake formation (behind the black holes) might be prevented which will lower the impact of dynamical friction on the black holes.As we discuss below, dynamical friction has only a sub-leading effect on the BBH dynamics, and thus, an exact treatment is not necessary.

BINARY MERGER TIME
Given the various forms of orbital decay discussed in the previous section, we now calculate the resulting BBH merger time.By differentiating Eq. 25 w.r.t.time, we find the total orbital decay rate as The individual components of the above equation are plotted in Fig. 2. In particular, it shows that at all binary separations, the decay term due to dynamical friction is dominated by at least one of the other two mechanisms 4 .As such, for ease of calculation, we safely ignore the decay resulting from dynamical friction.Moreover, qualitatively, one may assume that over long timescales, for BBH that comoves with the ISM, any momentum imparted to the gas by the BBH will eventually be accreted by them, hence negating the effect of dynamical friction.Defining where the subscript  indicate value at   then Eq. 31 yields 4 Fig. 2 only considers a subset of black hole masses and a fixed ⟨  ⟩.
However, one can check that the decay due to dynamical friction should remain subdominant for other relevant values of these parameters.Also, we note that in calculating the decay due to dynamical friction, we have generously set Δ = 10 9 yr in the expression for  max in Eq. 30.
This can be rewritten as This is a first-order linear ordinary differential equation, which on substituting () =  4 (), takes the form This equation has the exact solution where the integration factor takes the form Above, by definition, we have set (  ) = 1, yielding The merger time of an ISM accreting BBH is then the timescale  tot (the subscript implies net decay due to the various mechanisms detailed in Eq. 31) such that lim while the merger time of a BBH living in vacuum -i.e., purely due to emission of gravitational radiation -is (Peters 1964) Combining Eq. 38 and 39 then gives the relation between  tot and  GW as The above integral can be evaluated analytically using symbolic algebra tools but is not presented here due to its lengthy size5 .For comparative purposes, we instead provide the solution for the special case when  1, =  2, =   /2 (and set   = 0) yielding Since  GW can be readily calculated for a binary in the weak field limit, thus knowing the values of the initial mass of the black hole   , the ISM temperature at large distance from the BBH  ∞ and time-averaged particle density ⟨  ⟩, Eq. 42 makes the calculation of  tot straightforward.The solution to Eq. 42 is illustrated in Fig. 3, showing the quantitative effect of an ISM environment on the BBH merger time.
ratio of the binary  :=  2 / 1 ≳ 0.85.Thus, in Section 5, we interpolate between this general solution and the solution presented in Eq. 42 for  ∈ (0.8, 1].The general solution can be found in the jupyter notebook URL included in the Data Availability Statement.

IMPACT ON LVK OBSERVABLES
To estimate the impact of the ISM environment on the BBH merger properties, we apply the previous calculation for mass growth (Eq.10) and accelerated orbital decay (Eq.41) to a realistic population of stellar black holes.For the latter, we use the Bpass v2.2.1 detailed binary stellar models (Eldridge et al. 2017; Stanway & Eldridge 2018), which are then distributed over a range of formation redshifts.This distribution is a function of the redshift-dependent star formation rate (weighted with an initial mass function) and metallicity evolution of the Universe.Over the course of the subsequent stellar evolution, episodes of mass transfer and CEE may cause the orbit of these binaries to harden, resulting in BBH systems that could be potential sources for the LVK detectors (Acernese et al. 2015;LIGO Scientific Collaboration et al. 2015;Akutsu et al. 2021).
For more details on the implemented stellar population synthesis, we refer the reader to Appendix A in Ghodla et al. (2023) while the BBH merger rate calculation is discussed in Appendix A of the current work.

Stable mass transfer and common envelope evolution in stellar binaries
In separate studies (Briel et al. 2022(Briel et al. , 2023;;Ghodla et al. 2023), it has been observed that Bpass models underpredict the BBH merger rate.This can be attributed to the treatment of mass transfer episodes between progenitors of BBHs.Firstly, in the detailed treatment of the Bpass models, mass transfer is relatively more stable than in other rapid population synthesis codes that treat the interaction phase between the stellar binaries parametrically (e.g., Breivik et al. 2020;Riley et al. 2022).However, this is not just limited to Bpass, as similar effect has also been observed in the detailed models of Mesa (Gallegos-Garcia et al. 2021;Klencki et al. 2021;Marchant et al. 2021).Additionally, in Bpass super-Eddington accretion onto the black hole (from the secondary stellar companion) is allowed.
As such, when the mass transfer is stable, it also largely remains conservative, which reduces the scope for orbital hardening (Briel et al. 2023).Secondly, the Bpass models undergo more efficient6 CEE episodes as compared to the rapid population study codes.
Since CEE plays a key role in hardening the orbit of the progenitors of BBHs, one can expect a reduced BBH merger rate estimate.
For the current work, we take the Bpass models as a template and consider the impact of an ISM environment on its BBHs.

Redshift evolution of ⟨𝑛 𝐻 ⟩
As the galaxies in the early universe were presumably more gasrich, thus the supply of gas surrounding the BBHs is expected to last longer for those formed at higher redshifts .This implies that ⟨  ⟩ should be a monotonically increasing function of .Here, we assume this to take the form of a power-law as where  represents the redshift at BBH formation.The term ⟨  ⟩ 0 can be seen as ⟨  ⟩ evaluated at  = 0 or equivalently as ⟨  ⟩ when there is no redshift dependence.
In the simulation result shown in the following sections, we adopt two values for , namely  = 1 and 3 and set  ∞ = 100 K assuming a CNM surrounding.Moreover, only some BBHs may be found in an ISM environment that matches the requirements of our toy model.Thus, we present results assuming that of all BBHs formed in our simulation, only O (1)% of them are embedded in a favorable ISM surrounding while the rest either live in a vacuum or do not match the requirement of our toy model to experience a noticeable impact on their dynamics.The choice of this subset of O (1)% BBHs in the simulation is made on a random basis and contains no redshift bias.
Finally, although we do not consider the following scenario here, however instead of assuming that ⟨  ⟩ 0 takes the same value for all BBHs born in the simulation irrespective of their initial mass   , one may assume that due to a limited supply of matter, the surroundings of more massive black holes become ISM-poor at a faster pace (as massive binaries accrete at a faster rate, see Eq. 7 and Fig. 1) and thus will have a lower value of ⟨  ⟩ 0 .Additionally, for more massive black holes, a larger ⟨  ⟩ 0 also would raise the value of  at a relatively faster rate, thus causing ISM heating that will limit the accretion rate (see Section 2.3).Thus, to keep   low, one could employ a mass weighting in Eq. 43 and reformulate it as (cf.Eq. 12) where  is now the BBH mass-corrected variant of ⟨  ⟩ 0 that determines the time-averaged particle number density.

Source-frame BBH merger rate
More details on the source frame BBH merger rate calculation are presented in Appendix A1 and A2, while the resulting rates are shown in Fig. 4. Here, we assume that only 1.5% (top figure) and 3% (bottom figure) of the BBHs reside in an optimal ISM environment for their properties to be affected.These choices have been made to yield merger rates within the green band shown in the figure .A higher percentage would result in a merger rate that starts to diverge from the LVK inferred range7 .This implies that only a small fraction of the BBHs within our simulation should satisfy the condition of a favorable ISM environment.Separately, differences in implemented physics, such as the treatment of CEE or mass transfer stability, can affect the number of BBHs that can survive in the Bpass code compared to other population synthesis models.As more binaries would lead to more mergers, thus a variation in stellar evolution physics can also affect the choice of the aforementioned percentage between population studies.
Comparing the results presented in Fig. 4 to the case where BBHs evolve in a vacuum (black curve) shows the effectiveness of an ISM surrounding in hardening the BBH orbit.At the same time, the merger rate curve also approximately follows the expected trend.The shape of the curve and the strength of the merger rate are dependent on the value of ⟨  ⟩.Depending on the redshift evolution of the latter quantity, we find that an ⟨  ⟩ 0 ≈ 0.1−2 cm −3 results in a fair approximation to the LVK inferred rate when  = 1.However, a smaller ⟨  ⟩ 0 , e.g., ⟨  ⟩ 0 ≲ 0.1 leads to deviation from the LVK inferred trend.Similarly for  = 3, ⟨  ⟩ 0 ≈ 0.01 − 1 cm −3 also results in a fair approximation to the LVK inferred rate with ⟨  ⟩ 0 ≲ 0.01 resulting in deviation.On the other hand, large values of ⟨  ⟩ 0 (including some shown in Fig. 4) may be disfavored on physical grounds as that would require a strong and continuous supply of ISM and such ⟨  ⟩ 0 could lead to significant ISM heating (see Section 2.3, especially for massive black holes.

Detectable BBH mergers
The detectable number of BBH mergers as a function of their source frame mass is shown in Fig. 5 (for more details on the employed methodology, see Appendix A3).To assist in readability, only a subset of the ⟨  ⟩ 0 values shown in Fig. 4 are plotted in Fig. 5.The figure shows the number of detections that one may expect to observe on Earth in a year, given all three LVK detectors work in quadrature for a full year at the third observing run (O3) sensitivity with the noise curves taken from LVK-Collaboration (2020).We note that in reality, the KAGRA detector only became operational in the last phase of the third observation run.
We find that there is a tendency for the ISM-accreting BBHs to be relatively more massive than the preexisting observations.One way to prevent an excess of such systems could be by assuming a lower value of ⟨  ⟩ 0 for more massive black holes -or employing the mass-weighting approach presented in Eq. 44.This is because, due to a limited supply of matter, the surroundings of more massive black holes become ISM-poor at a faster pace (as massive binaries accrete at a faster rate, see Eq. 7).Moreover, for large ⟨  ⟩ 0 values, more massive black holes would generate stronger radiation feedback on their surrounding ISM.Thus, this should result in smaller ⟨  ⟩ 0 values being more favorable, which in turn will lower the mass of such black holes.
A large fraction of these ISM-accreting BBHs are sufficiently massive that one or both of their components fall within the pairinstability mass-gap (Fowler & Hoyle 1964).This mass-gap arises when massive stars develop the right conditions for efficient  −  + pair-production, which leads to a complete disruption of the star.As such, no black hole remnant would be left behind.Finally, if the ISMaccreting black holes do not gain significant angular momentum, they would develop relatively small values of the dimensionless spin parameter as compared to black holes living in a vacuum.This is in line with a majority of LVK observations to date (Abbott et al. 2021(Abbott et al. , 2023)).However, this might not remain true if the infalling gas acquires angular momentum for the BBH orbital dynamics, which is subsequently accreted by the black holes.There are also suggestions that while the matter gets accreted from the disk, the presence of large-scale magnetic fields near the horizon of the black hole could result in the formation of bipolar jets.The jets can feed on the spin angular momentum of the black hole, thus resulting in a spin down (e.g., Tchekhovskoy et al. 2011).

What if the BBHs have a nonzero 𝑣 ∞ ?
Till now, we have assumed that, for the most part, the BBHs are at rest w.r.t.their surrounding ISM.Here, we qualitatively discuss the impact of a nonzero  ∞ on our analysis.

Impact on mass
From Eq. 7 it can be seen that a BBH moving with a nonzero  ∞ w.r.t. the surrounding medium would experience a reduction in the mass accretion rate.If we assume that  ∞ =  ∞ , where  ∈ R, then it follows that In the case where the motion is subsonic (i.e.,  < 1), Eq. 45 suggests an utmost reduction in  by a factor of ≈ 2.8.As such, there is a fair possibility that the masses of subsonically moving BBHs could be influenced by the ISM 8 .On the other hand, even for a mild supersonic velocity of  ∞ = 2 ∞ ,  falls by a factor of ≈ 11.As such, BBHs moving at relatively higher velocities experience a diminishing impact from the ISM.

Impact on the BBH's orbital decay
In addition to the intensification of gravitational wave emission resulting from an increasing mass, e.g., see Section 3.1, earlier, we observed that isotropic accretion at  ∞ = 0 further accelerates the orbital decay process (Fig. 2).This is due to the accreted matter not contributing towards the BBH's angular momentum.
In the case where  ∞ ≠ 0, the rate of change of angular momentum  of a mass accreting binary (outside of gravitational wave emission) can be estimated as where  is the distance from the binary's center of mass.Now, depending on the magnitude and sign of the cross-product term, the accretion of mass will either contribute to or remove angular momentum from the binary.However, as the orbital velocity of LVK black holes is much larger than a subsonic  ∞ , thus averaged over an orbit, the contribution from ISM accretion to  would be negligible.Similarly, as we have already discussed in Section 6.1.1,since the mass accretion rate drastically drops for supersonic motion, as such, we expect a negligible change in  for large values of  ∞ as well.

What if 𝑣 ∞ = 0 but the accretion flow contains angular momentum?
It is possible that even though a BBH is comoving w.r.t.its surrounding ISM, the medium still possesses some angular momentum w.r.t. the center of mass of the BBH (in contrast to the disk around individual black holes which can extend to very small radii, e.g, Section 2.3, the circumbinary disk will likely disrupt as it reaches either of the black holes in the binary.Due to its large truncation radius, such a disk will not have significant radiation feedback but may heat the environment if it forms mini-disks around the individual black holes).Such a medium could comprise a generic ISM and/or the circumstellar medium, the latter being the leftover matter from the earlier episode of stellar evolution (e.g., sourced from stellar winds or due to Roche Lobe overflow from the L2 Lagrangian point).
For any such generic angular momentum possessing material to be accreted, it would need to be lowered to the location of the black hole.At this point, the specific angular momentum of the accreted fluid element would be the same as the specific angular momentum of the black hole w.r.t. the center of mass of the binary.As such, this will not affect the separation of the binary.On the other hand, it will make the black holes more massive, thus accelerating the gravitational radiation emission-induced inspiral.
The orbital decay of a BBH undergoing such an accretion process (assuming that the material is efficiently lowered to the location of the black hole such that Eq. 7 still remains a fair estimate) can be described as (cf.Eq. 34) The resulting solution is illustrated in Fig. B2 and shows that merely the growth in black hole mass (with   AM = 0) can considerably impact the merger time scale of black holes.

CONCLUSION
Many studies attempting to reproduce the LVK inferred BBH properties assume that the black holes live in a vacuum.While this assumption is justified for a short evolutionary timescale, for longer timescales of O (10 9 ) yr (typical for LVK BBHs), one might expect some deviations.This is because the BBHs are embedded in an ISM environment whose impact on the BBH could accumulate over a span of few billion years.
Assuming that for the most part, at least a few percent of BBHs in our study are immersed in a a favorable cold ISM environment, here we discussed the impact of the latter on the BBH dynamics.We showed that if the binary as a whole comoves w.r.t.its surrounding matter, then such a medium could act as an important driver of a subset of LVK BBH mergers.Additionally, it could also push some high-mass BBHs in the pair-instability mass gap.Apart from the above assumptions, another free parameter in the present work is the time-averaged particle number density ⟨  ⟩ of the surrounding ISM.While we approximately reproduce some of the LVK observables by fine-tuning ⟨  ⟩, it is unlikely that this choice would hold in a realistic setting.However, it could be possible that a certain subset of LVK's BBHs are impacted by their environment (including their dark matter surroundings), and we hope this work will shed light on the effect of such a medium on the BBH dynamics.

Figure 1 .
Figure 1.The mass evolution of a black hole embedded in a cold neutral medium of  ∞ = 100 K (see Eq. 10).The line styles represent the timeaveraged particle number density of the CNM.The line color distinguishes black holes with different initial masses   .See Fig. B1 for the case of a larger  ∞ .The case with   = 2 can be thought to instead represent a neutron star that accretes in the same fashion as a (horizon-possessing) black hole.

Figure 2 .
Figure 2. The absolute rate of orbital decay of equal component mass BBHs at different separations.The line styles represent the two masses considered.The line color represents the physical mechanism driving the orbital decay.Here, we consider a particle density of ⟨  ⟩ = 1 cm −3 .The two line styles for the blue curve (i.e., dynamical friction) apparently overlap due to the low resolution on the y-axis.

Figure 3 .
Figure 3.The merger time of equal component mass ISM accreting BBHs as a function of their (purely) gravitational wave emission induced merger time (see Eq. 42).The colors represent the different initial black hole masses in the binary (with mass-ratio = 1).The line styles represent the time-averaged particle density of the surrounding ISM.The ISM temperature is set to  ∞ = 100 K.

Figure 4 .
Figure 4.The source frame merger rate for the various employed ⟨  ⟩ 0 values.The top figure is for a linear redshift dependence (seeEq.43), while the bottom is for a cubic dependence.The top figure assumes only 1.5% of the BBHs have an optimal ISM surrounding, while the bottom assumes this to be 3%.These choices have been made to yield merger rates within the green band.The green curve shows the expected BBH merger rate as inferred by LVK, while the shaded region shows its 90% credible interval (data taken from TheLIGO Scientific Collaboration et al. 2021).The black curve is the BBH merger rate calculated assuming all BBHs live in a vacuum.

Figure 5 .
Figure 5.The detectable number of mergers on Earth in a year (as a function of the BBH mass) assuming that all three LVK detectors work in quadrature for a full year at the third observing run (O3) sensitivity.The legend labels are similar to those in Fig. 4. The grey region represents the BBH mergers observed by the LVK in O3.The top figure is for linear, while the bottom is for cubic redshift dependence of ⟨  ⟩.Similar to Fig. 4, we assume that only 1.5% (top figure) and 3% (bottom figure) BBHs have an optimal ISM environment.

Figure B2 .
Figure B2.Same as Fig. 3 but now assuming that the accreted mass contains the same amount of specific angular momentum as possessed by the black hole w.r.t. the center of mass of the binary.
47) which gives us the merger time relation as ∫  tot +