Modeling Dense Star Clusters in the Milky Way and beyond with the Cluster Monte Carlo Code

In the Hénon method, the cluster is assumed to be spherically symmetric and in dynamical (virial) equilibrium. The full 6D positions and velocities are reduced to 3D: the radial position, r, of each particle, and its radial and tangential velocities, v_r and v_t . Of course, these position and velocities change on a dynamical timescale as the particle orbits through the cluster potential. What the Monte Carlo method tracks over time is the energy (E) and angular momentum (J) of each particle's orbit, quantities that are conserved on dynamical timescales. Combined with the cluster potential, Φ(r), these quantities fully define an orbit for every particle in the cluster. New r, v_r , and v_t are randomly sampled from these orbits at every Monte Carlo time step.

While E and J do not change on dynamical timescales, the cumulative effect of many weak encounters drives their evolution on a relaxation timescale. During a Monte Carlo time step, we model these weak encounters as perturbations ΔE and ΔJ to the particle's energy and angular momentum. These perturbations are computed as a single effective encounter between each particle and its nearest neighbor in radius. In the center-of-mass frame of the two-body encounter, the magnitude of the particles' velocities is unchanged, but the velocity vectors are deflected by some angle β. In this frame, the energy of the two particles is conserved; however, in the cluster frame, such encounters exchange energy and angular momentum. To correctly capture two-body relaxation throughout the cluster, each of these effective encounters should give the correct value for the mean change in kinetic energy at each particle's position during the time step ΔT. We satisfy this constraint by computing an effective deflection angle β_e for each particle that results in the correct mean ΔE and ΔJ. See Figure 1.

Zoom In Zoom Out Reset image size

2.1.1. Calculating the Effective Scattering Angle

We now proceed to derive the effective scattering angle, β_e , following the logic of Stodółkiewicz (1982) and Joshi et al. (2000). We begin with an array of particles sorted by increasing radius from the cluster center such that r_i < r_i+1. Fundamentally we are interested in the average change in energy and angular momentum, ΔE and ΔJ, particle i experiences during some time interval ΔT. We begin by considering the change in velocity squared, ${\left({\rm{\Delta }}{v}_{i}\right)}^{2}$, for a single two-body interaction between the star with index i and a field star with mass m_f . This expression can be written in terms of β, the angle of deflection of particle i, or in terms of b, the impact parameter of the encounter at infinity, as (e.g., Binney & Tremaine 2008, Equation (3.54))

$$\begin{eqnarray}&&{\left({\rm{\Delta }}{v}_{i}\right)}^{2}=\displaystyle \frac{4\,{m}_{f}^{2}}{{\left({m}_{i}+{m}_{f}\right)}^{2}}{w}^{2}{\sin }^{2}\left(\displaystyle \frac{\beta }{2}\right)\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{4\,{G}^{2}\,{m}_{f}^{2}}{{w}^{2}\,{b}_{90}^{2}}\displaystyle \frac{1}{1+{\left(b/{b}_{90}\right)}^{2}},\end{eqnarray} \tag{ 3 }$$

where w = ∣ v _i − v _f ∣ is the relative speed of the particles at infinity and b₉₀ ≡ G(m_i + m_f )/w² is the impact parameter that would deflect the particles by β = 90°.

It is straightforward to extend Equation (3) to describe the effect of many encounters. The change in velocity over time due to many weak encounters can be calculated by considering the field particles encountered by particle i as it travels through the cluster. The number of particles encountered at an impact parameter b during an interval ΔT is simply the product of the local number density of particles, n, and the volume swept out at that impact parameter. For an infinitesimal annulus with inner radius b and outer radius b+db, this product is N_enc = 2π bnwΔTdb, where the length of the cylinder is the relative velocity of the particles, w, times ΔT. To compute the contribution of this result over all impact parameters, we multiply Equation (3) by N_enc and integrate over b from zero to some upper limit ${b}_{\max }$:

$$\begin{eqnarray}&&\left\langle {\left({\rm{\Delta }}{v}_{i}\right)}^{2}\right\rangle =\displaystyle \frac{8\pi {G}^{2}n{\rm{\Delta }}T\,{m}_{f}^{2}}{{{wb}}_{90}^{2}}{\int }_{0}^{{b}_{\max }}\displaystyle \frac{b\cdot {db}}{1+{\left(b/{b}_{90}\right)}^{2}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&=\ 4\pi {G}^{2}n{\rm{\Delta }}T\,{m}_{f}^{2}{w}^{-1}\mathrm{ln}\left[1+{\left({b}_{\max }/{b}_{90}\right)}^{2}\right]\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\simeq \ 8\pi {G}^{2}n{\rm{\Delta }}T\,{m}_{f}^{2}{w}^{-1}\mathrm{ln}{\rm{\Lambda }}\,.\end{eqnarray} \tag{ 6 }$$

The value of ${b}_{\max }$ to use in Equation (4) is a well-studied complication to the theory: on the one hand, in a system with uniform density, equal octaves of b contribute equally to $\left\langle {\left({\rm{\Delta }}{v}_{i}\right)}^{2}\right\rangle$, meaning that for sufficiently large systems ${b}_{\max }\gg {b}_{90};$ on the other hand, star clusters are finite systems whose densities eventually fall to zero at some finite radii. In going from Equation (5) to Equation (6), we have introduced the well-known Coulomb logarithm, defined as

$$\begin{eqnarray}&&{\rm{\Lambda }}\equiv \displaystyle \frac{{b}_{\max }}{{b}_{90}}=\displaystyle \frac{{w}^{2}{b}_{\max }}{G({m}_{i}+{m}_{f})},\end{eqnarray} \tag{ 7 }$$

and used our assumption of the weak-encounter regime to set $\mathrm{ln}\left[1+{\left({b}_{\max }/{b}_{90}\right)}^{2}\right]\simeq 2\mathrm{ln}\left[{b}_{\max }/{b}_{90}\right]$. Λ can be approximately computed at the half-mass radius of the cluster by equating ${b}_{\max }={r}_{h}$ and replacing w² with the mean-squared velocity from the virial theorem, $\left\langle {v}^{2}\right\rangle \simeq 0.45{{Gm}}_{f}/{r}_{h}$ (Binney & Tremaine 2008, p. 361). For a cluster with N stars of mass m_f , this yields

$$\begin{eqnarray}&&\mathrm{ln}{\rm{\Lambda }}=\mathrm{ln}\left[\gamma N\right],\end{eqnarray} \tag{ 8 }$$

with γ ≃ 0.2. However, comparisons to direct N-body simulations have suggested a γ of 0.11 (for equal-mass clusters; e.g., Giersz & Heggie 1994), 0.02 (for clusters with realistic mass functions; Giersz & Heggie 1996), or even 0.01 (the current default in CMC; Freitag et al. 2006; Rodriguez et al. 2018c), to be more realistic.

Equations (4)–(6) were computed assuming that every field star has the same mass and relative velocity, m_f and w, with respect to our test star, m_i . In realistic clusters, however, nearby particles will have a range of masses and velocities, determined by the distribution function, F( r , v , m). Therefore, to accurately reproduce the cumulative shift $\left\langle {\left({\rm{\Delta }}{v}_{i}\right)}^{2}\right\rangle$ from Equation (6), we must replace the ${m}_{f}^{2}{w}^{-1}$ factor with

$$\begin{eqnarray}&&\left\langle {\left({\rm{\Delta }}{v}_{i}\right)}^{2}\right\rangle =8\pi {G}^{2}n{\rm{\Delta }}T\mathrm{ln}{\rm{\Lambda }}{\left\langle {m}_{f}^{2}{w}^{-1}\right\rangle }_{F},\end{eqnarray} \tag{ 9 }$$

where we use ${\left\langle ...\right\rangle }_{F}$ to indicate an average over the local phase-space distribution function

$$\begin{eqnarray}&&{\left\langle {m}_{f}^{2}{w}^{-1}\right\rangle }_{F}\equiv \int {F}_{i}{F}_{f}{m}_{f}^{2}{w}^{-1}{d}^{3}{{\boldsymbol{v}}}_{i}{d}^{3}{{\boldsymbol{v}}}_{f}{{dm}}_{i}{{dm}}_{f},\end{eqnarray} \tag{ 10 }$$

with F_i = F( r _i , v _i , m_i ) and F_f = F( r _f , v _f , m_f ).

We now reach the crux of the Monte Carlo method: to get the angle of deflection β in terms of these quantities—and the correct mean value of ${m}_{i}{\left({\rm{\Delta }}{v}_{i}\right)}^{2}$—we equate the right-hand side of Equation (9) to the right-hand side of Equation (2):

$$\begin{eqnarray}&&\begin{array}{l}{\left\langle \displaystyle \frac{4\,{m}_{f}^{2}\,{w}^{2}}{{\left({m}_{i}+{m}_{f}\right)}^{2}}\right\rangle }_{F}{\sin }^{2}\left(\displaystyle \frac{{\beta }_{e}}{2}\right)\\ \ \ \ =\ 8\pi {G}^{2}n{\rm{\Delta }}T\mathrm{ln}{\rm{\Lambda }}{\left\langle {m}_{f}^{2}{w}^{-1}\right\rangle }_{F},\end{array}\end{eqnarray} \tag{ 11 }$$

where we have replaced a single deflection angle, β, with an effective deflection angle, β_e , and replaced the mass and velocity terms with their averaged equivalents. We now solve for β_e , which involves averaging the mass and velocity quantities over the relevant distribution functions. However, as Hénon (1971a) pointed out, the most straightforward way to sample the distribution function near star m_i is to pick the mass and velocity of its nearest neighbor, m_i+1. While this is not the same as a true average, the nearest neighbor represents a fair draw from the mass and velocity distributions and would reproduce the true distribution function after a sufficient number of time steps. Dispensing with the distribution function averages and substituting m_f → m_i+1, Equation (11) becomes

$$\begin{eqnarray}&&{\sin }^{2}\left(\displaystyle \frac{{\beta }_{e}}{2}\right)=\displaystyle \frac{2\pi {G}^{2}{\left({m}_{i}+{m}_{i+1}\right)}^{2}}{{w}^{3}}n\mathrm{ln}{\rm{\Lambda }}{\rm{\Delta }}T.\end{eqnarray} \tag{ 12 }$$

The only quantity that must be explicitly averaged is the local volumetric density of particles, n. In CMC, this is done with a moving average over 41 stars, from i − 20 to i + 20 for each star i.

For computational convenience, Equation (12) is more easily expressed in terms of "relaxation timescale" between particles i and i + 1:

$$\begin{eqnarray}&&\sin \left(\displaystyle \frac{{\beta }_{e}}{2}\right)=\displaystyle \frac{\pi }{4}\sqrt{\displaystyle \frac{{\rm{\Delta }}T}{{T}_{\mathrm{rel},i}}},\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&{T}_{\mathrm{rel},i}\equiv \displaystyle \frac{\pi }{32}\displaystyle \frac{{w}^{3}}{{G}^{2}n{\left({m}_{i}+{m}_{i+1}\right)}^{2}\mathrm{ln}{\rm{\Lambda }}}\end{eqnarray} \tag{ 14 }$$

is the typical time required for a particle to be deflected by 90° (Freitag & Benz 2001). Note that, in practice, we clamp the maximum β_e between any two particles at 90°.

2.1.2. Calculating New E and J

With β_e computed for each pair of particles, we now convert these pairwise deflection angles into ΔE and ΔJ for each orbit. We follow Hénon's original approach, placing our two-body encounter into 3D space. We denote the phase-space coordinates of the two interacting particles by (r_i , v_r,i , v_t,i ) and (r_i+1, v_r,i+1, v_t,i+1), with masses m_i and m_i+1, respectively. We choose our reference frame such that the z-axis is parallel to r _i and the (x, z)-plane contains v _i . In Cartesian coordinates, the two particle velocities are then

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{v}}}_{i} & = & ({v}_{t,i},0,{v}_{r,i}),\\ {{\boldsymbol{v}}}_{i+1} & = & ({v}_{t,i+1}\cos \phi ,{v}_{t,i+1}\sin \phi ,{v}_{r,i+1}),\end{array}\end{eqnarray} \tag{ 15 }$$

