Issue 
A&A
Volume 683, March 2024



Article Number  A186  
Number of page(s)  14  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/202348745  
Published online  18 March 2024 
Eccentric black hole mergers via threebody interactions in young, globular, and nuclear star clusters
^{1}
Physics and Astronomy Department Galileo Galilei, University of Padova, Vicolo dell’Osservatorio 3, 35122 Padova, Italy
email: marco.dallamico@phd.unipd.it
^{2}
INFNPadova, Via Marzolo 8, 35131 Padova, Italy
^{3}
Institut für Theoretische Astrophysik, ZAH, Universität Heidelberg, AlbertUeberleStr. 2, 69120 Heidelberg, Germany
email: mapelli@uniheidelberg.de
^{4}
INAFOsservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy
^{5}
Gran Sasso Science Institute (GSSI), 67100 L’Aquila, Italy
Received:
27
November
2023
Accepted:
22
January
2024
Eccentric mergers are a signature of the dynamical formation channel of binary black holes (BBHs) in dense stellar environments and hierarchical triple systems. Here, we investigate the formation of eccentric mergers via binarysingle interactions by means of 2.5 × 10^{5} direct Nbody simulations. Our simulations include postNewtonian terms up to the 2.5th order and model the typical environment of young (YSCs), globular (GCs), and nuclear star clusters (NSCs). Around 0.6% (1%) of our mergers in NSCs (GCs) have an eccentricity > 0.1 when the emitted gravitational wave frequency is 10 Hz in the source frame, while in YSCs this fraction rises to 1.6%. Approximately ∼63% of these mergers are produced by chaotic, resonant interactions where temporary binaries are continuously formed and destroyed, while ∼31% arise from an almost direct collision of two black holes (BHs). Lastly, ∼6% of these eccentric mergers occur in temporary hierarchical triples. We find that binaries undergoing a flyby generally develop smaller tilt angles with respect to exchanges. This result challenges the idea that perfectly isotropic spin orientations are produced by dynamics. The environment dramatically affects BH retention: 0%, 3.1%, and 19.9% of all the remnant BHs remain in YSCs, GCs, and NSCs, respectively. The fraction of massive BHs also depends on the host cluster properties, with pairinstability (60 ≤ M_{BH}/M_{⊙} ≤ 100) and intermediatemass (M_{BH} ≥ 100 M_{⊙}) BHs accounting for approximately ∼44% and 1.6% of the mergers in YSCs, ∼33% and 0.7% in GCs, and ∼28% and 0.4% in NSCs, respectively.
Key words: black hole physics / gravitational waves / methods: numerical / stars: black holes / stars: kinematics and dynamics / galaxies: star clusters: general
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Binarysingle encounters dominate the interactions between black holes (BHs) in the core of star clusters (Heggie 1975; Hut & Bahcall 1983; Hut 1983, 1993; Banerjee et al. 2010). In this region, BHs form a dynamically decoupled subcore where they can mostly interact via binarysingle scatter due to the large crosssection of this process (e.g. Breen & Heggie 2013a,b; Samsing et al. 2014). This springs from star clusters’ tendency to evolve towards energy equipartition (Spitzer 1969; Trenti & van der Marel 2013; Spera et al. 2016; Bianchini et al. 2016), combined with dynamical friction acting on the most massive bodies of the cluster (Meylan & Heggie 1997; Fregeau et al. 2002; Gürkan et al. 2004).
Most binary black holes (BBHs) in star clusters belong to the family of hard binaries; that is, binary systems with a binding energy larger than the average star kinetic energy of the cluster. Since hard binaries statistically tend to get harder during binarysingle encounters (Heggie 1975), BBHs tend to progressively decrease their semimajor axis, or even increase their total mass if a dynamical exchange with a single BH takes place (Hills & Fullerton 1980). Threebody interactions^{1} between BHs are therefore a key mechanism to speed up the merger before gravitational wave (GW) emission becomes efficient. These dynamical encounters not only efficiently produce BBH mergers (e.g. Sigurdsson & Phinney 1995; Portegies Zwart & McMillan 2000; Banerjee et al. 2010; Tanikawa 2013; Ziosi et al. 2014; Morscher et al. 2015; Rodriguez et al. 2015, 2016a,b; Mapelli 2016; Samsing et al. 2017, 2018a; Samsing & Ilan 2018; Trani et al. 2019, 2021; Dall’Amico et al. 2021) and BHneutron star mergers (e.g. Clausen et al. 2013; Arca Sedda 2020, 2021), but they may also cause eccentric mergers (Gültekin et al. 2006; Samsing et al. 2014, 2018b; Samsing & RamirezRuiz 2017; Samsing 2018; Rodriguez et al. 2018; Zevin et al. 2019; Arca Sedda et al. 2021; Trani et al. 2022; Codazzo et al. 2023). These are mergers in which the coalescence time of the binary is shorter than the timescale it takes for GW emission to circularise the orbit, such that the binary can merge with a nonzero eccentricity in the LIGO–Virgo sensitivity band (Abbott et al. 2019). In dynamical interactions, the energy exchange between the bodies can excite the eccentricity of a BBH and even induce it to merge. This is particularly true for threebody interactions, in which the system can evolve into a chaotic regime with temporary binaries that are continuously created and destroyed. In this unstable triple configuration, the single BH can perturb the temporary binary and induce it to merge rapidly enough that GWs are not sufficient to completely circularise the orbit (Samsing et al. 2014).
Isolated binaries, on the other hand, struggle to produce BBHs with a nonnegligible eccentricity at merger. Tidal effects, mass transfer episodes, and common envelope events usually circularise the orbit of a binary star even before it evolves into a BBH (Hurley et al. 2002). Even if supernova kicks can increase the eccentricity of the system, GWs efficiently circularise the orbits by the merger time (Peters 1964). Eccentric mergers are therefore commonly associated with BBHs formed in a dynamically active environment. Eccentricity, if detected in the waveform of a merger, might be used as a tool to infer the dynamical origin of a BBH (AmaroSeoane & Chen 2016; Chen & AmaroSeoane 2017; Gayathri et al. 2022; RomeroShaw et al. 2020, 2021; Zevin et al. 2021).
How often are these eccentric mergers produced by dynamical interactions? And in which environment should we expect them to be more frequent? Here, we aim to address these questions via direct Nbody simulations of threebody encounters between BBHs and BHs. We performed three different sets of binarysingle scattering experiments, each with different initial conditions appositely set to reproduce the properties of a class of star clusters: young star clusters (YSCs), globular clusters (GCs), and nuclear star clusters (NSCs). Our goal is to investigate the effect of the cluster properties on the interactions and to derive the influence that the hosting environment has on the outcomes and production of eccentric BBH mergers^{2}.
Threebody interactions are the fundamental mechanism at the base of hierarchical mergers; in other words, the process in which two BHs merge and their merger remnant collides with other BHs of the cluster, giving rise to multiple generations of BBHs (Miller & Hamilton 2002; Gerosa & Berti 2017; Fishbach et al. 2017; Rodriguez et al. 2019; Antonini et al. 2019; Doctor et al. 2020; Arca Sedda et al. 2023; Mapelli et al. 2021; Gerosa & Fishbach 2021; Atallah et al. 2023). Here, we discuss the impact of threebody recoil velocities on hierarchical mergers, and how this effect, combined with relativistic kicks and star cluster evaporation, could dynamically eject the BHs from the cluster.
2. Methods
2.1. Direct Nbody simulations
Threebody encounters are chaotic processes (Poincaré 1892; Valtonen & Karttunen 2006). Due to this nature, the orbits of the interacting bodies are highly unpredictable, and no general analytical solutions are known to exist. Potentially, even a perturbation of the Planckscale order applied to the initial conditions can exponentially grow and lead to different final configurations of the system (Samsing & Ilan 2018; Manwadkar et al. 2020; Boekholt et al. 2020, 2021; Portegies Zwart et al. 2022; Parischewsky et al. 2023). Therefore, the most convenient approach to studying the threebody problem from the perspective of BBH mergers is to use a numerical integrator over a large set of interactions and derive the statistical properties of the encounters.
Here, we used the direct Nbody code ARWV to simulate three different sets of threebody interactions, for a total of 2.5 × 10^{5} simulations between a BBH and a single BH. Each of our three sets of simulations was initialised with different initial conditions, designed to reproduce the properties of threebody encounters that take place inside YSCs, GCs, and NSCs. ARWV is an algorithmic regularisation direct Nbody code (Mikkola & Aarseth 1989, 1993; ArcaSedda & CapuzzoDolcetta 2019; Chassonnery et al. 2019; Chassonnery & CapuzzoDolcetta 2021) that solves the equations of motion of the interacting bodies with postNewtonian corrections up to the 2.5 order (Mikkola & Merritt 2008; Memmesheimer et al. 2004).
We integrated each system with ARWV for a minimum time of 10^{5} yr. We stopped the integration at this time only if at least one merger took place, or if the threebody encounter was over. If none of these two conditions was satisfied (e.g. the three bodies were still interacting), we carried on the integration with ARWV for a longer time. We stopped the simulation if at least one of the two aforementioned conditions was fulfilled, or if the time reached a maximum of 1 Myr. Figure 1 shows three examples of our threebody simulations computed with ARWV.
Fig. 1. Trajectories of three eccentric mergers triggered by a threebody interaction in NSCs (left and righthand panels) and GCs (central panel). In the three panels, the systems are centred in the centre of mass of the initial BBH, while the position at which the merger takes place is marked with a white star. The trajectories of the primary and secondary BH of the initial BBH (with mass m_{1} and m_{2}), and the intruder (m_{3}) are reported in red, blue, and grey, respectively. The incoming direction of the intruder is shown with an arrow. The insets in the lefthand and central panels display a closeup view of the merger region. In the lefthand and central panels, the initial binary is shown edgeon, while in the righthand panel the view is faceon. The initial coalescence time of these BBHs, i.e. the coalescence time of the initial binary at the beginning of the interaction, is longer than the Hubble time. 
Most BBH mergers did not take place during the simulation with ARWV. Therefore, we evolved the remaining binary population according to Peters (1964):
where
Here, G is the gravity constant, c the speed of light, m_{i} the primary mass, m_{j} the secondary mass, a the semimajor axis, and e the orbital eccentricity. We assumed that two BHs merge when their mutual distance is less than the sum of their innermost stable circular orbits for nonspinning BHs; that is, when r ≤ 6 G(m_{i} + m_{j})/c^{2}. All the binaries that survive after the end of the threebody integration with ARW are longlived and circularise their orbit through GW emission, unaffected by external perturbations. Hence, we can treat them with the simplified formalism described in Eq. (1), without the need for a more computationally expensive postNewtonian formalism.
When two BHs merged, we implemented the relativistic fitting equations reported in Lousto et al. (2012) to compute the relativistic kick exerted on the BH remnant by the anisotropic GW emission at merger. These are:
where
Here, q = m_{2}/m_{1} with m_{2} < m_{1}, η = q (1 + q)^{−2}, A = 1.2 × 10^{4} km s^{−1}, B = −0.93, H = 6.9 × 10^{3} km s^{−1}, V_{1, 1} = 3678 km s^{−1}, V_{A} = 2481 km s^{−1}, V_{B} = 1792 km s^{−1}, and V_{C} = 1506 km s^{−1}. The equations incorporate the components of the dimensionless spin parameter vectors, χ_{1} and χ_{2}, relative to the primary and secondary BHs. Specifically, χ_{1∥} and χ_{2∥} are the primary and secondary dimensionless spin components parallel to the orbital angular momentum of the binary, while χ_{1⊥} and χ_{2⊥} are the perpendicular spin components lying in the orbital plane. Furthermore, S_{∥} is the component parallel to the orbital angular momentum of the vector S = 2(χ_{1} + q_{2} χ_{2})/(1 + q)^{2}. Moreover, ϕ_{Δ} is the angle between the direction of the infall at merger and the inplane component of the vector Δ = (m_{1} + m_{2})^{2} (χ_{1} + q χ_{2})/(1+ q), while ϕ is the phase of the binary. For further details, we refer to Lousto et al. (2012). Finally, we computed the mass of the BH remnant and the magnitude of its spin with the fitting equations presented by JiménezForteza et al. (2017).
2.2. Initial conditions
Simulating a threebody interaction with spinning BHs requires 21 initial parameters for each simulation. Since covering a 21dimension parameter space with direct Nbody simulations is computationally prohibitive, we initialised three different sets of 5 × 10^{4}, 10^{5}, and 10^{5} simulations for YSCs, GCs, and NSCs, respectively. Table 1 reports all the initial conditions, the distribution used to generate each of them, and the interval in which the parameters are sampled. We initialised the properties of each encounter with the same procedure as in Dall’Amico et al. (2021), except for the masses, semimajor axis, and initial intruder velocity. Here, we summarise the main features of our initial conditions.
Initial conditions.
We assumed that the initial BH mass distribution does not change if the cluster is young, globular, or nuclear (Mapelli et al. 2021), and sampled the mass of all our single and binary BHs from a catalogue of synthetic BHs generated with the population synthesis code MOBSE (Mapelli et al. 2017; Giacobbo et al. 2018; Giacobbo & Mapelli 2018). With this method, our population is composed of firstgeneration BHs produced by the evolution of binary stars at metallicity Z = 0.1 Z_{⊙}, with Z_{⊙} = 0.02. MOBSE implements uptodate wind models for massive stars (Vink et al. 2001; Chen et al. 2015), the corecollapse supernova models by Fryer et al. (2012) and the (pulsational) pairinstability supernova treatment presented in Mapelli et al. (2020). We adopted the rapid corecollapse supernova model by Fryer et al. (2012). With this choice, our BH mass spectrum ranges between 5 and 60 M_{⊙} BHs.
The semimajor axis, a, of the initial binaries was sampled from a lognormal distribution as
with limits [max(a_{ej}, a_{gw}),a_{hard}], where
Here, a_{hard} is the limit for a binary to be considered hard (Heggie 1975), a_{ej} is the maximum semimajor axis for ejection by threebody encounters, and a_{gw} is the limit below which the semimajor axis shrinking by the emission of GWs becomes dominant with respect to dynamical hardening. In the above equations, m_{*} is the average mass of a star in the cluster, σ the typical 3D velocity dispersion of the cluster, v_{esc} ∼ 2σ the escape velocity, ρ_{c} the star cluster core density, and ξ ∼ 3 a numerically calibrated constant (Hills 1983; Quinlan 1996). The mean of the lognormal distribution of the semimajor axes was computed as the average of the logarithmic limits, such that it resulted in μ_{log(a/AU)} = 2.42, 1.22, 0.42 in the cases of YSCs, GCs, and NSCs, respectively. The dispersion was derived from the simulations of Di Carlo et al. (2019, 2020), and set as σ_{log(a/AU)} = 0.92 in all three samples of simulations, as was already done in Dall’Amico et al. (2021). With this prescription, all our BBHs are hard binaries for which GW emission is negligible if compared to hardening, but at the same time their semimajor axis is large enough such that previous interactions did not lead to a dynamical ejection of the binary from the cluster.
We assumed that the BHs were in thermal equilibrium with the rest of the cluster core so that the initial velocity at infinity, v_{∞}, of the single BH with respect to the centre of mass of the binary could be sampled from a MaxwellBoltzmann distribution (Heggie 1975). For the three sets of simulations, we assumed a 3D velocity dispersion of 5, 20, and 50 km s^{−1} in the cases of YSCs (Portegies Zwart et al. 2010), GCs (Pryor & Meylan 1993), and NSCs (Neumayer et al. 2020), respectively.
We set the initial distance of the single BH from the centre of mass of the binary as D = 10^{3} a. This guaranteed that the BBH had not been perturbed by the intruder before the beginning of the integration.
The impact parameter, b, of the interaction was sampled from a uniform probability distribution proportional to b^{2} (Hut & Bahcall 1983) in the range [0, b_{max}], with b_{max} defined as
This is the maximum impact parameter for a strong threebody interaction with a hard binary derived by Sigurdsson & Phinney (1993). As v_{∞} and a change between YSCs, GCs, and NSCs, the impact parameter in these three environments will also be different. Equation (12) assumes strong gravitational focusing, that is, , and a minimum intruderbinary star distance, r_{p} = a. In Appendix A, we discuss the impact of this assumption on our results.
LIGO–Virgo observations favour a population of lowspin BHs (Abbott et al. 2023). Therefore, we extracted the magnitude of the dimensionless spin, χ, of each BH from a MaxwellBoltzmann distribution with root mean square σ_{χ} = 0.1 and truncated to χ = 1. Furthermore, we isotropically sampled the spin directions, accounting for the fact that dynamical encounters randomise them (Rodriguez et al. 2016c).
3. Results
3.1. Outcomes
We divided the outcomes of the interactions into five classes, as a function of the system configuration at 1 Myr since the beginning of the simulation. We classified each simulation as follows:

Flyby: the final binary has the same components as the initial binary. If the encounter hardens it enough, this binary may merge during or after the simulation.

Exch13: the threebody interaction ends with a binary composed of the primary BH of the initial binary and the intruder (m_{1} − m_{3}). If this exchanged binary merges during the simulation, we still classify the interaction as an exch13 event.

Exch23: the threebody interaction ends with a binary composed of the secondary BH of the initial binary and the intruder (m_{2} − m_{3}). If this exchanged binary merges during the simulation, we still classify the interaction as an exch23 event.

Ionisation: the encounter splits the initial binary, resulting in three single BHs.

Triple: the system is still in an unstable triple configuration at the end of the simulation.
The upper panel of Fig. 2 classifies the end states of our binarysingle scattering experiments. Flybys are the most frequent end state in all three cluster types, followed by exchanges. This result is expected, since these encounters generally have a larger impact parameter, b, than the semimajor axis of the initial binary, a: if b ≫ a the intruder sees the binary as a pointlike object, and the interaction evolves into a flyby.
Fig. 2. Upper panel: percentage of different interaction outcomes for each cluster type. From left to right: (i) flyby events, (ii) exchanges in which the intruder replaces the secondary BH, (iii) exchanges in which the intruder replaces the primary BH, (iv) ionisations, and (v) unstable triples. Lower panel: percentage of BBH mergers. From left to right: (i) BBH mergers occurring after a flyby, (ii) an exchange interaction with the secondary BH replaced by the intruder, (iii) an exchange with the primary BH replaced by the intruder, and (iv) secondgeneration BBH mergers. In both the upper and lower panels, the colours mark the cluster type in which the interaction takes place. 
In YSCs, the number of ionisations is lower with respect to both GCs and NSCs. Vice versa, exchanges are more common in YSCs than in both GCs and NSCs. This happens because the typical dispersion velocities of YSCs are around 5 km s^{−1}, much smaller than GCs and NSCs, where the intruder can likely have a velocity higher than the critical velocity required to break up the binary system (Hut & Bahcall 1983). At 1 Myr, unstable triples are much more numerous in YSCs than in GCs, while we find no triple systems in NSCs at the end of our simulations. In YSCs, the encounter takes place later than in more dense clusters, since the interparticle distance is much larger while the dispersion velocity is lower. The intruder takes more time to reach the binary and the interaction begins at later times in the simulation, resulting in several systems that at 1 Myr are still in an unstable triple configuration.
3.2. BBH mergers
We find that 0.2%, 2.4%, and 11.8% of the simulations produce BBHs that merge within a Hubble time (13.8 Gyr) in YSCs, GCs, and NSCs, respectively. Of these mergers, 14.8%, 4.6%, and 4% take place in the first Myr of integration with ARWV. These results, which are also reported in Table 2, imply that the denser and more massive the cluster is, the larger the fraction of mergers produced via threebody interactions. This happens because the minimum binding energy of a hard binary is higher in more massive clusters. As a result, BBHs in NSCs are consistently closer to the GW regime than in YSCs, making it more likely for a single interaction to push them into the GW emission regime. The lower panel of Fig. 2 classifies these mergers by their formation channel: flyby, exch13, or exch23, and secondgeneration mergers (i.e., mergers that occur between the remnant of a previous merger and the third BH). The bar plot shows that flybys account for ∼53% of all our encounters, but ∼61 − 64% of the BBH mergers. Hence, flybys are more efficient in inducing mergers than, for example, exch23 events, which in turn are ∼11% of the outcomes but account only for ∼7% of the mergers. This happens because exchanges usually produce new BBHs with a larger total mass but also with a larger semimajor axis than BBHs involved in flyby events.
Percentage of peculiar events in YSCs, GCs, and NSCs.
For some systems that produce a BBH merger, we find that the initial binary would have merged within a Hubble time even without the interaction. These are binaries that, if evolved as an unperturbed binary, would have produced a BBH merger without undergoing any threebody encounter. These systems are 1.6%, 11.4%, and 45.0% of the BBH mergers in YSCs, GCs, and NSCs respectively. Even for these systems, the interaction has efficiently sped up the merger. This can be seen in Fig. 3, which shows the coalescence time of the unperturbed binaries (if the threebody interaction would not have taken place), and the coalescence time resulting from the threebody simulation. For instance, all the BBH mergers that take place during the simulation with ARWV within the first 1 Myr rapidly merge solely as a consequence of the threebody interaction.
Fig. 3. Distribution of coalescence times for unperturbed initial binaries (filled histograms), and for the same initial binaries perturbed by the threebody interaction (unfilled histograms). The upper, central, and lower panels show the cases of YSCs, GCs, and NSCs, respectively. 
Finally, in a few cases, we find the formation of secondgeneration BBH mergers. In these simulations, two BHs merge during the threebody interaction, and their remnant forms a new BBH with the remaining BH, which in turn is able to reach coalescence in less than a Hubble time. We find the same percentage of secondgeneration mergers in GCs and NSCs, while no secondgeneration systems form in YSCs via threebody interactions (Fig. 2).
Despite the three simulation sets being initialised with the same BH mass spectrum, the total mass of the BBH mergers produced in YSCs, GCs, and NSCs differ, as is shown by Fig. 4. This difference in the mass of BBH mergers is clearly an effect of the environment.
Fig. 4. Total mass of the BBH mergers in YSCs (blue), GCs (orange), and NSCs (purple). 
Mergers by dynamical exchange are slightly favoured in GCs and NSCs with respect to YSCs (Fig. 2). Since exchanges typically take place when the intruder is more massive than at least one of the two binary components, we should expect more massive BBH mergers in GCs and NSCs. Nevertheless, Fig. 4 shows that YSCs are more likely to produce massive mergers via threebody encounters than the other two sets of simulations. This happens because in YSCs only the most massive systems are able to merge within a Hubble time. In GCs and NSCs, given the larger velocity dispersions, the hardsoft boundary is shifted towards the lower semimajor axis (see Eq. (9)). Therefore, GW emission can also be efficient for relatively lowmass BBHs (Eq. (1)). In YSCs, instead, given the larger semimajor axis at formation, mergers are more biased towards the most massive BBHs. In this way, NSCs and GCs are more efficient in the production of BBH mergers with lowmass components than YSCs. This is further confirmed if we look at the percentage of massive BH remnants produced by these mergers in Table 2: pairinstability BHs, defined as BHs with mass in the range 60 − 100 M_{⊙}, are 44.3%, 33.3%, and 27.7% over all the mergers produced in YSCs, GCs, and NSCs, respectively. Intermediatemass BHs (IMBHs, i.e. BHs with a mass ≥10^{2} M_{⊙}) are 1.6%, 0.7%, and 0.4% over all the mergers in YSCs, GCs, and NSCs, respectively. The most massive remnant produced by a merger is ∼110 M_{⊙} in all three sets. In YSCs and GCs, the most massive remnants are produced by a flyby event, while in the NSC case the most massive remnant is produced by the merger of a secondgeneration BBH.
Pairinstability BHs and IMBHs form mainly via flybys, because these are the most frequent type of interaction among our simulations. Exch13 are the events most likely to create a BBH more massive than the initial binary. On the other hand, exch23 tend to produce less massive binaries than the initial ones. Nevertheless, the latter have a minor impact on the total mass distribution in Fig. 4 because they are rare events, representing only ≈7% of all mergers (Fig. 2). Finally, mergers of secondgeneration binaries are rare even among pairinstability BHs and IMBHs: they represent 0.25% and 0% of pairinstability BHs and IMBHs in GCs, and 0.12% and 4% in NSCs, respectively.
3.3. Eccentric mergers
Figure 5 shows the eccentricity of the BBH when the emitted GW frequency is 10 Hz in the source frame (hereafter, e_{10 Hz}) as a function of the coalescence time from the beginning of the simulation. All three sets present two distinct families of mergers regarding their eccentricity. In the first one, and also the most common one, the mergers follow the relation, t_{coal}(e), of Eq. (1), such that e_{10 Hz} decreases as the coalescence time of the binary increases. In these systems, the dynamical encounter ends before the merger takes place. This allows the BBH resulting from the interaction to evolve unperturbed up to the merger, such that the eccentricity evolution is ruled only by the angular momentum loss by GW emission^{3}.
Fig. 5. Orbital eccentricity at 10 Hz (e_{10 Hz}) as a function of the coalescence time (t_{coal}) of BBH mergers in YSCs (left), GCs (centre), and NSCs (right). Different colours indicate mergers through flybys and the two families of exchanges. The yellow stars show secondgeneration BBH mergers. The empty black dots highlight the BBH progenitors of secondgeneration mergers. The dashed vertical line divides the plot into two regions: up to 1 Myr we evolve the systems with ARWV, above 1 Myr all the dynamical interactions are concluded, and with Eq. (1) we stop the ARWV integration and evolve the remaining BBHs. The three insets show a zoom for mergers with e_{10 Hz} > 0.1. We note that some of the points at e_{10 Hz} ∼ 1 overlap. 
The second family of mergers is composed by systems that do not follow the relation of Eq. (1), but rather reach the merger with e_{10 Hz} > 10^{−3} in a relatively short time during the simulation with ARWV, while the dynamical interaction is still ongoing. Some of these systems even have e_{10 Hz} > 0.1 (insets of Fig. 5). In the following, we refer to the mergers with e_{10 Hz} > 0.1 simply as eccentric mergers. These are 1.6%, 1%, and 0.6% of all the mergers in YSCs, GCs, and NSCs, respectively (Table 2). Although YSCs exhibit a larger fraction of eccentric mergers, the only two eccentric events occurring in YSCs both have e_{10 Hz} ∼ 0.1. Conversely, eccentric events in GCs and NSCs span eccentricities up to and beyond e_{10 Hz} = 0.9. We can divide these eccentric mergers as a function of the type of interaction that triggered the coalescence:

Chaotic mergers are the product of threebody interactions in which temporary binaries with brandnew orbital parameters are continuously formed and destroyed. If the eccentricity of these temporary binaries is sufficiently high, and their lifetime is longer than the perturbation timescale of the outer body, a nearly radial merger can be triggered. These are the most common interactions to produce eccentric mergers, accounting for ∼63% of all the eccentric events with e_{10 Hz} > 0.1 in GCs and NSCs. The lefthand panel of Fig. 1 shows one of these interactions and the subsequent eccentric merger.

Prompt mergers are the second most common event to cause eccentric mergers, representing ∼31% of all these events in GCs and NSCs. The righthand panel of Fig. 1 shows an example of an eccentric merger that follows a prompt interaction. These mergers typically follow flyby events in which the intruder significantly extracts angular momentum from the binary, driving the two components into a nearly radial orbit and inducing a prompt merger. This is for example the case of the two eccentric mergers in the YSC set. We also call prompt mergers simulations in which the intruder tangentially intersects the orbital plane of the binary and approaches one of the two components with an almost antiparallel velocity, such that they rapidly merge in a nearly headon collision (as the simulation in the righthand panel of Fig. 1 demonstrates). This is the case of the most eccentric mergers in our simulations. In NSCs, five headon collisions trigger a merger, with e_{10 Hz} ∼ 1 (righthand panel of Fig. 5), while in GCs the maximum value of e_{10 Hz} is 0.87 (central panel Fig. 5).

