Orbital evolution of LIGO/Virgo binaries in stellar clusters driven by cluster tides, stellar encounters and general relativity

Origin of LIGO/Virgo gravitational wave events may involve production of binaries with relativistic components in dense stellar systems - globular or nuclear star clusters - and their subsequent evolution towards merger. Orbital parameters of these binaries (the inner orbit) and their motion inside the cluster (the outer orbit) evolve due to both external agents - random encounters with cluster stars and cluster tides due to the smooth cluster potential - and the internal ones - various sources of dissipation and precession within the binary. We present a numerical framework - Binary Evolution in Stellar Clusters (BESC) - that follows the evolution of the binary inner and outer orbits accounting for all these effects simultaneously, enabling efficient Monte Carlo studies. The secular effect of cluster tides is computed in the singly-averaged approximation, without averaging over the outer binary orbit. As to stellar encounters, we include the effects of both close and distant flybys on the inner and outer orbits of the binary, respectively. In particular, this allows us to explicitly account for the dynamical friction sinking the binary towards the cluster centre. Also, given our focus on the LIGO/Virgo sources, we include the general relativistic precession (which suppresses cluster tides at high eccentricities) and the gravitational wave emission (shrinking the binary orbit). We use BESC to illustrate a number of characteristic binary evolutionary outcomes and discuss relative contributions of different physical processes. BESC can also be used to study other objects in clusters, e.g. blue stragglers, hot Jupiters, X-ray binaries, etc.


INTRODUCTION
At the moment of writing, 76 compact binary mergers have been detected by the LIGO/VIRGO/KAGRA collaboration (The LIGO Scientific Collaboration et al. 2021), with 70 of them classified as black hole-black 1 hole (BH-BH) mergers.For these mergers to occur in the first place, some dynamical mechanism must exist that could bring the binary components close enough for the merger timescale due to the gravitational wave (GW) emission to be at least shorter than the Hubble time.In particular, for a circular 10M ⊙ + 10M ⊙ binary the semimajor axis would have to be  ≲ 0.09 AU.
The specific mechanism responsible for reducing the compact binary semi-major axis depends on the environment in which the binary resides.In the galactic field, the compact binaries can be naturally produced as a result of stellar evolution in isolated massive binaries via the common envelope evolution (Belczynski et al. 2016); the compact binary then exists in isolation and can only shrink its orbit via the GW emission.Another scenario in the field involves hierarchical triples that may shrink their orbits via the Lidov-Kozai eccentricity oscillations accompanied by the GW emission at the pericentre (Silsbee & Tremaine 2017).
★ E-mail: rrr@damtp.cam.ac.uk 1 In this work we will be using "BH binary" and "compact object binary" interchangeably.
In denser environments, i.e. in star clusters, compact binaries can either be primordial or form in close triple or quadruple stellar encounters.Their mergers in the cluster could be facilitated by the occasional close flybys of stars or other compact objects, which may harden the binary, i.e. reduce its semi-major axis, until the GW emission takes over (Rodriguez et al. 2016;Fragione & Kocsis 2018).Such close encounters with other stars in a cluster can be interpreted as resulting from the granularity of the gravitational potential inside the cluster.
Recently, Hamilton & Rafikov (2019a,b,c) have proposed another mechanism for reducing the semi-major axis of compact binaries in star clusters.They pointed out that under certain conditions the tidal forces arising due to the large-scale inhomogeneity of the smooth component of the cluster potential are able to drive binary eccentricity to high values, at which point the GW emission could become effective at reducing the binary semi-major axis.This mechanism is similar to the Lidov-Kozai effect with the tidal field of the whole cluster playing a role of the gravitational perturbation due to the third body.Using this idea, Hamilton & Rafikov (2019c) have carried out a simplified Monte Carlo calculation of the merger rate of compact binaries due to the cluster tides.However, this work ignored all dynamical effects of encounters with other stars in the cluster that must always be present in some level.
Given the number and complexity of physical processes affecting the orbital evolution of binaries (not necessarily consisting of compact objects) in clusters, their dynamics has often been explored using numerical models with different levels of approximation.For example, the Cluster Monte Carlo Code (Rodriguez et al. 2021) is a statistical framework that considers only the effects of stellar flybys by randomly changing the energy and angular momentum of every member of a cluster at every timestep.While that method accounts for effects such as (nonresonant) relaxation and close encounters, it ignores the effect of distant encounters (with pericentre distances much higher than the binary semimajor axis) or the cluster tidal field on the inner orbital parameters of the binary.Many other studies also consider only the effects of encounters with cluster stars on the binary orbital elements, for example Hamers & Tremaine (2017)  2  and Samsing et al. (2019).At the other extreme, Hamilton & Rafikov (2019c) and Bub & Petrovich (2020) modeled exclusively the effects of cluster tides on binary orbits.Full N-body simulations of stellar clusters naturally take into account both the smooth and spatiallyfluctuating components of the cluster potential, but they are usually computationally expensive (Li et al. 2023).
The goal of this work is to present a novel numerical framework -Binary Evolution in Stellar Clusters 3 (BESC) -for efficiently following the evolution of the orbital elements of an individual binary on time intervals as long as the Hubble time, by simultaneously accounting for the effects of the smooth and fluctuating components of the cluster potential, as well as the general relativistic effects.We describe the details of the physical and numerical implementations of this method and illustrate them with several individual examples, deferring to future work the use of BESC for statistical studies of populations of binaries in stellar clusters.Also, one of the goals of BESC is to help us understand the relative importance of cluster tides and stellar encounters for the orbital evolution of a binary and to determine the conditions under which it can become very eccentric.Our paper is organized as follows.After discussing in §2 the relative roles of different cluster-specific processes in binary evolution, we describe the implementation of the key ingredients of BESC in §3: encounters ( §3.1), cluster tides ( §3.2), and GR effects ( §3.3).We then provide a comparison of our singly-averaged implementation of cluster tides in BESC with the existing doubly-averaged results ( §4.1), and then describe some representative outcomes of binary evolution in clusters ( §4.2).We discuss our results in §5 (including the preliminary statistics of outcomes in §5.1) and conclude in §6.Finally, Appendix A is devoted to motivating the parameter choices used in BESC.

PHYSICAL PROCESSES AFFECTING BINARY EVOLUTION
Binaries orbiting in dense stellar systems are subject to a range of physical processes that may affect the dynamical characteristics of the binaries.Some of these processes affect only the inner orbit of the binary, i.e. the motion of the binary components around its centre of mass (CoM), while others affect the outer orbit of the binary, i.e. motion of the CoM within the cluster, or both.The outer orbit of the binary within the cluster is fully determined by the cluster potential Φ(r) and the initial conditions -position and velocity of the CoM in the cluster.However, because of the random flybys of cluster stars past the binary its outer orbit also evolves, 2 That paper actually considers a hot Jupiter orbiting around a star rather than a BH binary. 3BESC and the relevant supporting documentation are available at https: //github.com/RZCas/binary-evolution-in-a-cluster.
see §3.1.The effect of these stellar encounters can be qualitatively decomposed into • The dynamical friction due to the smooth density field of the cluster -a secular effect ( §3.1.6), • Random changes of the binary velocity in its orbit, equivalent to stochastic changes of the initial conditions for the subsequent integration of the outer orbit.
At the same time, the inner orbit of the binary evolves due to the following effects: • Close encounters with cluster stars (as for the outer orbit), which can change not only the eccentricity and inclination but also the semimajor axis of the binary ( §3.1), • Secular evolution of the inner orbit due to cluster tides -inhomogeneity of the smooth component of the cluster potentialleaving the binary semi-major axis unchanged ( §3.2), • General relativistic (GR) effects: orbital precession and GW emission ( §3.3).
Additional possible evolutionary factor is the variation of the cluster properties on long timescales due to non-resonant and resonant relaxation.However, in this work we will neglect overall cluster evolution for simplicity, to better focus on other aspects of binary dynamics.
Table 1 provides a comparison of the physical processes accounted for in our present work and in other existing studies of the binary evolution in stellar clusters.It will be further discussed in §5.2.The implementation of the various aforementioned physical effects will be described in §3.4,while a brief discussion of their relative importance for binary evolution is provided next.

