Binary Evolution, Gravitational-wave Mergers, and Explosive Transients in Multiple-population Gas-enriched Globular Clusters

Gas replenishment leading to the formation of 2P stars in GCs might form a disk-like structure in the cluster nuclei (e.g., Bekki 2010; Mastrobuono-Battisti & Perets 2013).

Following Bekki (2010), we consider a flat disk, i.e., with a constant aspect ratio. We estimate the aspect ratio by h/r ∼ c_s /v_K , where ${v}_{K}=\sqrt{G({M}_{\mathrm{gas}}+{M}_{\star })/{R}_{\mathrm{core}}}$ is the typical velocity in the central parsec. The speed of sound ${c}_{s}=\sqrt{{k}_{B}{T}_{\mathrm{gas}}/\mu {m}_{p}}$ ranges between 0.1 and 10 km s⁻¹ (e.g., Bekki 2010; Leigh et al. 2013), in correspondence with the gas temperature T_gas, such that c_s ≈ 0.6 km s⁻¹ corresponds to a temperature of 100 K, which is the typical temperature in star formation areas, where μ = 2.3, and m_p is the proton mass. Exponential disk models were also considered (Hénault-Brunet et al. 2015), but here we focus on simple models.

In our fiducial model, we consider c_s = 10 km s⁻¹, unless stated otherwise. Then, the aspect ratio h/r ≈ 0.23. Following Bekki (2010), we consider a velocity dispersion of σ_disk = 10 km s⁻¹ for stars embedded in the disk. As a conservative assumption, we consider the stellar/massive objects' density in the disk to be the same as in the core, i.e., n_⋆,disk ≈ n_⋆ = 10⁵ pc⁻³. However, it should be noted that, due to GDF, stars will migrate and experience inclination damping, and the effective density in the disk is expected to be higher (e.g., Artymowicz et al. 1993; Leigh et al. 2014; Grishin & Perets 2016).

We can estimate the volume ratio between the disk and the core volume by $\pi {R}_{\mathrm{core}}^{2}h/(4\pi {R}_{\mathrm{core}}^{3}/3)\sim 0.75h\,{r}^{-1}$. Then, under the assumption that all the second-generation gas is concentrated in the disk, we get a typical gas density of ρ_g,disk ∼ 1.74 × 10⁶ M_⊙ pc⁻³. The fraction of stars in the disk will change for thinner/thicker disks correspondingly.

The evolution of binaries in disks differs in several aspects from the evolution in a spherical configuration. For our discussion, the major ones are as follows: the velocity dispersion decreases, the gas density increases, and the total number of stars contained in the disk is only the volumetric fraction of the disk compared with the volume of the spherical core. The fraction might change with time owing to the interaction with gas.

Due to interactions with other stars, hard binaries tend to get harder, while soft binaries tend to get softer (Heggie 1975); see updated discussion and overview of these issues in Ginat & Perets (2021a, 2021b). Hence, in the absence of a gaseous environment stellar dynamical hardening plays an important role in binary evolution and in catalyzing binary mergers.

3.1.1. Hard Binaries

For hard binaries, the dynamical hardening rate (up to order-unity corrections calibrated usually from numerical simulations) is given by (Spitzer 1987)

$$\begin{eqnarray}&&{\left.\frac{{{da}}_{\mathrm{bin}}}{{dt}}\right|}_{3-\mathrm{body}}=-\frac{2\pi {{Gn}}_{\star }{m}_{\mathrm{pert}}(2m+{m}_{\mathrm{pert}}){a}_{\mathrm{bin}}^{2}}{{{mv}}_{\infty }},\end{eqnarray} \tag{ 2 }$$

where we consider a binary with equal-mass components, m = m₁ = m₂, and an external perturber with mass m_pert. For interactions with other massive objects only, n_⋆ and m_pert should be taken as n_• and ${\bar{m}}_{\bullet }$, respectively.

3.1.2. Soft Binaries

A binary is called a soft binary if its energy is lower than $\bar{m}{\sigma }^{2}$. This condition sets a critical SMA,

$$\begin{eqnarray}\begin{array}{rcl}{a}_{\mathrm{SH}} & = & \displaystyle \frac{2{{Gm}}^{2}}{\bar{m}{\sigma }^{2}}\\ & \approx & 200.53\mathrm{au}{\left(\displaystyle \frac{m}{10\,{M}_{\odot }}\right)}^{2}{\left(\displaystyle \frac{43.2\,\mathrm{km}\,{{\rm{s}}}^{-1}}{\sigma }\right)}^{2}\left(\displaystyle \frac{0.5\,{M}_{\odot }}{\bar{m}}\right).\end{array}\end{eqnarray} \tag{ 3 }$$

As can be seen, massive stars tend to be hard relative to the background stars in the cluster, due to the scaling ${a}_{\mathrm{SH}}\propto {m}^{2}/\bar{m}$. Hence, one should define the hardness of massive binaries relative to both low- and high-mass stars; in particular, the latter will give rise to softer binaries. We then get the following modified expression (Quinlan 1996; Kritos & Cholis 2020):

$$\begin{eqnarray}&&{a}_{\mathrm{SH},\bullet }\approx \displaystyle \frac{{Gm}}{4{\sigma }^{2}}\approx 1.25\mathrm{au}\left(\displaystyle \frac{m}{10\,{M}_{\odot }}\right){\left(\displaystyle \frac{43.2\,\mathrm{km}\,{{\rm{s}}}^{-1}}{\sigma }\right)}^{2}.\end{eqnarray} \tag{ 4 }$$

Soft wide binaries are prone to destruction owing to encounters with other stars. The dynamical evolution of massive binaries is dominated by interactions with other massive stars, and their number density in the core is elevated owing to mass segregation (Sigurdsson & Phinney 1995).

So as to bracket the effect of softening, we consider two possibilities: (1) Softening is dominated by encounters with stellar BHs, where we assume the number density of such objects to be n _b = n_• = 10³ pc⁻³, due to mass segregation to the core, where ${\bar{m}}_{\bullet }=10\,{M}_{\odot }$ (see discussion in Miller & Hamilton 2002). (2) Softening is dominated by low-mass (0.5 M_⊙) stars, if the cluster is not well segregated, and n _b = n_⋆ = 10⁵ pc⁻³.

