Binary mergers in the centers of galaxies: synergy between stellar flybys and tidal fields

Abstract

Galactic centers are very dynamically active environments, often harbouring a nuclear star cluster and supermassive black hole at their cores. Binaries in these environments are subject to strong tidal fields that can efficiently torque its orbit, exciting near unity eccentricities that ultimately lead to their merger. In turn, frequent close interactions with passing stars impulsively perturb the orbit of the binary, generally softening their orbit until their evaporation, potentially hindering the role of tides to drive these mergers. In this work, we study the evolution of compact object binaries in the galactic center and their merger rates, focusing for the first time on the combined effect of the cluster’s tidal field and flyby interactions. We find a significant synergy between both processes, where merger rates increase by a factor of $\sim 10-30$ compared to models in which only flybys or tides are taken into account. This synergy is a consequence of the persistent tides-driven eccentricity excitation that is enhanced by the gradual diffusion of $j_{z}$ driven by flybys. The merger efficiency peaks when the diffusion rate is $\sim 10-100$ slower than the tides-driven torquing. Added to this synergy, we also find that the gradual softening of the binary can lift the relativistic quenching of initially tight binaries, otherwise unable to reach extreme eccentricities, and thus expanding the available phase-space for mergers. Cumulatively, we conclude that despite the gradual softening of binaries due to flybys, these greatly enhance their merger rates in galaxy centers by promoting the tidal field-driven eccentricity excitation.

1 Introduction

Most nearby galaxies are known to harbor very dense stellar clusters at their core (Neumayer et al., 2011; Turner et al., 2012; Georgiev & Böker, 2014). These nuclear star clusters (NSCs) tend to have masses of roughly $10^{5}-10^{8}M_{\odot}$, a small effective radii of only a few parsecs, and usually have a supermassive black hole (SMBH) residing at their center. This is the case for our own galaxy, containing an NSC with a mass of $\sim 3\times 10^{7}M_{\odot}$ (Schödel et al., 2014, 2018) and the central SMBH SgrA^∗ with a mass of about $4\times 10^{6}M_{\odot}$ (Ghez et al., 2005; Gillessen et al., 2009). This makes the galactic center (GC) a very dense and dynamical environment, which makes it an ideal laboratory to test several astrophysical phenomena involving stellar dynamics, general relativity (GR), and gravitational waves (GWs).

In the last decade there has been an increasing interest in gravitational-wave astrophysics, due to GW detections from the LIGO and VIRGO interferometers. Compact object binary mergers are known to be one of the main GW emitting events (Abbott et al., 2016, 2017), therefore recent studies have focused on the dynamical mechanisms that lead to these mergers to better understand the sources of GWs. Particularly in GCs, one of the most important physical processes leading up to mergers is the torquing of binary’s orbits produced by the SMBH, process known as the von Zeipel-Lidov-Kozai (ZLK) mechanism (Antonini & Perets, 2012; Stephan et al., 2019; Hoang et al., 2019; Petrovich & Antonini, 2017). More generally, we can consider the tidal field from both the SMBH and the cluster which produces a behaviour similar to the ZLK mechanism (Petrovich & Antonini, 2017; Hamilton & Rafikov, 2019; Bub & Petrovich, 2020; Hamilton & Rafikov, 2021).

The strong tidal field due to the SMBH and NSC potentials has demonstrated to play a key role in the evolution of binaries in the GC, as it efficiently torques the orbits of these binaries resulting in near unity eccentricities that lead to mergers. Recent efforts have been made in accounting for different tidal field-driven physical processes in binary evolution such as cluster triaxiality (Bub & Petrovich, 2020) and relativistic phase space diffusion (Hamilton & Rafikov, 2023), however none has focused on combining these processes with others such as close encounters with field stars (flybys). Although the dynamics of wide binaries when under the influence of tidal fields (Petrovich & Antonini, 2017; Hamilton & Rafikov, 2019a; Rasskazov & Rafikov, 2023; Modak & Hamilton, 2023) and flybys (Collins & Sari, 2008; Leigh et al., 2017; Michaely & Perets, 2020; Hamilton & Modak, 2023) have been studied separately, the potential synergy between both effects has yet to be investigated in more depth.

Due to the high velocity dispersion $\sigma$ in the GC, most binaries in the GC are soft (e.g., $\sigma$ is larger than their Keplerian velocity) (Stephan et al., 2016; Alexander & Pfuhl, 2014; Rose et al., 2020; Gautam et al., 2024). One single encounter can potentially reset a binary’s orbit, changing its orientation and even pushing it towards extreme eccentricities. The semi-major axis of these binaries will also tend to increase due to these encounters (Heggie, 1975; Hills, 1975a, b). This orbit widening limits the efficiency for tidal field driven mergers, and previous works have simply prescribed it as a limiting evaporation timescale $t_{\text{evap}}$ (e.g., Stephan et al. 2016; Petrovich & Antonini 2017; Hoang et al. 2018), using $t_{\text{evap}}$ from Binney & Tremaine (2008). However, beyond the orbit softening, flyby interactions can produce a random walk in the binary’s angular momentum $\vec{j}$ which may lead to mergers as studied by Michaely & Perets (2020). A study looking into the effect from flybys in the cluster-driven dynamics is largely missing in this context (GCs) and the focus of our work.

The presented dynamics is reminiscent of that from wide binaries in general, including comets in the Oort Cloud. In particular, (Heisler & Tremaine, 1986) studied how the galactic tidal field may torque the orbit of binaries and found that flybys can contribute to the diffusion in $j$, becoming a considerable contribution to the loss rate of comets in the Oort Cloud. In contrast, as we show in this paper, binaries in GCs have slow diffusion rates compared to the torquing rates that drastically promote the merger fraction.

In this paper we study the binary evolution due to both stellar flybys and tidal fields, proving that by combining both mechanisms can drive these binaries in galactic centers to extreme eccentricities much more efficiently than when these processes are considered separately. The paper is organised as follows. We introduce the basic dynamics for binaries in tidal fields in section §2 and discuss the specific case of the GC including flybys in section §3. We present our results from numerical simulations in section §4 and we discuss the synergy found between tidal fields and flybys in §5, along with the effects on relativistic quenching in section §6. Finally, we present an overall discussion in section §7 and summarize and present the main conclusions in section §8.

2 Dynamical evolution due to tidal fields

As binaries are usually immersed in stellar clusters, an important aspect of their dynamical evolution will be driven by the smooth tidal field of said cluster. These tidal fields drive secular oscillations in a binary’s eccentricity and inclination, similar to the ZLK mechanism. In this section we model this evolution mathematically, following the work done by Hamilton & Rafikov (2019a, b). We will refer to the compact object binary system as the inner binary, while the orbit of this system around the central SMBH will be referred to as the outer binary.

2.1 Equations of motion

When studying orbits in axisymmetric potentials, the evolution of the inner binary can be approximated by averaging the tidal tensors over several periods of the outer binary. This is known as “torus-averaging”; when the outer orbit of the binary is not closed, it will trace a non-repeating path around the cluster which will, after many orbits, densely fill an axisymmetric torus. For any function that follows the outer orbit of the binary, a weighted volume average over this torus can be used instead of a time average.

As shown by Hamilton & Rafikov (2019a), in this torus-averaged limit the tidal tensor of a given potential $\Phi$ only has diagonal terms, and $\langle\Phi_{xx}\rangle=\langle\Phi_{yy}\rangle$. Here it is useful to introduce Milankovitch’s eccentricity and dimensionless angular momentum vectors $\{\vec{e},\vec{j}\}$, Milankovitch’s eccentricity and dimensionless angular momentum vectors (e.g., see Tremaine 2023), in order to simplify notation. For a smooth potential $\Phi$ that changes over scales that are much larger than that of the inner binary’s semi-major axis (expanded up to the quadrupole tidal field), this results in the torus-averaged potential