Relative importance of different physical effects
Efficiency of the physical processes listed above depends on the binary semimajor axis  and binary mass  b .To illustrate this, in Fig. 1 we show some characteristic values of the inner semi-major axes  separating different physical regimes (explained below), as a function of the binary mass  b .Our calculation assumes a cluster with a total mass  cl and characteristic radius .In all the estimates below we assume the radius of the binary outer orbit to be ∼  (which essentially makes the exact form of the cluster potential irrelevant).Different panels in Fig. 1 correspond to different values of  cl and .The meaning of the different lines is as follows: (i) Black solid: "Hard binary separation"  ℎ =  2 /(4 2 ) (Quinlan 1996) for a binary with component masses  2 <  1 ( 1 +  2 =  b ) and the local velocity dispersion  of cluster stars with mass  3 ≲  2 .This semi-major axis provides a separation between the 'hard' binaries that tend to harden (shrink their semi-major axes) in stellar encounters and 'soft' binaries, which are softened (driven to expand) by stellar flybys (Heggie 1975).As demonstrated by Quinlan (1996), binary hardening due to stellar ejections becomes efficient only at  <  ℎ , when the 3-body hardening rate d( −1 )/d becomes approximately constant.Thus,  ℎ represents an important scale separating two distinct regimes of the possible binary semimajor axis evolution due to stellar encounters.
(ii) Black dashed: the semimajor axis  3body at which the hardening timescale (Quinlan 1996) where  ≈ 15 is a constant factor and  is the local density of  cluster stars, becomes equal to the Hubble time,  ℎ =   .With this definition, 3-body hardening is an important actor and should be taken into account when following binary evolution whenever  3body <  <  ℎ .
(iii) Red: the semimajor axis  GW at which the merger timescale due to GW emission becomes equal to the Hubble time,  GW =   , assuming a circular binary.GW emission is efficient (on timescales ∼   ) for circular binaries when  <  GW ; for eccentric binaries,  GW can be much higher (see §3.3 for more details).
(iv) Green dashed: the semimajor axis  sec above which the period of secular oscillations of the orbital elements due to the cluster tides (see §3.2) becomes shorter than 0.1  ; we choose this somewhat arbitrary value to allow for at least a few oscillations to occur in the Hubble time, letting GW emission to shrink the binary orbit during the high- episodes.With this logic, secular evolution driven by cluster tides is efficient when  >  sec .
(v) Green solid: the semimajor axis   below which the binary pericentre precession due to GR effects becomes considerably faster than the cluster tide-driven precession.When this happens, the GR precession effectively suppresses eccentricity growth near its peak in the course of secular oscillations, severely reducing GW emission and semi-major axis evolution (Hamilton & Rafikov 2021).We calculate   using the constraint  GR = 20 on the dimensionless parameter  GR characterizing the strength of the GR precession, see §3.3 for details (we also set  = 0.5  tot / 3 in Eq. ( 9)).With this constraint in mind, cluster tides coupled with the GW emission would be efficient at shrinking the orbit of the binary when both  >  sec (to enable secular oscillations in the first place) and  >   .
Careful examination of the Fig. 1 allows us to make the following observations: • Cluster tides are important only at sufficiently high semimajor axes, tens to hundreds of AU, especially for higher  b .Only in densest (more massive and more compacts) clusters does  sec decrease to several AU (making cluster tides more important) for  b ∼  ⊙ , typical for neutron star-neutron star (NS-NS) binaries, see panel (c).
• The 3-body hardening is more efficient than softening at  ≲  sec , when the cluster tide effects are usually not so impor-tant.However, for lower  cl and smaller  there is a small range of  at high  b in which both hardening and cluster tides are important, i.e.where  sec <  <  ℎ .
• Hardening is much more efficient for heavier binaries, as they enter the "hard regime" earlier (since  ℎ ∝  2 ), and once  ≲  ℎ the hardening timescale  ℎ becomes independent of the binary mass, see equation ( 1).The semimajor axis range in which the 3-body hardening is important also shifts downwards in denser clusters; e.g. for • Unless the binary has very high eccentricity, 1 −  ≪ 1, GW emission is able to drive the binary to merger only at very low semimajor axes (0.01-0.1 AU) when the 3-body hardening is typically inefficient (except for the densest clusters, see panel (c)).Thus, the only two ways for the binary to achieve merger is to either (1) reach extreme eccentricity (via either cluster tides or 3-body interactions) or (2) have its semimajor axis substantially reduced as a result of a strong 3-body scattering event, such that the GW emission becomes important for shrinking its orbit.The results of both Monte-Carlo (Rodriguez & Loeb 2018) and population synthesis (Antonini & Gieles 2020) codes confirm that heavier and/or denser clusters have higher BH merger rates.
One important caveat to keep in mind is that in making this plot we neglected the stochastic component of the 3-body encounters and only took into account the average 3-body hardening rate.As we will show in §4, the real effect of random encounters is greatly enhanced by their stochastic nature, for which  ℎ and  ℎ are not very representative.Investigating the true role of stellar flybys on the binary orbital evolution will be the subject of Rasskazov & Rafikov (in prep.).

INGREDIENTS OF OUR NUMERICAL METHOD
Our numerical framework BESC is designed to follow the evolution of the outer and inner orbits of an individual binary in a stellar cluster, provided a set of binary and cluster parameters and initial conditions -masses of the components  1 and  2 , semimajor axis , eccentricity , inclination , argument of periapsis  and longitude of ascending node Ω, as well as the initial position and velocity of the binary within the cluster.BESC does not handle the production (although it accounts for their destruction, see §4.2) of binaries -they can be either primordial or formed as a result of dynamical interactions in a cluster.Note that all Keplerian elements of the inner orbit are defined with respect to the reference plane and direction inside the cluster, which are specified in the beginning of the calculation and remain fixed throughout the evolution.
We also specify the density profile of the cluster () as a function of distance from its centre , which allows us to obtain the potential Φ() by solving the Poisson equation.In this work we focus on spherically-symmetric clusters; however, our treatment of cluster tides can be trivially extended to axisymmetric (Hamilton & Rafikov 2019a) and fully triaxial (Bub & Petrovich 2020) systems.Our framework allows for a possibility of a mass spectrum of cluster stars, specified via the mass distribution function  () (normalized such that ∫  ()d = 1), but in the examples shown in this paper we will focus on a single-mass case for simplicity.
In this section we provide the details of implementation of each individual physical process accounted for by our framework ( §3.1-3.3), and then describe their synthesis into a framework for following the evolution of both the outer and the inner ( §3.4) orbits of the binary.

Random encounters with the cluster stars
Encounters with cluster stars play a very important role in binary evolution.Unlike cluster tides, they can change not only the eccentricity and inclination but also the semi-major axis of the binary.Moreover, they affect not only the inner but also the outer orbit of the binary, as well as its overall composition and fate ( §3.1.5).We will now provide a description of how we compute the rate of stellar encounters ( §3.1.1),model each individual encounter ( §3.1.2),and calculate the changes of the binary inner ( § §3.1.3, 3.1.4)and outer ( §3.1.6)orbits resulting from a stellar flyby.

Encounter rate
To determine when an encounter with a cluster star happens, we use the result of Hamers & Tremaine (2017) for the mean rate of stellar encounters with pericenter distances (relative to the binary CoM) below some chosen4  max : where  is the stellar number density,  b and  3 are the binary and the perturber masses, and  is the stellar mass distribution in the cluster; the factor in brackets accounts for gravitational focusing.Also,  rel = √  2 +  2 is an approximation for the relative velocity dispersion, where  is the cluster velocity dispersion,  is the (instantaneous) binary CoM velocity, and for simplicity an isotropic velocity distribution in the cluster frame is assumed.The isotropic velocity dispersion profile () is determined from the Jeans equation: where Φ() is the cluster potential.
If R is approximately constant along the outer orbit (which is the case when the stellar density and velocity dispersion vary little along the outer orbit or when the mean time interval between encounters R −1 is short compared to the period of the outer orbit), we can pick the time interval until the next 3-body interaction assuming the probability ( > ) for the time between the encounters  to exceed some  to be ( > ) = exp(−R), as in Hamers & Tremaine (2017).However, in many cases, we have to account for the variation of R along the trajectory, and we do this by adopting the following probability density distribution: reducing to the Hamers & Tremaine (2017) assumption for constant R.Here d is the probability that the encounter happens between  and  + d and  is the dimensionless auxiliary variable.Therefore, after each encounter, we determine the time  nxt until the next encounter in the following way: (i) Sample the value of  =  nxt from the exponential distribution, see equation (4a).
(ii) Integrate in time the outer orbit, the inner orbit (as described in §3.4), and the equation where (0) = 0 and R () is computed via the equation (2) at the current location of the binary CoM.
(iii) The moment of time when  =  nxt is  =  nxt and we assume that the next encounter has taken place.This procedure allows us to naturally account for the spatial inhomogeneities of  and  inside the cluster, which are sampled by the binary along its outer orbit.

Modeling of encounters
Once an encounter happens, we model it in the following way.
First, we randomly sample its initial parameters: (i) The perturber mass  3 from the distribution function  ( 3 ).In this paper we ignore the stellar mass spectrum and consider all perturbers to have the same mass.More realistically,  3 could be drawn from e.g the Salpeter distribution  ∝  −2.35 3 (Salpeter 1955) modified to account for the finite lifetime of stars (as was done in Hamers & Tremaine 2017).
(ii) The perturber initial velocity v rel (at infinity).To get it we first sample a velocity v in the cluster frame from the local (isotropic) velocity distribution  (), and then subtract the binary CoM velocity V from it.The local velocity distribution is considered to be Maxwellian, where  is the local velocity dispersion determined from the equation (3).The relative encounter velocity is then asymmetrically distributed, giving rise to the dynamical friction, which is discussed in more detail in §3.1.6.
(iii) The impact parameter  of the perturbing star (or any other passing object).For an adopted maximum encounter pericenter distance  max and  rel we calculate the maximum impact parameter  max and then randomly select  such that  2 is uniformly distributed between 0 and  2 max .(iv) The (uniformly distributed) orientation of the perturber orbital plane relative to its approach velocity at infinity (determined earlier).
(v) The initial mean anomaly of the binary (distributed uniformly).
In this work we adopt  max = 50 for the maximum pericenter distance of the perturbing star, see Fig. 2. Appendix A provides justification for this particular value of  max as well as for our choices of other BESC parameters discussed below - hyb ,  3body and  max .Our decision to include flybys as distant as 50 is motivated by the fact that the cumulative effect of distant encounters can have an important effect on the binary eccentricity, especially when  is high, and, as a consequence, on the overall binary merger rate (Samsing et al. 2019).However, that also means we have to process a very high number of flybys, scaling roughly as ∝  2 max and increasing as  grows.Luckily, the most distant encounters have a very small effect on the binary (individually), and in their course the orbit of the perturber does not deviate much from the hyperbolic trajectory.Because of that, for the flybys above a certain value of  =  hyb it is possible to reduce the computational cost by using some approximations instead of a full 3-body integration, as we describe next.We account for all flybys with pericentre distances of the third body  <  max = 50, but ignore the ones with  >  max (red trajectory).When  >  hyb = 10, we use the semi-analytical orbit-averaged approximation of Hamers & Samsing (2019b) to integrate the encounter (black trajectory).When  <  hyb (solid-to-dashed curve), we still use that approximation whenever the third body is more than  3body = 100 away from the binary centre of mass (black solid segment), but switch to the full 3-body integration whenever the incoming star approaches closer than  3body (blue dashed segment).See §3.1.2-3.1.4for more details.