where ϕ is a uniform random variate in the range [0, 2π], since spherical symmetry guarantees isotropic transverse velocities. The relative velocity is then

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{w}} & \equiv & {{\boldsymbol{v}}}_{i+1}-{{\boldsymbol{v}}}_{i}\\ & = & ({v}_{t,i+1}\cos \phi -{v}_{i,t},{v}_{t,i+1}\sin \phi ,{v}_{r,i+1}-{v}_{r,i}).\end{array}\end{eqnarray} \tag{ 16 }$$

We now define vectors w ₁ and w ₂ with magnitude equal to ∣ w ∣, such that w ₁ × w ₂ = w . In this right-handed coordinate system, we have

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{w}}}_{1} & \equiv & \left(\displaystyle \frac{-{w}_{y}w}{{w}_{p}},\displaystyle \frac{{w}_{x}w}{{w}_{p}},0\right),\\ {{\boldsymbol{w}}}_{2} & \equiv & \left(\displaystyle \frac{-{w}_{x}{w}_{z}}{{w}_{p}},\displaystyle \frac{-{w}_{y}{w}_{z}}{{w}_{p}},{w}_{p}\right),\end{array}\end{eqnarray} \tag{ 17 }$$

where ${w}_{p}=\sqrt{{w}_{x}^{2}+{w}_{y}^{2}}$. Since the distribution of field stars is assumed to be spherically symmetric, we randomly select the angle ψ ∈ [0, 2π] between the plane of relative motion, defined by ( r _i+1 − r _i , v _i+1 − v _i ), and the plane containing w and w ₁. The relative velocity after the encounter, w ^new, is then

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{w}}}^{\mathrm{new}} & = & {\boldsymbol{w}}\cos {\beta }_{e}+{{\boldsymbol{w}}}_{1}\sin {\beta }_{e}\cos \psi \\ & & +{{\boldsymbol{w}}}_{2}\sin {\beta }_{e}\sin \psi ,\end{array}\end{eqnarray} \tag{ 18 }$$

where β_e is the effective scattering angle of Section 2.1. The post-interaction particle velocities in the cluster frame, ${{\boldsymbol{v}}}_{i}^{\mathrm{new}}$ and ${{\boldsymbol{v}}}_{i+1}^{\mathrm{new}}$, become

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{v}}}_{i}^{\mathrm{new}} & = & {{\boldsymbol{v}}}_{i}-\left(\displaystyle \frac{{m}_{i+1}}{{m}_{i}+{m}_{i+1}}\right)({{\boldsymbol{w}}}^{\mathrm{new}}-{\boldsymbol{w}}),\\ {{\boldsymbol{v}}}_{i+1}^{\mathrm{new}} & = & {{\boldsymbol{v}}}_{i+1}+\left(\displaystyle \frac{{m}_{i}}{{m}_{i}+{m}_{i+1}}\right)({{\boldsymbol{w}}}^{\mathrm{new}}-{\boldsymbol{w}}).\end{array}\end{eqnarray} \tag{ 19 }$$

Since we set the z-axis parallel to r ⁱ , the new radial and transverse velocities for the first particle are ${v}_{r,i}^{\mathrm{new}}={v}_{z,i}^{\mathrm{new}}$ and ${v}_{t,i}^{\mathrm{new}}=\sqrt{{\left({v}_{x,i}^{\mathrm{new}}\right)}^{2}+{\left({v}_{y,i}^{\mathrm{new}}\right)}^{2}}$. The new specific orbital energy and specific angular momentum are then ${E}_{i}={\rm{\Phi }}({r}_{i})+\tfrac{1}{2}\left({\left({v}_{r,i}^{\mathrm{new}}\right)}^{2}+{\left({v}_{t,i}^{\mathrm{new}}\right)}^{2}\right)$ and ${J}_{i}={r}_{i}{v}_{t,i}^{\mathrm{new}}$, respectively, and similarly for the i + 1 particle.

Layered onto the basic two-body relaxation calculation is the possibility for strong interactions to occur between neighboring particles at each time step. Strong interactions can be separated into two groups: single–single interactions (which include physical collisions, tidal disruptions, tidal captures, and gravitational-wave captures of pairs of BHs) and binary interactions involving one or two binaries (integrated directly using fewbody). Strong interactions are sampled following the method of Freitag & Benz (2002; see also Fregeau & Rasio 2007), where we evaluate the quantity

$$\begin{eqnarray}&&{P}_{\mathrm{strong}}={nw}{{\rm{\Sigma }}}_{\mathrm{strong}}{\rm{\Delta }}T\end{eqnarray} \tag{ 20 }$$

between each pair of neighboring stars and binaries. Here P_strong is the probability for a strong interaction with cross section Σ_strong to occur, n is the local number density of stars, w is the relative velocity of the pair of objects at infinity (computed in the same coordinate system as Section 2.1.2), and ΔT is the time step. Accounting for gravitational focusing, Σ_strong can be written as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{strong}}=\pi {r}_{p}^{2}\left(1+\displaystyle \frac{2{GM}}{{r}_{p}{w}^{2}}\right),\end{eqnarray} \tag{ 21 }$$

where M is the total mass of the pair of interacting objects and r_p is the maximum distance of pericenter along a hyperbolic for each strong encounter to occur. At each time step, a random number X is drawn from a uniform distribution between 0 and 1. If X < P_strong, a strong interaction occurs; otherwise, the particles undergo standard two-body relaxation.

The value of r_p is determined by the type of interaction, as described in the following subsections. We note that both the strong encounters described here and the binary formation considered in Section 2.3 are actually performed before two-body relaxation is considered in CMC. Any particles that participate in strong encounters do not participate in two-body relaxation during the same time step. This is because, on a relaxation timescale, the change in a particle's orbit due to a strong encounter occurs instantaneously and is often significantly greater than the gradual changes due to two-body encounters. These large changes in velocity render the calculation in Section 2.1.1 invalid, since the relative velocity between field and test stars is no longer constant over a given time interval. However, because our time step is chosen as the minimum of the relaxation and strong encounter timescales (the latter of which is typically much smaller than the former), the deflection by two-body relaxation during a time step where a strong encounter occurs is typically minimal. See Section 2.8.

2.2.1. Single–Single Collisions

For collisions between single stars, CMC operates under the sticky-sphere approximation, where any two particles that touch radii are assumed to have collided. The cross section for a pair of single stars with masses m_i and m_i+1 and radii R_i and R_i+1 can be expressed as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{coll}}=\pi {\left({R}_{i}+{R}_{i+1}\right)}^{2}\left(1+\displaystyle \frac{2{GM}}{({R}_{i}+{R}_{i+1}){w}^{2}}\right),\end{eqnarray} \tag{ 22 }$$

where M = m_i + m_i+1.

For main-sequence and/or giant stars, we adopt the classic sticky-sphere approximation, where zero mass is lost during the collision. The mass of the collision product, m_new, is simply computed as m_i + m_i+1, with the new evolutionary phase of the star determined self-consistently using the stellar evolution prescriptions in COSMIC. In GCs, where the relative velocity of a pair of stars is typically much less than the escape velocity at the surface of the star, this approximation is reasonable (e.g., Lombardi et al. 2002). For compact objects, however, we employ more conservative assumptions for the mass retained by merger products. For collisions of white dwarfs (WDs) with other WDs, main-sequence stars, or giants, we follow the prescriptions described in Hurley et al. (2002, Section 2.7). As described in Kremer et al. (2019b), when a main-sequence star or giant collides with a BH/neutron star (NS), we assume that the accretion feedback supplies sufficient energy to unbind the stellar material from the system completely. In this case, negligible mass is accreted and the mass of the BH/NS is unchanged after the collision. We note that this is an extremely conservative assumption, particularly if the star is significantly more massive than the BH or NS. However, in practice such collisions are rare in most typical GC models; see Section 4.3.