Hence, the typical lifetime of a soft massive binary is given by (e.g., Binney & Tremaine 2008)

$$\begin{eqnarray}&&{\tau }_{\mathrm{evap},\mathrm{massive}}\approx \frac{({m}_{1}+{m}_{2})\sigma }{16\sqrt{\pi }{n}_{b}{\bar{m}}_{b}^{2}{Ga}\mathrm{ln}{\rm{\Lambda }}},\end{eqnarray} \tag{ 5 }$$

where $\mathrm{ln}{\rm{\Lambda }}$ is the Coulomb logarithm and n_b and ${\bar{m}}_{b}$ are the number density and the mass of the background stars, respectively, and change according to our choice between (1) and (2). The separation of the widest binaries that survives evaporation until the formation of second-generation stars, signed as τ_SG, taken here to be 100 Myr, is then given by

$$\begin{eqnarray}&&{a}_{\mathrm{widest}}=\max \left\{{a}_{\mathrm{SH},\bullet },\displaystyle \frac{({m}_{1}+{m}_{2})\sigma }{16\sqrt{\pi }{n}_{b}{\bar{m}}_{b}^{2}G{\tau }_{\mathrm{SG}}\mathrm{ln}{\rm{\Lambda }}}\right\}.\end{eqnarray} \tag{ 6 }$$

For our fiducial parameters, a_widest = 24.9 au for the segregated case and 200.53 au for the non-segregated case. In principle, binaries could soften and be disrupted via encounters during the gas replenishment episode; however, the GDF hardening described in the following is more efficient at this stage. Therefore, binary evaporation due to encounters sets the stage and determines the SMA of the widest binaries at the beginning of the gas enrichment stage, but it can be neglected during the time in which binaries are embedded in gas.

In gas-rich environments, such as the 2P gas environment of multiple-population GCs/YMCs (and AGN disks), GDF can play a major role in hardening. The evolution of binaries in gaseous media has been studied over a wide range of astrophysical scales from asteroids to MBHs (as discussed in the introduction).

The effect of gas was suggested to be modeled mainly via several approaches. One suggestion is that the accretion of gas onto a binary forms a circumbinary minidisk, due to accretion to the Hill sphere. In such disks, torques similar to the ones described by type I/II migration of planets in protoplanetary disks could lead to the shrinkage of the binary SMA (e.g., Artymowicz et al. 1991; McKernan et al. 2012; Stone et al. 2017; Tagawa et al. 2020). Such migration leads to very efficient mergers, far more efficient than the case of interaction dominated by GDF, as we discuss below. However, these issues are still debated, and some hydrodynamical simulations show that such torques might lead to outward migration (e.g., Moody et al. 2019; Duffell et al. 2020; Muñoz et al. 2020), while other hydrodynamical studies indicate that in thin disks one should have inward migration (Duffell et al. 2020; Tiede et al. 2020). We do note that most studies consider initially circular orbits and generally follow circular orbits, while eccentric orbits could evolve differently, with their orbital eccentricity possibly excited into very high eccentricities, as we discuss below in the context of modeling the evolution through GDF.

Therefore, the approach on which we focus here considers the effects of GDF (Ostriker 1999). When an object has a nonzero velocity relative to the background gas, the interaction with the gas reduces the relative velocity and therefore hardens binaries (e.g., Escala et al. 2004; Baruteau et al. 2011). The binary hardening induced by GDF for the circular case, with binary components with the same mass m₁ = m₂ = m, is given by Grishin & Perets (2016),

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{{da}}_{\mathrm{bin}}}{{dt}}\right|}_{\mathrm{GDF}}=-\displaystyle \frac{8\pi {G}^{3/2}{a}_{\mathrm{bin}}^{3/2}}{\sqrt{{m}_{1}+{m}_{2}}}{\rho }_{g}(t)\displaystyle \frac{m}{{v}_{\mathrm{rel}}^{2}}f\left(\displaystyle \frac{{v}_{\mathrm{rel}}}{{c}_{s}}\right);\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&f(x)=\left\{\begin{array}{l}\displaystyle \frac{1}{2}\mathrm{log}\displaystyle \frac{1+x}{1-x}-x,0\lt x\lt 1,\\ \displaystyle \frac{1}{2}\mathrm{log}\left({x}^{2}-1\right)+\mathrm{log}{{\rm{\Lambda }}}_{g},x\gt 1\end{array}\right.,\end{eqnarray} \tag{ 8 }$$

where f is a dimensionless function derived in Ostriker (1999), and v_rel is the velocity of the binary relative to the gas, taken as the Keplerian velocity of the binary, i.e., ${v}_{K}=\sqrt{G({m}_{1}+{m}_{2})/{a}_{\mathrm{bin}}}$, which dominates the relative velocity throughout most of the evolution.

Under this assumption, Equation (7) could be written as

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{{da}}_{\mathrm{bin}}}{{dt}}\right|}_{\mathrm{GDF}}=-8\pi \sqrt{\displaystyle \frac{{{Ga}}_{\mathrm{bin}}^{5}}{2m}}{\rho }_{g}(t)f\left(\displaystyle \frac{{v}_{K}}{{c}_{s}}\right).\end{eqnarray} \tag{ 9 }$$

For massive binaries, the effect of stellar hardening will be weaker than the effect on less massive stars, as can be seen directly from Equation (2). In contrast, the effect of gas hardening increases with mass (Equation (7)). Comparison of the two shows that hardening is dominated by gas hardening rather than stellar hardening. Moreover, although the effect of GDF decreases as the binary hardens, it decays more slowly than the three-body hardening, as could be seen from the scaling ${\dot{a}}_{\mathrm{hard},\star }\propto {a}^{2}$ and ${\dot{a}}_{\mathrm{GDF}}\propto {a}^{3/2}$, and therefore GDF dominates the evolution over stellar hardening throughout the evolution. After gas depletion, three-body hardening becomes the dominant dynamical process for wide binaries, while for sufficiently small separations the evolution is GW dominated.

For stellar-mass objects GW inspiral becomes important only at very small separations and can be neglected with regard to MS (or evolved) stellar binaries that merge before GW emission becomes important. However, GW inspiral plays a key role in the evolution of binaries composed of compact objects.