Changes of orbital elements: distant encounters
Keeping track of the binary orbital elements is extremely important for enabling GW-assisted mergers of compact binaries, especially when the binary can be driven to high eccentricity.For that reason, we approach calculation of the binary inner orbital element changes with great care, in particular, making sure that we can accurately compute binary evolution at high eccentricities.Since at  → 1 the angular momentum of the binary becomes very small, even the weak kicks experienced by the binary during its encounters with rather distant perturbers can have a significant effect on e.g. the binary eccentricity (Hamers & Samsing 2019a,b).Given the large number of such distant encounters, the direct 3-body integration of each encounter is numerically prohibitive, so we used a different method based on the work of Hamers & Samsing (2019b).
For distant perturbers, which we define as those with  > 10, we use the orbit-averaged (OA) approximation from Hamers & Samsing (2019b) 5 to calculate the changes of the orbital elements.This semianalytic approximation is based on the evolution equations (Eqs. (2) of that work) for the eccentricity and (dimensionless) angular momentum vectors e and j, accurate to octupole order (i.e. the perturbation on the binary from the third body is expanded in the series of (/)  up to  = 3, where  is the distance between the binary and the perturber) and assuming averaging over the binary orbital phase (i.e. over the inner orbit of the binary).Assuming also that the orbit of the perturbing star relative to the binary CoM is fixed (i.e. a Keplerian hyperbolic trajectory, considering the binary potential to be that of a point mass  b ), these evolutionary equations for e and j are then integrated numerically over the entirety of the perturber's orbit.This procedure was implemented by Hamers & Samsing (2019b) as a Python script6 , which we use in our calculations.
In our case, we approximate this calculation by integrating over the perturber's hyperbolic orbit from the moment it enters the sphere of radius  max = 10 4  (not shown in Fig. 2 for simplicity) to the moment it leaves that sphere.This implementation of the OA approximation is tested against the direct 3-body integrations in Appendix A, where we show, in particular, that adopting  max = 10 4  provides an accuracy sufficient for our purposes.At the same time, in most cases the computational expense of numerical integration over this radial range is acceptable, even for large number of distant encounters (but see §4.2).
In principle, one may adopt other approaches to computing the changes of the binary orbital elements.For example, Hamers & Samsing (2019a) considered an approximate analytical solution of the OA equations in the 'first-order approximation' -when equations are integrated by assuming that all binary parameters stay constant in the course of an encounter; additionally these formulae assume a potential truncation at the quadrupole order.The first order approximation does not work well when the binary eccentricity is very high and the angular momentum |j| of the binary is low, so even a small perturbation can easily change |j| significantly during the encounter.Since we are often interested in the case of extremely eccentric binaries, that approximation is not suitable for us.
This issue can be circumvented by switching to the 'second-order' (Hamers & Samsing 2019a) or 'third-order' (Hamers & Samsing 2019b) approximations that take into account the variation of the binary orbital parameters in the course of an encounter.However, that approach is not practically applicable to the octupole approximation for the potential (which needs to be resorted to for moderately close encounters), as the number of terms in the ensuing analytical expression exceeds 10,000 (Hamers & Samsing 2019b, Table 1).Such higher-order calculation is feasible for the quadrupole-order evolution equations considered in (Hamers & Samsing 2019a, Eqs. 5), however, the quadrupole approximation is inaccurate in the case of an unequal-mass binary, which is also one of the possibilities we consider here.Because of that, we do not use any of the analytical approximations and just integrate the OA octupole equations instead.
Averaging over the (fast) binary orbital motion, which is the backbone of the OA approximation, is only justified when the third body moves sufficiently slowly, i.e. when the mean motion of the binary is much faster than the angular speed of the perturber at periapsis of its hyperbolic orbit (so called secular encounter).As shown by Hamers & Samsing (2019a), that condition always applies for hard binaries.However, as clear from Figure 1, the binaries where the effects of cluster tides are important, which are of great interest for us, are soft ( ≳  ℎ ) in most cases.For such binaries, one can show using the results of Hamers & Samsing (2019a, see their Eqs.1-3) that the OA approximation is valid only when As will be discussed later, in our calculations we typically consider the binary to be ionized (disrupted) and stop the simulation when  exceeds 10 3 AU.For a typical initial configuration considered in this paper ( 1 =  2 = 10M ⊙ ,  tot ∼ 10 6 M ⊙ ,  ∼ 1 pc) the hard binary separation  ℎ ∼ 10 AU.Therefore, / ℎ ≲ 100 and for  > 10 -the condition we adopted for considering an encounter as distant and for adopting the OA approximation -the secular constraint ( 7) is fulfilled.

Changes of orbital elements: close encounters
Some of the simplifying assumptions listed above (e.g. that of a secular encounter, purely Keplerian hyperbolic trajectory of the perturber, etc.) get gradually violated as the pericenter distance of an encounter becomes smaller than  hyb = 10.In that case, switching to a full 3-body integrations, as has been done by many in the past (Hamers & Tremaine 2017;Rodriguez et al. 2021), becomes inevitable.
To treat such situations in a computationally efficient way, we developed a 'hybrid' approach: we still use the OA approximation when the distance between the binary and the third body exceeds  3body = 100, switching to a precise 3-body integration only when the third body approaches the binary closer than  3body .This scheme is illustrated in Fig. 2. The 3-body integrations are carried out using ARCHAIN, an implementation of algorithmic regularization developed specifically to treat small- systems (Mikkola & Merritt 2008).During some particularly chaotic encounters, this switching between 3-body and orbit-averaged can happen multiple times as the perturber leaves and returns to  3body sphere, sometimes with an exchange interaction between the binary and the third body (when the third body becomes bound to one of the binary components, and the other one is ejected).
While switching from the OA approximation to a full 3-body integration at  3body , we must provide the 3-body integrator with the initial condition for the mean anomaly of the binary, which is averaged over in the OA approximation.This initial phase of the binary is chosen randomly every time.The numerical errors resulting from the 'hybrid' approximation and from choosing randomly the initial binary phase during the switchovers are discussed in Appendix A.
To summarize, we find this approach to provide a sufficient accuracy for our purposes.

Encounter outcomes
We stop the encounter calculation when one of the following outcomes takes place: (i) Two bodies are gravitationally bound to each other and form a binary, while the third one is at least  3,max away and is gravitationally unbound from the binary.
(ii) All three bodies become unbound from each other and are at a distance of at least  3body from each other -in that case, we consider the binary to be destroyed (ionized).
(iii) The number of times when we switch to precise 3-body integration routine in our hybrid procedure (see §3.1.3)exceeds 20.In those highly chaotic cases (which constitute less than 10 −5 of all encounters) we simply ignore the effect of that flyby.
The flowchart illustrating these outcomes as part of our numerical framework is shown in Fig. 3.

Changes of the outer orbit and dynamical friction
Encounters with the cluster stars change not only the inner but also the outer orbit of the binary.After every encounter, we calculate the change in the binary CoM velocity, which modifies its outer orbit around the cluster.During the distant encounters ( >  hyb = 10) the binary is assumed to act on the incoming star as a point mass  b (see §3.1.3),thus, the change of the binary CoM motion is given by the outcome of a two-body encounter between the masses  b and  3 .For closer encounters ( <  hyb ) the change of the outer orbit depends on the outcome of the 3-body integration which activates as a part of our hybrid method (see §3.1.4).
At the same time, we assume encounters to be local, i.e. after an encounter the binary re-emerges at the same point in space inside the cluster.New CoM velocity after each encounter is used as an initial condition for integration of the post-encounter outer orbit of the binary, see §3.4.
Given the randomness of the encounter parameters (see §3.1.2),the outcome of each individual encounter is stochastic and leads to dynamical (non-resonant) relaxation.However, because the binary sees a greater flux of cluster stars arriving from the direction of its orbital motion, a large number of encounters leads to a net effectdynamical friction (DF) -which causes the outer orbit of the binary to gradually decay towards the cluster centre (see §4.2).This effect is directly captured by our treatment of encounters in BESC, without any semi-analytical modeling (cf.Rodriguez et al. 2021).
At the same time, our constraint  <  max = 50 for an encounter to be considered may cause us to underestimate the magnitude of the DF.This is because DF is contributed not only by close but also by very distant encounters, with impact parameters comparable to the size of the cluster .By lowering the maximum impact parameter to ∼  max we are effectively reducing the Coulomb logarithm in the expression for the DF (Binney & Tremaine 2008) from ln[ 2 /( b )] ∼ 8 to ln[ max  2 /( b )] ∼ 5, where we adopted typical values  = 1 pc,  = 10 2 AU,  b = 20 ⊙ , and  = 20 km s −1 (typical for  cl = 10 6  ⊙ ).Thus, we expect our neglect of encounters with  >  max to result in us underestimating the magnitude of DF by ≲ 50%, which we consider acceptable given the intrinsic uncertainties of the model.