$$\langle\Phi\rangle=\frac{3a_{\text{in}}^{2}}{2}\langle\Phi_{zz}+\Phi_{xx}\rangle\left[\frac{1}{2}\Gamma(5e_{z}^{2}-j_{z}^{2})+\frac{1}{4}e^{2}(1-5\Gamma)\right],$$ (1)

where

$$\Gamma\equiv\frac{\langle\Phi_{zz}-\Phi_{xx}\rangle}{3\langle\Phi_{zz}+\Phi_{xx}\rangle}.$$ (2)

From this potential, a doubly-averaged Hamiltonian can be written in terms of $\Gamma$ as $H=H_{K}+H_{1}$ where $H_{K}$ is the Keplerian term and the dimensionless perturbing term is given by (Hamilton & Rafikov, 2019b):

$$\begin{split}H_{1}=2\Bigl{[}1-\frac{3}{2}e^{2}(5\Gamma-1)-3\Gamma(\vec{j}\cdot\hat{j}_{\text{out}})^{2}+15\Gamma(\vec{e}\cdot\hat{j}_{\text{out}})^{2}\Bigr{]}.\end{split}$$ (3)

The associated dimensionless time is $\tau=t/\tau_{\text{sec}}$, where $\tau_{\text{sec}}$ is the secular timescale given by

$$\tau_{\text{sec}}^{-1}=\frac{3a_{\text{in}}^{3/2}}{2\sqrt{GM_{\text{bin}}}}\langle\Phi_{zz}+\Phi_{xx}\rangle.$$ (4)

The equations of motion for the binary that we use throughout this work can be written in a vectorial form as:

$$\begin{split}\frac{d\vec{e}}{d\tau}&=\frac{1}{2}(5\Gamma-1)(\vec{j}\times\vec{e})-5\Gamma(\vec{e}\cdot\hat{j}_{\text{out}})(\vec{j}\times\hat{j}_{\text{out}})\\ &\qquad+\Gamma(\vec{j}\cdot\hat{j}_{\text{out}})(\vec{e}\times\hat{j}_{\text{out}}),\\ \end{split}$$ (5)

$$\begin{split}\frac{d\vec{j}}{d\tau}&=\Gamma(\vec{j}\cdot\hat{j}_{\text{out}})(\vec{j}\times\hat{j}_{\text{out}})-5\Gamma(\vec{e}\cdot\hat{j}_{\text{out}})(\vec{e}\times\hat{j}_{\text{out}}),\end{split}$$ (6)

where we shall assume that $\hat{j}_{\text{out}}$ remains fixed throughout the evolution. The latter approximation assumes that either the tidal field is spherically symmetric or the outer orbit lies in the meridional plane of an axisymmetric field (Petrovich & Antonini, 2017).

As discussed by Hamilton & Rafikov (2019b, a), $\Gamma$ determines the phase-space structure of the Hamiltonian, and hence is a useful parameter to classify different dynamical regimes for binary evolution. For example, $\Gamma=1$ when the potential is Keplerian (only a SMBH at the center), and hence the Hamiltonian in equation (3) reduces to the quadrupole ZLK Hamiltonian. The secular timescale reduces to the well-known ZLK timescale (e.g., Holman et al. 1997):

$$\tau_{\rm ZLK} = \frac{4a^{3}_{\text{out}}(1-e_{\text{out}}^{2})^{3/2}M^{1/2}_{\text{bin}}}{3a_{\text{in}}^{3/2}M_{\text{BH}}G^{1/2}}$$ (7)

$$\simeq 3\times 10^{5}\left(\frac{M_{\rm bin}}{10M_{\odot}}\right)\left(\frac{M_{\text{BH}}}{4\times 10^{6}M_{\odot}}\right)^{-1}$$

$$\times \left(\frac{a_{\text{in}}}{25\rm au}\right)^{-3/2}\left(\frac{a_{\text{out}}(1-e^{2}_{\text{out}})^{1/2}}{0.3\rm pc}\right)^{3}\text{year}.$$

We note that the ZLK Hamiltonian assumes the test particle approximation and sets $e_{\text{out}}=0$, so that the quadrupole level expansion is a good approximation of the secular dynamics. If, instead, the outer orbit is eccentric, the $z$-component of the angular momentum is no longer conserved, and hence the octupole expansion is required (Naoz, 2016). This expansion allows for qualitatively different behavior, where extremely high eccentricities for the inner orbit can be achieved along with orbital flipping (Naoz et al., 2011, 2013; Hoang et al., 2018).

Similarly, we can obtain the Hamiltonian for a binary evolving in a thin Galactic disk by replacing $\Gamma=1/3$ (see discussion in §7.1), while realistic spherical potentials can be reproduced with $0\leqslant\Gamma\leqslant 1$ as shown by Hamilton & Rafikov (2019a).

An important property of the secular Hamiltonian in Equation (3) is its axissymmetry (i.e., invariability to rotations around $\hat{j}_{\text{out}}$ ), which implies the $z$-component of the angular momentum $\propto\sqrt{a_{\text{in}}(1-e_{\text{in}}^{2})}\cos i$ with $\cos i_{\text{in}}=\hat{j}_{\text{in}}\cdot\hat{j}_{\text{out}}$ is conserved. Thus, the Hamiltonian is integrable and we can express the maximum eccentricity due to tidal fields as (Hamilton & Rafikov, 2019a):

$$e_{\text{max}}=\sqrt{1-\frac{10\Gamma\cos^{2}{i_{0}}}{5\Gamma+1}},$$ (8)

where $i_{0}$ is the initial initial when $e_{0}^{2}$ is nearly 0.

Finally, in addition to the equations of motion (5) and (6), we add an extra term in Equation (5) to account for the relativistic precession given by

$$\frac{d\vec{e}}{dt}\Big{|}_{\text{GR}}=\frac{\dot{\omega}_{\text{GR}}}{(1-e^{2})^{3/2}}\vec{j}\times\vec{e},$$ (9)

where

$$\dot{\omega}_{\text{GR}}=\frac{3G^{3/2}M_{\text{bin}}^{3/2}}{a_{\text{in}}^{5/2}c^{2}}.$$ (10)

2.2 Merger conditions

A major application we consider in the galactic center is the evolution and potential merger of compact object binaries. When the two components of these binaries come very close together, they lose energy by GW radiation. The merger timescale for this gravitational radiation is defined by Peters (1964) as

$$\tau_{\text{GW}}=\frac{3}{85}\left(\frac{a_{\text{in}}^{4}c^{5}}{G^{3}m_{1}m_{2}M_{\text{bin}}}\right)(1-e_{\text{in}})^{7/2},$$ (11)

where $m_{1}$ and $m_{2}$ If a binary’s orbit shrinks significantly in one secular cycle, then a rapid merger will take place. For this to happen, the timescale for the change in the semi-major axis of the binary due to GW losses ($\tau_{\text{GW}}$) must be shorter than the timescale for the change in the pericentric distance $a_{\text{in}}(1-e_{\text{in}})$ due to tidal fields, which, at high eccentricities, is $\simeq\sqrt{1-e^{2}}\tau_{\rm sec}$. This results in the condition

$$1-e_{\text{merger}}\lesssim 3\times 10^{-5}\left(\frac{a_{\text{out}}(1-e_{\text{out}}^{2})^{1/2}}{0.1\text{pc}}\right),$$ (12)

for binaries with $M_{\text{bin}}=10M_{\odot}$ initially at $a_{\text{in}}=10$au when only ZLK cycles are considered (only SMBH). A similar expression can be obtained for a more general cluster potential (arbitrary $\Gamma$; see Equation 15 for a secular timescale adding a cluster with a Hernquist potential).