For a circular binary in the quadruple approximation, the GWs' inspiral rate is given by Peters (1964),

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{da}}{{dt}}\right|}_{\mathrm{GW}}=-\displaystyle \frac{64{G}^{3}{m}_{1}{m}_{2}({m}_{1}+{m}_{2})}{5{c}^{5}{a}^{3}},\end{eqnarray} \tag{ 10 }$$

where G is the gravitational constant and c is the speed of light.

Without gas dissipation, the maximal SMA for GW merger within a Hubble time is given by

$$\begin{eqnarray}\begin{array}{rcl}{a}_{\max ,\mathrm{GW}} & = & {\left(\frac{64{\tau }_{\mathrm{Hubble}}{G}^{3}{m}_{1}{m}_{2}({m}_{1}+{m}_{2})}{5{c}^{5}}\right)}^{1/4}\\ & \approx & 0.07\,\mathrm{au}{\left(\frac{m}{10\,{M}_{\odot }}\right)}^{3/4}.\end{array}\end{eqnarray} \tag{ 11 }$$

A compact binary that is driven by GDF to separations below ${a}_{\max ,\mathrm{GW}}$ would eventually inspiral and merge, even if it survived the gas replenishment stage, and would produce a GW source.

In Figure 1 we compare the different hardening processes of binaries in gas-embedded regions. As can be seen, for large separations the evolution is dominated by the gas hardening, while for smaller separations (at late times after the gas depletion) three-body hardening and finally GWs dominate the evolution. The transition between the different regimes is determined by the gas density in the cluster, as well as stellar density. Unless stated otherwise, we consider for our fiducial model a background of stars with typical masses of $\bar{m}=0.5\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

In Figure 2 we present the evolution of binaries with an initial separation of a₀ = 1 au, due to GDF, for different ambient gas densities. The gas hardening mechanism is generally very effective and leads to binary migration to small separations within short timescales, given a sufficiently dense gaseous environment. As we discuss below, such gas-assisted evolution would then give rise to high rates of GW mergers of BH binaries, comparable to the BH merger rates inferred from the aLIGO-VIRGO-KAGRA (LVK) collaboration (Abbott et al. 2016, 2021).

Zoom In Zoom Out Reset image size

It should be noted that the gas could still dominate the evolution even after reaching a_GW, as long as the gas was not depleted and the timescale for GWs mergers is larger than the GDF-induced merger timescale. In principle, GDF-dominated evolution might even be identified in the GW inspiral (in future space missions) before the merger, under appropriate conditions, if GDF still dominates the evolution in LISA frequencies.

We find that circular binaries shrink and reach final small separations, dictated by the initial conditions, which are not sufficiently small to allow for GW emission alone to drive the binaries to merger even after a Hubble time. Nevertheless, at such a short period, these very hard binaries are more likely to merge owing to dynamical encounters in the long term compared with the primordial population of binary BHs, and they should be appropriately accounted for in simulations of GC stellar populations.

In Figure 3, we introduce the evolution of binaries with different initial separations under the combined effect of GDF, three-body hardening, and GWs. It could be seen that although the merger timescales of wider binaries are slightly larger, all the binaries are expected to merge within a Hubble time. Hence, the effect of the presence of gas in the initial stages is robust across all separations and will modify the binary population. For wide enough binaries, we enter the subsonic range. In order to avoid the discontinuity in Equation (8), we take it as a constant in a small environment around Mach 1—for ${ \mathcal M }\lt 1.01$, we consider $f({ \mathcal M })\equiv f(1.01)$, where the widest binary we consider corresponds to ${ \mathcal M }\approx 0.97$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 4 we consider different sound speeds, all of them in the supersonic regime. Higher sound speed leads to larger merger timescales, although the results are robust and do not change steeply between the different choices of sound speed in this regime.

In Figure 5, we demonstrate the dependence of gas hardening on different binary masses. As can be seen from Equation (7), lower-mass binaries harden over longer timescales, due to the dependence on the mass that scales as $\propto \sqrt{m}$, for an equal-mass binary with companions m₁ = m₂ = m. The final SMA of the binary also depends on the mass of the binary, such that more massive binaries will attain smaller final SMAs.

Zoom In Zoom Out Reset image size

Heretofore we have considered gas hardening induced by GDF. However, there are other approaches to model gas hardening.

In AGN disks, gas hardening is also modeled using processes similar to migration models in protoplanetary disks (as was suggested in the context of AGN disks; McKernan et al. 2012; Stone et al. 2017; Tagawa et al. 2020). Gas is captured in the Hill sphere of a binary and leads to the formation of a circumbinary minidisk. The disk applies a torque on the binary that leads to separation decay similar to migration type I/II in protoplanetary disks, although there were studies that pointed out that this torque could lead to a softening rather than hardening (Moody et al. 2019). Notwithstanding, we will assume that the formation of a minidisk can take place in GCs and compare the resulting hardening with our GDF model. The typical timescale for hardening due to migration torques is given by (e.g., McKernan et al. 2012), with the parameters relevant to our system,

$$\begin{eqnarray}&&{\tau }_{\mathrm{type}\ \mathrm{II}}\sim 46\,\mathrm{yr}\left(\displaystyle \frac{0.01}{\alpha }\right){\left(\displaystyle \frac{0.23}{h/r}\right)}^{2}\left(\displaystyle \frac{40{\mathrm{yr}}^{-1}}{{{\rm{\Omega }}}_{\mathrm{bin}}({a}_{\mathrm{bin}})}\right),\end{eqnarray} \tag{ 12 }$$

where α is the Shakura−Sunyaev parameter, h/r is the aspect ratio, and ${{\rm{\Omega }}}_{\mathrm{bin}}=\sqrt{G({m}_{1}+{m}_{2})/{a}_{\mathrm{bin}}^{3}}$ is the angular frequency of the binary. We adopt h/r = 0.23 and α = 0.01 as a conservative value for the viscosity parameter of the disk. We substitute the Ω_bin that corresponds to a binary with a separation of 1 au. Under these assumptions, the migration timescales, which could be used to approximate the hardening timescales, are shorter than the typical migration timescales we derived using the GDF model. These timescales are also shorter than the ones obtained in AGN disks (e.g., Stone et al. 2017; Tagawa et al. 2020), as expected. We therefore expect the merger rates we derived to be similar in this case, and even higher for the lowest gas densities models, where the rates were limited by slower hardening. There were more recent studies that suggested modified migration timescales, here taken for an equal-mass binary

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{type}\ \mathrm{II},{\rm{K}}} & = & \displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{disk}}}{{{\rm{\Sigma }}}_{\mathrm{disk},\min }}{\tau }_{\mathrm{type}\ \mathrm{II}},\\ {{\rm{\Sigma }}}_{\mathrm{disk},\min } & = & \displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{disk}}}{1+0.04K},\\ K & = & {\left(\displaystyle \frac{{m}_{1}}{{m}_{1}+{m}_{2}}\right)}^{2}{\left(\displaystyle \frac{h}{r}\right)}^{-5}{\alpha }^{-1}.\end{array}\end{eqnarray} \tag{ 13 }$$