Secular evolution due to cluster tides
A distinct feature of BESC, present in only a handful of other studies (Hamilton & Rafikov 2019c; Bub & Petrovich 2020), is the inclusion of the effects of cluster tides on the evolution of the inner orbit of the binary.
Cluster tides arise because of spatial inhomogeneity of the cluster potential and cause binary orbital elements to vary in a secular fashion, without changing the binary semi-major axis.Their importance for the dynamics of relativistic (and other) binaries in stellar clusters has been pointed out by Hamilton & Rafikov (2019a,b,c), who developed an analytical formalism for treating their effects; even earlier their significance was recognized for the dynamics of the Oort Cloud comets (Heisler & Tremaine 1986;Brasser et al. 2006).Hamilton & Rafikov (2019b) explored the cluster tide-driven dynamics of a binary, showing that under certain circumstances a binary may be driven to very high eccentricities, similar to the Lidov-Kozai effect.As the intensity of the GW emission can rise dramatically during these high-eccentricity episodes, proper treatment of cluster tides is important for accurately following evolution of compact object binaries in a cluster.
Effects of cluster tides can be treated using two different levels of approximation.In a singly-averaged (SA) approximation, the evolution equations are averaged over the inner orbit of the binary only, and the local tidal field tensor (which serves as an input for the evolution) is sampled along the outer orbit of the binary in a cluster.In a doubly-averaged (DA) approach one also performs averaging of the tidal tensor over the outer orbit of the binary.As expected, the SA approach is more accurate and can reveal important effects which are hidden in the DA calculations, see Hamilton & Rafikov (2023) and §4.1 of this work.This is why here we use the full SA evolution equations to model cluster tides.
In BESC, we account for the tidal forces from the cluster potential using the numerical framework7 based on Eqs. ( 4)-( 5) from Bub & Petrovich (2020), formulated in the SA approximation.In this framework, the outer orbit of the binary is first evolved using the galpy package (Bovy 2015).Then we calculate the tidal tensor from the cluster potential at  points in this orbit, so that we could interpolate it later.We have established that 30 point per outer period is enough for almost all practical applications (apart from some rare cases where the eccentricity reaches extreme values, see §4.1).Finally, the differential equations for the evolution of the inner binary parameters are solved, with the tidal tensor interpolated in space at every timestep as mentioned above.

General relativistic effects
GR effects are the final key ingredients of our modeling framework.Emission of the GW is what causes the binary semi-major axis  to shrink and ultimately drives the binary to merge.GR precession changes only the binary periastron, without affecting other orbital elements, but it still plays a key role in the dynamics.
We include the GR effects by adding the following terms to the evolution equations for the semimajor axis  and eccentricity  of the inner orbit and the argument of periapsis  (Peters 1964): where  is the speed of light and  =  2 / 1 is the binary mass ratio.Contributions given by the equations ( 8b), (8c) are added to the corresponding evolution equations for  and  from §3.2.As discussed in Hamilton & Rafikov (2019b, 2021), when the semimajor axis  is low enough and/or the eccentricity approaches unity, pericentre precession due to the GR becomes so fast (faster than the rate of cluster tide-driven evolution) that it suppresses the eccentricity oscillations arising from the cluster tides.This has important physical implications, reducing the maximum value that  could reach and lowering the intensity of the GW emission, which considerably slows down binary evolution towards merger.The role of GR precession in counteracting the effects of cluster tides is characterized by the dimensionless parameter where  is the parameter (defined in Hamilton & Rafikov 2019a) characterizing the strength of the cluster tide in the DA approximation;  is of the order of the (squared) dynamical frequency of the cluster.One can think of  GR as the ratio of the GR precession rate (8c) for a circular binary to the rate of secular evolution due to cluster tides  −1 sec ∼  b ( b is the binary period), down to factors of order unity.As demonstrated by Hamilton & Rafikov (2021), for  GR ≳ 20 GR precession effectively suppresses the effects of cluster tides on binary evolution.
Fast GR precession at low  leads to certain numerical issues: in this regime the GR precession can be so fast that accounting for it in  evolution numerically (necessary since the cluster tide contribution depends on ) requires a very small timestep, grinding the calculation to a halt.At the same time, in this regime the cluster tides are no longer important and accounting for them is not really necessary.
For that reason we have adopted the following strategy: if at some point in the calculation (usually, after a close encounter that strongly reduces ) we find that then we (1) stop calculating the contribution of cluster tides to the evolution of binary orbital elements and (2) evolve the binary apsidal angle analytically, using () = (d/d) GR  +  0 ( 0 is the initial value of ).A threshold value of 20 for  GR was motivated by Fig. 9 of Hamilton & Rafikov (2021) and serves as a good compromise between the accuracy and numerical efficiency of the computation.This procedure allows us to avoid the aforementioned bottleneck due to the fast GR precession, dramatically speeding up the calculation.
If at some later time the binary semi-major axis is increased (as a result of a close encounter) and  GR drops below 20, we restart the numerical calculation of the orbital evolution accounting for the cluster tides.
To calculate  GR , we average  over the outer orbit (Hamilton & Rafikov 2019a) until the next encounter.In many cases, the corresponding segment of the outer orbit covers only a few orbital periods around the cluster centre and  may not have fully converged to its long-term average.That method, however, is good enough for our purposes, and quick enough to not slow down the computations.The  GR () plots shown later in Figs.5-11 do not use these values of  but rather the ones obtained using Eqs.( 39) and (D7) from Hamilton & Rafikov (2019a).We have checked that those agree on average with the values used in our numerical evolution.

Summary: evolution of the outer and inner orbits
The three physical processes described above are combined together in a unified framework for integrating the evolution of the outer and inner orbits of the binary, which conceptually operates as follows (see Figure 3 for an illustration of various steps).
(i) After each encounter with a cluster star we sample the random variable  nxt , which eventually determines the time when the next encounter occurs ( §3.1.1).
(ii) We then integrate the orbit of the binary CoM in the cluster potential using the galpy package (Bovy 2015).In this work, since we limit ourselves to the spherical potentials (a simplification that can be easily relaxed, see §5.3) the turning points (inner   and outer   ) and the plane of the outer orbit are preserved between encounters.
(iii) We simultaneously evolve the orbital elements of the inner orbit of the binary due to the (smooth) effects of the cluster tides ( §3.2) and general relativity ( §3.3).
(iv) Also, knowing the encounter rate at any point of the outer orbit (see §3.1.1)we integrate the equation ( 5) until  reaches  nxt .
(v) At this point we declare the next encounter to take place, stop the integration of the inner and outer orbits and process the encounter as described in §3.1; the (impulsive) changes of the inner orbital elements ( § §3.1.3,3.1.4)and of the outer orbit ( §3.1.6)are used to update the binary parameters.
(vi) We then determine the outcome of an encounter ( §3.1.5)and, if the binary survives, repeat the cycle.We consider the encounters to be instantaneous and local, thus we restart orbit integration from the time and point in space, at which the last encounter has taken place.This way, BESC follows the binary evolution until it either merges or is destroyed, or until a Hubble time has passed at which point we stop the calculation.

RESULTS
We now show the results of some tests of our framework, first focusing on the cluster tide-driven evolution and turning off the effects of encounters ( §4.1).This allows us to compare our framework against some known results and to demonstrate the importance of using the SA approximation for cluster tides.We then present some illustrative examples of binary evolution including both cluster tides and encounters with cluster stars ( §4.2).