2.2.2. Binary Interactions

2.2.3. Collisions during Binary Interactions

While star clusters are likely born with a nonnegligible fraction of their stars in binaries (e.g., Ivanova et al. 2005; Milone et al. 2012), they can also dynamically produce binaries many Myr after their initial phase of star formation is complete. Modeling this physics correctly is critical to understanding the long-term evolution of DSCs.

2.3.1. Three-body Binary Formation

Two-body encounters cannot create bound binaries from single stars (at least for Newtonian point masses) without some way to dissipate energy from the encounter. But if a two-body encounter occurs in the presence of a third body, excess kinetic energy can be carried away by the interloping particle, leaving behind a gravitationally bound system. While this process can involve more than three particles (e.g., Tanikawa et al. 2013), CMC uses a probabilistic treatment of three-body binary formation (Morscher et al. 2013) based on the formalism from Ivanova et al. (2005, 2010) and O'Leary et al. (2006).

Before either two-body relaxation or strong encounters are considered, the sorted list of stars is parsed three at a time, and the rate of binary formation for two objects of masses m₁ and m₂ from interactions with a third star (m₃) is calculated as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}_{3\mathrm{bb}}(\eta \geqslant {\eta }_{\min }) & = & \sqrt{2}{\pi }^{2}{G}^{5}{n}^{2}{\left\langle w\right\rangle }^{-9}\\ & & \times {\left({m}_{i}+{m}_{i+1}\right)}^{5}{\eta }_{\min }^{-5.5}(1+2{\eta }_{\min })\\ & & \times \left[1+2{\eta }_{\min }\left(\displaystyle \frac{{m}_{i}+{m}_{i+1}+{m}_{i+2}}{{m}_{i}+{m}_{i+1}}\right)\right],\end{array}\end{eqnarray} \tag{ 23 }$$

where n is the local stellar density of the cluster at that point, $\left\langle w\right\rangle =4\sigma /\sqrt{3\pi }$ is the average relative velocity at infinity for two particles in a Maxwellian distribution, and σ is the local 3D velocity dispersion. The rate is computed as a function of the binary hardness ratio, defined as

$$\begin{eqnarray}&&\eta \equiv \displaystyle \frac{{{Gm}}_{i}{m}_{i+1}}{a\left\langle m\right\rangle {\sigma }^{2}},\end{eqnarray} \tag{ 24 }$$

where a is the semimajor axis of the hypothetical binary, $\left\langle m\right\rangle$ is the average local mass, and σ is the average local velocity dispersion, both computed over the closest 20 particles to the encounter. Note that this is half the number of particles used for computing the averages in Section 2.1.1, since three-body binary interactions are a fundamentally local process. Since only hard binaries are expected to survive long enough to influence the long-term evolution of the cluster, we only form binaries with a minimum hardness of ${\eta }_{\min }=5$. By default, CMC only considers three-body binary formation for neighboring triplets of BHs in the cluster. Both of these choices are available as user parameters in CMC.

We note that Equation (23) differs from the original presentation in Ivanova et al. (2005, 2010) in two key ways. In Ivanova et al. (2005) the rate of three-body binary formation depends on the distinct number densities of the three masses. We simplify the expression here by replacing n₂, n₃, and n_c —the number densities of m₂, m₃, and the overall cluster core—with n, the local number density of the cluster. But for a given BH, we only consider Γ_3bb whenever both adjacent field stars are also BHs, meaning that in Monte Carlo fashion the probability calculated from Equation (23) is implicitly multiplied by ${\left({n}_{\mathrm{BH}}/n\right)}^{2}$ (or n₂ n₃/n² for particles with different masses); see also the discussion regarding binary–binary encounters in Section 2.2.2. In practice, the majority of three-body binary formation occurs in the BH-dominated region of the core (even when three-body binaries are allowed to form from any stars), meaning that n_c ≈ n₂ ≈ n₃ ≈ n_BH for most triple encounters.

The second key difference concerns the relative velocities: when computing Equation (23), we use $\left\langle w\right\rangle$, the averaged relative velocity between two particles. But the original expression from Ivanova et al. (2005) also contains the terms ${v}_{12}/\left\langle w\right\rangle$, where v₁₂ is the relative velocity between m₁ and m₂, and ${v}_{3}/\left\langle w\right\rangle$, where v₃ is the relative velocity between m₃ and the m₁/m₂ center of mass. For computing the probability of binary formation, we assume $\left\langle w\right\rangle ={v}_{12}={v}_{3}$ for computational convenience. However, when a binary is formed, we explicitly calculate v₁₂ and v₃ using the same coordinate system presented in Section 2.1.2. This is necessary to determine the relevant energies of the encounter and to sample the correct distribution of η for that binary. We draw a random η between ${\eta }_{\min }$ and 50 using the normalized distribution given by $\tfrac{d{{\rm{\Gamma }}}_{3\mathrm{bb}}}{d\eta }$, and we use this to determine the semimajor axis of the new binary (the eccentricities are assumed to be thermal). The relative velocities of the binary's center of mass and the third particle are then adjusted to conserve energy.

2.3.2. Two-body Binary Formation

While two-body binary formation is impossible when the two masses are Newtonian point particles, relaxing either of those assumptions can allow for binaries to form during close two-body encounters. If at least one of the masses has a nonnegligible radius (i.e., is not a WD, NS, or BH), close encounters can raise tides on one (or both) of the objects, allowing kinetic energy to be deposited into the structure of the stars themselves (e.g., Fabian et al. 1975). If this energy exceeds the overall kinetic energy of the encounter, the two stars can form a binary through tidal capture. We assume polytropic stellar models for the main-sequence stars undergoing tidal capture processes and use polytropic index n = 1.5 for low-mass main-sequence stars (<0.7 M_⊙) and n = 3 for higher-mass main-sequence stars. The calculations of the induced oscillation energy of the main-sequence stars and the cross sections during tidal capture follow prescriptions in Portegies Zwart & Meinen (1993) and Kim & Lee (1999), respectively. We track the first passage of a single–single close encounter. If a star is tidally captured after this passage, we set the binary semimajor axis to be twice the pericenter distance immediately after tidal capture, assuming angular momentum conservation. Note that giant stars are not included for tidal capture because of the highly uncertain reaction of their envelopes to the tidal force. In addition to tidal capture, if a compact object or a main-sequence star collides with a giant star during single–single interactions, the compact object/main-sequence star can form a binary with the giant core by transferring the orbital energy to the giant envelope and ejecting it. We treat main-sequence stars as point mass during this process and calculate the cross sections and final binary properties from the giant collisions using Equations (4)–(6) in Ivanova et al. (2006). See C. S. Ye et al. (2021, in preparation) for a complete description of these processes.

If instead of extended point particles we relax the Newtonian assumption, it becomes possible to form binaries through GW emission during close encounters (as occurs during BH mergers in fewbody). In CMC, the maximum pericenter distance where the GW emission of the two BHs carries away sufficient energy to take the two particles from unbound (positive total energy) to bound (negative total energy) is calculated as (Quinlan & Shapiro 1987)

$$\begin{eqnarray}&&{r}_{p}={\left(\displaystyle \frac{85\pi }{6}\right)}^{2/7}\displaystyle \frac{G({m}_{i}+{m}_{i+1})}{{c}^{10/7}{w}^{4/7}}.\end{eqnarray} \tag{ 25 }$$

BH binary formation via GW emission is decided probabilistically following a similar procedure to the collisions described in Section 2.3.1, using the above value of r_p in Equation (21). When a binary is formed, its impact parameter at infinity is sampled uniformly within the area of its capture cross section. The semimajor axis and eccentricity are calculated from conservation of energy and angular momentum, accounting for the loss due to GW emission. Many of these binaries are extremely wide (with semimajor axes of hundreds of AU or more), but their similarly extreme eccentricities (e ≳ 0.99) mean that they may still merge before they interact with neighboring particles. For these BBHs, we compare the timescale for GW emission to drive the binary to merge. We then compute the rate for this binary to encounter other stars using Equations (36) and (38) (albeit with the actual mass and semimajor axis of this binary, instead of the averaged quantities), and in true Monte Carlo fashion, we draw a random encounter time from an exponential distribution with that rate. If the encounter occurs before the BBHs would merge and the binary is pathologically wide (as described in Section 2.2.2), then the binary is disrupted. Otherwise, its semimajor axis and eccentricity are evolved by integrating the orbit-averaged equations for $\tfrac{{da}}{{dt}}$ and $\tfrac{{de}}{{dt}}$ from Peters (1964).

CMC evolves every star and binary in the cluster forward by time step ΔT using detailed prescriptions for stellar evolution. These prescriptions, originally from Hurley et al. (2000, 2002), have been significantly updated over many years and form the basis for the COSMIC population synthesis package, which is now directly integrated into CMC. See Section 3.1, Chatterjee et al. (2010), and Breivik et al. (2020a) for more details.

We perform stellar evolution in CMC after the dynamical encounters have taken place, but before the new cluster potential and orbits are computed. This is to ensure that any dynamical changes to individual stars and binaries (particularly collisions and mass transfer in post-encounter binaries) can be accounted for using the same stellar physics, but before the new cluster orbits are calculated. Collisions in fewbody, where multiple stars can potentially collide during a single encounter, are passed to the stellar evolution module in the same order that they occur after the encounter is complete; see also Section 2.2.3.