These factors lengthen significantly the typical migration timescales, such that for our fiducial model we expect τ_typeII,K ≈ 71,515 yr. This timescale is still much shorter than the expected timescale calculated via the GDF model.

Another approach to modeling gas-induced inspirals is discussed in Antoni et al. (2019). They simulate Bondi–Hoyle–Lyttelton (BHL) supersonic flows and derive the corresponding energy dissipation, fitted to an analytical theory. While the overall gas hardening timescales could be comparable or shorter for the parameters that are in our major interest, there are significant differences in the scaling. The typical inspiral timescale is given by (Equation (52) in Antoni et al. (2019))

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{BHL}} & = & 61\,\mathrm{Myr}{\left(\displaystyle \frac{{a}_{0}}{\mathrm{au}}\right)}^{0.19}{\left(\displaystyle \frac{{v}_{\mathrm{rel}}}{100\,\mathrm{km}\times {{\rm{s}}}^{-1}}\right)}^{3.38}\\ & & \times {\left(\displaystyle \frac{20\,{M}_{\odot }}{{m}_{1}+{m}_{2}}\right)}^{1.19}\left(\displaystyle \frac{7.72\times {10}^{7}\,{\mathrm{cm}}^{-3}}{{n}_{\mathrm{gas}}}\right),\end{array}\end{eqnarray} \tag{ 14 }$$

where a₀ is the initial separation of the binary and n_gas is the number density of the gas, such that ρ_gas = n_gas m_p , where m_p is the proton mass.

Each model for gas hardening sets a different critical initial separation from which the binary will merge within a Hubble time. The timescales dictated from both the type II migration and BHL mechanism are even shorter than the ones expected by our fiducial model.

Hence, we will conclude that in all the approaches that we considered to model gas hardening the process is very efficient and leads to a robust rate of mergers, which modifies significantly the binaries' population, while the major difference between them is the time of the merger, dictated by the different gas hardening timescales.

The evolution of binaries in a gaseous medium is significantly different for noncircular binaries. Here we derive and solve the equations for an orbit-averaged eccentric evolution of an initially eccentric binary embedded in gas, but we leave a more detailed discussion on the implications for the dynamical three-body hardening of eccentric binaries to future studies.

For simplicity, we will assume that the Keplerian velocity of the binary components dominates the relative velocity to the gas, and that the gas velocity is zero relative to the center of mass of the binary. Hence, the relative velocity between the binary and the gas in the center-of-mass frame is given by

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{\mathrm{rel}}=\displaystyle \frac{{\rm{\Omega }}a}{2\sqrt{1-{e}^{2}}}\left[e\sin f\hat{r}+(1+e\cos f)\hat{\varphi }\right].\end{eqnarray} \tag{ 15 }$$

The orbit equations for the GDF for a binary with two equal masses are then given by

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\overline{{da}}}{{dt}}\right|}_{\mathrm{GDF}}=\displaystyle \frac{4{a}^{3/2}}{m\sqrt{2{Gm}(1-{e}^{2})}}\left[{F}_{r}e\sin f+{F}_{\varphi }(1+e\cos f)\right],\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\overline{{de}}}{{dt}}\right|}_{\mathrm{GDF}}=\displaystyle \frac{1}{2m}\sqrt{\displaystyle \frac{a(1-{e}^{2})}{2{Gm}}}\left[{F}_{r}\sin f+{F}_{\varphi }(\cos f+\cos E)\right],\end{eqnarray} \tag{ 17 }$$

where ${{\boldsymbol{F}}}_{\mathrm{drag}}={F}_{r}\hat{r}+{F}_{\varphi }\hat{\varphi }$, f is the true anomaly, and E is the eccentric anomaly. The orbit-averaged equations are given by

$$\begin{eqnarray}\begin{array}{cll}{\left.\frac{\bar{{da}}}{{dt}}\right|}_{\mathrm{GDF}} & = & \frac{4{F}_{0}{\left(1-{e}^{2}\right)}^{2}}{\pi {m}_{\mathrm{bin}}{{\rm{\Omega }}}^{3}{a}^{2}}{\int }_{0}^{2\pi }\\ & & \times \frac{I{df}}{{\left(1+e\cos f\right)}^{2}\sqrt{1+2e\cos f+{e}^{2}}},\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}\begin{array}{ccl}{\left.\frac{\bar{{de}}}{{dt}}\right|}_{\mathrm{GDF}} & = & \frac{{F}_{0}{\left(1-{e}^{2}\right)}^{3}}{\pi {m}_{\mathrm{bin}}{{\rm{\Omega }}}^{3}{a}^{3}}{\int }_{0}^{2\pi }\\ & & \times \frac{I(e+\cos f){df}}{{\left(1+e\cos f\right)}^{2}{\left(1+2e\cos f+{e}^{2}\right)}^{3/2}},\end{array}\end{eqnarray} \tag{ 19 }$$

where F₀ is given by ${{\boldsymbol{F}}}_{\mathrm{drag}}={F}_{0}{{\boldsymbol{v}}}_{\mathrm{rel}}I/{v}_{\mathrm{rel}}^{3}$. The orbit-averaged equations for GWs are given by

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\overline{{da}}}{{dt}}\right|}_{\mathrm{GW}}=-\displaystyle \frac{64{G}^{3}{m}_{1}{m}_{2}({m}_{1}+{m}_{2})}{5{c}^{5}{a}^{3}{\left(1-{e}^{2}\right)}^{7/2}}\left(1+\displaystyle \frac{73}{24}{e}^{2}+\displaystyle \frac{37}{96}{e}^{4}\right),\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\overline{{de}}}{{dt}}\right|}_{\mathrm{GW}}=-\displaystyle \frac{304{G}^{3}{{em}}_{1}{m}_{2}({m}_{1}+{m}_{2})}{15{c}^{5}{a}^{4}{\left(1-{e}^{2}\right)}^{5/2}}\left(1+\displaystyle \frac{121}{304}{e}^{2}\right).\end{eqnarray} \tag{ 21 }$$