In order for BH (stellar) binaries to reach extremely high eccentricities of $1-e_{\text{max}}\lesssim 10^{-5}$ ($1-e_{\text{max}}\lesssim 10^{-3}$), we require that $\dot{\omega}_{\text{GR}}\tau_{\mathrm{sec}}\lesssim 0.01$ ($\dot{\omega}_{\text{GR}}\tau_{\mathrm{sec}}\lesssim 0.1$) as shown by Petrovich & Antonini (2017). In general, we will define extreme eccentricities as $1-e_{\text{max}}\lesssim 10^{-4}$ following equation (12), hence binaries that achieve these values will be considered as possible mergers.

3 Dynamics in the Galactic center

In the case of a spherical $\gamma$-family cluster with a distance scale $s$ and with a central SMBH in which we have a circular orbit, $\Gamma$ becomes

$$\Gamma=\frac{3M_{\text{BH}}+M_{\text{cl}}a_{\text{out}}^{3-\gamma}(a_{\text{out}}+s)^{\gamma-4}(3a_{\text{out}}+s\gamma)}{3[M_{\text{BH}}+M_{\text{cl}}a_{\text{out}}^{3-\gamma}(a_{\text{out}}+s)^{\gamma-4}(a_{\text{out}}+s(4-\gamma))]},$$ (13)

where we have considered the tidal averaged tidal tensors $\langle\Phi_{xx}\rangle$ and $\langle\Phi_{zz}\rangle$ as derived by Bub & Petrovich (2020).

In the simplified case that we consider a spherical Hernquist potential with $\gamma=1$ to model the center of our galaxy, this reduces to

$$\Gamma=\frac{3M_{\text{BH}}+M_{\text{cl}}a_{\text{out}}^{2}(a_{\text{out}}+s)^{-3}(3a_{\text{out}}+s)}{3[M_{\text{BH}}+M_{\text{cl}}a_{\text{out}}^{2}(a_{\text{out}}+s)^{-3}(a_{\text{out}}+3s)]},$$ (14)

where $s$ is the scale radius of the potential.

For small values of $a_{\text{out}}/s$ we have $\Gamma\approx 1$, where the SMBH dominates over the cluster and hence the dynamical evolution of the binary reduces to the ZLK mechanism as mentioned in section §2.1. The same happens when $a_{\text{out}}/s$ becomes very large; the cluster and SMBH in this limit behave as a single body so that as we get further away from the central SMBH, $\Gamma\rightarrow 1$.

The secular timescale in these kinds of environments then becomes (Bub & Petrovich, 2020):

$$\begin{split}\tau_{\text{sec}}^{-1}=&\frac{3a_{\text{in}}^{3/2}}{2\sqrt{GM_{\text{bin}}}}\Bigg{(}\frac{GM_{\text{BH}}}{2a_{\text{out}}^{3}}\\ &+\frac{GM_{\text{cl}}}{a_{\text{out}}^{\gamma}(a_{\text{out}}+s)^{3-\gamma}}\left[2-\frac{3a_{\text{out}}+s\gamma}{2(a_{\text{out}}+s)}\right]\Bigg{)}.\end{split}$$ (15)

Figure 1: Full evolution of a binary due to both tidal fields and flyby effects. The initial parameters are $a_{\text{out}}=0.3$ pc, $a_{\text{in}}=25$ AU, $e_{\text{in}}=0.1$, $\omega_{\text{in}}=\pi/4$, $\Omega_{\text{in}}=0$, $i_{\text{in}}=80$°. Each stellar component of the binary has a mass of $5M_{\odot}$, and we consider the fiducial scenario of $M_{\text{BH}}=4\times 10^{6}M_{\odot}$. The red curves represent the evolution in a flyby-free scenario, simply by integrating Equations (5) and (6). In panels (d) and (e), the shaded region contains the median and both upper and lower quartiles.

3.1 Fly-bys

The main focus of this work is to study the additional effect produced by flybys on the evolution of binaries in tidal fields. Due to the high velocity dispersion in the GC, the binaries residing there are mostly soft ($v_{\text{orb}}<\sigma$, with $\sigma$ the velocity dispersion of the cluster), hence we expect flyby encounters to produce a systematic drift in the semi-major axis of the binary. This is a consequence of the Heggie-Hills law, derived using thermodynamic arguments: "hard binaries tend to harden, whereas soft binaries tend to soften" (Hills, 1975a, b; Heggie, 1975).

In order to classify binaries between hard and soft, we can write the velocity dispersion in the galactic center as according to Kocsis & Tremaine (2011)

$$\sigma(r)=\begin{cases}280\text{kms}^{-1}\sqrt{0.1\text{pc}/r}\sqrt{1-0.035(r/0.1\text{pc})^{2.2}},\\ \phantom{280\text{kms}^{-1}\sqrt{0.1\text{pc}/r}\sqrt{1-0.03}}\text{if }r<0.22\text{pc}\\ 250\text{kms}^{-1}\sqrt{0.1\text{pc}/r},\phantom{\sqrt{01\text{pc}/r}}\text{if }r>0.22\text{pc}\end{cases}$$ (16)

The hard-soft boundary then is given by $v_{\rm{kep}}=\sigma$. For example, if we consider a twin binary with two $5M_{\odot}$ components at $a_{\text{out}}=0.3$pc , this boundary would be given by $a_{\text{in}}=0.2$au. If we increase the mass of the binary component this limit will occur at a larger value; for $10M_{\odot}$ components it occurs at $a_{\text{in}}=0.4$au, and for $50M_{\odot}$ at $a_{\text{in}}=2.1$au. These continue to be relatively small values, further reinforcing that most binaries in the galactic center can be classified as soft.

The time between flyby encounters

$$t_{\rm enc}=(n_{\star}\sigma_{p}v_{p})^{-1},$$ (17)

where $\sigma_{p}=\pi b_{\text{max}}^{2}$ is the geometric cross section considered for the flybys and $n_{\star}$ is the stellar number density of the environment. For the galactic center it is considered as $n_{\star}\approx 10^{6}\text{pc}^{-3}$ (Hamers & Tremaine, 2017) and $n_{\star}\approx 0.1\text{pc}^{-3}$ for the solar neighbourhood. It can be more accuratelly modelled by (Kocsis & Tremaine, 2011)

$$n_{\star}(r)=\frac{2.8\times 10^{6}\text{M}_{\odot}\text{pc}^{-3}}{\langle m_{\star}\rangle}\left(\frac{r}{0.22\text{pc}}\right)^{-\gamma},$$ (18)

with $\gamma=1.2$ inside 0.22pc and $\gamma=1.75$ outside 0.22pc. We take $\langle m_{\star}\rangle\approx 1M_{\odot}$ as the average stellar mass in the GC; Alexander & Pfuhl (2013) presents equations that result in $\langle m_{\star}\rangle$ ranging from $1M_{\odot}$ to $2M_{\odot}$, therefore to simplify calculations we take all flybys to have a solar mass.

4 Numerical simulations

Because the interaction time due to a flyby (order $a_{\text{in}}/\sigma$, $\sim 0.2$ years) is much shorter than the secular timescale, we evolve an individual binary using the secular equation of motions in (5) and (6) subject to discrete instances where a flyby takes place, modifying the orbital elements of the binary, and continue the secular evolution until the next flyby takes place. In other words, we switch back and forth between tidal fields and flybys.

The flybys are modelled using direct N-body integrations with REBOUND (Rein & Liu, 2012), with an encounter rate as specified in equation (17) defining $b_{\text{max}}=5a_{\text{in}}$ to avoid considering flybys that are too distant. This definition allows us to stay close to the regime in which the impulse approximation is valid, which will prove to be relevant in section §5. We repeated our fiducial numerical experiments in which we doubled $b_{\text{max}}$ to $10a_{\text{in}}$, and found that the outcome was statically similar to those with $b_{\text{max}}=5a_{\text{in}}$, plus the fact that a higher value for $b_{\text{max}}$ implies a higher frequency for flyby encounters, which is much more computationally expensive to simulate.