At the end of every Monte Carlo time step, a new position and corresponding velocity are chosen for each particle by randomly sampling a point on its orbit, weighted by the amount of time the particle spends at that point along its orbit. We do so by first computing the orbital pericenter ${r}_{\min }$ and apocenter ${r}_{\max }$ as the two roots of the energy equation:

$$\begin{eqnarray}&&Q(r)=2E-2{\rm{\Phi }}(r)-{J}^{2}/{r}^{2}=0.\end{eqnarray} \tag{ 26 }$$

The probability P(r) of finding the particle at a position r within interval dr is proportional to the time the particle spends traversing dr. Therefore,

$$\begin{eqnarray}&&P(r){dr}=\displaystyle \frac{{dt}}{T}=\displaystyle \frac{{dr}/| {v}_{r}| }{{\int }_{{r}_{\min }}^{{r}_{\max }}{dr}/| {v}_{r}| },\end{eqnarray} \tag{ 27 }$$

where T is half the orbital period and $| {v}_{r}| =\sqrt{Q(r)}$. We draw a random sample from P(r) and use the resultant r_new as the particle's position for the next time step. Note that P(r) ∝ 1/∣v_r ∣ becomes infinite at the turning points of the orbit, where ∣v_r ∣ = 0. Changing coordinates by defining a suitable function r = r(s) solves this problem (see Section 2.6 of Joshi et al. 2000). We must also set a new velocity consistent with the new orbital position: the new radial velocity is given magnitude $| {v}_{r}^{\mathrm{new}}| =\sqrt{Q({r}_{\mathrm{new}})}$ with a randomly selected sign, while the new transverse velocity is given by ${V}_{t}^{\mathrm{new}}=J/{r}_{\mathrm{new}}$. Once the new positions and velocities have been calculated, the particles are then resorted in order of radial distance.

After new positions and velocities have been sampled, we remove any particles that are no longer bound to the cluster. For clusters in isolation, this is accomplished by removing any particles with positive energy after the dynamics and sorting steps are complete. Such objects are routinely formed during strong encounters (which can eject particles from the central regions of the cluster) or through standard two-body relaxation, as many particles will slowly diffuse to positive energies.

CMC can also model clusters that are tidally limited, assuming the clusters to be on circular orbits within their galaxies. For a cluster of mass M_C on a circular orbit about the center of a point-mass galaxy, the tidal boundary of a cluster, r_t , can be calculated as

$$\begin{eqnarray}&&{r}_{t}^{3}={\left(\displaystyle \frac{{M}_{C}}{3{M}_{G}}\right)}^{1/3}{R}_{G},\end{eqnarray} \tag{ 28 }$$

where M_G is the mass of the galaxy and R_G is the distance from the cluster to that galaxy (e.g., Spitzer 1987). Although sometimes called the Jacobi radius in stellar dynamics, Equation (28) is identical to the calculation of the Roche surface in binary stellar physics. As the cluster loses mass, the initial tidal boundary (specified by the user) decreases at a rate $\propto \,{M}_{C}^{1/3}$, following Equation (28). ¹⁵

Of course, the Jacobi surface is not spherical; r_t actually represents the location of the Lagrange points along the line from the cluster to the galactic center. Since the Hénon method assumes the cluster to be spherically symmetric, approximate methods must be used to determine whether a star is stripped from the cluster by the galactic potential. CMC uses an energy-based criterion (Giersz et al. 2008) that removes any star with an energy above the potential at the tidal radius

$$\begin{eqnarray}&&E\gt \alpha {\rm{\Phi }}({r}_{t}),\end{eqnarray} \tag{ 29 }$$

where the parameter α is given by

$$\begin{eqnarray}&&\alpha =1.5-3{\left(\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{N}\right)}^{1/4}.\end{eqnarray} \tag{ 30 }$$

As previously, N is the initial number of stars in the cluster and $\mathrm{ln}{\rm{\Lambda }}$ is the Coulomb logarithm from Equation (8). Numerical testing (Giersz et al. 2008; Chatterjee et al. 2010) has shown that Equation (29) produces significantly better agreement with direct N-body simulations than simply stripping stars with apocenters beyond the tidal boundary.

Because we assume the cluster potential to be spherically symmetric, computing it is straightforward. The potential at point r, located between two particles at r_i and r_i+1, is given by

$$\begin{eqnarray}&&{\rm{\Phi }}(r)=-G\left(\displaystyle \frac{{M}_{i}}{r}+\sum _{i=j+1}^{N}\displaystyle \frac{{m}_{i}}{{r}_{i}}\right),\end{eqnarray} \tag{ 31 }$$

where ${M}_{i}\equiv {\sum }_{j=1}^{i}{m}_{j}$ is the total mass interior to star i.

While Equation (31) is formally correct, it is computationally expensive to evaluate whenever the potential at an arbitrary radius is required (such as when computing new orbital positions). However, we can accelerate this process by evaluating the potential at every particle, Φ_i , and interpolating the potential between that particle and its nearest neighbor, Φ_i+1. Because the stars are already sorted by radius, the potential at each particle can be computed recursively from the outermost star inward with

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{N+1}=0,\\ {{\rm{\Phi }}}_{i}={{\rm{\Phi }}}_{i+1}-{{GM}}_{i}\left(\displaystyle \frac{1}{{r}_{i}}-\displaystyle \frac{1}{{r}_{i+1}}\right),\end{array}\end{eqnarray} \tag{ 32 }$$

where ${M}_{N}={\sum }_{i=1}^{N}{m}_{i}$ is the total mass and M_i−1 = M_i − m_i . With Φ_i evaluated for every particle in the cluster, the potential at any radius r_i < r < r_i+1 can be expressed as

$$\begin{eqnarray}&&{\rm{\Phi }}(r)={{\rm{\Phi }}}_{i}+\left(\displaystyle \frac{1/{r}_{i}-1/r}{1/{r}_{i}-1/{r}_{i+1}}\right)\left({{\rm{\Phi }}}_{i+1}-{{\rm{\Phi }}}_{i}\right).\end{eqnarray} \tag{ 33 }$$

The potential at any r is then computed by finding the index i such that r_i ≤ r ≤ r_i+1 and applying Equation (32).

In CMC, the time step ΔT is computed as the minimum of several important physical timescales. The first such timescale is the local average relaxation timescale—the average time to deflect a particle by angle ${\beta }_{e}^{\max }$ within a certain localized region of the cluster. Following Equation (14), we define this as

$$\begin{eqnarray}&&{\widehat{T}}_{\mathrm{rel},i}={\left(\displaystyle \frac{{\beta }_{e}^{\max }}{\pi /2}\right)}^{2}\displaystyle \frac{\pi }{32}\displaystyle \frac{{\left\langle w\right\rangle }^{3}}{4{G}^{2}n\langle {m}^{2}\rangle \mathrm{ln}{\rm{\Lambda }}},\end{eqnarray} \tag{ 34 }$$

where ${\beta }_{e}^{\max }$ can be set by the user, with larger ${\beta }_{e}^{\max }$ corresponding to larger relaxation time steps. By default, CMC sets ${\beta }_{e}^{\max }=\sqrt{2}$ (though we typically use ${\beta }_{e}^{\max }=1$ for systems of point masses). Whereas T_rel,i is defined for pairs of particles i and i + 1, ${\widehat{T}}_{\mathrm{rel},i}$ is a broader rolling average over particles centered on particle i. The local number density n is calculated using the same rolling average as in Equation (12), with similar rolling averages used to compute $\left\langle {m}^{2}\right\rangle$ and the relative velocity between two particles, $\left\langle w\right\rangle$ (calculated from the averaged 3D velocity dispersion, assuming a locally Maxwellian velocity distribution). To accurately capture the overall relaxation dynamics, we require that ${\rm{\Delta }}T\lt {\widehat{T}}_{\mathrm{rel},i}$ for all i particles. In Freitag & Benz (2001, Equation (11)), this was accomplished by selecting a time step in each radial bin that was smaller than ${\widehat{T}}_{\mathrm{rel},i}$ by a factor of 0.005–0.05. In CMC, we instead choose a single relaxation time step for the entire cluster, taken to be the minimum ${\widehat{T}}_{\mathrm{rel},i}$ for all particles i. This choice ensures that our global time step is sufficiently short for all particles and to resolve the complicated dynamics in the cluster core, while making the parallelization of the code significantly easier; see Pattabiraman et al. (2013).

In addition to the relaxation timescale, we also compute characteristic timescales for collisions (T_coll), binary–binary encounters (T_bb), and binary–single encounters (T_bs). In general, the rate of collisions can be found by integrating the cross section for strong encounters over the distribution function of velocities for the two particles (Freitag & Benz 2002; Binney & Tremaine 2008):

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{strong}}=\displaystyle \frac{1}{n}\int F({{\boldsymbol{v}}}_{i})F({{\boldsymbol{v}}}_{i+1})w{{\rm{\Sigma }}}_{\mathrm{strong}}{d}^{3}{{\boldsymbol{v}}}_{i}{d}^{3}{{\boldsymbol{v}}}_{i+1},\end{eqnarray} \tag{ 35 }$$

where, as before, F is the distribution function (ignoring the mass and position), $w=\left|{{\boldsymbol{v}}}_{i+1}-{{\boldsymbol{v}}}_{i}\right|$ is the relative velocity, and Σ_strong is taken from Equation (21). Assuming a Maxwellian distribution for the velocity distribution, Equation (35) can be evaluated as (e.g., Binney & Tremaine 2008, p. 626)

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{strong}}=4\sqrt{\pi }{{nr}}_{p}^{2}\sigma \left(1+\displaystyle \frac{{GM}}{2{r}_{p}{\sigma }^{2}}\right),\end{eqnarray} \tag{ 36 }$$