Evolution in the absence of encounters: comparison with the DA approximation
In Fig. 4 we show some examples of binary evolution calculated using BESC with the stochastic stellar encounters completely turned (C) Illustration of the sensitivity of singly-averaged calculation to numerical parameters.Same initial conditions as in (A), with different lines corresponding to SA integration with the different number of interpolation points (per outer orbit) used to compute the potential derivatives when calculating the binary evolution: 30 (solid black), 100 (dashed blue), 300 (dotted red), 1000 (dot-dashed green).See §4.1 for more details on the chaotic nature of the singly-averaged evolution.
off, allowing us to focus on the effect of cluster tides and to compare the performance of the DA and SA approximations ( §3.2).Top row shows the evolution of the semi-major axis  while the bottom panels illustrate the behavior of 1 −  for binaries that reach very high eccentricities as a result of the specially chosen initial conditionsthe initial binary inclination (relative to the plane of the outer orbit) is very close to 90 • .We first focus on panels (A) and (B), in which we compare the evolutionary tracks of the binaries computed using the doubly-averaged (red, reproduced from Figs. 4 and 5 of Hamilton & Rafikov 2022, which should be consulted for parameters of these calculations) and singly-averaged (blue and black, computed using BESC) approximations with the same initial conditions, correspondingly.
The DA track in Fig. 4A shows an initially librating binary in a weak GR regime (see Hamilton & Rafikov (2022) for a detailed discussion of the different evolutionary phases), capable of reaching an extreme eccentricity  max ≈ 1 − 10 −5 in the course of cluster tide-driven secular oscillations.Its phase space trajectory eventually crosses the separatrix and becomes circulating (around  ≈ 4.8 Gyr), then the binary merges at  ≈ 7.4 Gyr.Fig. 4B shows another example from Hamilton & Rafikov (2022): a lower-mass binary with smaller initial semimajor axis , but with the same outer orbit in the same potential.It also starts with a librating phase space trajectory (although in the moderate GR regime) and goes through the same phase space transformations as the one in the example A before merging (see Hamilton & Rafikov (2022) for details).
As to the SA results, they are initially similar to the DA tracks, but start to diverge from them after a few secular cycles.In addition to that, evolution in SA approximation sensitively depends on the initial longitude of ascending node Ω 0 of the binary, which is the parameter distinguishing the black and blue tracks.The reason for this is that SA integration accounts for additional eccentricity fluctuations on the outer orbital period timescale, on top of the large-amplitude, smooth DA evolution of  on a much longer secular timescale.These eccentricity fluctuations are typically negligible, but they can become very significant for  → 1, as even a small change in  can strongly affect the angular momentum of the binary and amplify its GW emission.
In the example B those SA oscillations are small (peak eccentricity is almost the same for DA and SA tracks), and the DA-SA discrepancy increases gradually over many secular cycles, with both SA examples going through the same phase space transformations as DA, but with a delay (and as a result, merging much later than the DA track shows).
On the other hand, in the example A the eccentricity oscillations of the SA tracks are much larger, which is obvious from their deviations from the DA track at highest .As a result, even though the timing of the first 6-7 eccentricity maxima is approximately the same, an abrupt divergence between the three evolutionary trajectories in both  and  happens soon after.As a consequence of the much stronger GW emission during the high- episodes on the SA tracks, both SA integrations in panels (A) result in an order of magnitude shorter binary merger times than predicted by the DA calculation.
Interestingly, even in the absence of GW emission (which we can easily switch off in BESC, not shown here) but with the GR precession still present, the SA trajectories diverge from the DA ones much faster in the example A than in the example B. This behavior is a manifestation of the so-called relativistic phase space diffusion (RPSD), which was studied by Hamilton & Rafikov (2023).It takes place when, in the course of its SA evolution, the binary can reach eccentricity so high that the time it spends around the peak of  (over which GR precession drives a large-amplitude swing of the apsidal angle) becomes shorter than the outer orbital period of the binary.This causes sharp jumps of the (otherwise well conserved) secular integrals of motion of the binary (i.e. its perturbation Hamiltonian) and leads to a noticeable evolution of the period of secular oscillations (in the absence of both flybys and GW emission).
Hamilton & Rafikov (2023) have derived the following criterion for the SA eccentricity oscillations to be significant and, at the same time, for RPSD to occur (see their Eq.37): i.e. that the binary inner orbital plane must be very close to orthogonal to the plane of the outer orbit (a necessary condition for reaching a very high ).One can easily check that for the initial parameters specified in the captions of Fig. 4 the example A satisfies this condition, while the example B does not, which explains their different SA-DA behavior even in the absence of GW emission.Since Hamilton & Rafikov (2023) did not include GW emission when comparing DA and SA approximations for cluster tides, the examples shown in Fig. 4 provide the first illustration of the coupling between the RPSD and the GW emission.Their synergy makes the binary merge much sooner, as even a small RPSD-related increase of  around the eccentricity maximum can significantly increase the rate of the semimajor axis decay.
Although only two SA trajectories (for two values of Ω 0 ) are shown in Fig. 4A,B for every set of initial conditions, SA trajectories calculated for other initial Ω 0 show similarly erratic divergence from the DA track and from each other.This illustrates the essentially chaotic behaviour of the SA oscillations, in a sense that even a tiny change in the input parameters of the calculation can lead to a significant change in the binary evolution.
This statement is true not only for physical inputs (e.g.varying Ω 0 in the examples A and B), but also for numerical parameters of the calculation, as we illustrate in Fig. 4C.There we show (for the same initial conditions as in the example A) the result of varying the number of points  on the outer orbit that are used to approximate the derivatives of the gravitational potential when computing the tidal tensor for the SA calculation (see §3.2).One can see that when  periodically reaches extreme values, results of the calculation become very sensitive to the value of , and the SA tracks show no obvious convergence even as we increase  to  = 1000.This means that the long-term binary evolution in the SA approximation is realistically hardly predictable, and that one can draw meaningful conclusions only upon examining the evolution of large samples of binaries with similar initial conditions.
As we will see next, inclusion of encounters adds so much additional stochasticity to the binary evolution that this completely obviates the need to worry about the sensitivity to various input parameters.For practical applications -calculating the evolution including encounters -evaluating cluster potential at  = 30 points turns out to be fully adequate.

Binary evolution with all physical effects included: examples of different outcomes
We now illustrate the performance of BESC with all the physics (described in §3) fully incorporated, including stellar encounters.Figs.5-10 show several examples of individual evolutionary tracks of BH-BH binaries illustrating a variety of possible outcomes: a merger (Fig. 5), an impulsive binary disruption (ionization, Fig. 6), an exchange interaction (Fig. 7), a diffusive disruption of a binary ( exceeding 10 3 AU, Fig. 8), an ejection from the cluster (Fig. 9), a binary that survives for a Hubble time (Fig. 10).Each of these cases is discussed in §4.2.1-4.2.6.Unless specified otherwise, we have chosen the following initial parameters for these evolutionary tracks: • Spherical Hernquist potential of the cluster Φ() = −  cl /(+ ) with  cl = 10 6 M ⊙ and  = 2 pc.Velocity dispersion is assumed to be isotropic and given by the Eq. ( 10) from Hernquist (1990).
The initial parameters, which are different for the individual tracks that we present are as follows: • Initial inclination  0 is drawn from a uniform distribution of cos  0 .
• Apsidal longitude  0 is uniformly distributed (although the longitude of ascending node is set to Ω 0 = 0).This ensures that some differences between the evolutionary tracks exist from the start, however, the key reason for their subsequent divergence is the stochasticity due to stellar encounters.
Our choice of a (cusped) Hernquist modes may be more suitable for nuclear star clusters than for globular clusters with their cored density profiles.However, it allows us to illustrate some outcomes that are unlikely in the cored models, e.g.see §4.2.5.
Our choice of the initial semimajor axis  is motivated by the following considerations.The initial value of the GR parameter  GR (0) = 0.26 satisfies the constraint (10), i.e. cluster tides are not suppressed by the GR precession initially.However, given a very steep dependence of  GR on ,  GR ∼  −4 , this would no longer be the case if  were much smaller than 100 AU; the threshold ( 10) is reached at (from Eq. 9) Thus, in order to explore the effects of the cluster tides we should start with  ≳ 50 AU.
On the other hand, the higher the semimajor axis, the less probable it is for the binary to eventually merge, since softer binaries harden due to encounters less efficiently.In addition to that, when dealing with encounters in BESC the encounter rate grows approximately as ∼  2 (Eq.2), and for very large  the number of encounters that needs to be followed becomes so high that the computational cost of the calculation becomes prohibitive.To avoid that, we stop the simulation whenever the binary semimajor axis reaches 10 3 AU.We have not observed any cases of a binary going above that limit but then shrinking its orbit and merging.
As a result, starting evolution at  = 10 2 AU, roughly at the hard/soft boundary, appears to be a reasonable compromise between the two aforementioned limiting factors.shows  GR ( ) defined by the equation ( 9) and illustrating the importance of GR precession (red line is the constraint ( 10)).See §4.2.1 for details.

Example 1: binary merger
of the outer orbit.Note that  io is different from the inclination  relative to the fixed cluster frame (set in the beginning of integration, when  =  io ), since the plane of the outer orbit changes as a result of stellar encounters; see below for the significance of  io .In panel (d) we show the outer   and inner   turning points of the outer orbit of the binary.In panel (e) we show the evolution of an important dimensionless parameter Θ io = (1 −  2 ) cos 2  io , an essential integral of motion in both the DA (Hamilton & Rafikov 2019b) and SA (Hamilton & Rafikov 2023) approximations.As shown in Hamilton & Rafikov (2019b), Θ io ≪ 1 is a necessary condition for the cluster tides to be able to increase  close to 1 for a binary that was not very eccentric initially, which means that initially  io must be very close to 90 • .Finally, panel (f) illustrates the behavior of  GR (see Eq. 9), which characterizes the importance of the GR precession in suppressing the effect of cluster tides.All these variables are shown as functions of time.
One can see that the outer orbit of the binary gradually decays over time (panel (d)) as a result of dynamical friction, which is a consequence of multiple stellar encounters, see §3.1.6.Encounters also result in rapid fluctuations of   and   but this diffusive effect is overwhelmed by the overall orbital decay.The outer orbit never becomes too radial as   and   remain similar to each other.Eventually, around  = 0.25 Gyr the binary sinks to the centre of the cluster, where the encounter rate rapidly grows because of increased stellar density  and velocity dispersion .This sinking time is in general agreement with the Chandrasekhar's formula with  and  evaluated at the characteristic radius of the cluster .
During this time, the inner orbit evolves as well (panels (a)-( c)).In this example it stays moderately eccentric through most of the evolution until eventually one of the stellar encounters, that are very frequent in the cluster center, sends  to a very high value, 1 −  ≲ 10 −3 .With the periastron distance below 10 −2 AU (panel (a)), the binary rapidly shrinks its orbit due to the GW emission and merges (the merger timescale for an  b = 20 ⊙ ,  = 1 binary with  = 100 AU and   = 0.01 AU is around 4 Myr, see Peters 1964).
Panel (f) shows that until the very last moment  GR stays below the critical value of 20.This implies that cluster tides are not suppressed by the GR precession and are being accounted for by BESC.Given that the initial inclination of the binary is not close to 90 • , cluster tides acting alone are unable to increase  to values close to unity.But they are still capable of changing  by ∼ 0.1 in the course of the evolution.Nevertheless, the final increase of , which ultimately results in a merger, is due to stellar encounters, which are abundant in the cluster centre -this manifests in the orbital parameters starting to change rapidly at the end of the binary lifetime.We can thus conclude that in this example it is the encounters and not the cluster tides that drive the binary to merger.This can also be inferred from the behavior of Θ io (panel (e)) that becomes very small (while  is not close to unity) only for a short period of time in the end of the evolution.This does not give cluster tides a chance to increase  of this binary to values close to unity (Hamilton & Rafikov 2019b).