In Figures 6 and 7 we introduce the evolution of eccentric binaries. In Figure 6, we present the evolution due only to GDF, and in Figure 7, we also introduce the effect of GW emission. As can be seen, the eccentricities become extremely high within short timescales, indicating that the pericenter shrinks significantly. Once the pericenters are sufficiently small, the effect of GWs becomes more significant, the orbit shrinkage is accompanied by eccentricity damping, and the binaries are driven into approximately circular orbit when entering the LIGO-Virgo-KAGRA (LVK) GW bands.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Such eccentric evolution could play a key role in the evolution of the binary populations, as eccentric binaries merge within potentially far shorter timescales than circular binaries. We note, however, that some studies of a circumbinary gas-disk evolution of binaries suggest that they are only excited to moderate eccentricities ∼0.45 (Tiede et al. 2020). Nevertheless, if binary migration occurs through such processes, the overall shrinkage is rapid irrespective of the eccentricity, leading to a fast migration timescale (see previous subsection).

We further discuss these issues, and in particular the implications for the delay time distribution of GW sources from this channel, in Section 4.4. We note that the consideration of eccentric binaries' gas hardening, little studied before, should play a similarly important role in binary evolution in AGN disks, possibly in a different manner than in cases where circumbinary disk evolution is assumed (Samsing et al. 2020; Tagawa et al. 2021).

In the following we estimate the GW merger rate of binary BHs from the gas-catalyzed channel studied here. We will consider old-formed GCs and YMCs separately, given their different formation history.

In all the models we considered for gas hardening, all the binaries are expected to merge within a Hubble time. However, different gas hardening models suggest different merger timescales. As discussed above, our GDF models suggest that eccentric binaries merge rapidly, and some of the hydrodynamical studies discussed above suggest that even circular binaries merge during the early gas phase. Since most GCs formed very early, such mergers would not be detected by VLK, given the effectively limited look-back time. However, the younger equivalents of GCs, so-called YMCs, continue to form and generally follow the star formation history in the universe. Hence, mergers in such YMCs could occur sufficiently late (and hence closer by) and be detected by VLK, and the contribution of YMCs to the total VLK rate will be the dominant one for the eccentric cases (or for all binaries, according to, e.g., the circumbinary disk migration models). It should be noted that there is observational evidence for gas replenishment also in YMCs (e.g., Li et al. 2016). If, however, gas densities are lower or the binaries are initially circular/in low eccentricity, the final SMA of the binaries could be larger, leading to longer GW merger time catalyzed by three-body hardening (driving the delay time distribution to longer timescales), in which case the contribution from old GCs would be the dominant one.

The rates as a function of the redshift change according to the geometric structure of the 2P stars. Formation of 2P stars in disks is characterized by lower velocity dispersions, which lead to earlier mergers, where for the case of spherical constellation the higher velocity dispersion leads to later mergers.

We will start by estimating the number of mergers per cluster,

$$\begin{eqnarray}&&{N}_{\mathrm{merge}}\sim {f}_{\mathrm{disk}}{f}_{\mathrm{bin},\mathrm{surv}}{f}_{\geqslant 20\,{M}_{\odot }}{f}_{\mathrm{ret}}{f}_{\mathrm{merge}}{N}_{\star },\end{eqnarray} \tag{ 22 }$$

where f_disk is the fraction of stars that reside in the disk, f_bin,surv is the fraction of binaries among massive stars that will survive stellar evolution (i.e., SNe), ${f}_{\geqslant 20\,{M}_{\odot }}$ is the fraction of stars with masses that exceed 20 M_⊙, f_ret is the retention fraction of BHs in the cluster, f_merge is the fraction of binaries that merge among the surviving binaries embedded in the disk, and N_⋆ is the number of stars in the cluster.

Following our geometrical considerations in Section 2.1, we set f_disk in the range [2%, 20%]. However, even large fractions could be taken into account if there is a significant capture of objects to the disk.

The binarity fraction of massive BHs is ∼0.7, although even higher values are quite plausible for the massive-star progenitors of BHs (e.g., Sana et al. 2012); stellar evolution may reduce this fraction to a typical value of f_bin,surv = 0.1 (e.g., Antonini & Perets 2012). We use a Kroupa mass function for the cluster, such that the fraction of stars with masses larger than 20 M_⊙ is 2 × 10⁻³ for a non-segregated environment; for segregated ones we take a fraction of 0.01. The retention fraction from the cluster is taken to be 10% (e.g., Kritos & Cholis 2020 and references therein). Taking into account the initial survival fraction of wide binaries, we consider f_merge ≈ 0.49 − 0.61 for our fiducial model. The lower value corresponds to massive background stars and the upper limit to low-mass background stars ($\bar{m}=0.5\,{M}_{\odot }$); see the discussion below Equation (6).

Following Rodriguez et al. (2016), we consider logarithmically flat distribution of initial SMA in the range [10⁻², 10⁵] au, where the lower limit is close to the point of stellar contact and the upper one to the Hill radius. It should be noted that although the choice of logarithmically flat is common, there were other choices of distribution considered, based on observational data (see Antonini & Perets 2012 for further discussion).

For our fiducial model, N_⋆ = 10⁵ and M_cluster = 10⁵ M_⊙.

In order to calculate the GW merger from old GCs, we follow the calculation of Rodriguez et al. (2016) and Kritos & Cholis (2020),

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{old}}(z)=\displaystyle \frac{1}{{V}_{c}(z)}{\int }_{{z}_{\min }}^{z}{{\rm{\Gamma }}}_{\mathrm{old}}(z^{\prime} ){n}_{\mathrm{old}}(z^{\prime} )\displaystyle \frac{{{dV}}_{c}}{{dz}^{\prime} }{\left(1+z^{\prime} \right)}^{-1}{dz}^{\prime} ,\end{eqnarray} \tag{ 23 }$$