and the relevant timescale between encounters is simply ${T}_{\mathrm{strong}}={{\rm{\Gamma }}}_{\mathrm{strong}}^{-1}$. The timescales of relevance to CMC are then (see Fregeau & Rasio 2007, for details)

$$\begin{eqnarray}&&{T}_{\mathrm{coll}}^{-1}=16\sqrt{\pi }{n}_{s}\left\langle {R}^{2}\right\rangle \sigma \left(1+\displaystyle \frac{G\left\langle {MR}\right\rangle }{2{\sigma }^{2}\left\langle {R}^{2}\right\rangle }\right),\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{bs}}^{-1}=4\sqrt{\pi }{n}_{s}{X}_{\mathrm{bin}}^{2}\langle {a}^{2}\rangle \sigma \left(1+\displaystyle \frac{G\langle M\rangle \langle a\rangle }{2{\sigma }^{2}{X}_{\mathrm{bin}}\langle {a}^{2}\rangle }\right),\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{bb}}^{-1}=16\sqrt{\pi }{n}_{\mathrm{bin}}{X}_{\mathrm{bin}}^{2}\left\langle {a}^{2}\right\rangle \sigma \left(1+\displaystyle \frac{G\left\langle {Ma}\right\rangle }{2{\sigma }^{2}{X}_{\mathrm{bin}}\left\langle {a}^{2}\right\rangle }\right),\end{eqnarray} \tag{ 39 }$$

where all of the averaged quantities above, including the number densities of single and binary stars, n_s and n_bin, and the velocity dispersion σ, are computed over the innermost 300 particles in the cluster (i.e., the central regions of the core). Note that if no binaries are present within these 300 stars, T_bs and T_bb are set to infinity for that time step. Conversely, if only binaries are present, then T_bb is computed normally, while T_bs is set to infinity.

Finally, while the stellar evolution package has its own internal time step, the mass loss from stars can drive significant dynamical changes in the cluster, particularly in the early stages of evolution. We compute a relevant timescale for stellar evolution based on the total fraction of mass lost in the previous time step:

$$\begin{eqnarray}&&{T}_{\mathrm{se}}=0.001\left(\displaystyle \frac{{M}_{\mathrm{cl}}}{{\rm{\Delta }}{M}_{\mathrm{se}}}\right){\rm{\Delta }}{T}_{\mathrm{prev}},\end{eqnarray} \tag{ 40 }$$

where M_cl is the total cluster mass, ΔT_prev is the previous time step, and ΔM_se is the mass lost during the current time step. CMC uses the minimum individual time step from Equations (34)–(40) as the global time step for the cluster. However, we note that for a typical cluster model (such as those presented in Section 4.2), the central relaxation timescale is the smallest of the relevant timescales for more than 96% of the time steps within the first 100 Myr.

Finally, to ensure energy conservation, a correction must be made to the energies of the stellar orbits. While the two-body relaxation of the Hénon method intrinsically conserves energy, the process of calculating new orbits and sampling positions and velocities is done after dynamical encounters have altered the energy and angular momenta of each orbit, but before the potential has been updated to account for the changes in the orbits. In other words, the Φ in Equation (31) is perpetually one time step behind the new E and J, allowing this time-dependent shift in the potential to do work on the particle. While the energy error during a single time step is negligible, this effect can become significant during long-term integrations, leading to a spurious drift in the energy of the cluster.

To counter this, we adopt the energy conservation scheme first developed by Stodółkiewicz (1982) in his update of the method. First, note that the change in energy for a single particle is expressible as

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{i}^{\mathrm{corr}}=\int \displaystyle \frac{\partial U({r}_{i})}{\partial t}{dt}.\end{eqnarray} \tag{ 41 }$$

Since we are interested in the work done by the change in potential between the previous and current time steps, Equation (41) can be approximated as the average of the change in potential energy between the positions of the particle at the previous and current time steps:

$$\begin{eqnarray}&&\int \displaystyle \frac{\partial U({r}_{i})}{\partial t}{dt}=\left[{\rm{\Delta }}{\rm{\Phi }}({r}_{i}^{\mathrm{prev}})+{\rm{\Delta }}{\rm{\Phi }}({r}_{i}^{\mathrm{curr}})\right]/2,\end{eqnarray} \tag{ 42 }$$

where ΔΦ ≡ Φ^curr − Φ^prev. To implement this conservation scheme, CMC records the four relevant energies for each particle every time step—${{\rm{\Phi }}}^{\mathrm{prev}}({r}_{i}^{\mathrm{prev}})$, ${{\rm{\Phi }}}^{\mathrm{prev}}({r}_{i}^{\mathrm{curr}})$, ${{\rm{\Phi }}}^{\mathrm{curr}}({r}_{i}^{\mathrm{prev}})$, and ${{\rm{\Phi }}}^{\mathrm{curr}}({r}_{i}^{\mathrm{curr}})$—and uses them to compute ${\rm{\Delta }}{E}_{i}^{\mathrm{corr}}$ for each particle. The velocity magnitudes for each particle are then scaled to increase or decrease the particle's kinetic energy by ${\rm{\Delta }}{E}_{i}^{\mathrm{corr}}$ (while maintaining the ratio of transverse to radial velocities as determined by the Hénon method), ensuring that the cluster conserves energy over many relaxation times.

In addition to the energy drift, a second issue arises from calculating stellar orbits for the current time step using the potential from the previous time step. For any given particle, the potential includes the contribution of that same particle at its last position in the cluster, and when the new orbit is calculated, it includes this spurious "self-force" from the previous time step. While this effect is relatively minor for clusters with equal-mass particles, it can be significant in clusters with a realistic mass spectrum, particularly when the most massive particles (e.g., BHs) are concentrated in the cluster center, where only a handful of particles determine the potential (Freitag & Benz 2001; Fregeau & Rasio 2007).

To solve this, we add a correction term to remove the particle's own contribution to the previous potential at its previous position:

$$\begin{eqnarray}{{\rm{\Phi }}}_{\mathrm{sf}}(r)=\left\{\begin{array}{cl}{{GM}}_{i}/r & r\geqslant {r}_{i}^{\mathrm{prev}},\\ {{GM}}_{i}/{r}_{i}^{\mathrm{prev}} & r\lt {r}_{i}^{\mathrm{prev}}.\end{array}\right.\end{eqnarray} \tag{ 43 }$$

This self-force correction, Φ_sf(r), is added to the energy Equation (26) when the new radial positions are sampled from the particle orbits (Section 2.5).

The Monte Carlo method, and Hénon's method in particular, has been used for 50 yr for the study of DSCs, starting with the initial papers by Hénon (1971a, 1971b). The technique was reintroduced with improvements by Stodółkiewicz (1982), and it is this version that most modern implementations (Giersz 1998; Joshi et al. 2000; Freitag & Benz 2001; Sollima & Mastrobuono Battisti 2014) are based on. While Monte Carlo methods have enjoyed great success in the modeling of DSCs, that success relies on a series of assumptions that we have implicitly or explicitly made in Section 2. We now briefly describe each of these assumptions and the limitations they impose on the method.

Underlying nearly all of the formal averages and encounter rates presented in Section 2 is Boltzmann's molecular chaos assumption, where we have assumed that the distribution functions of our test stars and field stars are independent and separable. This assumption is explicitly used in Equations (10) and (35), where we assume that the joint distribution function for particles i and i + 1 can be written as

$$\begin{eqnarray}&&F({{\boldsymbol{v}}}_{i},{{\boldsymbol{v}}}_{i+1})=F({{\boldsymbol{v}}}_{i})F({{\boldsymbol{v}}}_{i+1}).\end{eqnarray} \tag{ 44 }$$

Equation (44) can be justified in the following way: since equal octaves of impact parameter contribute equally to Equation (4), the majority of two-body encounters will occur at ${b}_{90}\ll b\ll {b}_{\max }$. Because of this, it is highly unlikely for two particles in a sufficiently large cluster to encounter each other more than once within their lifetimes, or at least before other encounters have "reset" the velocities of the two particles. While this assumption is valid for large-N systems, it does mean that the Monte Carlo method cannot resolve resonant effects between particles, such as scalar resonant relaxation (Rauch & Tremaine 1996) around central-massive BHs or resonant behavior between heavy objects in the centers of GCs (Meiron & Kocsis 2019).

The assumption of large N is also employed when computing the two-body diffusion of ${\left({\rm{\Delta }}v\right)}^{2}$. When computing Equation (6), we assumed ${\rm{\Lambda }}={b}_{\max }/{b}_{90}\gg 1$, which allowed us to substitute $\mathrm{ln}\left(1+{{\rm{\Lambda }}}^{2}\right)\simeq 2\mathrm{ln}({\rm{\Lambda }})$ and neglect any higher-order contributions to the change in velocity, e.g., $\left\langle {\left({\rm{\Delta }}{v}_{i}\right)}^{4}\right\rangle$. This neglect of the higher-order terms is known as the Fokker–Planck approximation and is largely justified because any terms beyond second order are at least a factor of $\mathrm{ln}{\left({\rm{\Lambda }}\right)}^{-1}$ smaller than the $\left\langle {\left({\rm{\Delta }}v\right)}^{2}\right\rangle$ term (Henon 1973; Binney & Tremaine 2008). This is also related to our assumption of orbit averaging, where we assumed that N is sufficiently large that T_relax ≫ T_dyn in Equation (1), allowing us to treat the orbits as essentially fixed on dynamical timescales.

Naturally, one may ask how large must N be for the above assumptions to be valid. While hard to derive explicitly,

$$\begin{eqnarray}&&N\gtrsim {10}^{3}\displaystyle \frac{{m}_{\max }}{\left\langle m\right\rangle }\end{eqnarray} \tag{ 45 }$$

has been shown to be a decent working criterion (Freitag 2008). Note that this criterion depends on the ratio between the largest mass in the cluster, ${m}_{\max }$, and the average stellar mass, $\left\langle m\right\rangle$. For clusters with realistic initial mass functions (IMFs), the ratio can be as large as 100:1 (once stellar evolution has reduced the mass of the most massive stars). In our experience, N must be ≳10⁵ initially for such clusters to show acceptable energy conservation and good agreement with direct N-body results (Chatterjee et al. 2010).

Hénon's method also assumes that the cluster is spherically symmetric. This assumption is not explicitly related to the size of the cluster or the relaxation timescale, but it is necessary for computing new particle orbits, determining nearest neighbors, and computing the potential. While this is an appropriate assumption for many GCs and some NSCs, this does mean that clusters with net angular momentum (e.g., rotating GCs, if the rotation creates a significant departure from spherical symmetry) or triaxial stellar distributions (e.g., NSCs that have recently undergone a merger) cannot be modeled with CMC. To address this, similar techniques, referred to as orbit-following Monte Carlo methods, have been developed that explicitly follow the orbits of stars in arbitrary potentials on a dynamical timescale, either by directly integrating the orbits in a given potential (the so-called Princeton Monte Carlo, of Spitzer & Hart 1971; Spitzer & Thuan 1972) or by tracking the diffusion of orbits over an integer number of orbital periods (the so-called Cornell Monte Carlo of Marchant & Shapiro 1980). These techniques have been used extensively for the study of triaxial NSCs and central-massive BH dynamics (see, e.g., Hopman 2009; Madigan et al. 2011; Vasiliev 2015, for recent examples). We note that the assumption of spherical symmetry does not mean that we have assumed the velocities to be isotropic. CMC can model systems with significant velocity anisotropy, since the ratio of v_r to v_t is calculated self-consistently every time step.

Finally, because the Hénon method operates on a relaxation timescale, we are implicitly assuming that GCs can be modeled as quasi-adiabatic transitions between different cluster models in dynamical (virial) equilibrium. While this assumption is largely valid for the long-term evolution of clusters in most cases, there are notable exceptions, such as tidal shocking (Ostriker et al. 1972), or the initial "violent relaxation" of the cluster after its (potentially highly nonvirialized) birth (Lynden-Bell 1967). These processes involve clusters that are sometimes far from equilibrium and evolve to a virialized state on a dynamical timescale. For these cases, more computationally expensive orbit-following techniques, such as the aforementioned Princeton/Cornell methods or direct N-body integrations, must be employed.

While COSMIC was originally based on the version of BSE incorporated into CMC, including its updates to compact object physics (Chatterjee et al. 2010; Rodriguez et al. 2016c), the two codes have diverged over the past several years, and the use of COSMIC as a community-driven population synthesis code has kept it up to date with the latest developments in binary stellar evolution. Here we review the changes to binary stellar evolution (and in particular binary mass transfer) in the latest version of COSMIC (v3.4), which is paired with the release of CMC. For a complete description of the physics included in COSMIC and the various options, see Hurley et al. (2000, 2002) and Breivik et al. (2020a).

Since the release of COSMIC v3.3, the treatment of Roche overflow mass transfer from nondegenerate stars has been expanded to include multiple assumptions for both mass-loss rates from the donor and mass accretion rates onto the accretor. The default mass-loss rate is determined using the expression given by Hurley et al. (2002), which steeply increases with the ratio of the donor radius to its Roche lobe radius

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{don}}=F({M}_{\mathrm{don}}){\left[\mathrm{ln}\displaystyle \frac{{R}_{\mathrm{don}}}{{R}_{\mathrm{rl},\mathrm{don}}}\right]}^{3}\,{M}_{\odot }\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 46 }$$