Example 2: binary ionized impulsively
Fig. 6 shows a binary that ends up being destroyed (ionized) by experiencing a close flyby that causes all three bodies to become unbound from each other.It starts from the same initial conditions as in the previous example (except for the initial orientation of the inner orbit), and the only reason for the different outcome is the stochastic nature of the flybys.As in the previous example, the outer orbit of the binary decays with time due to the DF, such that the rate of stellar encounters steadily increases.Simultaneously with that, the inner semi-major axis steadily grows to ∼ 200 AU, while the periastron distance decays to ∼ 10 AU.The corresponding eccentricity growth is accomplished predominantly through stellar flybys; even though cluster tides are operating effectively ( GR ≪ 20), they change  only at the level of ∼ 0.1 since Θ io does not get low enough, see panel (e).
At  ≈ 0.08 Gyr, still far from the centre of the cluster, the binary experiences a strong penetrative encounter with a cluster star ( = 0.014,  rel,∞ = 32 km s −1 = 3.3 √︁  b /), which completely unbinds it.This happens when  is still substantially below 10 3 AU, which is our threshold for following its evolution.Such impulsive ionization is a rather frequent evolutionary outcome given the large initial  that we adopt in these examples (our reasoning behind that is explained above in §4.2), making the binary 'soft' from the start.

Example 3: binary experiencing an exchange interaction and then merging
A binary can also experience an exchange interaction whereby, in the course of a close encounter with a cluster star, one of the binary components gets ejected and the perturber takes its place, see Fig. 7.
In this example the binary again sinks to the cluster centre due to the DF within 0.2 Gyr.In the process, its eccentricity shows a substantial evolution (panel (b)) as it gets increased to 1 −  ∼ 10 −4 around  = 0.08 Gyr by encounters with cluster stars (we verified that).Once the binary has reached a high eccentricity (i.e.low angular momentum), it becomes easy for any given encounter to increase the eccentricity even further, although this process is stochstic.After reaching this peak,  then decreases to ∼ 0.1, again, due to encounters.Cluster tides play some, but not decisive, role in this eccentricity evolution.
Upon reaching the cluster center, around 0.18 Gyr, a close encounter with a star ( = 0.088,  rel,∞ = 4.8 km/s = 0.11 √︁  ( 1 +  2 )/  ) takes place, unbinding one of the black holes and binding another one to the incoming star, thereby forming a new 10 ⊙ + 1 ⊙ binary.This new binary has a rather small semi-major axis < 0.1 AU (so that the GR precession completely suppresses secular effect of cluster tides), which is necessary to unbind a heavy BH companion, and it eventually merges due to the GW emission at  ≈ 0.7 Gyr.
Even though this new binary sits right at the cluster centre where the stellar density is high, the rate of close stellar encounters capable of changing its  drops dramatically.This is because R ∝  2 and the close stellar flybys are very rare, see panel (a).At the same time, encounters are still capable of temporarily dislodging the new binary from the cluster centre by ∼ 0.2 pc from time to time, see panel (f).

Example 4: binary ionized diffusively
Fig. 8 shows an example of a binary that is disrupted (ionized) in a diffusive fashion.As the binary sinks towards the cluster centre due to dynamical friction, it gets softened by encounters that become more frequent as  decreases.Eventually,  crosses 10 3 AU at which point we stop the calculation since it would get fully unbound soon after anyway.Here binary ionization occurs as a result of a large number of small changes of  that accumulate in a diffusive fashion, unlike the case shown in Fig. 6, in which a single close encounter disintegrated the binary when its  was still not too far from the hard-soft boundary.
We also note that in this example, due to the large , cluster tides enable somewhat stronger  evolution than in the previous examples, with the corresponding Δ ∼ 0.2.

Example 5: binary ejected
One of the possible consequences of the close encounters is that the binary is sometimes ejected from the cluster.In our simulations that almost always happens in Hernquist potentials with low mass  cl = 10 5 M ⊙ .One of those cases is shown in Fig. 9.Here the binary sinks to the cluster centre rather quickly, within 0.1 Gyr, where its  get rapidly reduced by encounters to below 1 AU (effectively turning off cluster tides, see panel (e)).A strong encounter around 0.16 Gyr then kicks the binary into a very radial orbit inside the cluster with   ≈ 3 pc but   ≈ 0 pc (and also increases its ).The outer radius   then decays by the DF and by 0.4 Gyr the binary orbit is again fully confined to the central regions of the cluster.This cycle repeats around 0.42 Gyr and 0.7 Gyr, until a particularly strong encounter completely ejects the binary from the cluster around  = 1 Gyr.
The exact value of the semimajor axis at which the binary becomes prone to ejections can be estimated as follows.For  1 ≈  2 ,  3 ≲  b the typical velocity of the cluster star  3 after a close flyby ( ∼ ) is  ej,3 ≈ √︁  b / (Merritt 2013, Eq. 8.21) and therefore the typical binary recoil velocity is The escape velocity from the centre of a Hernquist potential is

√︁
2  cl / (Hernquist 1990, Eq. 15).Thus,  ej,bin >  esc and the binary can get ejected from the cluster centre when its  drops below This estimate is consistent with Fig. 9 since  drops to around 0.1 AU right before the ejection, while  is somewhat larger than  ej during the earlier sharp increases of   that did not result in an ejection.Equation ( 14) also agrees with the results on binary ejections obtained using CMC (Weatherford et al. 2023, see their Figs. 6).
The ejection probability is higher in a Hernquist potential than in a cored potential (e.g.Plummer), which is not surprising: in a Hernquist potential  esc from the cluster centre is finite, while the central stellar density diverges (unlike the Plummer case), so the rate of encounters able to eject the binary increases dramatically (even for small ) as the binary sinks into the centre.However, if  is too small, an encounter with  ∼  capable of ejecting the binary may not happen at all before the binary merges due to the GW emission.This is what we observe in the example shown in Fig. 7, in which after the exchange  is significantly below  ej (in this example  esc is also higher due to higher  cl ).Thus, ejections are more likely for  just below  ej , but not much below.

Example 6: binary survives for a Hubble time (with kicks disabled)
To demonstrate the impact of DF on the fate of the binary, we performed a set of simulations where we do not take into account the velocity kicks received by the binary CoM as a result of stellar flybys.This artificially eliminates DF, but all the inner orbital element changes due to encounters are still fully accounted for.In such runs the outer orbit stays unchanged at  = 2 pc resulting in much lower encounter rates than in the runs including DF (see previous examples).
In Fig. 10 we show one such run in which a binary gets gradually hardened over time, but its  evolves slowly due to the infrequent flybys, so that the binary survives for 10 Gyr.This almost never happens when the binary CoM recoil is properly accounted for giving rise to DF, emphasizing the importance of encounter-driven evolution of the outer orbit for determinig the overall evolutionary outcome.

Example 7: cluster tides-dominated evolution (with kicks disabled)
In the previous examples, the encounters affected the binary evolution stronger than the cluster tides -i.e. the eccentricity oscillations in the absence of encounters described in §4.1 are severely washed out by the stochastic changes of  due to encounters.In part, this is because in none of these examples the binary inner and outer orbital planes were kept orthogonal, which is necessary for large amplitude oscillations of  (Hamilton & Rafikov 2019b).Also, encounters very easily conceal the signs of smooth cluster tide-driven evolution in time series of .We will investigate the relative contribution of the two effects in a more quantitative manner in a subsequent work (Rasskazov & Rafikov 2023, in preparation).
For now, we just show in Fig. 11 a particular example with the specially designed initial conditions, which allow cluster tides to reveal themselves rather clearly.To achieve that, we put the binary on a rather distant orbit in the cluster (  =   = 4) and artificially maintain this outer orbit by turning off DF (like in §4.2.6).We have also assumed a higher initial semimajor axis  = 300 AU and the initial inclination of  0 = 89.9• to enable significant eccentricity oscillations.This high  prevents GR precession from suppressing the cluster tides, see panel (f).
One can see that 1 −  gradually decreases to ≲ 10 −2 , a behavior that reverses around 0.75 Gyr in an almost symmetric fashion.The () dependence is rather smooth and is only weakly perturbed by encounters, which are rather weak and infrequent so far from the cluster centre (but are still capable of gradually softening the binary).This behavior is strongly indicative of the secular cluster tide-dominated evolution, similar to Lidov-Kozai cycles, which is triggered in this case by the low value of Θ io ≲ 10 −2 throughout the evolution.This conjecture is also confirmed by the estimate of the secular timescale expected8 for such a binary,  sec ≈ 1.6 Gyr, which is in good agreement with the approximate oscillation half-period of 0.7 Gyr in Fig. 11a,b.This binary eventually crosses  = 10 3 AU threshold at which point we stop our calculation.