where Γ_old is the rate of mergers in old GCs; n_old is the GC number density, which is taken to be in the range [0.33, 2.57]E³(z) Mpc⁻³ (Portegies Zwart & McMillan 2000; Rodriguez et al. 2016; Kritos & Cholis 2020); dV_c /dz is the comoving volume; and (1 + z)⁻¹ accounts for the time dilation. The comoving volume is given by Hogg (1999),

$$\begin{eqnarray}&&\displaystyle \frac{{{dV}}_{c}}{{dz}}=\displaystyle \frac{4\pi {c}^{3}}{{H}_{0}^{3}E(z)}{\left({\int }_{0}^{z}\displaystyle \frac{{dz}^{\prime} }{E(z^{\prime} )}\right)}^{2},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&E(z)=\sqrt{{{\rm{\Omega }}}_{M}{\left(1+z^{\prime} \right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}},\end{eqnarray} \tag{ 25 }$$

where Ω_K = 0, Ω_M = 0.3, and Ω_Λ = 0.7 (Planck Collaboration et al. 2016).

As a conservative estimate, we take the merger rate Γ_old to be Γ_old ∼ N_merge/τ_GC, where τ_GC is taken to be 10 Gyr. In Figure 8 we present the cumulative rate of expected mergers in old GCs (in blue). There are two types of contributions to the rate: The first are eccentric binaries, such as those with initial eccentricity of 2/3, which corresponds to the mean value of a thermal eccentricity distribution, that will merge within short timescales, i.e., with negligible delay time. These practically follow the star formation rate (SFR). In this case observed contributions are likely to rise from YMCs. The second case corresponds to low-eccentricity/circular binaries, in which there will be a delay time that corresponds to a typical time of ∼ 10⁴ Myr. These contributions will be observed in old GCs.

Zoom In Zoom Out Reset image size

In this case, the major contribution from our channel to currently observable GW sources would originate not from old GCs but from YMCs. We define a YMC as a cluster formed later than redshift 2 and mass > 10⁴ M_⊙ such that we assume for that case that the 2P formation already occurred. The formation rate of YMCs follows the SFR, which enables us to write the merger rate from YMC as (Banerjee 2021)

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal R }}_{\mathrm{young}}(z) & = & \displaystyle \frac{{N}_{\mathrm{mrg}}}{{N}_{\mathrm{samp}}}\displaystyle \frac{1}{2{\rm{\Delta }}{t}_{\mathrm{obs}}}\displaystyle \frac{{\displaystyle \int }_{{M}_{\mathrm{cl},\mathrm{low}}}^{{M}_{\mathrm{cl},\mathrm{high}}}{{\rm{\Phi }}}_{\mathrm{CLMF}}(M){dM}}{{\displaystyle \int }_{{M}_{\mathrm{GC},\mathrm{low}}}^{{M}_{\mathrm{GC},\mathrm{high}}}{{\rm{\Phi }}}_{\mathrm{CLMF}}(M){dM}}\\ & & \times \displaystyle \frac{{\displaystyle \int }_{0}^{z}{{\rm{\Psi }}}_{\mathrm{SFR}}({z}_{f}){{dz}}_{f}}{{\displaystyle \int }_{3}^{6}{{\rm{\Psi }}}_{\mathrm{SFR}}({z}_{f}){{dz}}_{f}}{\rho }_{\mathrm{GC}}.\end{array}\end{eqnarray} \tag{ 26 }$$

N_mrg is the number of mergers expected in N_samp clusters, Δt_obs = 0.15 Gyr (Banerjee 2021) is the uncertainty in the cluster formation epoch, Φ_CLMF ∝ M⁻² (e.g., Portegies Zwart et al. 2010) is the cluster mass function, and we consider $[{M}_{\mathrm{cl},\mathrm{low}},{M}_{\mathrm{cl},\max }]=[{10}^{4},{10}^{5}]\,{M}_{\odot }$ as the available mass range for YMCs and [M_GC,low, M_GC,high] = [10⁵, 10⁶] M_⊙ as the typical present-day masses for GCs. ρ_GC is the observed number density of GCs per unit comoving volume. Ψ_SFR(z) is the cosmic SFR, which is given by Madau & Dickinson (2014),

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{SFR}}(z)=0.01\frac{{\left(1+z\right)}^{2.6}}{1+{\left[(1+z)/3.2\right]}^{6.2}}{M}_{\odot }\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 27 }$$

We consider N_mrg/N_samp = N_merge and spatial densities in the range [0.33, 2.57] Mpc⁻³, following Banerjee (2021 and references therein). In Figure 8, we present the cumulative rate of expected mergers in YMCs and GCs. For YMCs, the rate follows the SFR (in general, with a small correction due to the delay time—which is short) and hence peaks in relatively low redshifts. For the eccentric case, the dominant contribution will rise from YMC, while for circular ones the dominant contribution is from GCs.

It should be noted that, in general, there could be a nonnegligible delay time for the binary merger. However, for all the parameters we checked for the disk configuration, the merger timescales are extremely short and are negligible in terms of redshifts.

The total contribution to the GW merger rate from YMCs is in the range ${{ \mathcal R }}_{\mathrm{young}}\approx [0.08,25.51]\,{\mathrm{Gpc}}^{-3}{\mathrm{yr}}^{-1}$, which intersects the expected range of LVK, i.e., ${23.9}_{-8.6}^{+14.3}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ (Abbott et al. 2021), where the range is bracketed by the models with lowest and highest rates (see Table 1).

Table 1. Rates from YMCs for Redshifts z ≤ 1, for Different Choices of Parameters

Model	${{ \mathcal R }}_{\mathrm{YMC}}(z\leqslant 1)\,({\mathrm{Gpc}}^{-3\ }\,{\mathrm{yr}}^{-1})$	Model	${{ \mathcal R }}_{\mathrm{YMC}}(z\leqslant 1)\,({\mathrm{Gpc}}^{-3\ }\,{\mathrm{yr}}^{-1})$
ρ− cs− n+	0.32	ρ+ cs− n+	2.55
ρ− cs− n−	0.08	ρ+ cs− n−	0.64
ρ− cs+ n+	3.28	ρ+ cs+ n+	25.51
ρ− cs+ n−	0.82	ρ+ cs+ n−	6.35