Temporary triple mergers take place when the system evolves as a hierarchical triple. Stable triple systems cannot form from three isolated unbound bodies (Naoz 2016); however, temporary stable hierarchical triples can be created via threebody interactions. In this configuration, the intruder sets in an outer orbit, perturbing the initial binary and causing it to merge rapidly enough that GW emission does not efficiently circularise the binary’s orbit. This is the case for the system shown in the central panel of Fig. 1. The merger occurs only if the system remains stable for a sufficient period of time for the perturbations to be effective. Binary systems involved in chaotic mergers exhibit considerably shorter lifespans than those within temporary triple systems. In contrast, temporary triple systems can survive for several orbital periods of the third body, which has enough time to orbit around the inner binary perturbing it. We define an eccentric merger as a temporary triple merger if the merger event involves the inner binary and if the outer orbiting third body is able to carry out at least two complete stable orbits with approximately the same period around the inner binary. If these conditions are not met, we define the merger as a chaotic merger. Due to the low stability of these systems, temporary triple mergers are responsible only for ∼6% of the eccentric mergers in GCs and NSCs.
Even if flybys are the most common formation path of BBH mergers in our simulations (lower panel of Fig. 2), eccentric mergers come in almost the same proportion from flybys and exchanges in both GCs and NSCs. Hence, eccentric mergers are more likely to form from exchanges. No secondgeneration BBH mergers belong to eccentric mergers. With our assumption of no further dynamical interactions beyond 1 Myr, secondgeneration binaries, after their formation, circularise their orbits before reaching coalescence. On the other hand, all the mergers that give birth to one of the two BHs that compose secondgeneration binaries are eccentric mergers. We refer to these systems as BBH progenitors of secondgeneration mergers. In the GC scenario, the progenitor systems of the two secondgeneration BBHs are both eccentric mergers with an eccentricity of 0.82 and 0.67, respectively, with the first being an almost headon collision in a prompt event, and the latter coming from a chaotic merger. This is also true for the NSC case, in which all the nine progenitors have e_{10 Hz} > 0.1, four of which merge in a headon collision with e_{10 Hz} ∼ 1.
Finally, Fig. 6 shows the same systems as Fig. 5 but with the eccentricity computed at the GW frequency of 1 Hz in the source frame. This frequency is the predicted lower limit in the sensitivity band of the Einstein Telescope (Punturo et al. 2010; Maggiore et al. 2020). In this regime, eccentric mergers constitute 0.9%, 1.8%, and 9% of all mergers in NSCs, GCs, and YSCs, respectively (Table 2). This is because, at lower frequencies, the two BHs are at a larger orbital separation and GWs still have to efficiently circularise the orbit, such that more systems have higher eccentricities.
3.4. Dynamical ejections
Our BBHs can escape from the cluster as a consequence of dynamical recoil. Threebody interactions with hard binaries tend to reduce the binary internal energy and convert it into kinetic energy of the system. This induces a recoil velocity, v_{rec}, both on the binary and the single object. If this velocity is larger than the escape velocity, v_{esc}, of the cluster, both the binary and the single object are dynamically ejected from the cluster (Heggie 1975; Hills 1975; Sigurdsson & Phinney 1993).
This has major implications for BBH mergers: on the one hand, the energy exchange speeds the merger up, reducing the semimajor axis and thus the coalescence time, t_{coal}; on the other hand, it might also kick out the binary from the cluster, preventing further encounters. If the ejection happens before the binary merges, the BH remnant produced by the coalescence will not be able to dynamically interact with other bodies of the cluster and produce new binaries. This may strongly affect the efficiency of the hierarchical merger mechanism by which repeated BH mergers produce massive BHs (e.g. Miller & Hamilton 2002; Mapelli et al. 2021). Figure 7 shows the impact of these dynamical ejections on BBH mergers. Here, we assume that all the mergers that take place during the threebody encounter merge inside the cluster. Therefore, the plots report only the BBH mergers that take place after 10^{5} yr and for which the threebody interaction is concluded. These recoil velocities span from less than 1 km s^{−1} up to ∼40 km s^{−1} in YSCs, ∼200 km s^{−1} in GCs, and ∼400 km s^{−1} in NSCs. These differences are, once again, explained by the large binding energy of the binaries in GCs and NSCs, which translates into larger recoil kicks. In YSCs, for example, the BBHs have a larger semimajor axis than binaries in more dense clusters and the recoil velocities are lower.
Fig. 7. BBH mergers that take place after the threebody interactions in YSCs, GCs, and NSCs (from top to bottom). All the mergers that take place during the threebody interaction, i.e. < 10^{6} yr, are assumed to merge inside the cluster. The x axis reports the recoil velocity caused by the energy exchange in the threebody encounter. The y axis shows the coalescence time of the binary from the beginning of the simulation. Different outcomes are shown with different markers and colours. In the YSCs case, the horizontal dashed line reports the typical evaporation time of a YSC. The vertical dashed line in all the plots shows the escape velocity of the cluster. Distributions of the recoil velocities are displayed as marginal histograms. The colour legend is the same as the scatter plot. 
In the three simulation sets, the most likely population of ejected mergers is the one produced by exch23 events. The marginal plots of Fig. 7 and Table 3 show that these BBHs generally have larger recoil velocities with respect to the other two families of mergers. The fraction of BBH mergers ejected after the threebody interaction is 6% in NSCs, 18% in GCs, and 18% in YSCs. In YSCs, we must also keep into account the evaporation of the cluster; that is, when the cluster dissolves because of stellar mass loss and tidal stripping (Spitzer 1987; Heggie & Hut 2003; Binney & Tremaine 2008). Due to evaporation, most of the BBH mergers in YSCs happen in the field even without being dynamically ejected (e.g. Rastello et al. 2021; Torniamenti et al. 2022). If we assume a typical evaporation time, t_{evap} ∼ 1 Gyr for a YSC^{4} (Torniamenti et al. 2022), 48% of the BBH mergers happen inside the cluster and 5% take place in the field because of dynamical recoil, while the remaining 47% occur in the field because of cluster evaporation. Finally, all the eccentric mergers in GCs and NSCs merge inside the cluster. They are mostly the product of chaotic interactions that take place over a short timescale and rapidly lead to the merger of two of the three BHs.
Threebody recoil velocities.
GWs emitted from a spiraling BBH are generally irradiated anisotropically due to the asymmetry of the system. This induces linear momentum transfer on the remnant BH produced by the merger, which translates into a relativistic kick that might accelerate the remnant even up to a few thousand of km s^{−1} (Fitchett 1983; Maggiore 2018). Figure 8 shows this relativistic kick produced by all the BBH mergers in our three sets of simulations. All three distributions peak at ∼200 km s^{−1}, with velocities that are approximately one order of magnitude higher than the ones reported in Fig. 7. Hence, GW recoils are more efficient in ejecting BHs from their parent cluster than threebody recoils. Due to GW kicks, only 3.6% and 21.2% of the BH remnants are retained in GCs and NSCs, respectively. This fraction drops to 0 in YSCs.
Fig. 8. Distribution of the relativistic recoil kick of the remnants produced by the BBH mergers in YSCs, GCs, and NSCs. The dashed vertical lines represent the escape velocities of the three types of star clusters. 
We can now count the overall fraction of BH remnants retained by the cluster for which (1) the relativistic recoil is below the escape velocity of the cluster, (2) the BBH progenitor is not ejected after a threebody interaction, and (3) the binary is able to merge inside the cluster before its evaporation. This fraction is 0% for YSCs, 3.1% for GCs, and 19.9% for NSCs (Table 2). Finally, all the secondgeneration BH remnants in GCs and NSCs are kicked out of the cluster due to GW recoils.
3.5. Orbital plane tilt
Figure 9 shows the tilt angle, i, defined as the angle between the orbital plane of the initial binary and the orbital plane of the final binary left after the threebody interaction. These distributions show that threebody interactions are capable of inducing large tilt angles in the BBH population, in good agreement with what was found by Trani et al. (2021, see also Banerjee et al. 2023). The magnitude of the tilt depends on the outcome of the interaction. Binaries that underwent a flyby generally experience smaller tilt angles with respect to exchanges, with the distribution peaking at ∼15° for the former and ∼90° for the latter. Since there is no strong correlation between the initial angle, θ (Sect. 2.2), of the interaction and the outcome, the distributions of Fig. 9 are a direct product of the interactions. This means that flybys statistically induce small perturbations in the orbital plane of the initial binary if compared to exchanges. Exchanges, on the other hand, favour the production of more massive newborn binaries with an orbital plane likely tilted with respect to the original one. This implies that the orientation of spins in dynamically assembled BBHs might not be perfectly isotropic. Our finding is consistent with Bouffanais et al. (2019), who assumes an isotropic spin orientation for exchanges, but nearly aligned spins for flybys, which are less perturbed by the encounter.
Fig. 9. Tilt angle distributions of the BBH orbital plane in NSCs at the end of the simulation with respect to the initial binary orbital plane. Unfilled histograms are flybys (red), exch13 (blue), and exch23 (grey). The filled histogram (yellow) shows the total distribution. We do not show GCs and YSCs because they have exactly the same behaviour. 
4. Discussion and caveats
4.1. Merger rate density of eccentric mergers
In Sect. 3.3, we derived the fraction of eccentric mergers in YSCs, GCs, and NSCs; that is, those that have e_{10 Hz} > 0.1. Through semianalytical simulations, Mapelli et al. (2022) were able to derive the merger rate density of BBH mergers in YSCs, GCs, and NSCs. They find ℛ_{YSC} = 3.0 Gpc^{−3} yr^{−1}, ℛ_{GC} = 4.4 Gpc^{−3} yr^{−1}, and ℛ_{NSC} = 1.3 Gpc^{−3} yr^{−1} in their fiducial model. These rates refer to firstgeneration BBH mergers. Since our firstgeneration BHs are drawn from the same population as Mapelli et al. (2022), we combined these rates with the fraction of eccentric mergers produced by threebody interactions in our three simulation sets. In this way, we find the following rates of eccentric mergers as a function of the type of cluster: ∼2 × 10^{−3} Gpc^{−3} yr^{−1} in YSCs, ∼4 × 10^{−2} Gpc^{−3} yr^{−1} in GCs, and ∼8 × 10^{−3} Gpc^{−3} yr^{−1} in NSCs. These rates represent a rough estimate, since a more realistic evaluation requires taking into account the possible differences between the coalescence times of our eccentric mergers and the BBH merger population of Mapelli et al. (2022).
4.2. Single versus multiple threebody interactions
Our eccentric mergers (Fig. 5) qualitatively match the properties of eccentric mergers simulated by Samsing & RamirezRuiz (2017) and Zevin et al. (2019). For example, Samsing & RamirezRuiz (2017) and Samsing et al. (2018a) find that ∼1% of all the mergers by binarysingle encounters in GCs are eccentric mergers, which is in excellent agreement with the fraction of eccentric events that we find in this work. On the other hand, Samsing (2018) and Rodriguez et al. (2018) find a fraction of eccentric mergers ∼5 times higher than the one we find in our simulations. This difference is a consequence of our choice to simulate one single interaction for each binary.
Assuming a single dynamical encounter per binary is particularly well founded in YSCs, since threebody interactions between BHs are expected to happen with a frequency of ∼1 per cluster (Di Carlo et al. 2019). In more dense clusters, it is more likely for a BBH to experience multiple interactions during its lifetime. However, as we presented in Sect. 3.3 and in Fig. 5, all the eccentric mergers take place while the dynamical interaction is still ongoing, meaning that in our simulations the typical timescale for an eccentric threebody merger is smaller compared to the timescale of a threebody interaction. Our results imply that even one single threebody encounter in a BBH lifetime is sufficient to speed up the merger and to produce an eccentric event in the LIGO–Virgo band. Nevertheless, the masses, eccentricities, and fraction of ejected BHs reported in this work must be interpreted as lower limits, since further interactions might produce more massive and eccentric mergers and more dynamical ejections. This is especially true for all the mergers that take place after 1 Myr in our simulations. The effects of multiple threebody interactions on the BBH population are going to be presented in a future paper.
Finally, our work does not consider the special case of an active galactic nucleus (AGN) disc. In AGN discs, eccentric mergers can be significantly boosted if binarysingle scatters occur with small mutual inclinations (of less than a few degrees). This ultimately leads to a relatively flat distribution of the spintilt angle in eccentric mergers (Samsing et al. 2022), which is very different from the one we find here for NSCs. The reason for this difference is that threebody interactions in star clusters require completely distinct initial conditions with respect to the AGN case. We will include AGN discs in forthcoming studies.
5. Summary
In this work, we have presented the results of 2.5 × 10^{5} threebody simulations performed via direct Nbody integration with the ARWV code (Mikkola & Aarseth 1989; Chassonnery et al. 2019). Our simulations incorporate postNewtonian corrections up to the 2.5th order and adopt initial conditions that mimic the properties of YSCs, GCs, and NSCs. With this approach, we aim to investigate the influence of the host environment on: (1) the outcomes of threebody encounters, (2) the populations of BBH mergers produced through interactions, and (3) the production of BBH mergers with nonnegligible eccentricities in the LIGO–Virgo frequency range. Our results can be summarised as follows:

We divide the outcomes into flybys, exchanges in which the primary or the secondary BH component is replaced by the intruder, ionisations, and triples. Flybys dominate the interactions in all the simulation sets, accounting for approximately ∼53% of all the outcomes. YSCs differ from GCs and NSCs, with fewer ionisations (around ∼5% compared to ∼14% for GCs and NSCs) but more exchanges (about ∼41% compared to ∼33% for GCs and NSCs), and also a nonzero number of systems that are in an unstable triple configuration at the end of the simulation.

Threebody interactions in GCs and NSCs produce a higher number of mergers compared to YSCs. Approximately 2.4% and 11.8% of the simulations in GCs and NSCs, respectively, lead to a BBH merger within a Hubble time, compared to the 0.2% of the simulations in YSCs. Flybys are the most effective pathway to produce mergers as they significantly decrease the coalescence time. Of the three types of clusters we considered, YSCs are less efficient at producing lowmass BBH mergers than both GCs and NSCs.

Pairinstability BH remnants (60 − 100 M_{⊙}) are ∼44%, ∼33%, and ∼28% of all the mergers produced in YSCs, GCs, and NSCs, while IMBHs (> 100 M_{⊙}) are 1.6%, 0.7%, and 0.4%, respectively. Finally, we find secondgeneration BBH mergers only in GCs and NSCs, accounting for 0.08% of all mergers in both sets.

The percentage of BBH mergers with an orbital eccentricity higher than 0.1 at a GW frequency of 10 Hz in the source frame (e_{10 Hz} > 0.1) is 1.6%, 1%, and 0.6% in YSCs, GCs, and NSCs, respectively. We additionally investigated the percentage of BBH mergers with an eccentricity exceeding 0.1 at a GW frequency of 1 Hz in the source frame. This frequency corresponds to the lower boundary of the sensitivity band of the Einstein Telescope. We find that 9% of all the BBH mergers in YSCs are eccentric in this frequency range, while this fraction drops to 1.8% for mergers in GCs and 0.9% for mergers in NSCs.

In both GCs and NSCs, the most frequent interactions leading to eccentric mergers are chaotic exchange events, accounting for approximately 63% of all eccentric mergers. These involve the creation and destruction of several temporary binaries before the merger takes place. Prompt interactions, including flybys in which the intruder extracts enough angular momentum from the binary to cause a radial merger, and headon collisions between the intruder and one of the binary components, account for approximately 31% of eccentric events. Finally, BBH mergers in temporary stable hierarchical triples contribute to approximately 6% of all eccentric mergers in GCs and NSCs. In our simulations, all the progenitors of secondgeneration BBHs are eccentric mergers in both GCs and NSCs.