DISCUSSION
This study is primarily devoted to the description of a newly developed code BESC and to illustrating its performance via some illustrative examples.For that reason, here we are not providing quantitative results on binary dynamics except for the preliminary statistics of the outcomes in §5.1 below.
BESC enables a statistical analysis of the evolution of large samples of binaries in stellar clusters by following individual trajectories of each binary in great detail.The initial parameters of these samples should be provided to BESC as inputs since it does not follow the production of binaries in clusters.As we showed already in §4.1, the evolution of an individual binary is typically very chaotic, especially in presence of encounters, see §4.2.As a consequence, meaningful conclusions about binary evolution in stellar clusters are possible only in a statistical sense, based on large samples of simulated binaries.This is the approach we will adopt in our followup work.
Although the examples that we provided in this study are limited to BH-BH binaries, BESC can certainly be applied to studying other exotic binaries encountered in stellar clusters -blue stragglers, X-ray binaries, etc.One just needs to adopt the relevant physical ingredients, as discussed in §5.3 below.

Preliminary statistics of the outcomes
A detailed analysis of the statistics of binary evolution outcomes based on BESC calculations with realistic inputs (stellar mass spectrum, distribution of initial semimajor axes and eccentricities, characteristics of the outer orbits) is left for future work.Nevertheless, we can provide a preliminary (rather low number) statistics of the outcomes at least for the initial conditions used in this work, see §4.2.To this goal, we have performed 8 sets of simulations with these initial conditions, about 100 realizations for every combination of the Hernquist and Plummer (not shown earlier) potentials, cluster mass  cl = 10 5 M ⊙ , 10 6 M ⊙ , with and without the effects of encounters on the outer orbit (i.e. with or without the DF, see § §4.2.6, 4.2.7).
Based on this sample, we found that the fraction of binaries that merge within a Hubble time is higher in Hernquist models than in the Plummer ones: 76% vs 34% for  cl = 10 5 M ⊙ , and 45% vs 19% for  cl = 10 5 M ⊙ .Turning off flyby kicks on the outer orbit (i.e.preventing the outer orbit from evolving and reaching high stellar density in the cluster center) reduces the merged fraction to around 6% regardless of other parameters.This shows that DF is an extremely important factor in determining the fate of the binary: it allows the systems starting far from the cluster centre to eventually merge as the DF brings them into the much denser central regions where they can experience 3-body hardening more efficiently.None of the binaries end up being disrupted either impulsively ( §4.2.1) or diffusively (via exceeding  in = 10 3 AU, §4.2.4) in  cl = 10 5 M ⊙ clusters.However, in  cl = 10 6 M ⊙ clusters more than half of the binaries get ionized either way.This is because for our choice of the initial  = 10 2 AU binaries are hard in the former and soft in the latter cases, see Fig. 1b,c.Finally, as mentioned in §4.2.5, we observe binaries being ejected from the cluster in significant numbers (18%) only in a Hernquist model with  cl = 10 5 M ⊙ .

Comparison with existing results
To the best of our knowledge, BESC is the first code that evolves the binary inner orbital elements taking into account not only the stochastic effects of close stellar encounters but also the secular effect of cluster tides, while also self-consistently evolving the outer orbit of the binary CoM in the cluster (directly capturing the effect of dynamical friction).Earlier studies focused on just one physical process, e.g.Hamers & Tremaine (2017), Samsing et al. (2019) and Rodriguez et al. (2021) include only encounters, while Hamilton & Rafikov (2019c) and Bub & Petrovich (2020) considered only cluster tides.BESC incorporates these effects in a semi-analytical fashion, making it more computationally efficient than the direct N-body calculations (e.g.Li et al. 2023), which naturally account for both effects (as long as they are able to reliably treat close encounters).
Until now, one of the most comprehensive treatments of the binary evolution in stellar clusters has been provided by the Cluster Monte-Carlo (CMC) code developed by Rodriguez et al. (2021).We summarize the differences between the CMC and BESC as follows: • Unlike BESC, CMC does not account for cluster tides in the evolution of the inner orbit.
• CMC ignores the effect of distant flybys ( > 2) on the inner binary parameters, while BESC explicitly accounts for them in a semi-analytical fashion.
• CMC directly includes the effect of very close flybys ( < 2) on both the inner and outer binary orbits, but for the more distant encounters ( > 2) it accounts for their effect only on the outer orbit, in a prescribed fashion.
• Therefore, CMC should account for the dynamical friction, although Rodriguez et al. (2021) do not explicitly mention this; BESC accounts for the dynamical friction directly, by explicitly calculating the outcome (change of the CoM motion) of each encounter.
• CMC follows only the energy and angular momentum of the outer orbit but not its actual shape; the phase of the outer orbit is randomly sampled at each time step.BESC fully traces the outer orbit as a function of time, which is important in aspherical clusters.
• Since CMC describes every member of the cluster, it accounts for the self-consistent evolution of the cluster as a whole, whereas BESC assumes a prescribed model for the cluster density (which may however be made time-dependent).
These differences of BESC with CMC, as well as with other existing numerical approaches to binary evolution in clusters, are also summarized in Table 1.

Future extensions
The implementation of BESC described in this paper makes a number of approximations regarding the setup and relevant physics, many of which can be easily relaxed in the future.
For example, in this paper we have not considered non-spherical potentials, but the cluster tide module implemented by Bub & Petrovich (2020) which is used in BESC works also in fully triaxial potentials.We assumed that all perturbers have mass  3 = 1M ⊙ but the code can naturally account for a mass spectrum of perturbers (the effect of varying  3 will be explored in Rasskazov & Rafikov, in prep.).We assume the cluster potential/density distribution to be constant in time, not accounting for cluster relaxation, core collapse, mass segregation, etc.In the future, these effects can be modeled by evolving the cluster mass distribution and potential in a prescribed manner.
Since our focus in this paper is on BH-BH binaries, we accounted for the GR precession and the dissipative effect of the GW emission.However, BESC can be easily extended to treat the phenomena including less relativistic binaries in clusters, for example normal stellar binaries, stars with planets, blue stragglers, X-ray binaries, etc.In such cases one would need to include additional sources of precession -rotationally-and tidally-induced quadrupoles -and consider the possibility of tidal dissipation inside the binary components instead of the GW emission.

SUMMARY
We developed a numerical framework -BESC -that allows one, for the first time, to simulate evolution of a binary in a stellar cluster while simultaneously including a variety of physical effects: stellar flybys, tidal forces from the smooth cluster potential, and GR effects.Effect of cluster tides are accounted for at the singly-averaged level, without averaging over the outer orbit of the binary (which we compute explicitly).Moreover, when accounting for stellar flybys we incorporate the effects (on the inner and outer orbit of the binary) of both the close stellar encounters (via the direct 3-body integration of an encounter) as well as the distant flybys (using a semi-analytic secular approximation).
We then carried out a number of binary evolution calculations with various parameters, and our preliminary results can be summarized as follows.
• In the absence of encounters, we performed a comparison of the singly-averaged (SA) framework for the cluster tide-driven evolution with the existing results Hamilton & Rafikov (2022) based on the doubly-averaged (DA) approximation.We found that the binary eccentricity oscillations on the outer orbital timescale can cause binary evolution to deviate from the DA approximation in an unpredictable manner at very high eccentricities  → 1.As a result, in most cases the binary reaches higher eccentricities and merges sooner than in the DA calculations.
• Encounters with cluster stars play a very important role in a variety of ways.In particular, they tend to disrupt the smooth secular evolution driven by the cluster tides adding a strong chaotic component to the variation of binary orbital elements.
• Another key outcome of stellar encounters is the dynamical friction that causes the outer orbit of the binary to decay towards the cluster center.This effect is directly accounted for in BESC thanks to its inclusion of the distant encounters.
• Cluster tides are most effective when the binary semi-major axis is in the range  ∼ 10 2 − 10 3 AU, however, these binaries can be easily ionized by passing stars.Also, for lower  the GR precession presents an important obstacle to reaching  → 1 in the course of cluster tide-driven evolution.
• Cluster tides can dominate binary eccentricity evolution (over the encounters) only when the binary is very far from the cluster centre and evolution of its outer orbit due to the dynamical friction is suppressed.Nevertheless, as we demonstrate more quantitatively in the future (Rasskazov & Rafikov 2023, in prep.), the effect of cluster tides on the binary eccentricity evolution cannot be ignored.
• In agreement with other studies (e.g.Weatherford et al. 2023), we sometimes observe exchange interactions or the binaries escaping the cluster.
• The BH binary merger rate seems to be mainly determined by the mean density at the outer orbit of the binary, which is higher in the Hernquist (cusped) potential compared to the Plummer (cored) potential.Stellar encounters play the key role in promoting mergers in part via the encounter-driven dynamical friction that sinks the binary towards the cluster center, where the stellar density is high.
We also note that BESC can be used to follow the evolution of other binary systems within a cluster, such as hot Jupiters (improving on the numerical method of Hamers & Tremaine 2017), blue stragglers, X-ray binaries, etc.