where ${\dot{M}}_{\mathrm{don}}$ is the donor mass-loss rate, M_don is the donor mass in units of M_⊙, R_don is the donor radius, R_rl,don is the donor's Roche lobe radius, and

$$\begin{eqnarray}&&F({M}_{\mathrm{don}})=3\times {10}^{-6}\,{[\min \left({M}_{\mathrm{don}},5.0\right)]}^{2}.\end{eqnarray} \tag{ 47 }$$

Version 3.4 of COSMIC also includes the Roche lobe overflow prescription of Claeys et al. (2014), which is calibrated to binary systems with zero-age main-sequence component masses less than 10 M_⊙. In this case, the mass-loss rate of the donor is determined similarly to Equation (46), but with an overall multiplicative factor, f, given by

$$\begin{eqnarray}f=\left\{\begin{array}{ll}1000, & Q\lt 1\\ \max \left(1,\displaystyle \frac{1000}{Q}\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{\mathrm{ln}Q}{0.15}\right)}^{2}\right]\right) & Q\geqslant 1,\end{array}\right.\end{eqnarray} \tag{ 48 }$$

where Q ≡ M_acc/M_don. In both cases, mass loss is always limited by the thermal timescale such that

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{don}}^{\max }=\max \left({\dot{M}}_{\mathrm{don}},\displaystyle \frac{{M}_{\mathrm{don}}}{{\tau }_{\mathrm{KH}}}{M}_{\odot }{\mathrm{yr}}^{-1}\right),\end{eqnarray} \tag{ 49 }$$

where

$$\begin{eqnarray}&&{\tau }_{\mathrm{KH}}=3.1\times {10}^{7}{\left(\displaystyle \frac{M}{{M}_{\odot }}\right)}^{2}\left(\displaystyle \frac{{R}_{\odot }}{R}\right)\left(\displaystyle \frac{{L}_{\odot }}{L}\right)\,\mathrm{yr}\end{eqnarray} \tag{ 50 }$$

is the Kelvin–Helmholtz time.

Four prescriptions are now included to determine the amount of mass that the accretor is allowed to gain, in addition to the original one in Hurley et al. (2002). The original BSE prescription assumes that stars with radiative envelopes (i.e., main sequence, Hertzsprung gap, and core helium burning) can accrete mass at a rate limited to 10 times the thermal rate of the star, while stars with convective envelopes (i.e., giant branch, asymptotic giant branch) can accrete mass at an unlimited rate. Three of the newly added prescriptions consider modifications to this prescription: one that limits stars with radiative envelopes to accrete mass at the thermal rate, one that applies accretion limits of 10 times the thermal rate for all stars—not just those with radiative envelopes—and one that applies accretion limits of the thermal rate for all stars. We also include another choice where accretion is a fixed fraction of the mass lost by the donor, as used in Belczynski et al. (2008).

In COSMIC, initial conditions are generated by sampling particle positions and velocities from either a Plummer (1911), Elson et al. (1987), or King (1966) profile. Both the Elson and Plummer profiles are sampled piecewise, with the particle positions drawn first from the cumulative mass distribution for each given profile. To self-consistently sample velocities for each particle at its given position, we draw velocities from the distribution function via rejection sampling up to the local escape speed (Aarseth 2003). For the Plummer profile we use the analytic form of the distribution function, while for the Elson profile we numerically solve for the distribution function in energy space (see, e.g., Grudić et al. 2018, Appendix B). ¹⁶ The King profile is generated in a similar fashion: the differential equations for cluster density, potential, and enclosed mass are solved numerically, with the latter being used to sample random positions from the cluster center. The velocities are then sampled from the analytic distribution function via the aforementioned rejection sampling technique.

By default, COSMIC draws radii and velocities from the available distributions before converting them into Hénon units (Heggie & Mathieu 1986; Heggie 2014), where the gravitational constant and initial cluster mass are both set to unity (G = M₀ = 1) and the sum of the initial cluster kinetic and potential energies is E_0,kin + E_0,pot = −1/4. (This is not the same as the total initial cluster energy E₀ = E_0,kin + E_0,pot +E_0,bin, where E_0,bin is the initial negative-valued total binding energy in primordial binaries, if any.) The units of mass, length, and time are then

$$\begin{eqnarray}&&\begin{array}{l}{U}_{m}={M}_{0},\\ {U}_{l}=\displaystyle \frac{-{{GM}}_{0}^{2}}{4{E}_{0}},\\ {U}_{t}=\displaystyle \frac{{{GM}}_{0}^{5/2}}{{\left(-4{E}_{0}\right)}^{3/2}},\end{array}\end{eqnarray} \tag{ 51 }$$

respectively. The above length units are such that the initial virial radius of the cluster is unity, while the units of time roughly correspond to the crossing time of a particle at the virial radius. However, because CMC operates on a relaxation timescale, we instead use a more appropriate unit of time,

$$\begin{eqnarray}&&{U}_{t}^{\mathrm{rel}}\equiv \displaystyle \frac{{N}_{0}}{\mathrm{log}\gamma {N}_{0}}{U}_{t},\end{eqnarray} \tag{ 52 }$$

where N₀ is the initial number of particles (single stars or binaries). These units are the default used throughout CMC.

Once the positions and velocities have been selected, the user can optionally add additional stellar physics to each of the particles, such as stellar masses, radii, and binary companions. All of the options available in COSMIC, including IMFs (e.g., Kroupa 2001; Salpeter 1955), binary initial conditions (e.g., following Sana et al. 2012; Moe & Di Stefano 2017), and stellar metallicities, can be used to generate realistic cluster initial conditions for CMC. When generating binary initial conditions for cluster simulations, we truncate the upper limit of the P_orb distribution at the hard−soft boundary for each binary at its respective radius in the cluster (see, e.g., Heggie 1975).

Input in CMC is handled at the command-line level by passing an initialization file containing all the parameters controlling the physics, parallelization, and runtime options (unspecified parameters are set to default values). The parameter files contain many of the same options (particularly related to stellar evolution) as the parameter files for COSMIC, allowing initial conditions, population synthesis, and cluster dynamics simulations to be performed with identical options. The initial conditions are saved as tables in the Hierarchical Data Format (HDF5) and passed as input to the initial parameters file.