We establish the velocity of our flybys to be constant at $v_{p}=200$ kms^-1, considering that within 1 pc of the center of our galaxy $\sigma$ varies roughly between 100 and 300 kms^-1 as according to equation (16). These simplifications allow for controlled simulations, which makes comparisons with analytical estimates such as the ones made in section §5 easier.

Figure 1 shows the complete dynamical evolution of an individual binary with $a_{\text{out}}=0.3$pc. The red curves show the evolution of an identical binary in a flyby-free scenario, where we only get tidal field driven evolution. It is quite clear by comparing these curves with our results that the inclusion of flybys in our model results in a more dramatic evolution which involves extreme eccentricities, orbital flips and jumps between the $e_{\text{in}}-\omega_{\text{in}}$ phase-space curves.

The individual changes in $a_{\text{in}}$ presented in panel (a) follow a random walk, but show an overall tendency of making the binary wider. This is expected due to the Heggie-Hills law, as mentioned earlier. We can also easily see in panel (b) that by only taking into account tidal field driven evolution, such high eccentricities would never have been achieved by this particular binary. However, this becomes possible when including flybys in the model, which allows the binary to reach values of $1-e_{\text{in}}\approx 10^{-5}$. As discussed in section §2.2, with these extreme eccentricities a rapid merger can be produced due to GW tidal captures.

It is also interesting to note that in panel (c) we can see orbital flips, in which the orbit switches between prograde and retrograde. Considering that the Hamiltonian shown in Equation (3) is only of quadrupole order and assumes $e_{\text{out}}=0$, these flips should not be possible as $j_{z}=\sqrt{1-e_{\text{in}}^{2}}\cos i_{\text{in}}$ is fixed (e.g., Naoz 2016). However, as flybys modify both $e_{\text{in}}$ and $i_{\text{in}}$ driving a slow random walk in $j_{z}$ (see also the examples in Figure 5), the orbit is allowed to flip, a phenomenon observed during the high eccentricity episodes.

Panels (d) and (e) show the variation in $\omega_{\text{in}}$ and $i_{\text{in}}$ product of each flyby encounter. We can see that both these variations in general are very small, however, along with the variations in eccentricity due to each encounter, they are enough to make the binary’s orbit jump from one phase-space curve to another as is shown in panel (f). Since flybys not only affect the magnitude of the eccentricity but also the orientation of the orbit, we get changes both in $\omega_{\text{in}}$ and in $e_{\text{in}}$ which produce these ‘jumps’ in the phase-space portrait. In particular, in a flyby-free scenario the orbit would librate around the fixed point $\omega_{\text{in}}=\pi/2$, however by including flybys in our model the orbit generally circulates one as shown in panel (f).

4.1 Population statistics

Figure 2: CDFs of the maximum eccentricities reached in the galactic center after 10, 25, 50 and 100 initial secular timescales $\tau_{0}$, as labeled. The dashed histogram represents the CDF for a flyby free population, evolved for $100\tau_{0}$.

In order to quantify the overall effects of tidal fields and flybys on binary evolution, we generated a series of simplified population synthesis models. It is important to note that all binaries within the population are independent from one another, hence we are not considering possible effects due to binary-binary interactions. We also consider a dynamical stability criterion, so as to not consider binaries that become extremely wide and hence unstable when (Grishin et al., 2016)

$$a_{\text{in}}>0.25a_{\text{out}}\left(\frac{M_{\text{bin}}}{M_{\text{BH}}}\right)^{1/3}.$$ (19)

In Figure 2 we show the CDF of the maximum eccentricities achieved by a population of $\sim 700$ binaries in the galactic center at different stages of their evolution. The initial eccentricities of these binaries follow a thermal distribution such that $f(e)de=2ede$, and their initial inclinations are uniform such that $\cos{i_{0}}$ is also uniform between -1 and 1. The eccentricity of the outer orbits is set to 0. All binaries are twins with a total mass of $10M_{\odot}$ and are evolved for $100\tau_{0}$, $\tau_{0}$ being the initial secular timescale of each individual binary according to equation (15). We consider the fiducial scenario of our GC, where $M_{\text{BH}}=4\times 10^{6}M_{\odot}$, $M_{\text{cl}}\approx 8M_{\text{BH}}$, and $s=4$ pc. The initial values for both $a_{\text{in}}$ and $a_{\text{out}}$ are drawn from a uniform distribution, ranging between $15-50$au for $a_{\text{in}}$ and $0.1-1$pc for $a_{\text{out}}$. This means that the total population of binaries is initially soft, and that the velocity dispersion does not exceed $275$kms^-1 making our assumption of $v_{p}=200$kms^-1 appropriate.

After only a few secular timescales, the initial thermal eccentricity distribution is modified into a wider distribution, achieving higher eccentricities as expected from the tidal fields. We can also observe the cumulative nature of this type of dynamical evolution; the longer we let the binaries evolve, the higher the probability of them reaching extreme eccentricities. This is expected from a chaotic nature of the binary evolution depicted in Figure 1.

Figure 3: CDFs of the maximum eccentricities and minimum pericenter distances $q_{\text{in}}=a_{\text{in}}(1-e_{\text{in}})$ reached in the galactic center after 100 initial secular timescales $\tau_{0}$ for different dynamical evolution mechanisms: tidal fields, flybys and a combination of both.

Figure 4: CDFs of the maximum eccentricities achieved by diffusion in different orbital parameters, evolved for $50\tau_{0}$. The black dashed curve represents evolution considering diffusion in all the orbital parameters due to flybys; it is the same as the yellow curve shown in figure 2.

Figure 5: Evolution of orbital parameters for three binaries with different values of $\mathcal{R}$, all with the same initial conditions of $a_{\text{in}}=25$au, $a_{\text{out}}=0.3$pc, $e_{\text{in}}$=0.1 and $e_{\text{out}}=0$. All were evolved for $50\tau_{0}$, $\tau_{0}$ being the initial secular timescale for each binary. Note the difference in the scale of the vertical axis in the top row; as $\mathcal{R}$ increases, $a_{\text{in}}$ grows more drastically.

In Figure 3 we can see the CDFs of the maximum eccentricities achieved after $100\tau_{0}$ by three different populations of binaries. One population was evolved considering only tidal field effects, the second considering only flybys and the third taking into account both effects (the same population as shown in figure 2). The figure shows that tidal fields are the least efficient dynamical mechanism in producing extreme eccentricities; only about $\sim 1\%$ of binaries evolved with tidal fields alone reach $1-e_{\text{max}}<10^{-4}$. Flybys on the other hand prove to be slightly more effective, driving about $3\%$ of binaries in the population to said values. However, when combining both mechanisms a clear synergy is found, increasing this probability up to $\sim 34\%$. For reference, the bottom panel shows the minimum pericenter distances for the same population of binaries. Although several binaries reach very high eccentricities, as was mentioned in section §3.1 these binaries will tend to get softer, which means that $a_{\text{in}}$ and hence the pericentric distance $q=a_{\text{in}}(1-e_{\text{in}})$ will grow as well. This is where the stability criterion mentioned previously in equation (19) becomes important. As can be seen in the CDF, the extreme values of the eccentricities achieved seem to dominate over the increase in $a_{\text{in}}$, producing a considerable amount of binaries with pericentric distances smaller than $\sim 10^{-3}$au, equivalent to a merger timescale $\tau_{\rm GW}<3\times 10^{5}$ yrs (Eq. [11]). The latter may be shorter than the ZLK timescale in Equation (7) and, therefore, expected to lead to a merger within one secular cycle.

5 Tides vs. flybys

Figure 6: Merger fraction as a function of the diffusion parameter $\mathcal{R}$. The dashed lines indicate the merger fractions for the limits $\mathcal{R}=0$ (tidal fields alone) and $\mathcal{R}=\infty$ (flybys alone). The error bars are evaluated assuming a Poisson distribution as $\sqrt{N}$, $N$ being the amount of mergers for each bin.