APPENDIX A: CHOICES OF BESC PARAMETERS AND THE ACCURACY OF HYBRID INTEGRATION
We start by providing a more quantitative justification for our choice of the upper limit  max = 50 for an encounter to be considered at all. Figure A1 shows, for several typical binary evolutionary tracks (in different panels), the total change in eccentricity (Δ) total over the entire binary lifetime due to the flybys with pericentre distances  <  max only (i.e., with more distant flybys excluded), as a function of  max /.We see that (Δ) total typically saturates and becomes essentially independent of  max at  max / ∼ 20.In other words, accounting for the encounters more distant than  max = 50 would have almost no effect on binary eccentricity evolution, with good margin.
Next we describe some tests of accuracy of the orbitaveraged/hybrid approach (described in §3.1.3)for calculating the outcome of a stellar flyby and motivate our choices of the parameters  hyb ,  3body and  max employed by BESC.We do this by comparing the result of a hybrid calculation with a direct 3-body integration of an encounter and analyzing the resultant deviations.For every value of the encounter pericenter distance  we sample 100 initial conditions, with all the orbital angles randomized and the rest of parameters set as follows: • The masses of all three bodies are  1 =  2 =  3 = 5M ⊙ • Binary semimajor axis is  = 1 AU • Binary eccentricity  = 0.1 or  = 0.999 • Initial velocity of the third body (at infinity)  = 3 km/s Cumulative change of the binary eccentricity (Δ) total that it has experienced in its lifetime due to stellar encounters with periastron distances  <  max , plotted as a function of  max / for several representative binary evolutionary tracks.In these tracks the effect of encounters on the outer orbit has been turned off (i.e.no dynamical friction), the outer orbit lies between   = 1.9 pc and   = 2.1 pc.Calculation assumes Hernquist potential with total mass  cl = 10 6 M ⊙ and scale radius  = 2 pc.The initial binary inclination and eccentricity are  = 89.8• and  = 0.5, respectively.One can see that (Δ) total saturates for  max / ≳ 20, motivating our choice of  max = 50 in BESC, see §3.1.3.
For every such initial condition we simulate an encounter and measure the change of the binary orbital elements (Δ, Δ etc.) using all 3 approaches: 3-body, orbit-averaged and hybrid.In Fig. A2 we show the fraction of the simulated encounters where the relative error in Δ, Δ etc. for the orbit-averaged/hybrid approach (assuming 3-body approach is precise) exceeds 20% (e.g. when Δ hybrid /Δ 3body − 1 > 0.2).
As expected, the orbit-averaged approximation breaks down at / ≲ 5 where the perturbations to both the binary and the perturber's orbit become too strong.For that reason we set  hyb = 10, so that all encounters with  <  hyb are processed using the hybrid approach, including the direct 3-body integration.
In some cases, there is also a small fraction of significant errors at large /.These deviations at large / are due to the Δ(/) function going through zero at a certain /, so that even very small differences in Δ (between the 3-body integration and other methods) lead to large relative errors.This phenomenon is illustrated in Fig. A3, where we plot |Δ| for two sets of initial conditions (described in the caption), both for relatively low  = 0.1; in both cases Δ passes through zero at / ≳ 15.
Even with the hybrid approach, the fraction of discrepancies with the 3-body simulations starts to rise when / ≲ 6.However, that happens not because of hybrid approach being inaccurate compared to 3-body, but rather due to the chaotic nature of close interactions, i.e. the strong dependence of the flyby outcome on the initial binary mean anomaly (which is picked randomly in both hybrid and 3-body methods).For the same reason, 3-body integrations with different initial mean anomalies would also disagree with each other.
At higher / the number of discrepancies goes down, but only when  3body -the distance where we switch from the orbit-averaged to 3-body integration -is large enough, see Fig. A2.For high enough  3body ≳ 100, the fraction of discrepancies falls to zero at  ≳ 8 while it sometimes stays at a nonzero level for lower  3body .Fig. A4 illustrates why that happens.It shows the evolution of (a) perturber separation  and (b) binary eccentricity during an encounter (computed with the 3-body integrator).Zooming in on the moment of the closest approach in panel (c), one can see the eccentricity oscillations on the orbital period of the binary with the amplitude comparable with (or even larger than) the resultant total eccentricity change.In the 3-body approach those oscillations are naturally accounted for, while in the orbit-averaged approach they are averaged out.In the hybrid approach, however, neither of those two things happen, and the initial binary phase (randomly chosen when the perturber separation crosses  3,max ) appears to affect the result more than it actually does.While the Δ values obtained using the hybrid approach are still precise on average, their dispersion is artificially increased.This effects weakens as we choose larger  3body .Based on these considerations, in our calculations we set  3body = 100 and use the orbit-averaged approximation whenever  >  hyb = 10, switching to the hybrid one for  <  hyb , as mentioned in §3.1.3.That way, the fraction of simulations with > 20% errors in the eccentricity calculation stays below 2%.
We start/stop the encounter processing when the distance from the third body to the binary is  =  max = 10 4 .While this appears as a very distant threshold, in some cases we really need to start the hybrid or purely orbit-averaged integration that early to fully account for the effect of a flyby.This is illustrated in Fig. A5 where we show the dependence of Δ on  max for two sets of initial conditions (with pericentre distances  ∼ 10 in both cases).One can see that in these cases Δ converges on its true value only for  max ≳ 10 3 .To be fair, in both these cases Δ is quite small (these parameter sets are close to the Δ = 0 points discussed above and illustrated in Fig. A3) and for most other choices of the encounter initial conditions such high  max is not required for an accurate Δ calculation.Nevertheless, since increasing  max above 10 4  does not increase the orbit-averaged/hybrid calculation time significantly, we set  max = 10 4  in BESC.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure1.Characteristic values of the inner semimajor axis  delineating different physical regimes for a binary evolving in a stellar cluster.Different curves correspond to the following: (black solid) hard/soft binary separation  ℎ , (black dashed) the semimajor axis  3body at which 3-body hardening timescale  3body becomes shorter than Hubble time, (red solid)  =  GW at which the GW decay timescale becomes shorter than Hubble time, (green solid)  =   at which GR precession suppresses cluster tides ( GR = 20), (green dashed)  =  sec at which the secular timescale is equal to 0.1 Hubble time.The cluster has characteristic radius  = 1 pc (top) or 3 pc (bottom) and total mass  cl = 10 4 M ⊙ (left), 10 5 M ⊙ (middle) or 10 6 M ⊙ (right).See §2.1 for additional details.

Figure 2 .
Figure 2. Schematic illustration of the different regimes of encounter integration used in BESC.We account for all flybys with pericentre distances of the third body  <  max = 50, but ignore the ones with  >  max (red trajectory).When  >  hyb = 10, we use the semi-analytical orbit-averaged approximation ofHamers & Samsing (2019b) to integrate the encounter (black trajectory).When  <  hyb (solid-to-dashed curve), we still use that approximation whenever the third body is more than  3body = 100 away from the binary centre of mass (black solid segment), but switch to the full 3-body integration whenever the incoming star approaches closer than  3body (blue dashed segment).See §3.1.2-3.1.4for more details.

Figure 3 .
Figure 3.The flowchart illustrating the operational logic and major modules of BESC, as well as the possible evolutionary outcomes.See §3.1.5,3.4 for details.

Fig. 5 Figure 5 .
Fig.5illustrates the orbital evolution of a binary that ends up merging, and we use it to describe the presentation of binary characteristics in Figs.5-10, which all have the same structure.In panels (a)-(c) we show the 'raw' orbital parameters of the inner orbit -(a) the semimajor axis  and the periastron distance  p , which is indicative of the GW emission strength, (b) the deviation of the binary eccentricity  from unity, and (c) the cosine of the inclination  io of the inner orbit relative to the (instantaneous) plane

Figure 6 .
Figure6.Same as Fig.5except the binary is disrupted (ionized) by a 3-body encounter (see the legend at the top for the parameters of this evolutionary track).See §4.2.2 for a discussion.

Figure 7 .
Figure 7. Same as Fig. 5 except the binary experiences an exchange interaction at the moment marked by the blue dashed vertical line.See §4.2.3 for a discussion.

Figure 8 .
Figure8.Same as Fig.5except the binary inner semimajor axis  eventually exceeds 10 3 AU, at which point we abandon the calculation and consider the binary to be effectively disrupted.See §4.2.4 for a discussion.

Figure 9 .
Figure9.Same as Fig.5, except the total cluster mass is  cl = 10 5 M ⊙ and the binary ends up being ejected from the cluster.See §4.2.5 for a discussion.

Figure 10 .
Figure 10.Same as Fig. 5 except the binary baricenter velocity kicks due to encounters are disabled (i.e.no dynamical friction, the (circular) outer orbit of the binary is unchanged, see panel (d)) and the binary survives for 10 Gyr.See §4.2.6 for a discussion.

Figure 11 .
Figure11.Same as Fig.5, but now showing the evolution of an initially wide ( = 300 AU) and distant (  =   = 4 pc) binary in a Hernquist potential with total mass  cl = 10 5 M ⊙ and radius  = 1 pc with the velocity kicks due to encounters disabled.This example illustrates the possibility of cluster tides dominating the eccentricity evolution (eventually the binary is ionized).See §4.2.7 for a discussion.

Figure A2 .Figure A3 .
Figure A2.Fraction  of stellar flybys (for initial conditions randomly chosen as described in Appendix A) for which the hybrid (black) and orbit-averaged (red) calculations of the orbital elements changes in an encounter -(top) eccentricity, (middle) inclination, (bottom) longitude of the ascending node and the argument of periapsis -result in > 20% deviations relative to the exact 3-body integration.Initial binary eccentricity is  = 0.1 (left) and  = 0.999 (right).Hybrid calculations in panels (a)-(d) use different values of  3body indicated in panels; the bottom row assumes  3body = 100.

Figure A4 .Figure A5 .
Figure A4.(a)Evolution of the perturber-binary distance  in units of  and (b) binary eccentricity change Δ during an encounter, calculated using the direct 3-body integration ( is time in units such that the binary orbital period is 2 ).Grey bands show an interval, which is zoomed over in panel (c).Black and red curves use the same parameters except for the different binary initial mean anomaly (0 and 85.9 • correspondingly).The red line is not shown in the top two panels for the sake of clarity as it almost coincides with the black one.

Table 1 .
Comparison of the physical processes accounted for in different models of a binary evolution in a stellar cluster Paper Outer orbit Cluster tides Cluster evolution Dynamical friction Strong encounters Weak encounters