The percentage of remnant BHs that are not expelled from the cluster is 0% for YSCs, 3.1% for GCs, and 19.9% for NSCs. These are BHs that are not dynamically ejected from the cluster by the threebody and GW relativistic recoil kicks, and for which the progenitor BBHs merge before the evaporation of the star cluster. In YSCs, ∼50% of the BBH mergers take place in the field after the cluster has evaporated. Relativistic recoil kicks due to anisotropic GW emission are the primary cause of dynamical ejections, with typical velocities that exceed 100 km s^{−1}. This strongly affects hierarchical mergers.

Threebody interactions alter the inclination of the original orbital plane, causing tilt angles that are not isotropically distributed, but that rather depend on the interaction outcome. Exchanges tend to produce new binary systems that have an isotropically oriented orbital plane with respect to the original one, while flybys usually result in relatively minor perturbations of ∼15° on the orbital plane. This result challenges the idea that dynamics produces perfectly isotropic spin orientations.

We estimate the merger rate density of eccentric BBH mergers to be ∼2 × 10^{−3} Gpc^{−3} yr^{−1} for YSCs, ∼4 × 10^{−2} Gpc^{−3} yr^{−1} for GCs, and ∼8 × 10^{−3} Gpc^{−3} yr^{−1} for NSCs. These rates must be regarded as lower limits, as we only considered a single threebody interaction per binary in our simulations. Additional dynamical interactions during the lifetime of these binaries may lead to an increase in the number of eccentric mergers.
The data underlying this article are available at the following Zenodo link: https://zenodo.org/records/7684085
This evaporation timescale refers to YSCs with mass ∼10^{4} M_{⊙} and must be considered as an upper limit. Processes like galactic perturbations and encounters with giant molecular clouds might, in principle, accelerate the disruption of the cluster (e.g. Gieles et al. 2006).
Acknowledgments
We are grateful to Roberto Capuzzo Dolcetta, Pauline Chassonnery, and Seppo Mikkola for making the ARWV code available to us. We thank Alessandro Trani for the useful comments. We also thank the members of the DEMOBLACK team for the helpful discussions. M.D. acknowledges financial support from Cariparo Foundation under grant 55440. M.M. and S.T. acknowledge financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK under contract no. 770017. M.A.S. acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie SkłodowskaCurie grant agreement No. 101025436 (project GRACEBH, PI: Manuel Arca Sedda).
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, ApJ, 883, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, Phys. Rev. X, 13, 041039 [Google Scholar]
 AmaroSeoane, P., & Chen, X. 2016, MNRAS, 458, 3075 [NASA ADS] [CrossRef] [Google Scholar]
 Antonini, F., Gieles, M., & Gualandris, A. 2019, MNRAS, 486, 5008 [NASA ADS] [CrossRef] [Google Scholar]
 Arca Sedda, M. 2020, Commun. Phys., 3, 43 [CrossRef] [Google Scholar]
 Arca Sedda, M. 2021, ApJ, 908, L38 [Google Scholar]
 ArcaSedda, M., & CapuzzoDolcetta, R. 2019, MNRAS, 483, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Arca Sedda, M., Li, G., & Kocsis, B. 2021, A&A, 650, A189 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arca Sedda, M., Mapelli, M., Benacquista, M., & Spera, M. 2023, MNRAS, 520, 5259 [NASA ADS] [CrossRef] [Google Scholar]
 Atallah, D., Trani, A. A., Kremer, K., et al. 2023, MNRAS, 523, 4227 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S., Baumgardt, H., & Kroupa, P. 2010, MNRAS, 402, 371 [Google Scholar]
 Banerjee, S., Olejak, A., & Belczynski, K. 2023, ApJ, 953, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Bianchini, P., van de Ven, G., Norris, M. A., Schinnerer, E., & Varri, A. L. 2016, MNRAS, 458, 3644 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton: Princeton University Press) [Google Scholar]
 Boekholt, T. C. N., Portegies Zwart, S. F., & Valtonen, M. 2020, MNRAS, 493, 3932 [NASA ADS] [CrossRef] [Google Scholar]
 Boekholt, T. C. N., Moerman, A., & Portegies Zwart, S. F. 2021, Phys. Rev. D, 104, 083020 [CrossRef] [Google Scholar]
 Bouffanais, Y., Mapelli, M., Gerosa, D., et al. 2019, ApJ, 886, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Breen, P. G., & Heggie, D. C. 2013a, MNRAS, 432, 2779 [NASA ADS] [CrossRef] [Google Scholar]
 Breen, P. G., & Heggie, D. C. 2013b, MNRAS, 436, 584 [CrossRef] [Google Scholar]
 Chassonnery, P., & CapuzzoDolcetta, R. 2021, MNRAS, 504, 3909 [NASA ADS] [CrossRef] [Google Scholar]
 Chassonnery, P., CapuzzoDolcetta, R., & Mikkola, S. 2019, arXiv eprints [arXiv:1910.05202] [Google Scholar]
 Chen, X., & AmaroSeoane, P. 2017, ApJ, 842, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068 [Google Scholar]
 Clausen, D., Sigurdsson, S., & Chernoff, D. F. 2013, MNRAS, 428, 3618 [CrossRef] [Google Scholar]
 Codazzo, E., Di Giovanni, M., Harms, J., Dall’Amico, M., & Mapelli, M. 2023, Phys. Rev. D, 107, 023023 [NASA ADS] [CrossRef] [Google Scholar]
 Dall’Amico, M., Mapelli, M., Di Carlo, U. N., et al. 2021, MNRAS, 508, 3045 [CrossRef] [Google Scholar]
 Di Carlo, U. N., Giacobbo, N., Mapelli, M., et al. 2019, MNRAS, 487, 2947 [NASA ADS] [CrossRef] [Google Scholar]
 Di Carlo, U. N., Mapelli, M., Bouffanais, Y., et al. 2020, MNRAS, 497, 1043 [NASA ADS] [CrossRef] [Google Scholar]
 Doctor, Z., Wysocki, D., O’Shaughnessy, R., Holz, D. E., & Farr, B. 2020, ApJ, 893, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Fishbach, M., Holz, D. E., & Farr, B. 2017, ApJ, 840, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Fitchett, M. J. 1983, MNRAS, 203, 1049 [NASA ADS] [CrossRef] [Google Scholar]
 Fregeau, J. M., Joshi, K. J., Portegies Zwart, S. F., & Rasio, F. A. 2002, ApJ, 570, 171 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [Google Scholar]
 Gayathri, V., Healy, J., Lange, J., et al. 2022, Nat. Astron., 6, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Gerosa, D., & Berti, E. 2017, Phys. Rev. D, 95, 124046 [NASA ADS] [CrossRef] [Google Scholar]
 Gerosa, D., & Fishbach, M. 2021, Nat. Astron., 5, 749 [NASA ADS] [CrossRef] [Google Scholar]
 Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011 [Google Scholar]
 Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959 [Google Scholar]
 Gieles, M., Portegies Zwart, S. F., Baumgardt, H., et al. 2006, MNRAS, 371, 793 [Google Scholar]
 Gültekin, K., Miller, M. C., & Hamilton, D. P. 2006, ApJ, 640, 156 [CrossRef] [Google Scholar]
 Gürkan, M. A., Freitag, M., & Rasio, F. A. 2004, ApJ, 604, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D. C. 1975, MNRAS, 173, 729 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D., & Hut, P. 2003, The Gravitational MillionBody Problem: A Multidisciplinary Approach to Star Cluster Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Hills, J. G. 1975, AJ, 80, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Hills, J. G. 1983, AJ, 88, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Hills, J. G., & Fullerton, L. W. 1980, AJ, 85, 1281 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
 Hut, P. 1983, ApJ, 268, 342 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1993, ApJ, 403, 256 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P., & Bahcall, J. N. 1983, ApJ, 268, 319 [NASA ADS] [CrossRef] [Google Scholar]
 JiménezForteza, X., Keitel, D., Husa, S., et al. 2017, Phys. Rev. D, 95, 064024 [CrossRef] [Google Scholar]
 Lousto, C. O., Zlochower, Y., Dotti, M., & Volonteri, M. 2012, Phys. Rev. D, 85, 084015 [NASA ADS] [CrossRef] [Google Scholar]
 Maggiore, M. 2018, Gravitational Waves: Volume 2: Astrophysics and Cosmology (Oxford: Oxford University Press) [CrossRef] [Google Scholar]
 Maggiore, M., Van Den Broeck, C., Bartolo, N., et al. 2020, JCAP, 2020, 050 [CrossRef] [Google Scholar]
 Manwadkar, V., Trani, A. A., & Leigh, N. W. C. 2020, MNRAS, 497, 3694 [CrossRef] [Google Scholar]
 Mapelli, M. 2016, MNRAS, 459, 3432 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M., Giacobbo, N., Ripamonti, E., & Spera, M. 2017, MNRAS, 472, 2422 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M., Spera, M., Montanari, E., et al. 2020, ApJ, 888, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M., Dall’Amico, M., Bouffanais, Y., et al. 2021, MNRAS, 505, 339 [CrossRef] [Google Scholar]
 Mapelli, M., Bouffanais, Y., Santoliquido, F., Arca Sedda, M., & Artale, M. C. 2022, MNRAS, 511, 5797 [NASA ADS] [CrossRef] [Google Scholar]
 Memmesheimer, R.M., Gopakumar, A., & Schäfer, G. 2004, Phys. Rev. D, 70, 104011 [NASA ADS] [CrossRef] [Google Scholar]
 Meylan, G., & Heggie, D. C. 1997, A&ARv, 8, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Mikkola, S., & Aarseth, S. J. 1989, Celest. Mech. Dyn. Astron., 47, 375 [Google Scholar]
 Mikkola, S., & Aarseth, S. J. 1993, Celest. Mech. Dyn. Astron., 57, 439 [NASA ADS] [CrossRef] [Google Scholar]
 Mikkola, S., & Merritt, D. 2008, AJ, 135, 2398 [Google Scholar]
 Miller, M. C., & Hamilton, D. P. 2002, MNRAS, 330, 232 [Google Scholar]
 Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 2015, ApJ, 800, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Naoz, S. 2016, ARA&A, 54, 441 [Google Scholar]
 Neumayer, N., Seth, A., & Böker, T. 2020, A&ARv, 28, 4 [Google Scholar]
 Parischewsky, H. D., Trani, A. A., & Leigh, N. W. C. 2023, SciPost Phys. Core, 6, 016 [NASA ADS] [CrossRef] [Google Scholar]
 Peters, P. C. 1964, Phys. Rev., 136, 1224 [Google Scholar]
 Poincaré, H. 1892, Les méthodes nouvelles de la mécanique céleste (Paris: GauthierVillars et fils) [Google Scholar]
 Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17 [Google Scholar]
 Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., Boekholt, T. C. N., Por, E. H., Hamers, A. S., & McMillan, S. L. W. 2022, A&A, 659, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pryor, C., & Meylan, G. 1993, in Structure and Dynamics of Globular Clusters, eds. S. G. Djorgovski, & G. Meylan, ASP Conf. Ser., 50, 357 [Google Scholar]
 Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Class. Quant. Grav., 27, 194002 [Google Scholar]
 Quinlan, G. D. 1996, New Astron., 1, 35 [Google Scholar]
 Rastello, S., Mapelli, M., Di Carlo, U. N., et al. 2021, MNRAS, 507, 3612 [CrossRef] [Google Scholar]
 Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, Phys. Rev. Lett., 115, 051101 [Google Scholar]
 Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016a, Phys. Rev. D, 93, 084029 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., Haster, C.J., Chatterjee, S., Kalogera, V., & Rasio, F. A. 2016b, ApJ, 824, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., Zevin, M., Pankow, C., Kalogera, V., & Rasio, F. A. 2016c, ApJ, 832, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., AmaroSeoane, P., Chatterjee, S., et al. 2018, Phys. Rev. D, 98, 123005 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., Zevin, M., AmaroSeoane, P., et al. 2019, Phys. Rev. D, 100, 043027 [Google Scholar]
 RomeroShaw, I., Lasky, P. D., Thrane, E., & Calderón Bustillo, J. 2020, ApJ, 903, L5 [NASA ADS] [CrossRef] [Google Scholar]
 RomeroShaw, I., Lasky, P. D., & Thrane, E. 2021, ApJ, 921, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Samsing, J. 2018, Phys. Rev. D, 97, 103014 [CrossRef] [Google Scholar]
 Samsing, J., & Ilan, T. 2018, MNRAS, 476, 1548 [NASA ADS] [CrossRef] [Google Scholar]
 Samsing, J., & RamirezRuiz, E. 2017, ApJ, 840, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Samsing, J., MacLeod, M., & RamirezRuiz, E. 2014, ApJ, 784, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Samsing, J., MacLeod, M., & RamirezRuiz, E. 2017, ApJ, 846, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Samsing, J., MacLeod, M., & RamirezRuiz, E. 2018a, ApJ, 853, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Samsing, J., Askar, A., & Giersz, M. 2018b, ApJ, 855, 124 [NASA ADS] [CrossRef] [Google Scholar]
 Samsing, J., Bartos, I., D’Orazio, D. J., et al. 2022, Nature, 603, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Sigurdsson, S., & Phinney, E. S. 1993, ApJ, 415, 631 [NASA ADS] [CrossRef] [Google Scholar]
 Sigurdsson, S., & Phinney, E. S. 1995, ApJS, 99, 609 [NASA ADS] [CrossRef] [Google Scholar]
 Spera, M., Mapelli, M., & Jeffries, R. D. 2016, MNRAS, 460, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Spitzer, L. 1969, ApJ, 158, L139 [NASA ADS] [CrossRef] [Google Scholar]
 Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton: Princeton University Press) [Google Scholar]
 Tanikawa, A. 2013, MNRAS, 435, 1358 [NASA ADS] [CrossRef] [Google Scholar]
 Torniamenti, S., Rastello, S., Mapelli, M., et al. 2022, MNRAS, 517, 2953 [NASA ADS] [CrossRef] [Google Scholar]
 Trani, A. A., Spera, M., Leigh, N. W. C., & Fujii, M. S. 2019, ApJ, 885, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Trani, A. A., Tanikawa, A., Fujii, M. S., Leigh, N. W. C., & Kumamoto, J. 2021, MNRAS, 504, 910 [NASA ADS] [CrossRef] [Google Scholar]
 Trani, A. A., Rastello, S., Di Carlo, U. N., et al. 2022, MNRAS, 511, 1362 [NASA ADS] [Google Scholar]
 Trenti, M., & van der Marel, R. 2013, MNRAS, 435, 3272 [NASA ADS] [CrossRef] [Google Scholar]
 Valtonen, M., & Karttunen, H. 2006, The ThreeBody Problem (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zevin, M., Samsing, J., Rodriguez, C., Haster, C.J., & RamirezRuiz, E. 2019, ApJ, 871, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Zevin, M., RomeroShaw, I. M., Kremer, K., Thrane, E., & Lasky, P. D. 2021, ApJ, 921, L43 [NASA ADS] [CrossRef] [Google Scholar]
 Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703 [Google Scholar]
Appendix A: Impact of the maximum impact parameter
Threebody scatterings are computationally advantageous with respect to full star cluster simulations. The price to pay is that the results depend on the choice of initial parameters, which must convey information about the properties of the host star clusters and their binary systems. For example, to compute the upper limit of the impact parameter, b_{max} (eq. 12), we assumed that the distance of closest approach, r_{p}, is equal to the semimajor axis, a, of the binary system and that the strong gravitational focusing limit is valid, or in other words that . These assumptions ensure that the majority of our threebody interactions are classified as hard encounters, where significant energy exchange occurs between the single and binary BHs. How critical is this assumption, and what are the effects on our results if we change this prescription? If we relax this hypothesis, the upper limit of the distribution takes the more general form
The distance of closest approach is typically defined as r_{p} = k a, with k assumed as an arbitrary constant. To study how the strong gravitational focus approximation impacts our results, we ran 5 × 10^{4} additional threebody interactions initialised as the main NSC set, but with impact parameters generated according to eq. A.1, with k = 1, 2, 3, 4, 5.
Figure A.1 shows the cumulative distribution function of the impact parameter if we consider b_{max} computed as in eq. 12 or eq. A.1, with various choices for k in the NSC case. By calculating b_{max} from eq. 12, it is more likely to generate interactions where the impact parameter favours closer encounters compared to the sample generated with eq. A.1. The difference between the distributions is marginal with k = 1, but it gradually grows with increasing k. For example, there is a probability of approximately 34% for a simulation that has an impact parameter ≤10 AU in our fiducial model, while this probability decreases to 23% for k = 1 and drops to 5% for k = 5.
Fig. A.1. Cumulative distribution function of the impact parameter computed assuming b_{max} from eq. 12, and from eq. A.1 with k = 1, 2, 3, 4, 5. 
The probability of drawing an impact parameter, b > 10^{3} AU is < 0.1% in our fiducial model, while this fraction rises to 3% and 4%, respectively, in the cases of k = 4 and k = 5 with eq. A.1. Assuming large values of k produces a nonnegligible fraction of simulations where the impact parameter is of the same order of magnitude as the interparticle distance in the core of a NSC. In these simulations, the single BH is likely to interact with other members of the star cluster before reaching the target binary, and the simulated interaction may not occur.
Figure A.2 shows the outcomes of threebody interactions and BBH mergers as a function of the assumed b_{max}. The simulation sets with larger values of b_{max} result in fewer exchanges and ionisations while favouring more flybys, with this effect becoming more pronounced with larger values of k. For instance, in the cases of k = 4 and k = 5, almost eight to nine simulations out of ten evolve as a flyby. In contrast, our fiducial model produces flybys in almost 50% of cases. This happens because, when the strong gravitational focus approximation is not assumed, interactions with larger impact parameters become more common. In such encounters, the single BH is more likely to perceive the binary system as a pointlike object, leading more frequently to weaker energy exchanges between the binary and the single object. It also disfavours exchanges, which require closer passages between the single BH and one of the two binary components.
Fig. A.2. Same as Fig. 2 but assuming eq. A.1 for the choice of b_{max}. The plot shows the outcome of the threebody simulations (upper panel) and the formation channel of the BBH mergers (lower panel) for the main NSC set initialised with eq. 12 (dark red, fiducial model) and the smaller NSC sets initialised with eq. A.1 assuming k = 1 (orange), 2 (yellow), 3 (light green), 4 (dark green), and 5 (blue). 
The lower panel of Fig. A.2 shows that the choice of b_{max} also affects the percentage of mergers, even if less dramatically than the percentage of flybys, exchanges, and ionisations. The dependence of the percentage of BBH mergers on the choice of b_{max} is further shown in Fig. A.3. The fraction of BBH mergers decreases only by 1.7% from our fiducial model to the simulation set with k = 5. This plot also shows that the merger fraction linearly decreases with k as P_{merg} = α k + β, with α = −0.38 ± 0.09 and β = 11.83 ± 0.31.
Fig. A.3. Percentage of BBH mergers over all the simulated encounters in NSCs, where b_{max} was computed with eq. 12 and eq. A.1, assuming k = 1, 2, 3, 4, 5. The bar colours are the same as in Figs. A.1 and A.2. The dashed grey line is the linear fit to the percentage of mergers as a function of the k parameter. 
BBH mergers are less affected by the choice of b_{max} because both Eqs. 12 and A.1 encode the distribution of our semimajor axis, which shapes the overall impact parameter distribution through the b_{max} dependence on a. In our approach, we sample a in a way that ensures that all our binaries are hard, while also being sufficiently large to prevent efficient GW emission.
Different approaches are assumed in the literature to sample the parameter space of close encounters. For example, Hut & Bahcall (1983) and Zevin et al. (2019) computed b_{max} as eq. A.1 and set k = 1 with arbitrary fixed values for a. Samsing et al. (2014) chose k = 5 to ensure the sampling of distant interactions, but at the same time they kept the semimajor axis fixed at 10^{−5} AU, which corresponds to very hard binaries, close to the merger. Quinlan (1996), on the other hand, fixed a but varied the k parameter.
Here, we adopted a physically motivated semimajor axis distribution, ranging from the limit between hard and soft binaries down to the threshold for efficient GW decay. Moreover, in our models we aim to fully explore the parameter space for hard encounters. The price that we pay is that we do not consider encounters with large k. Their effect is shown in figs. A.2 and A.3, for completeness.
Finally, fig. A.4 shows the relative total energy variation of the binaries, dE/E = E_{fin}−E_{in}/E_{fin}, where E_{fin} and E_{in} are the total binary energy at the beginning and end of the simulation, respectively. The peak of the distribution shifts from larger to lower values of the relative total binary energy variation, when k increases. For instance, the distributions for our fiducial model and for eq. A.1 with k = 1 peak more than one order of magnitude above the distributions obtained for k = 4 and k = 5. In the former, only 2.6% and 2.7% of the binary systems yield dE/E < 10^{−2}, while this number increases to 20% and 31% for the latter.
Fig. A.4. Histograms of the normalised variation in the binary total energy with respect to the final binary total energy as a function of b_{max}. The colours are the same as in Figs. A.1, A.2, and A.3. 
All Tables
All Figures
Fig. 1. Trajectories of three eccentric mergers triggered by a threebody interaction in NSCs (left and righthand panels) and GCs (central panel). In the three panels, the systems are centred in the centre of mass of the initial BBH, while the position at which the merger takes place is marked with a white star. The trajectories of the primary and secondary BH of the initial BBH (with mass m_{1} and m_{2}), and the intruder (m_{3}) are reported in red, blue, and grey, respectively. The incoming direction of the intruder is shown with an arrow. The insets in the lefthand and central panels display a closeup view of the merger region. In the lefthand and central panels, the initial binary is shown edgeon, while in the righthand panel the view is faceon. The initial coalescence time of these BBHs, i.e. the coalescence time of the initial binary at the beginning of the interaction, is longer than the Hubble time. 

In the text 
Fig. 2. Upper panel: percentage of different interaction outcomes for each cluster type. From left to right: (i) flyby events, (ii) exchanges in which the intruder replaces the secondary BH, (iii) exchanges in which the intruder replaces the primary BH, (iv) ionisations, and (v) unstable triples. Lower panel: percentage of BBH mergers. From left to right: (i) BBH mergers occurring after a flyby, (ii) an exchange interaction with the secondary BH replaced by the intruder, (iii) an exchange with the primary BH replaced by the intruder, and (iv) secondgeneration BBH mergers. In both the upper and lower panels, the colours mark the cluster type in which the interaction takes place. 

In the text 
Fig. 3. Distribution of coalescence times for unperturbed initial binaries (filled histograms), and for the same initial binaries perturbed by the threebody interaction (unfilled histograms). The upper, central, and lower panels show the cases of YSCs, GCs, and NSCs, respectively. 

In the text 
Fig. 4. Total mass of the BBH mergers in YSCs (blue), GCs (orange), and NSCs (purple). 

In the text 
Fig. 5. Orbital eccentricity at 10 Hz (e_{10 Hz}) as a function of the coalescence time (t_{coal}) of BBH mergers in YSCs (left), GCs (centre), and NSCs (right). Different colours indicate mergers through flybys and the two families of exchanges. The yellow stars show secondgeneration BBH mergers. The empty black dots highlight the BBH progenitors of secondgeneration mergers. The dashed vertical line divides the plot into two regions: up to 1 Myr we evolve the systems with ARWV, above 1 Myr all the dynamical interactions are concluded, and with Eq. (1) we stop the ARWV integration and evolve the remaining BBHs. The three insets show a zoom for mergers with e_{10 Hz} > 0.1. We note that some of the points at e_{10 Hz} ∼ 1 overlap. 

In the text 
Fig. 6. Same as Fig. 5 but with the orbital eccentricity computed at 1 Hz (e_{1 Hz}). 

In the text 
Fig. 7. BBH mergers that take place after the threebody interactions in YSCs, GCs, and NSCs (from top to bottom). All the mergers that take place during the threebody interaction, i.e. < 10^{6} yr, are assumed to merge inside the cluster. The x axis reports the recoil velocity caused by the energy exchange in the threebody encounter. The y axis shows the coalescence time of the binary from the beginning of the simulation. Different outcomes are shown with different markers and colours. In the YSCs case, the horizontal dashed line reports the typical evaporation time of a YSC. The vertical dashed line in all the plots shows the escape velocity of the cluster. Distributions of the recoil velocities are displayed as marginal histograms. The colour legend is the same as the scatter plot. 

In the text 
Fig. 8. Distribution of the relativistic recoil kick of the remnants produced by the BBH mergers in YSCs, GCs, and NSCs. The dashed vertical lines represent the escape velocities of the three types of star clusters. 

In the text 
Fig. 9. Tilt angle distributions of the BBH orbital plane in NSCs at the end of the simulation with respect to the initial binary orbital plane. Unfilled histograms are flybys (red), exch13 (blue), and exch23 (grey). The filled histogram (yellow) shows the total distribution. We do not show GCs and YSCs because they have exactly the same behaviour. 

In the text 
Fig. A.1. Cumulative distribution function of the impact parameter computed assuming b_{max} from eq. 12, and from eq. A.1 with k = 1, 2, 3, 4, 5. 

In the text 
Fig. A.2. Same as Fig. 2 but assuming eq. A.1 for the choice of b_{max}. The plot shows the outcome of the threebody simulations (upper panel) and the formation channel of the BBH mergers (lower panel) for the main NSC set initialised with eq. 12 (dark red, fiducial model) and the smaller NSC sets initialised with eq. A.1 assuming k = 1 (orange), 2 (yellow), 3 (light green), 4 (dark green), and 5 (blue). 

In the text 
Fig. A.3. Percentage of BBH mergers over all the simulated encounters in NSCs, where b_{max} was computed with eq. 12 and eq. A.1, assuming k = 1, 2, 3, 4, 5. The bar colours are the same as in Figs. A.1 and A.2. The dashed grey line is the linear fit to the percentage of mergers as a function of the k parameter. 

In the text 
Fig. A.4. Histograms of the normalised variation in the binary total energy with respect to the final binary total energy as a function of b_{max}. The colours are the same as in Figs. A.1, A.2, and A.3. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.