Figure 7: Evolution of a binary that begins in a GR quenched regime in which mergers are forbidden (colored region). Flyby interactions allow this binary to become softer and move out of this regime, hence becoming a merger candidate. The top panel shows the evolution of $\epsilon_{\text{GR}}$ as according to equation (9), and the bottom panel shows the evolution of $a_{\text{in}}$ and $q_{\text{in}}$.

Figure 8: Contours of the diffusion parameter $\mathcal{R}$ in Equation (22) as a function of $a_{\text{in}}$ and $a_{\text{out}}$, considering the fiducial scenario of $M_{\text{BH}}=4\times 10^{6}M_{\odot}$, $M_{p}=1M_{\odot}$ and $M_{\text{bin}}=10M_{\odot}$. The boundary for unstable binaries is plotted as according to Equation (19), while the limit for hard binaries is defined by $v_{\text{kep}}>\sigma$. We also plot the region in which relativistic quenching ($\epsilon_{\text{GR}}>0.1$) does not allow for mergers.

In this section we analyze in more depth the diffusive effects of flybys on the binary evolution.

As discussed in the previous section, the flybys impart random changes in all orbital elements ($e-i-\omega$), but it remains unclear whether the diffusion in $e_{\rm max}$ is dominated by the behavior of a single element or it is a collective process. As a controlled experiment, we show in Figure 4 the CDFs of the maximum eccentricities achieved by populations in which flyby encounters were replaced by random changes in individual orbital elements using the statistical distribution from N-body experiments. We clearly observe that diffusion in the inclination is the main effect driving binaries to high eccentricities. This translates into diffusion of $j_{z}$ as the dominant driver of maximum eccentricity grows. In turn, the diffusion of $\omega_{\text{in}}$ and $e_{\text{in}}$ adds little compared to the case where no diffusion occurs (i.e., tidal fields only).

Guided by these numerical experiments, we define a parameter $\mathcal{D}$ that quantifies the average diffusion rate of the specific angular momentum vector due to a single flyby encounter as

$$\mathcal{D}=\frac{\langle\|\Delta\vec{j}\|^{2}\rangle_{(f,b,i,\omega,\Omega)}}{t_{\text{enc}}}.$$ (20)

This coefficient not only includes changes in the specific angular momentum $j=\sqrt{1-e^{2}}$, but also changes in its orientation which is relevant to model the changes in $j_{z}$. As explained in Appendix A.1, the $\langle\cdot\rangle$ average is carried over an ensemble of binaries with random orientations and phases relative to the stellar perturber.

From this diffusion coefficient we define a dimensionless parameter $\mathcal{R}$ that quantifies the expected diffusion in $\vec{j}$ due to flybys after a secular timescale $\tau_{0}$ has occurred as

$$\mathcal{R}=\sqrt{\mathcal{D}\tau_{0}}.$$ (21)

Thus, the limits $\mathcal{R}\rightarrow 0$ and $\mathcal{R}\rightarrow\infty$ refer to scenarios completely dominated by tidal fields and by flybys, respectively. In Appendix A.1 we show the full derivation of $\mathcal{D}$, in which we have used the impulse approximation to find an analytical expression for $\Delta\vec{j}$. From this we obtain

$$\begin{split}\mathcal{R}\simeq&0.02\left(\frac{M_{\text{bin}}}{10M_{\odot}}\right)^{-\frac{1}{4}}\left(\frac{M_{p}}{M_{\odot}}\right)\left(\frac{M_{\text{BH}}}{4\times 10^{6}M_{\odot}}\right)^{-\frac{1}{2}}\\ &\times\left(\frac{v_{p}}{200\text{km/s}}\right)^{-\frac{1}{2}}\left(\frac{a_{\text{in}}}{25\text{au}}\right)^{\frac{3}{4}}\left(\frac{a_{\text{out}}}{0.3\text{pc}}\right)^{\frac{3}{2}}\\ &\times\left(\frac{n_{\star}}{10^{6}\text{pc}^{-3}}\right)^{\frac{1}{2}}\left(\frac{f(\beta)}{1.8}\right)^{\frac{1}{2}},\end{split}$$ (22)

where $\beta=s/a_{\text{out}}$ and $f(\beta)$ is defined in equation (A8). The accuracy of this expression was verified numerically, by simulating the diffusion in a binary’s angular momentum vector due to only flyby interactions using REBOUND. Figure A2 shows that the numerical equivalent for $\mathcal{R}$ also converges to $0.02$ when considering the fiducial scenario of $M_{\text{bin}}=10M_{\odot}$, $M_{p}=M_{\odot}$, $M_{\text{BH}}=4\times 10^{6}M_{\odot}$, $a_{\text{in}}=25$au and $a_{\text{out}}=0.3$pc.

In general, we expect different types of evolution when varying $\mathcal{R}$, as it should become much more diffusive and incoherent as this parameter grows. Figure 5 shows this behavior, where we have evolved three binaries with different values of $\mathcal{R}$ by varying $M_{\text{BH}}$ in each case. All three binaries have initial parameters of $a_{\text{in}}=25$au, $a_{\text{out}}=0.3$pc, $e_{\text{in}}=0.1$ and $e_{\text{out}}=0$, and were evolved for $50\tau_{0}$ where $\tau_{0}$ is the initial secular timescale for each case according to equation (15). It is important to note that $\mathcal{R}$ increases with $\tau_{\text{sec}}$, therefore in a more diffusive scenario we will have more flybys per secular cycle as the cycles grow longer. We can see in the rightmost column how coherent cycles due to tidal fields become much more difficult to identify when $\mathcal{R}=0.06$, as opposed to the case where flybys start to become less important as shown in the first column where $\mathcal{R}=0.006$. Panel (d) more specifically shows a much less chaotic evolution in which we can easily distinguish coherent cycles with some perturbations such as extreme eccentricities and variations in the secular timescale, as opposed to panel (f) where the cycles are greatly perturbed by the flybys and hence such large eccentricities are not achieved.

The bottom row shows the $z-$component of the angular momentum $j_{z}$ which should be a conserved quantity in a flyby-free scenario, however we can see how its value drifts due to the diffusion caused by flybys, especially as $\mathcal{R}$ increases. Due to the increasingly chaotic nature of the evolution of binaries in more diffusive regimes, it becomes much harder to reach the extreme eccentricities required to produce a merger, as can be seen in the case for $\mathcal{R}=0.06$. The large amount of flybys per secular timescale constantly resets the orbital parameters of the binary, preventing it from completing regular eccentricity cycles and hence decreasing the chances of achieving eccentricities of $1-e\lesssim 10^{-4}$.

5.1 Experiments changing $\mathcal{R}$

In Figure 6 we show the merger fractions of different binary populations as a function of $\mathcal{R}$, which was varied by changing the value of $M_{\text{BH}}$. These populations were all integrated for $\sim 15$ Myr, which is roughly $50\tau_{\text{sec}}$ for a binary starting with $a_{\text{in}}=25$au at $0.3$pc considering the fiducial case of $\mathcal{R}=0.02$. We consider any binary that reaches $1-e_{\text{max}}\leq 10^{-4}$ to be a merger, as was discussed in section §2.2. As mentioned previously, a less diffusive regime comes with shorter secular timescales; therefore, by integrating a population with $\mathcal{R}=0.006$ for $15$Myr we are actually letting them evolve for $\sim 500\tau_{\text{sec}}$. However, as we keep lowering $\mathcal{R}$ by increasing $M_{\text{BH}}$, binaries start to become unstable very fast according to equation (19) and hence the likelihood of them merging before they evaporate decreases. We also move into a regime in which the binaries only experience a few flyby interactions in each cycle, which could also contribute to decreasing the probability of reaching extreme eccentricity values as occurs in figure 6. Our experiments suggest that there is an optimal value for $\mathcal{R}\sim 0.02$ at which the merger fraction peaks.