Note. ρ_± correspond to ρ_GC = 0.33E³(z) Mpc⁻³ and ρ_GC = 2.57E³(z) Mpc⁻³, c_s± correspond to c_s− = 1 km s⁻¹ and c_s+ = 10 km s⁻¹, and n_± correspond to high density of progenitors and low fraction of hard binaries (n₊, segregated environment) and low density of progenitors and high fraction of hard binaries (n₋, non-segregated environment). These correspond also to different fractions of soft/hard binaries; see Section 3.1.2. Here we present the rates expected for initially eccentric binaries (e.g., e₀ = 0.66).

Download table as:  ASCIITypeset image

In Table 1 we present our calculated rates for different choices of parameters. As expected, higher gas densities lead to larger merger rates and higher sound speeds correspond to thicker disks that host more stars and hence yield more mergers.

Given the early epoch of gas replenishment, gas-catalyzed mergers operate on primordial binaries in the clusters. The merging components are therefore likely distributed similarly to the primordial distribution of binary components. However, even very wide binaries can merge in this channel compared with only relatively close binaries merging in, e.g., isolated binary evolution channels for GW mergers. This could give rise to significant differences in the expected masses and mass ratios of the merger objects.

Interaction with gaseous media could excite binaries to high eccentricities, due to the dependence of the drag force on the relative velocity between the gas and the binary, which changes along the orbit such that the effect is the strongest at the apocenter. Evolution of eccentric binaries hence shortens significantly the expected merger timescales, as larger separations correspond to small pericenters, in which GWs could dominate the evolution. In this case, eccentric merger could be observed in LVK, but this is not the case for slow mergers in the circular case.

We should remark in passing on the possibility of triples. In triples, the outer component migrates faster than the inner binary, potentially leading to an unstable configuration and effective chaotic three-body interaction (see, e.g., a reversed case of triples expanding owing to mass loss, leading to similar instability, in Perets & Kratter 2012); such chaotic encounters could give rise to eccentric mergers. This possibility and its potential contribution will be discussed elsewhere.

As we showed, the presence of gas modifies the binary population in GCs. It leads to an efficient merger of binaries, together with the formation of binaries via the L2 and L3 mechanisms (which were initially used to study the formation of Kuiper-belt binaries (Goldreich et al. 2002) and recently were applied to calculate the formation rate in AGN disks; Tagawa et al. 2020).

After the gas dissipation, the initial properties of the binary, as well as the gas, dictate the final separation, to which all the binaries with initial separations larger than the final separations will converge.

Therefore, gas hardening leaves a significant signature on the binary population and its properties, which sets the ground for further dynamical processes in general and specifically for later dynamical mergers.

In addition to the contribution of the channel to the total rate of GWs, the modification of the properties of binaries (mass, separation, etc.) caused by the gas hardening sets unique initial conditions for the other GW channels. This will induce an indirect signature of the gas hardening on the expected observed mergers. We introduced analytical results that could in principle be plugged in as initial conditions for the later evolution of GCs and the dynamical channels for GW production in such environments. The binary abundance changes owing to the gas hardening, since a significant fraction of binaries could merge, while others form. Furthermore, additional L2/L3-formed binaries could participate and produce GW sources, beyond the primordial binaries considered here. Nevertheless, since stars might be far more abundant than BHs, L2/L3 processes might mostly produce mixed BH−star binaries and may not contribute to the GW merger rate, but they may form other exotic binaries such as X-ray sources, etc., and/or produce microtidal disruption events (Perets et al. 2016; disruption of stars by stellar BHs).

The gas-catalyzed dynamics discussed here could take place in any other gas-rich environments, with the proper scaling. While enhanced GW merger rates were discussed in the context of AGN disks (McKernan et al. 2012; Stone et al. 2017; Tagawa et al. 2020, and references therein), they are usually discussed in the context of the evolution of a particular binary or the overall BH merger rate. However, in those cases too, the whole binary populations of both compact objects and stars will change their properties.

A very similar process could take place for young binaries embedded in star formation regions (Korntreff et al. 2012). In this case, the effect is limited to a shorter timescale and compact objects might not yet have formed and are therefore not directly affected (but their progenitor massive stars are).

YMCs are still relatively little studied in the context of the production of GW sources, although their contribution to the total estimated rate of GWs is potentially not negligible (Portegies Zwart & McMillan 2000; Banerjee 2021). In these clusters, gas can be present up to smaller redshifts, such that the effect from the channel we suggested for GWs could potentially be observed. Hence, their overall contribution to the currently observed merger rates in LVK will be more significant (as can be seen also in Figure 8). Our rate estimates, discussed below, account for both GCs and their younger counterparts, YMCs.

All the dynamical processes that take place in the early stages of GC evolution might be affected by the presence of gas, e.g., few-body dynamics.

One aspect is that wide binaries that formed during the gas epoch are protected from evaporation by the gas hardening, as they harden within timescales shorter than the typical evaporation/ionization timescales.

GDF could also enhance mass segregation (Indulekha 2013; Leigh et al. 2014). The energy dissipation leads to a change in the velocity dispersion in short timescales, such that massive objects will fall toward the center of the cluster. Moreover, since the more massive objects are prone to merge (as can be seen from Equation (7), or visually from Figure 5), the relaxation will be affected by the modified mass function induced by the gas hardening.

5.4.1. GW Recoils, Spins, and Mass-gap Objects

The focus of the current paper is the merger of BHs and the production of GW sources due to gas interactions in multiple-population clusters. However, the evolution of stars and other compact objects such as WDs and NSs could be significantly affected in similar ways. Though some of these aspects are discussed in a companion paper (Perets 2022), we postpone a detailed exploration of these objects to a later stage and only briefly mention qualitatively some potentially interesting implications.

A fraction of the gas could be accreted on objects in the cluster. Gas accretion changes the velocities of the accretors and the overall mass function of objects in the cluster, such that there is a shift toward higher masses (e.g., Leigh et al. 2014), which might affect the dynamical GW channels in clusters that operate after the gas replenishment epoch, since we enrich the abundance of massive objects that are likely to be the progenitors of GWs. Stars that accrete gas could evolve into compact objects that in turn might produce novae. Enhanced accretion in the early stages of the cluster evolution could potentially modify the nova rates and properties (Maccarone & Zurek 2012) and the production of accretion-induced collapse of WDs into NSs (Perets 2022).