Three forms of output are generated during a typical CMC run. First, a series of tab-separated data files containing specific parameters of the cluster (e.g., mass, half-mass radius, number of binaries) are printed every time step and can be easily read in any scripting language. Second, specific events, such as stellar collisions, BH mergers, binary–single encounters, etc., are recorded as they occur in human-readable log files (with a provided Python parser). Third, snapshots of the full cluster state are periodically written to several HDF5 tables at user-specified intervals. Separate files exist for snapshots containing all stars in the cluster and those containing only the BHs (allowing the latter to be resolved at shorter intervals). Each snapshot is saved in its respective file as a pandas-readable HDF5 table with keys corresponding to the simulation time when the snapshot was written. These can be directly imported into Python using the read_hdf command in pandas (McKinney 2010; Reback et al. 2021).

Finally, CMC periodically saves checkpoints by directly dumping the state of every processor to a binary file associated with that processor. These checkpoints allow bit-for-bit restarting of a CMC run should the run be either interrupted or terminated owing to computing cluster queue-time limits. While the checkpoints contain the state of the cluster, stars, and binaries up to that point, both the physical parameters of the run (set in the initialization file) and the current random number seed can be changed by the user at restart.

CMC is designed to be compiled against the standard libraries of the Message Passing Interface (MPI; Clarke et al. 1994), enabling distributed-memory parallelization across multiple cores and machines. Many of the physical processes in CMC, such as stellar evolution and nearest-neighbor encounters, are fundamentally local processes that are amenable to simple parallelization. However, the radial sorting of stars, the computation of new stellar orbits in the global cluster potential (Section 2.5), and even the computation of the potential itself (Section 2.7) are global processes that require a customized parallel approach.

Because the calculation of the cluster potential in Equation (33) and the nearest-neighbor dynamics require the particles to be sorted radially, the particles are resorted each time step after new positions are selected (Section 2.5) in order of increasing radius from the cluster center. CMC uses a custom parallel sorting implementation, in which the local data on each processor are sorted using a Quicksort algorithm (Hoare 1961) before being combined globally across processors using a parallel Sample Sort algorithm (Frazer & McKellar 1970). This implementation allows for rapid sorting and redistribution of particles across multiple distributed-memory processes. For more details on the parallel sorting method, see Section 3.4 of Pattabiraman et al. (2013). Once the sorting is complete, the particles are redistributed across all MPI processes with a series of all-to-all communications.

The calculation of the stellar orbits is more straightforward and only requires that every MPI process have access to the global potential of the cluster at each time step. To that end, each process maintains a global list of the radial positions and masses for every star in the cluster, allowing each process to compute Equation (33) independently. These lists are updated every time step. For a complete description of this and the other parallelization strategies employed in CMC, see Pattabiraman et al. (2013).

The Plummer profile was one of the first used to fit realistic observations of GCs to a theoretical model, and it remains popular since, unlike models based on more sophisticated distribution functions, most of its key features can be expressed analytically (Heggie & Hut 2003). This also makes it a favorite tool of the N-body simulator since the enclosed mass, velocity dispersion, and distribution function can all be expressed analytically and are easily sampled to generate cluster initial conditions. Here, we use COSMIC to generate a realization of a Plummer sphere with 10⁸ initial point-mass particles. We choose this particular setup because the evolution of a Plummer sphere to core collapse is a well-studied problem in the literature (e.g., Freitag & Benz 2001; Freitag et al. 2006; Binney & Tremaine 2008). We then integrate these initial conditions twice using CMC: once with no binary physics, and once with three-body binary formation and strong binary encounters enabled.

In Figure 2, we show the evolution of the two clusters from their initial state to core collapse, defined as having < 1000 stars within the theoretical core radius (see Casertano & Hut 1985, for the description of the latter). We also show the number of binaries as a function of time, since it is the formation and hardening of binaries in the core that provide an energy source to support the cluster against continued collapse. Since we stop both models at core collapse, we do not show the post-collapse evolution of either cluster; however, the cluster model with binary physics has its core collapse time delayed by 0.15 in units of the half-mass relaxation time:

$$\begin{eqnarray}&&{T}_{\mathrm{rlx}}=0.138\displaystyle \frac{N}{\mathrm{log}(\gamma N)}{\left(\displaystyle \frac{{r}_{h}^{3}}{{GM}}\right)}^{1/2},\end{eqnarray} \tag{ 53 }$$

where we have used the more precise definition from Spitzer (1987), rather than the approximate Equation (1). Similarly, the number of binaries in the core increases dramatically during the final moments of core collapse, as the central density increases sharply. In both cases the cluster collapses at about 16.8T_rlx, consistent with both theoretical predictions and previous numerical integrations (Cohn 1979, 1980; Spurzem & Aarseth 1996; Quinlan 1996; Freitag & Benz 2001).

Zoom In Zoom Out Reset image size

It is also well known that the collapse of a cluster proceeds homologously, with the central density and time to core collapse evolving at a self-similar rate throughout the cluster's central regions (Hénon 1961; Lynden-Bell 1975). Both homologous models and 1D Fokker–Planck models suggest that the central density profile of a cluster near collapse should approach a power law of the form ρ ∝ r^−2.23 (Cohn 1980; Lynden-Bell & Eggleton 1980; Heggie & Stevenson 1988). In Figure 3, we show the density of stars (averaged over 20 nearest neighbors) in our cluster model without binary formation as it approaches core collapse. The density of the central regions increases dramatically as the core approaches final collapse, with the central density increasing by nearly 4 orders of magnitude within the last 10⁻⁵ T_rlx of the evolution. At the moment of core collapse, CMC reproduces the r^−2.23 distribution over more than 15 orders of magnitude. To the best of our knowledge, this is the largest fully collisional star-by-star model of a Plummer sphere presented to date.

Zoom In Zoom Out Reset image size

Finally, a critical test of any N-body scheme is how well it conserves energy, since this serves as an independent benchmark of the self-consistency of the physics. By itself, two-body relaxation as implemented in the Hénon algorithm automatically conserves energy; in practice, with the energy-conserving corrections described in Section 2.9, a CMC run typically conserves energy to within one part in ∼ 10⁵ over many relaxation times. We show the differential change in energy for both models as a function of time in Figure 4.

Zoom In Zoom Out Reset image size

To demonstrate the capability of CMC to generate realistic cluster models, we generate 10 models of star clusters with random initial conditions drawn from the same distribution function. We begin by sampling positions and velocities for 10⁶ particles using a King (1966) profile (with W₀ = 6) and a virial radius of 1 pc. The properties for each star and binary are sampled using the same initial condition generators as COSMIC: each star is first assigned a mass from a Kroupa (2001) IMF. We randomly assign 10% of those stars to be binaries, keeping the initial IMF draw as the primary mass, and assigning a secondary mass from a uniform distribution between 0.1 and 1 times the primary mass. For isolated binaries, the "independent" sampler in COSMIC assigns binary semimajor axes from any number of user-specified distributions, e.g., flat in the log-period (Abt 1983), or the distributions from Sana et al. (2012) and Renzo et al. (2019). However, cluster environments come with an additional complication: many wide binaries that are observed in galactic fields (and are reproduced by the distributions in COSMIC) are significantly larger than the typical interparticle separation of the cluster and would be very rapidly destroyed by dynamical encounters (or simply not formed at all). To avoid this, we truncate the orbital periods at the local hard−soft boundary for a binary with the given mass at that position in the cluster (this approximately corresponds to limiting the minimum orbital velocity to the local three-dimensional velocity dispersion). For the binaries considered here, we assume a flat-in-log distribution of semimajor axes and a thermal distribution of eccentricities (Ambartsumian 1937; Heggie 1975). We set stellar metallicities equal to Z = 0.0002 (approximately 1% of the solar value) and do not place the clusters in an external galactic potential.

In Figure 5, we show a projected snapshot of one of the cluster models rendered using the fresco package (Rieder & Pelupessy 2019), which incorporates the stellar positions, effective temperatures, and radii to construct a mock HST observation. We also show the surface brightness and velocity dispersion profiles for that snapshot, calculated using the cmctoolkit (Rui et al. 2021a, 2021b). The toolkit can directly convert the HDF5 output snapshots of CMC into observational quantities, including surface brightness and velocity dispersion profiles, number density profiles, counts of blue stragglers and other stellar types, and so on. Combined with the public release of CMC, these two software packages contain all the necessary tools to create models of GCs and SSCs that can be directly compared to observations in the local universe.

Zoom In Zoom Out Reset image size

One of the main advantages of N-body modeling is the ability to study the internal dynamics of lower-luminosity components of the cluster. In Figure 6, we show the populations of WDs, NSs, and BHs that are present throughout the cluster model shown in Figure 5 (each overlaid over the synthetic observation), as well as the projected and smoothed 2D number densities as calculated with the cmctoolkit. As expected, the heavier components—the carbon–oxygen and oxygen–neon WDs, NSs, and BHs—have segregated into the central regions of the cluster. The BHs, in particular, are highly concentrated, with 50% of the BHs lying within 0.2 pc of the cluster center, where the number density of BHs is more than 7000 pc⁻³.

Zoom In Zoom Out Reset image size