5.2 Understanding the optimal value of $\mathcal{R}$ to increase the merger fractions

We provide an analytical explanation at the order-of-magnitude level of why there is a particular value for $\mathcal{R}$ at which the merger fraction peaks as shown in figure 6.

We can write our merger condition as $j_{\text{min}}=\sqrt{1-e_{\text{max}}^{2}}$ reaching values below a critical specific angular momentum $j_{\text{crit}}$. The latter is set to $j_{\text{crit}}\simeq 0.02$ in our simulations ($1-e_{\rm crit}=10^{-4}$). In the absence of flybys ($\mathcal{R}=0$), the fraction of mergers for initially circular and isotropic orbits is simply $\sqrt{3/5}j_{\rm crit}\simeq 0.015$ in our simulations. This changes slightly to $\simeq 0.52j_{\rm crit}\simeq 0.01$ for an initial thermal distribution.

As Figure 4 shows, the dominant effect of flybys on the merger rates comes from changes in inclinations, which corresponds to changes in $j_{z}$ at fixed $j$. Binaries that reach $|j_{z}|<j_{\text{crit}}$ can merge, provided that the binary lies in that range for at least a significant fraction of a secular timescale so it can be assisted by the tidal field. In a completely diffusive scenario, $j_{z}$ undergoes a random walk due to stellar encounters such that $\langle j_{z}^{2}\rangle=\mathcal{D}t$, with $\mathcal{D}$ a constant diffusive parameter. Thus, in a relaxation timescale $1/\mathcal{D}$ binaries outside the range $|j_{z}|<j_{\text{crit}}$ will be kicked in, monotonically increasing the fraction of mergers with $\mathcal{R}\propto\mathcal{D}$ consistent with the increase in Figure 6 below $\mathcal{R}\sim 0.02$.

As we continue to increase the diffusion rate due to stellar encounters $\mathcal{D}$, the merger fraction is expected to turn over when the binaries are kicked in and out of the $|j_{z}|<j_{\text{crit}}$ range in timescale shorter than the secular time. The critical diffusion timescale $t_{\text{crit}}$ to linger in this region is given by

$$t_{\text{crit}}\sim j_{\text{crit}}^{2}/\mathcal{D}.$$ (23)

Thus, we expect flybys to most efficiently boost the merger rate when $t_{\text{crit}}\sim\tau_{0}$, a condition that can be expressed in terms of our definition of $\mathcal{R}$ as

$$\mathcal{R}=(\tau_{0}\mathcal{D})^{1/2}\sim j_{\text{crit}},$$ (24)

which would be the value of $\mathcal{R}$ for which the merger fraction peaks. In our simulations, this estimate predicts a peak value at $\mathcal{R}\sim 0.02$, which perfectly matches the observed peak in Figure 6.

However, it is important to note that our simulations are not in a fully diffusive regime and we observe sudden large jumps in $a_{\text{in}}$ (see Figure 5) due to penetrative encounters ($b<a_{\text{in}}$). Therefore, Equation (24) serves as a general approximation to where the merger fraction should peak, but can be slightly off due to the non-diffusive nature of these penetrative flybys. In order to illustrate the approximate nature of our estimates, we repeat the numerical experiments for different values of $j_{\rm crit}$ in Figure B1. We observe that even though there is a value for $\mathcal{R}$ at which the merger fraction peaks, and it decreases for lower values of $j_{\text{crit}}$, the peaks do not occur exactly at $\mathcal{R}\sim j_{\text{crit}}$.

6 Lifting the relativistic quenching from flybys

Considering the properties of the GC, it is likely that a considerable amount of binaries residing in it will be subject to relativistic quenching. This presents an obstacle when producing mergers, as it does not allow binaries to reach extreme eccentricities. We can quantify the relative importance of relativistic precession using the following dimensionless parameter $\epsilon_{\text{GR}}$ as (Petrovich & Antonini, 2017):

$$\begin{split}\epsilon_{\text{GR}}&\equiv\frac{\tau_{\text{sec}}}{\tau_{\text{GR}}}\\ &=\frac{4GM_{\text{bin}}^{2}a_{\text{out}}^{3}(1-e_{\text{out}}^{2})^{3/2}}{c^{2}a_{\text{in}}^{4}M_{\text{BH}}}\\ &\times\left[1+\frac{M_{\text{cl}}}{M_{\text{BH}}}\frac{2a_{\text{out}}^{3-\gamma}}{(a_{\text{out}}+s)^{3-\gamma}}\left(2-\frac{3a_{\text{out}}+s\gamma}{2(a_{\text{out}}+s)}\right)\right]^{-1},\end{split}$$ (25)

where binaries with values of $\epsilon_{\text{GR}}>0.1$ are considered as quenched.

When this parameter becomes smaller relativistic quenching can be lifted, allowing for a potential merger driven by the tidal fields. This can be achieved through the softening of binaries, as this will increase their secular timescale $\tau_{\rm sec}$. Conveniently, this paper has focused on the effect of flybys on binaries in the GC, which tend to soften said binaries and hence can effectively lift the effects of relativistic quenching.

In figure 7 we can see a specific example where softening driven by flybys lifts relativistic quenching and allows for a merger. The initial conditions of the simulated binary are $a_{\text{in}}=1$au and $a_{\text{out}}=0.3$pc, which makes it a soft binary severely quenched by GR with an initial dimensionless parameter $\epsilon_{\text{GR}}\sim 200$ . However, due to flyby interactions increasing $a_{\text{in}}$, $\epsilon_{\text{GR}}$ reaches values smaller than 0.1 hence lifting relativistic quenching. It is clear that while subject to GR quenching the eccentricity of the binary’s orbit remains practically constant at a low value, such that $a_{\text{in}}\sim q_{\text{in}}$. However once it escapes the GR quenching regime, extreme eccentricity oscillations are induced such that the pericentric distance reaches values as low as $10^{-3}$au, allowing for a potential merger.

7 Discussion

We have shown that stellar flybys can greatly enhance the merger rates of binaries in the centers of galaxies by assisting the eccentricity growth to extreme values that are driven by the tidal torques exerted by the SMBH and the central cluster.

More specifically, our results demonstrate that the diffusive effects of flybys and the torquing of binary’s orbits due to tidal fields produce comparable merger rates (within a factor of 3; see Figure 3). In turn, for regimes in which the diffusion rate is much slower than the tidal field torquing rate ($\mathcal{R}\sim 0.01$; see §5) flybys added to the tidal fields are able to produce a merger rate $\sim 30$ times larger that the one produced by tidal fields alone. This is due to the chaotic behavior introduced by flyby interactions, as they constantly reset the binary’s orbital parameters and consequently can produce random and dramatic changes in the orbit. Most importantly, the diffusion of the $z-$component of the angular momentum $j_{z}$.

We also find that the softening of GR quenched binaries due to flyby interactions allows them to escape the relativistic quenching regime, thus expanding the available phase-space for mergers.

7.1 Dynamical regimes in the galactic center

In Figure 8 we show the different dynamical regimes as a function of $a_{\text{in}}$ and $a_{\text{out}}$, along with a contour plot of $\mathcal{R}$. It is important to note that in this plot, Equation (22) has been adapted such that $v_{p}$, $n_{\star}$ and $f(\beta)$ are written in terms of $a_{\text{out}}$. We can see that for binaries that are extremely tight ($a_{\text{in}}\lesssim 0.1$au), $\mathcal{R}$ reaches very low values which implies an almost insignificant diffusion in the angular momentum. If these binaries are subject to relativistic quenching, even if they are soft, it is unlikely that flyby interactions will have enough impact on their orbit to push them out of this regime. Considering this, the lifting of GR quenching due to flyby interactions is probably most efficient when $a_{\text{in}}\gtrsim 0.1$ au