We should point out that our scenario suggests a robust merger not only of BHs but also of NSs and WDs. These mergers might leave unique signatures. Besides their contribution to the production of short GRBs and GW sources, binary NS mergers are a promising channel to the production of heavy elements via r-process (e.g., Freiburghaus et al. 1999) and would affect the chemical evolution of the clusters.

Thermonuclear explosions of WDs could produce Type Ia SNe, whether via a single-degenerate channel (WD and a nondegenerate companion; Whelan & Iben 1973) or a double-degenerate channel (two WDs; Iben & Tutukov 1984). Both of these channels will be affected by the gas accretion. First, as we mentioned (Leigh et al. 2014 and references therein), the mass function will change. This in turn might change the characteristics of the SNe and their rate. Furthermore, regardless of the mass variation, a large fraction of the compact-object binaries are expected to merge within short timescales, which will also affect the SN rate.

Mergers of WDs could yield a remnant merged object with small or absent natal kick and hence constitute another channel for NS formation. Accretion could potentially change the retention fraction and potentially explain the retention problem in the formation of pulsars (Perets 2022).

The amount and origin of gas in GCs during the formation of 2P stars are still uncertain (Bekki 2017). In this channel, we suggest that the amount of gas dictates a final SMA, such that the separation distribution/GW rate could be used to constrain the gas abundance in the cluster and its lifetime.

For sufficiently low gas densities (or lower densities following gas depletion), gas hardening is not efficient enough to lead to a merger. In this case, the terminal SMA of the binary will exceed a_GW, such that GWs will not be emitted without a further dissipation process. However, if the gas remains for longer timescales, further hardening will occur. For the whole parameter space we considered, the early stages of the hardening process are very efficient, i.e., wide binaries harden and become hard binaries on short timescales.

This channel of production of GW sources could serve as a tracer to later star formation, as it is coupled to the gas that accompanies this formation. The amount of gas and its decay with time are determined by the star formation history. Since these parameters play a role in gas hardening and hence in the final separation distribution at the end of the gas epoch, they could potentially serve to constrain the 2P gas and star formation phase and may help explain some of the differences between 1P and 2P stellar populations.

For example, we might speculate that the inferred difference between the 1P and 2P binary fractions (e.g., Lucatello et al. 2015) could be explained by gas-catalyzed hardening and mergers of MS stars residing in the gaseous region. Such 1P binaries that also accrete a significant mass of 2P gas would appear and be part of the 2P populations, while outside the gas regions binaries are not affected. In this case some of the 2P binaries preferentially merge compared with 1P stars outside the 2P gas region, leading to an overall smaller binary fraction.

That being said, the many uncertainties and degeneracies involved might be challenging in directly connecting current populations with the early conditions directly.

In the following we discuss potential caveats of our model/scenario.

• 1.  
  The specific scenario for formation of 2P stars is still unknown/debated, and hence there are large uncertainties in the amount of gas in the cluster and its source during the different stages of evolution. Moreover, some explanations for the different chemical composition of the so-called 2P stars might require lower gas masses than the total mass of 2P stars. In these cases, the phenomena we described might be somewhat suppressed, though, as we have shown that even lower gas densities could be highly effective and will not qualitatively change the results.
• 2.  
  The expected production rates of GW sources depend on the initial parameters of the clusters we consider, including the gas densities, stellar and binary populations, star formation histories, etc. All of these contain many uncertainties, which we did not directly address in this initial study, limited to a small number of models so as to provide an overall estimate to bracket the expected GW rates from this channel. Nevertheless, all of our models show that gas-catalyzed mergers in multiple-population clusters could produce a significant and even major contribution to the GW merger rate and could play a key role in the general evolution of stars and binaries in such clusters.
• 3.  
  The interaction of gas with binaries is complex and includes many physical aspects. Here we assumed that the gas density in the cluster, or at least in the region in which the binaries evolve, is spatially constant. Most of the gas should be concentrated in the star-forming region, preferentially toward the inner parts of the cluster. Outer parts of the cluster might be more dilute. Future study could relax the simplified assumption of a constant spatial density and account for a more detailed distribution of gas, stars, and binaries.
• 4.  
  We assumed that the relative velocity between the objects and the gas is dominated by the Keplerian velocity of the binary dominant. A more realistic approach, but requiring a detailed Monte Carlo or N-body simulation, could account for the detailed velocity distribution of binaries in the cluster.
• 5.  
  As we mentioned in Section 5.5, objects embedded in gas could accrete from it and change their mass over time. As a result, their dynamics will change both in the cluster and as binaries (Roupas & Kazanas 2019). Here we considered constant masses throughout the evolution and neglected the effects of gas accretion. This is a somewhat conservative assumption, in regard to catalysis of mergers, as more massive objects are prone to merge even faster in gas (see Equation (7) and Figure 5).
• 6.  
  We considered several choices for the gas depletion, assuming an exponential decay, with a fiducial model of 50 Myr and a lifetime of 100 Myr. However, the formation epochs of stars could set different scenarios, e.g., in which gas is abundant in the cluster for longer timescales of ∼100 Myr, but only intermittently (Bekki 2017), which will change the picture, or when several wide-scale gas replenishment episodes occur over timescales of even many hundreds of Myr or even Gyr, as might be the case for nuclear clusters.
• 7.  
  In our analysis we considered for simplicity only equal-mass binaries. Though we do not expect a major change in the results, the generalization to binaries with different masses is more complex and requires more detailed population studies, beyond the scope of the current study.
• 8.  
  It should be noted that there were studies that suggested more limited efficiency of GDF (e.g., Li et al. 2020; Toyouchi et al. 2020) than considered here. A more detailed comparison is left for further studies.
• 9.  
  Although the initial parameters of our disk suggest a thick disk, in later stages the disk will be thinner and finally fragment to enable star formation. Hence, for these stages/initial thin disks, the gas hardening epoch should be limited to the regime in which the disk is stable.
• 10.  
  We restrict ourselves to binaries that are not likely to be disrupted by interactions with other stars. Further disruptions could take place and are encapsulated in f_bin,surv (see Equation (22)).