In a realistic GC, massive stars evolve in ∼10–100 Myr after cluster formation, leaving behind a large population of BH remnants. These BHs rapidly sink to the center, where they begin to participate in the three- and four-body interactions that form and dynamically harden binaries. This process continues without limit until all the BH binaries are either ejected from the cluster (having each also ejected several single BHs during three-body encounters) or merge owing to GW emission. The rate of BH depletion is determined by the energy required to support the cluster against continued collapse (e.g., Breen & Heggie 2013), and in clusters with sufficiently small initial virial radii, this process can completely remove the BH subsystem from the cluster, producing the core-collapsed clusters observed in the Milky Way (Kremer et al. 2019a). Overall, the dynamics of BH subsystems are expected to play a key role in shaping the overall evolution of DSCs (e.g., Merritt et al. 2004; Mackey et al. 2007, 2008; Breen & Heggie 2013; Wang et al. 2016; Chatterjee et al. 2017b; Arca Sedda et al. 2018; Kremer et al. 2018b, 2020b, 2020c; Antonini & Gieles 2020; Weatherford et al. 2021). The particular cluster simulation shown in Figure 5 still retains 125 BHs by the present day, producing the large core radius observed in the surface brightness profile.

In addition to the global dynamical and observable properties of present-day clusters, CMC is ideal for understanding the rate and properties of many high-energy transients and compact binary sources that originate in DSCs. As previously stated, the BHs in the central regions of a cluster will continually form binaries and undergo hardening encounters, providing a finite energy source to support the cluster against gravothermal collapse. While this behavior is critical to understanding the overall evolution of the cluster and its present-day appearance, it also produces a large number of BBHs and is thought to provide a key formation channel for the GW sources detected by LIGO, Virgo, and KAGRA (LVK). This process has been studied theoretically for many years (Kulkarni et al. 1993; Sigurdsson & Hernquist 1993; Portegies Zwart & Mcmillan 2000; O'Leary et al. 2007; Downing et al. 2010, 2011; Aarseth 2012; Tanikawa 2013; Bae et al. 2014; Askar et al. 2016; Fragione & Kocsis 2018), including in several papers using CMC (e.g., Rodriguez et al. 2015, 2016a, 2016b, 2018a, 2018b, 2019,2020, 2021b; Chatterjee et al. 2017a, 2017b; Kremer et al. 2019c, 2020c, 2020a; Martinez et al. 2020; González et al. 2021; Holgado et al. 2021; Martinez et al. 2021; Weatherford et al.2021).

With over 50 GW detections so far, much work has been done trying to determine the origin of the observed population of merging BBHs. One well-studied technique to discriminate between BBHs formed dynamically in clusters and those formed from isolated binary evolution in the field is to measure the spins of the BHs (e.g., Rodriguez et al. 2016d; Farr et al. 2017, 2018), since dynamically formed BBHs are expected to have isotropically distributed spin–orbit misalignments (vs. the relatively small spin–orbit misalignments expected from binary star evolution). Using the distribution of BBH spins, the latest catalog from the LVK suggests that between 25% and 93% of these systems may have been assembled dynamically (The LIGO Scientific Collaboration et al. 2021), with more sophisticated analyses of BBH formation models (including those created using CMC) reaching similar conclusions (Wong et al. 2021; Zevin et al. 2021).

In addition to statistical evidence, individual BBH events have also suggested a dynamical origin. One particular system, GW190521, was detected with at least one component mass inside the "upper mass gap," where individual BHs are prevented from forming by stellar collapse due to pulsational pair instabilities, which eject large amounts of material from (and potentially destroy) any star with a helium core mass between about 40 and 130 M_⊙ (e.g., Woosley 2017). While it is extremely difficult (though not impossible) to describe such systems as the product of isolated stellar evolution (e.g., Farmer et al. 2019; Belczynski 2020), these systems can be easily produced in DSCs, through either the repeated mergers of BHs (e.g., Rodriguez et al. 2018b, 2019; Fragione & Silk 2020; Fragione et al. 2020a) or massive star mergers occurring prior to BH formation (e.g., Di Carlo et al. 2019; Kremer et al. 2020a; González et al. 2021; Weatherford et al. 2021).

Figure 7 shows the total masses and spins of all BBHs that merge within 13.8 Gyr of cluster formation from our 10 GC models. We subdivide the population according to the "generation" of the BHs in the binary, with first generation, or 1G, BHs being those created from stellar collapse, and second generation, or 2G, referring to those created in a previous merger. The more massive BHs, particularly those within the mass range consistent with GW190521, are predominantly the result of previous mergers, with four of the six systems being composed of 2G+2G BBHs (that is, binaries whose components are both the product of previous mergers). The remaining two BBHs are created from previous stellar mergers of massive stars prior to their collapse into BHs. These stellar-merger products can contain unique envelope/core mass ratios, allowing them to bypass the typical pulsational pair-instability limits for single massive stars. While repeated BH mergers are the dominant mechanism for producing GW190521-like systems in these models, clusters with higher central densities, fractal initial conditions, or different stellar binary fractions can prefer formation of repeated mergers through stellar collisions (e.g., Di Carlo et al. 2019; Kremer et al. 2020a; González et al. 2021). We note that while our models also ignore the possibility of accretion from BH/star collisions, the majority of these occur after the most massive stars (including stellar-merger products) have expired, since the BHs typically require ∼100 Myr to dynamically segregate into the cluster center. As such, M_BH > M_star for the majority of BH/star collisions.

Zoom In Zoom Out Reset image size

We also show in Figure 7 the distribution of the effective spins of the cluster BBH population as a function of the mass. This mass-weighted projection of the BH spins onto the orbital angular momentum, ${\chi }_{\mathrm{eff}}\equiv \hat{L}\cdot \left[({m}_{1}{{\boldsymbol{\chi }}}_{1}+{m}_{2}{{\boldsymbol{\chi }}}_{2})/({m}_{1}+{m}_{2})\right]$, is the spin parameter most easily constrained by GW detectors (e.g., Ajith et al. 2011; Abbott et al. 2016). Since we have assumed that BHs from collapsing stars are born with zero spin, consistent with recent theoretical studies of angular momentum transport in massive stars (e.g., Fuller & Ma 2019; Fuller et al. 2019), BBHs whose components are formed from stellar collapse also have zero effective spins across all masses. On the other hand, when BHs are created from previous mergers, their spins inherit the angular momenta of their parent binaries, producing a population of spinning BBHs. Because the BH spins are assumed to be isotropically distributed on the sphere, the median effective spin for all cluster BBHs is still centered at χ_eff = 0. This is consistent with the observation of GW190521, which appears to be consistent with having zero effective spin. However, as consistent with our previous results (e.g., Rodriguez et al. 2018b, Figure 3), any population of 2G+2G BBHs (like the ones that predominately form GW190521-like binaries here) should have a range of effective spins centered at χ_eff = 0.

As a demonstration of the utility of the CMC/COSMIC integration, we directly compare the BBH population created through dynamical processes to that created from the evolution of isolated stellar binaries. We evolve a population of 5 × 10⁴ binaries with masses ≥18 M_⊙ using the same stellar physics, stellar metallicity, and binary initial conditions (the truncation in binary orbital period) with COSMIC. In Figure 7, we show one of the most distinct features of the BBHs produced by GCs: the production of repeated mergers of BHs in star clusters. When two BHs merge in a GC, their merger product receives a recoil kick due to the asymmetric emission of GWs, which depends on both the mass ratio of the system and the spins of the BHs (e.g., Merritt et al. 2004; Campanelli et al. 2007; Lousto & Zlochower 2008). Because three-body encounters preferentially form BBHs with near-equal-mass components (e.g., Sigurdsson & Phinney 1993), it is the BH spins that primarily determine merger retention in GCs. While large spins (χ ∼ 1) leave few second-generation BHs in the star cluster, mergers of BHs born with smaller spins (χ ≲ 0.2) can be retained in many GCs, where they will continue to participate in dynamical processes, form binaries, and merge again (Rodriguez et al. 2019). This process can produce a significant population of BBHs with large (χ ∼ 0.7) spins and masses in the pulsational pair-instability mass gap.

While an observation of a BBH with χ_eff < 0 is strongly suggestive of dynamical formation, isolated binary evolution can also produce BBHs with spin–orbit tilts greater than 90° given the right orientation and magnitude of the supernova kicks (Kalogera 2000; Rodriguez et al. 2016d). Therefore, any analysis of BBH spin–orbit misalignments must include both dynamical and binary evolutionary channels using identical stellar physics. The integration of COSMIC into CMC makes this trivial; in Figure 8, we show the masses and spin–orbit misalignments (θ_LS) from the 10 GC models and 5 × 10⁴ isolated massive binaries evolved using COSMIC. In the top panel, the limit imposed by pulsational pair-instability physics on the BH mass function is readily apparent, with no BHs forming with masses greater than 44.5 M_⊙, and no BBHs forming with total masses greater than 89 M_⊙ (a number that is largely insensitive to choices of accretion and other binary physics; van Son et al. 2020). Figure 8 also clearly demonstrates the difference in the spin–orbit misalignments between isolated binaries and cluster-formed BBHs. Given the supernova prescriptions in COSMIC (here we have used the "rapid" prescription from Fryer et al. 2012), BHs with final core masses ≳20 M_⊙ form via the "direct collapse" or failed supernova scenario, where the BHs form without any natal kick or change to the supernova orbital parameters. For BBHs with total masses ≲45 M_⊙, the supernova kicks can still produce significant spin–orbit misalignments, with the 99th percentile of allowed misalignments approaching complete anti-alignment (θ_LS = 180°). However, as the mass increases, the amount of spin–orbit misalignment decreases, until eventually all BBHs are aligned with their orbital angular momenta. Contrast this with BBHs formed in clusters, which are isotropically distributed on the sphere (i.e., $P({\theta }_{\mathrm{LS}})d{\theta }_{\mathrm{LS}}\propto \sin ({\theta }_{\mathrm{LS}})d{\theta }_{\mathrm{LS}}$). These results can be easily reproduced and expanded on using the apples-to-apples comparisons enabled by the CMC and COSMIC integration.

Zoom In Zoom Out Reset image size