Our results suggest that previous studies such as Petrovich & Antonini (2017) and Hoang et al. (2018) (who considered a regime similar to our own with $a_{\text{out}}\lesssim 0.5$pc and $\epsilon_{\text{GR}}\lesssim 0.1$) have underestimated the effect of interactions with field stars, presenting it as merely a limiting evaporation timescale that should inhibit merger rates. Our findings on the other hand, suggest that flyby interactions are actually quite efficient when producing mergers, much more so than tidal fields alone. This is due mainly to the diffusion in the angular momentum vector, quantified by our dimensionless parameter $\mathcal{R}$, that for these previous works would be of order $\mathcal{R}\sim 0.01$. This, combined with the ability to lift relativistic quenching and hence increase the phase-space for mergers (which was not taken into account in the mentioned studies either), provides strong evidence for the importance of flybys in producing mergers in the GC and shows that their effects must be accounted for.

7.2 Other astrophysical environments

Our results are applicable the evolution of wide binaries in various astrophysical environments. Next we discuss our results in the context of previous results.

Solar neighborhood

The dynamical evolution of wide binaries is equivalent to that of comets in the Oort Cloud studied by Heisler & Tremaine (1986). This is a different regime than the one considered in this paper, where the evolution should be flyby dominated (higher value for $\mathcal{R}$) due to the orbits of binaries being wider and having longer secular timescales. We can write an expression for $\mathcal{R}$ in this environment as follows

$$\begin{split}\mathcal{R}\simeq&9\left(\frac{M_{\text{bin}}}{M_{\odot}}\right)^{-\frac{1}{4}}\left(\frac{M_{p}}{0.162M_{\odot}}\right)\left(\frac{v_{p}}{40\text{km/s}}\right)^{-\frac{1}{2}}\\ &\times\left(\frac{a_{\text{in}}}{2.5\times 10^{4}\text{au}}\right)^{\frac{3}{4}}\left(\frac{n_{\star}}{1.14\text{pc}^{-3}}\right)^{\frac{1}{2}}\left(\frac{\rho_{0}}{0.185M_{\odot}\text{pc}^{-3}}\right)^{-\frac{1}{2}},\end{split}$$ (26)

where we have considered the secular timescale as

$$\tau_{g}=\frac{M_{\text{bin}}^{1/2}}{2\pi\rho_{0}G^{1/2}a_{\text{in}}^{3/2}},$$ (27)

with $\rho_{0}=0.185\pm 0.02M_{\odot}\text{pc}^{-3}$ being the local density as used by Heisler & Tremaine (1986).

We can see that for an average comet in the Oort Cloud we get $\mathcal{R}\simeq 9$ as opposed to the $\mathcal{R}\simeq 0.02$ obtained for compact object binaries in the GC. This implies that in the Oort Cloud flybys dominate over tidal fields, as the diffusion rate from these interactions is $\sim 9$ times faster than the torquing due to the tidal field. This is mainly due to the difference in the amount of flybys in a secular timescale, $\tau_{0}/t_{\mathrm{enc}}$. As shown in equation (A7), we get about 68 encounters for an average binary in the GC, while in the solar neighbourhood the average comet suffers $\sim 2\times 10^{4}$ flyby encounters (considering equation (27) as $\tau_{0}$ and $t_{\mathrm{enc}}\simeq 0.015\mathrm{Myr}$ as used by Heisler & Tremaine 1986). This increase in the amount of close encounters inevitably results in more diffusion of the angular momentum $\vec{j}$, hence the higher value of $\mathcal{R}$ obtained in this regime. Consistently, the numerical experiments by Heisler & Tremaine (1986) show the secular evolution is largely detuned due to the stellar flybys, making it hard to complete a clean secular cycle in the $\omega-e$ space (see figure 4 therein). This in turn makes it difficult to achieve the high eccentricities required by a merger, as is discussed by Stegmann et al. (2024).

Globular clusters

A recent study by Rasskazov & Rafikov (2023) looked at the evolution of compact object binaries when subject to tidal fields and stellar encounters in globular clusters, similar to the work done throughout this paper. The main difference is that these globular clusters are less dense and do not contain a central SMBH as the NSCs we have considered in this study, which means that they have a much lower velocity dispersion. This results in a higher $a_{\text{in}}$ hard-soft limit for binaries, implying that many of the binaries in these environments are actually hard, unlike our work. This is particularly important when speaking of GR quenching, as we have found that flybys are able to diminish this effect only when they are capable of softening binaries. In the case of hard binaries, flybys should tend to promote GR quenching rather than lift it, further inhibiting merger rates.

By looking at the behavior in the various evolution examples in Rasskazov & Rafikov (2023), we suggest that similar to the galactic field, the diffusion parameter $\mathcal{R}$ is large (dominated by flybys). We note that equation (22) is not applicable here because the encounters are generally far from being impulsive given the low dispersion velocities. More work would need to be done to extend the analysis to these environment, including flybys that can be secular in nature (e.g., Heggie & Rasio 1996; Hamers & Samsing 2019).

8 Conclusions

In this work we have studied the combined effect of cluster tidal fields and flyby interactions on the evolution of binaries in the galactic center. Our main result is the dramatic synergy between these two physical processes at driving binaries into extremely high eccentricities ($1-e_{\text{max}}\lesssim 10^{-4}$) and their subsequent merger.

Although stellar flybys tend to increase the semi-major axis of the binaries thus inhibiting the mergers, they also give rise an stochastic evolution of other orbital elements that can dramatically shrink the pericenter distance before the binary becomes unbound. This dynamics has been previously studied showing that it could lead to mergers of wide binaries in the galactic field, consistent with our results for the galactic center including only flybys that drive mergers in few a percent of the systems. We extend these results and show that this rate increases by an order of magnitude when flybys act in concert with the static tidal field from the SMBH and the central cluster.

We find that the main cause of this synergy is the persistent tidal field-driven eccentricity excitation that is enhanced by the diffusion of $j_{z}$ due to flyby encounters. We calculated the diffusion rate of the angular momentum vector $\vec{j}$ due to flybys using the impulse approximation and assess the various regimes with a dimensionless parameter $\mathcal{R}$ (Eq. 22) that measures the expected diffusion within a secular timesale (i.e, torquing timescale due to tidal fields). Our experiments suggest that the merger rates peak in the range $\mathcal{R}\sim 0.01-0.1$, meaning when the diffusion rate is $\sim 10-100$ times more slow than the torquing due to the tidal field. This regime is typical of the central parsec of our GC (see Figure 8), while other environments such the solar neighborhood and stellar clusters tend to have much higher values of $\mathcal{R}$.

We also show that the gradual softening of binaries in the GC can lift the relativistic quenching of initially tight binaries, thus expanding the available phase-space for mergers. This is also likely to contribute to the previously discussed synergy, given that by increasing the amount of binaries that can potentially merge in a population we can expect a higher merger rate.

In summary, we conclude that despite the gradual softening of binaries due to stellar encounters, these greatly enhance merger rates in GCs by promoting the tidal field driven eccentricity excitation This behavior has been ignored in previous works and further reinforces that galactic centers are ideal environments for the production of compact object binary mergers.

We would like to thank Antranik Sefilian, Carolina Charalambous, Diego Muñoz, Mark Dodici, and Scott Tremaine for useful discussions. We are also grateful to an anonymous referee for providing helpful comments on the maniscript. CP acknowledges support from CATA-Basal AFB-170002, ANID BASAL project FB210003, FONDECYT Regular grant 1210425, CASSACA grant CCJRF2105, and ANID+REC Convocatoria Nacional subvencion a la instalacion en la Academia convocatoria 2020 PAI77200076. CH is supported by the John N. Bahcall Fellowship Fund at the Institute for Advanced Study.

Appendix A Impulse approximation

From Collins & Sari (2008), we can write the variation in the velocity (impulse) $\Delta\vec{v_{i}}$ of one component of the binary with mass $m_{i}$ ($i=\{1,2\}$) considering an impact parameter $\vec{b_{i}}$ as follows:

$$\Delta\vec{v_{i}}=\frac{2Gm_{p}}{v_{p}b_{i}}\hat{b_{i}}.$$ (A1)

The impact parameter and the velocity of the perturber must be orthogonal such that $\vec{b}_{i}\cdot\vec{v}_{p}=0$. The trajectory of the perturber is approximated to a straight line, such that $\vec{r}(t)=\vec{b}_{i}+\vec{v}_{p}t$. The perturber’s velocity is relative to only one component of the binary; therefore, we can relate both impact parameters $b_{1}$ and $b_{2}$ through the equation $\vec{b_{2}}=\vec{b}_{1}-\vec{r}+\hat{v}_{p}(\vec{r}\cdot\hat{v}_{p})$ considering that $\vec{b_{2}}$ must also be perpendicular to $\vec{v}_{p}$. The overall change in relative velocity can then be expressed as $\Delta\vec{v}=\Delta\vec{v}_{1}-\Delta\vec{v}_{2}$, where both $\Delta\vec{v}_{1}$ and $\Delta\vec{v}_{2}$ obey the equation A1.

We can define $\Delta\vec{v}$ as the sum of the velocity components perpendicular and parallel to the instantaneous orbital velocity; $\Delta\vec{v}=\Delta v_{\perp}+\Delta v_{\parallel}$. We can rewrite it using equation A1 as

$$\Delta\vec{v}=\frac{2GM_{p}}{v_{p}}\left(\frac{\hat{b}_{1}}{b_{1}}-\frac{\hat{b}_{2}}{b_{2}}\right)$$ (A2)

Figure A1: Geometry used to calculate $\hat{r}\times\hat{b}_{1}$ and $\hat{r}\times\hat{b}_{2}$. The orbit lies on the equatorial plane of a sphere for reference. The flyby trajectory is approximated to a straight line, and the gray crosses indicate the points where it intersects either the sphere or its equator.

When the perturbing potential is axisymmetric (as is the case in this scenario), the inner angular momentum and eccentricity vectors $\vec{j}$ and $\vec{e}$ must always be orthogonal such that $\vec{j}\cdot\vec{e}=0$. Therefore, if there is a variation in the eccentricity vector due to fly-by encounters, there must also be a variation in $\vec{j}$. We can write the change in the angular momentum vector as

$$\Delta\vec{j}=\frac{\vec{r}\times\Delta\vec{v}}{\sqrt{GM_{\text{bin}}a_{\text{in}}}}$$ (A3)

which, using A1, results in

$$\Delta\vec{j}=\sqrt{\frac{G}{M_{\text{bin}}a_{\text{in}}}}\frac{2M_{p}}{v_{p}}\left(\frac{\vec{r}\times\hat{b}_{1}}{b_{1}}-\frac{\vec{r}\times\hat{b}_{2}}{b_{2}}\right).$$ (A4)

A.1 Diffusion coefficient analysis

It would be of interest to analyse the relative contributions of the tidal and the impulse elements on the evolution of the angular momentum of the inner binary, for it is an indicator of how likely the binary is to merge at a given moment. To this purpose, we define a diffusion coefficient $\mathcal{D}\equiv\langle\|\Delta j\|^{2}\rangle_{(f,b,i,\omega,\Omega)}/t_{\text{enc}}$. From equation A4, we can write $\|\Delta\vec{j}\|$ as

$$\|\Delta\vec{j}\|=\sqrt{\frac{G}{M_{\text{bin}}a_{\text{in}}}}\frac{2M_{p}}{v_{p}}r\left\|\frac{\hat{r}\times\hat{b}_{1}}{b_{1}}-\frac{\hat{r}\times\hat{b}_{2}}{b_{2}}\right\|,$$ (A5)

where the term involving the cross products can be evaluated numerically by averaging over $10^{6}$ different values calculated using random impact parameters $b_{i}$ and orbital parameters $f$, $\omega$, $i$ and $\Omega$ so as to cover all possible stages of the binary’s orbit. Using this numerical average, we can write $\langle\|\Delta j\|^{2}\rangle_{(f,b,i,\omega,\Omega)}$ as

$$\begin{split}\langle\|\Delta\vec{j}\|^{2}\rangle_{(f,b,i,\omega,\Omega)}\approx&6.33\times 10^{-6}\left(\frac{M_{\text{bin}}}{10M_{\odot}}\right)^{-1}\left(\frac{M_{p}}{M_{\odot}}\right)^{2}\\ &\left(\frac{v_{p}}{200\text{km/s}}\right)^{-2}\left(\frac{a_{\text{in}}}{25\text{au}}\right),\end{split}$$ (A6)

where we have used the average $\langle r^{2}\rangle=a_{\text{in}}^{2}(1+\frac{3}{2}e^{2})=\frac{7}{4}a_{\text{in}}^{2}$, considering that $\langle e^{2}\rangle=0.5$ due to the initial thermal distribution assumed.

A.2 Number of flybys

The number of interactions in a secular timescale can be calculated as $\tau_{0}/t_{\text{enc}}$, following equations (4) and (17)

$$\begin{split}\frac{\tau_{0}}{t_{\text{enc}}}&\approx 67.5\left(\frac{M_{\text{bin}}}{10M_{\odot}}\right)^{\frac{1}{2}}\left(\frac{M_{\text{BH}}}{4\times 10^{6}M_{\odot}}\right)^{-1}\\ &\times\left(\frac{a_{\text{out}}}{0.3\text{pc}}\right)^{3}\left(\frac{a_{\text{in}}}{25\text{au}}\right)^{\frac{1}{2}}\left(\frac{n_{\star}}{10^{6}\text{pc}^{-3}}\right)\\ &\times\left(\frac{v_{p}}{200\text{km/s}}\right)\left(\frac{f(\beta)}{1.8}\right),\end{split}$$ (A7)

where $\beta=s/a_{\text{out}}$, and

$$f(\beta)=\left[\frac{1}{2}+\frac{4(1+3\beta)}{(1+\beta)^{3}}\right]^{-1}.$$ (A8)

A.3 Numerical experiments for $\mathcal{R}$

In order to confirm that the expression derived for $\mathcal{R}$ using the impulse approximation is accurate, we simulated a binary that evolves only by flyby interactions using REBOUND. From this numerical experiment, in which we considered $M_{\text{bin}}=10M_{\odot}$, $M_{p}=M_{\odot}$, $M_{\text{BH}}=4\times 10^{6}M_{\odot}$, $a_{\text{in}}=25$au, $a_{\text{out}}=0.3$pc and $v_{p}=200$kms^-1, we calculated the mean displacement of $\|\vec{j}\|$ in different time windows until finding a value at which it converges after one secular timescale $\tau_{0}$. This is shown in Figure A2, where we can clearly see how $\mathcal{R}$ converges to $\approx 0.02$ after one secular timescale, which matches the value we retrieve using equation (22) when considering the same initial conditions. In order to keep this simulation as similar as possible to the derivations using the impulse approximation, we do not consider the changes in $a_{\text{in}}$ after each encounter, as this would change the amount of flybys per secular timescale.

Figure A2: Mean square displacement of $\sqrt{\langle\|\Delta\vec{j}\|^{2}\rangle}$ as a function of time window width $w/\tau_{0}$.

Appendix B Merger fraction and $j_{\text{crit}}$

To show that the optimal value of $\mathcal{R}$ at which the merger fraction peaks is proportional to the $j_{\text{crit}}$ chosen to classify mergers, a set of populations was simulated at different values for $\mathcal{R}$. Their merger fractions for three different values of $j_{\text{crit}}$ are plotted in figure B1 as a function of $\mathcal{R}$, and we can see that for smaller $j_{\text{crit}}$ the peak value tends towards a smaller $\mathcal{R}$ as well.

Figure B1: Merger fraction as defined by different values of $j_{\text{crit}}$ as a function of $\mathcal{R}$.