The Impact of Flyby Interactions on the Rotation of Open Clusters

Introduction

Coeval clusters of stars have long been utilized as fundamental laboratories for astrophysics. In particular, dynamical signatures of these clusters inform our knowledge across numerous domains; the rotation rates of constituent stars constrain gyrochronology [e.g., 26, 27, 28], their dissolution constrains the tidal field they reside within [e.g., 29, 30] and the inhomogeneities within it [e.g., 11, 31, 32], samples of globular clusters reveal the assembly history of our Milky Way [e.g., 51, 52], kinematics of young clusters reveal substructure in star forming environments [e.g., 54, 53], and their internal kinematics inform us of expectations for gravitational wave detections [e.g., 33].

One such dynamical signature, bulk rotation of the constituent stars about a central axis, has been studied extensively in clusters via numerical simulations. Orbit-averaged solutions to the 2D Fokker-Planck equations show that initially rotating clusters experience an accelerated mass loss and subsequent core-collapse [37], and that rotating clusters evolve more rapidly [38]. Simulations of clusters in a tidal field reveal the preferential escape of stars on prograde orbits relative to the clusters’ orbital velocities [40, 41, 39], the long-term evolution of bulk rotation towards an approximately solid-body limit equal to one-half the clusters’ orbital frequencies [39, 19], and the two-body driven degradation of primordial angular momentum and spin through the clusters’ lifetimes [20, 42, 84] which occurs on the order of the half-mass relaxation time [21].

These simulations, by and large, have been performed with globular clusters (GCs) in mind; rotational detections in GCs [e.g., 44, 35, 43] are more robust than those observed in open clusters (OCs). Their high star counts and advantageous location out of the Galactic plane allow for significant ($>3\sigma$) rotation detections [e.g., 36, 34].

The ever-increasing astrometric precision of Gaia[45], however, has finally shed light onto rotation in OCs by revealing correlated proper motions and/or radial velocities. [46] detect rotation in eight OCs, with an additional nine candidates. [48, 47] report $\geq 3\sigma$ rotation detections in four nearby open clusters: Alpha Persei, the Pleiades, Praseape, and the Hyades. [49] report rotation signatures for 10 OCs (with 16 additional candidates), albeit at a lower precision (typically $<2\sigma$). [63] recently published the largest rotating OC sample to date; they find significant rotation across a range of cluster ages, that the rotation signal weakens as a function of age, and that older OCs may be slightly more likely to display prograde rotation, though the overall alignment between the orbital and internal angular momentum vectors is still weak. These results are in broad agreement with the expectation from theory that older clusters should slowly lose their spin on a timescale comparable to the two-body relaxation time [e.g., 21].

The source of cluster rotation is often attributed to the primordial spin of their parent molecular clouds [e.g., 55, 34, 36]. Velocity gradients indicative of net rotation in giant molecular clouds have been characterized both by observation [e.g., 23, 50] and hydrodynamic simulations [e.g., 56, 24]. One such rotation detection in particular, however, challenges this paradigm; [16] report a mean rotational velocity in the old open cluster M67 [, 60] of $0.81 \pm 0.06~\mathrm{km ~s^{-1}}$ within its Jacobi radius. [22] find that M67 is dynamically old, having undergone $\approx 13$ relaxation times and having lost $\approx 90 \%$ of its initial mass. [39, 19] find that clusters orbiting in a tidal field born with primordial rotation evolve towards solid body rotation with an asymptotic angular velocity of $\Omega / 2$, where $\Omega$ is the angular velocity of the cluster about the host galaxy. Given the size and orbital velocity of M67, this corresponds to a rotational velocity of order $\sim 0.1 ~\mathrm{km ~s^{-1}}$, making the observed rotational velocity far larger than expected. M67 is not the only cluster challenging theoretical expectations for spin down; [63] find at least $10$ clusters classified as significant rotators with ages greater than $1 ~\mathrm{Gyr}$. Note that M67 does not appear in this rotation sample, but [63] make cuts to ignore clusters with prominent extra-tidal features, which M67 has [62, 16]. Additionally, both [46] and [49] classify Ruprecht 147 [old; 73] as a rotating OC. We choose to highlight M67 specifically in this work due to its large reported rotational velocity, its old age, and the extensive studies of its dynamical history in the literature [e.g., 22, 74].

Motivated both by these observations and the pending Gaia DR4 release, we present a theoretical framework, along with numerical validation, to describe another potentially significant source of cluster rotation: tidal shocking from flyby interactions with massive galactic over-densities. We tailor our analysis for OCs, and thus consider interactions with giant molecular clouds (GMCs) in this work, although our methodology could be abstracted to any star cluster interacting with a massive over-density. In § 2, we derive an analytic approximation for the angular momentum donated to stars via a GMC–cluster encounter and describe the initial conditions and methodology behind our $N$-body simulations of GMC–cluster encounters. We present the resultant changes in energy, mass, angular momentum, and rotational velocity of our simulated clusters in § 3. We explore the limitations of our analytic model in § 4 and speculate on the existence of observations of dynamically relaxed OCs as well. Finally, we provide our concluding thoughts in § 5.

Methods

Spin-Up via Individual Flyby Interactions

Consider a spherical star cluster of mass $M_c$ described by a Plummer density profile [6] with half-mass radius $r_h$, at rest at the origin with an oncoming GMC with mass $M_{\rm GMC}$, impact parameter at infinity $b$, and relative velocity at infinity $V$. To derive the expected change in internal cluster angular momentum for a given GMC–cluster interaction, we adopt the following simplifying assumptions to derive an analytic, closed-form expression:

1. The impulse approximation is valid, i.e., $t_{\rm enc} \ll t_{\rm dyn}$, where the encounter time is defined as $t_{\rm enc} \equiv b/V$. We determine the dynamical time of the cluster from [13] as

   $$\begin{equation} t_{\rm dyn} \sim 2 \times 10^{-2} ~\mathrm{Myr} \left( \frac{M_c}{10^6 ~\mathrm{M_{\odot}}} \right)^{-1/2} \left( \frac{r_{\rm vir}}{1 ~\mathrm{pc}} \right)^{3/2}, \tag*{(1)} \end{equation}$$

   where virial equilibrium is assumed. For a Plummer profile, the virial radius $r_{\rm vir}$ is equated to the half-mass radius by $r_{\rm vir}=1.3 r_h$ [e.g., 14]. We thus ignore any orbital motion of stars during the encounter.

2. The encounter is weak, such that $b \gg b_{90}$, where $b_{90} = G (M_c + M_{\rm GMC}) / V^2 \approx G M_{\rm GMC}/V^2$, and $M_{\rm GMC} \gg M_c$.

3. The impact parameter is much larger than the size of the system, i.e., $b \gg r_h$, such that the GMC can be treated as a point mass.

In this work, we focus solely on a cluster of mass $M_c = 1000 ~\mathrm{M_{\odot}}$, impact parameters $b < 100 ~\mathrm{pc}$, and a relative velocity at infinity of $V = 20 ~\mathrm{km ~s^{-1}}$. At these parameters, note that $t_{\rm enc} < 4.9 ~\mathrm{Myr}$ and $t_{\rm dyn} \approx 4.9 ~\mathrm{Myr}$, computed from Equation (1). Though assumption (1) is not always strictly satisfied at every simulated impact parameter, $N$-body simulations show that the impulse approximation is a fairly reasonable approximation of the relative energy gain for $t_{\rm enc} \leq t_{\rm dyn}$ for a point mass perturber [see the right panel of Figure 1 in 15]. Additionally, the stars that receive the largest kicks relative to the cluster center of mass are those at large radii, where the local dynamical time is larger. Thus, we consider the impulse approximation a reasonable assumption in computing the angular momentum change for all simulations explored hereinafter.

We orient the GMC’s velocity along the $z$-axis, and place the GMC’s orbit in the $yz$-plane, such that the GMC is at position $(0, b, 0)$ at closest approach. We assume the GMC travels on a straight-line trajectory, unaffected by the star cluster’s gravity, which is a reasonable assumption for a hyperbolic trajectory with $M_{\rm GMC} \gg M_c$. A star in the cluster with initial position $\mathbf{r}_i = (x_i, y_i, z_i)$ relative to the cluster origin will experience both a velocity kick $\Delta \mathbf{v}_i$ relative to the cluster center of mass velocity, as well as a displacement relative to the center of mass $\Delta \mathbf{r}_i$ caused by the acceleration from the GMC during the encounter. The change in angular momentum of the cluster as a whole, relative to its center of mass, is then

where $\Delta \mathbf{L}^{(v)} \equiv \sum_i m_i \mathbf{r}_i \times \Delta \mathbf{v}_i$ and $\Delta \mathbf{L}^{(r)} \equiv \sum_i m_i \Delta \mathbf{r}_i \times \Delta \mathbf{v}_i$.

Starting from Equations (3.54a) and (3.54b) in [12], which give the final velocity kick for an object in the center of mass frame from a hyperbolic encounter, we can approximate the change in angular momentum for the star cluster as

where $b_{\rm eff} \equiv \sqrt{x^2 + (b - y)^2}$, $\rho(r)$ is the Plummer density profile $\rho(r) = (3 M_c / 4\pi a^3 ) \left[ 1 + r^2 / a^2 \right]^{-5/2}$, and $a$ is the Plummer radius, related to the half-mass radius by $r_h = 1.305 a$. We derive Equations (2) and (3) in detail in Appendix § 6, and note that the only surviving term of order $b_{90}^2$ or lower is the $x$-component of $\Delta \mathbf{L}$ under our chosen coordinate system for the encounter.

Following assumptions (2) and (3), we can derive a closed form expression for the cluster’s change in internal angular momentum. Writing

where $\epsilon \ll 1$, and keeping all leading order terms following expansion about $b_{90} / b$ and $\epsilon$, reduces Equation (3) to

where $\langle r^2 \rangle$ is the mean squared radius of stars in the cluster. Note that, at given encounter parameters $M_{\rm GMC}$, $b$, and $V$, the change in angular momentum of the cluster is proportional to $M_c \langle r^2 \rangle$, which is some factor multiplied by the cluster’s moment of inertia.

In Figure 1, we show a schematic of a flyby encounter between a GMC and a star cluster in our chosen coordinate system to demonstrate how the angular momentum of the cluster changes. We show the cluster both before (gray points) and after (black points) the interaction, and focus on the position and velocity displacements relative to the center of mass motion of two stars in particular; one that sits between the center of mass and the GMC (star 1) and another that sits below the center of mass (star 2). Relative to star 2, star 1 is pulled towards the $-z$ direction during the encounter, and receives a velocity kick in the $y$ and $z$ directions. Both the displacement and the velocity kick of star 1 relative to star 2 lead to changes in angular momentum of the same order in the $x$ direction.

In Figure 2, we plot the change in angular momentum in the $x$ direction $\Delta L_x$ as a function of impact parameter $b$ for a cluster of mass $1000 ~\mathrm{M_{\odot}}$, half-mass radius $r_h = 3 ~\mathrm{pc}$, and relative encounter velocity $V = 20 ~\mathrm{km ~s^{-1}}$. Solid lines show the numerical solution to Equation (3); for the radial portion of the integral, we integrate from $r=0$ to $r=r_t = 14.75 ~\mathrm{pc}$, which is the truncation radius used when generating initial conditions for our $N$-body models (see § 2.2). Dashed lines refer to Equation (5), where we have solved for the mean squared radius as $\langle r^2 \rangle = M_c^{-1}\int_0^{r_t} r^2 \rho(r) d^3r$ using a Plummer density profile. The two formulations of $\Delta L_x$ agree well for $b \gtrsim 20 ~\mathrm{pc}$ at our chosen GMC masses of $10^4 ~\mathrm{M_{\odot}}$ and $10^5 ~\mathrm{M_{\odot}}$.

N-body Simulations

To simulate encounters between clusters and GMCs, we utilize REBOUND [1] to perform all $N$-body simulations, using the IAS15 integrator [3]; a non-symplectic, 15th order Gauss-Radau integrator with adaptive time-stepping. We store simulation outputs using the Simulationarchive [4].

We initialize star clusters using McLuster [5]. Each star cluster is generated from a Plummer profile [6] with a total mass $M_c = 1000 ~\mathrm{M_{\odot}}$ and half-mass radius $r_h = 3 ~\mathrm{pc}$, which corresponds to a Plummer radius $a_c = 2.3 ~\mathrm{pc}$. By default, these models are truncated at a tidal radius of $r_t = 14.75 ~\mathrm{pc}$. We generate stellar masses using the canonical Kroupa [7] initial mass function, with a minimum and maximum stellar mass of $0.08 ~\mathrm{M_{\odot}}$ and $100 ~\mathrm{M_{\odot}}$, respectively. We do not consider any primordial binaries or mass segregation, nor do we consider stellar evolution in this work.

Throughout all simulations, we treat the GMC as a point mass. Galactic GMCs are found to be extended objects with sizes roughly proportional to $M_{\rm GMC}^{1/2}$ [e.g., 9, 8, 10], such that a typical GMC of mass $10^5 ~\mathrm{M_{\odot}}$ is $\sim 30 ~\mathrm{pc}$ in extent [see 8]. For simulations with $b < 30 ~\mathrm{pc}$, the energy and angular momentum gain of the simulated cluster will be larger than it would have been had the GMC been ascribed a realistic size. Thus, we consider the results of our simulations to be an upper limit on the angular momentum gain of a cluster that undergoes an encounter with a given set of impact parameters, as exploring the effect of the GMC size on the imparted angular momentum is outside of the scope of this work.

We carry out the GMC–cluster encounter in the same manner as [11]. That is, we begin the simulation with the cluster centered at the origin, and the GMC placed at position $(0, b, -D)$ with velocity $(0, 0, V)$, where $D = 500 a_c$, corresponding to $\approx 1150 ~\mathrm{pc}$ given our chosen cluster initialization. Here, $b$ and $V$ are the impact parameter and relative velocity at infinity, respectively. The GMC is placed on a straight-line trajectory, such that $\mathbf{r}_{\rm GMC}(t) = (0, b, -D + V t)$. In all simulations, the GMC is at least $10$ times more massive than the cluster, making a straight-line path a reasonable approximation. Note that, while the GMC is on a fixed trajectory, the cluster’s center of mass is allowed to respond to the GMC, such that the closest approach distance and relative velocity at impact are smaller and larger than $b$ and $V$, respectively.

To avoid sudden jumps in the potential felt by simulated stars, we ramp up (and down) the mass of the GMC in a smooth fashion, based on the distance $r_{\rm GMC}$ from the GMC to the origin. The simulated GMC mass $M_{\rm GMC}^{\star}$ as a function of $r_{\rm GMC}$ is

where $x \equiv 10 (D - r_{\rm GMC})/D$, and $M_{\rm GMC}$ is the peak mass of the GMC. We incorporate the simulated GMC as an external force in our simulations, not an additional particle, utilizing REBOUNDx [2] to do so.

The total duration of each GMC–cluster simulation is three times the amount of time taken for the GMC to reach position $(0, b, 0)$, i.e., $T_{\rm total} = 3 T_{\rm enc}$, where $T_{\rm enc} \equiv D/V$ (not to be confused with the relevant encounter timescale $t_{\rm enc} \equiv b /V$). This gives stars with velocities greater than the escape velocity sufficient time to leave the cluster post-encounter.

In this work, we present GMC–cluster simulations with $M_{\rm GMC} = 10^4$ and $10^5 ~\mathrm{M_\odot}$, impact parameters $b = 5$, $10$, $20$, $40$, $60$, $80$, and $100 ~\mathrm{pc}$, and a constant relative velocity at infinity of $V = 20 ~\mathrm{km ~s^{-1}}$. For each of the possible $14$ combinations of $M_{\rm GMC}$, $b$, and $V$, we run $10$ realizations with different initial seeds for the cluster. At our chosen $V$, each simulation runs for a total of $T_{\rm total} = 168.6 ~\mathrm{Myr}$. For comparison, we also run $10$ simulations where a cluster with the same initial conditions evolves in isolation for the same time $T_{\rm total} = 168.6 ~\mathrm{Myr}$, which corresponds to $\approx 2$ initial half-mass relaxation times for a cluster with a Plummer density profile with mass $M_c = 1000 ~\mathrm{M_{\odot}}$ and $r_h = 3 ~\mathrm{pc}$.

Results

Energy Gain and Mass Loss from the Cluster

In the left panel of Figure 3 we compare the relative energy gain in each simulation with the analytic prediction for such a GMC–cluster encounter by [11] as a function of impact parameter $b$. The quantity $\Delta E_{\rm sim}$ is taken as the difference in total energy of all stars in the simulation at the beginning and end, i.e., $E_f - E_0$, computed internally via REBOUND, where the center of mass motion has been subtracted off. The analytic prediction $\Delta E_{\rm Gieles}$ is computed from their Equation (15) of [11] using the half-mass radius $r_h$ at the initial timestep, after converting our impact parameters, which are defined at infinity, to theirs, which are defined at closest approach between cluster and GMC. Each diamond represents the median of $10$ simulations at a given $b$ and $M_{\rm GMC}$, and the error bars represent the 25th–75th percentile spread, i.e., the inter-quartile range, of $\Delta E_{\rm sim} / \Delta E_{\rm Gieles}$ for those $10$ simulations. We see good agreement at low ($b < 30 ~\mathrm{pc}$) impact parameters (and thus, a large relative energy gain), which is the most interesting regime in terms of adding net angular momentum to stars in the cluster (see Figure 2). At larger impact parameters, while our simulated results still show consistency with the predictions of [11], the spread in the energy gain increases dramatically. This is expected; [11] predict only the second order change in the energy, as the first order change in internal energy $\Delta E^{(1)} = \sum_i m_i \mathbf{v}_i \cdot \Delta \mathbf{v}_i$ sums to zero for a perfectly isotropic velocity distribution. For a finite $N$ system, $\Delta E^{(1)}$ becomes comparable to $\Delta E_{\rm Gieles}$ for less energetic interactions (i.e., high $b$, low $M_{\rm GMC}$). Nonetheless, the median and spread of all simulations explored in this work see energy increases that broadly agree with theoretical expectation, validating our continued employment of the impulse approximation.

In the right panel of Figure 3, we show the fractional mass loss of the cluster, where $\Delta M \equiv M_{c,0} - M_{c,f}$, as a function of the predicted relative energy gain of the encounter, where $E_{0}$ is the total initial cluster energy in the center of mass frame, computed internally via REBOUND. $M_{c,0}$ and $M_{c,f}$ refer to the bound mass at the initial and final timestep. Boundedness at a given timestep is determined iteratively; we compute the relative energy of each star as $\epsilon_i = \tfrac{1}{2}|\mathbf{v}_i - \bar{\mathbf{v}}|^2 + \Phi_i (\mathbf{r}_i)$, where $\mathbf{v}_i$ is the velocity of the $i$-th star relative to the center of mass velocity of the current bound subset, and $\Phi_i (\mathbf{r}_i)$ is the potential of the current bound subset at position $\mathbf{r}_i$, computed via direct summation. Stars with $\epsilon > 0$ are removed from the bound subset, and the procedure is repeated until the bound membership has converged. Filled diamonds again represent the median of $10$ simulations at a given $b$ and $M_{\rm GMC}$, and error bars represent the 25th–75th percentile spread. The gray line represents the empirical fit found by [11] between the energy gain and mass loss, i.e., $\Delta M / M_0 \approx 0.25 \Delta E / |E_0|$.

Our simulations cover a large range in $\Delta E / |E_0|$, while the simulations explored by [11] probe relative energy gains from $10^{-3}$ to $10$. Our simulations with $\Delta E_{\rm Gieles} / |E_0| \approx 1$ agree well with their empirical fit, but deviate at larger and smaller relative energy gains. We also show, in orange, both the median and fractional spread in the fractional mass loss for our simulations where the cluster evolved in isolation, to control for the effects of mass lost via two-body relaxation. We see that, for weak encounters with $\Delta E_{\rm Gieles} / |E_0| < 10^{-2}$, the encounter is negligible in terms of causing mass loss, as the entirety of the mass lost from the cluster can be explained via two-body processes.

Angular Momentum

In Figure 4, we present a comparison of the change in angular momentum from our simulations with the analytic approximation for $\Delta L_x$ for weak encounters under the impulse approximation, i.e., Equations (3) (solid lines) and (5) (dashed lines). Filled diamonds represent the median change in internal angular momentum in the $x$-direction across a given subset of $10$ simulations, while error bars again represent the 25th–75th percentile spread. Both Equations (3) and (5) are computed using simulations snapshots two encounter times before the closest passage between cluster and GMC, i.e., $t_{\rm pre-enc} = T_{\rm enc} - 2 b /V$, where all positions and velocities have been shifted to the cluster’s center of mass frame. We choose this timestep because the cluster has only been minimally perturbed by the GMC [see Figure 6 in 11], while also having relaxed from its initial conditions. We compute Equation (3) exactly; for each simulation Equation (3) is evaluated via direct summation over the stellar positions at $t_{\rm pre-enc}$ over a grid of $b$ values from $1$ to $105 ~\mathrm{pc}$. The solid red and black lines show the median of these curves for simulations with $M_{\rm GMC} = 10^4 ~\mathrm{M_{\odot}}$ and $M_{\rm GMC} = 10^5 ~\mathrm{M_{\odot}}$, respectively. To compute Equation (5), we evaluate $\langle r^2 \rangle$ relative to the center of mass at $t_{\rm pre-enc}$ for every simulation explored in this work, including the control set with no GMC. We use the median $\langle r^2 \rangle$ at $t_{\rm pre-enc}$ across all simulations to evaluate Equation (5), as taking $\langle r^2 \rangle$ at $t=0$ ignores the pre-encounter radial evolution of the cluster, and would under-predict $\Delta L_x$. Note that the angular momenta of all stars relative to the center of mass are considered here, regardless of whether or not the stars are unbound to the cluster by the encounter or by two-body processes.

For all encounters tested, we see good agreement between simulation and theory. At impact parameters $b > 40~\mathrm{pc}$, the change in angular momentum for clusters that have encounters with $10^4 ~\mathrm{M_{\odot}}$ GMCs are dominated by stochastic variability due to the finite number of stars. The first order contribution to the change in angular momentum $\Delta L_{x}^{(1)}$, which vanishes under perfect symmetry, is comparable to the signal $\Delta L_x$ in this regime. Since $\Delta L_x \propto M_{\rm GMC}^2$, encounters with $M_{\rm GMC} = 10^5 ~\mathrm{M_{\odot}}$ produce a signal two orders of magnitude larger than those with $M_{\rm GMC} = 10^4 ~\mathrm{M_{\odot}}$, and thus, the signal-to-noise ratio is adequate out to values of $b = 100 ~\mathrm{pc}$. For very close encounters, i.e., $b=5$, Equation (5) over-predicts the angular momentum. At low impact parameters, the size of the GMC also becomes important; we expect our simulations and analytic prescription to provide an upper limit on $\Delta L_x$ in this regime, as a point mass approximation overestimates the energy gain in such an encounter [see Figure 7 in 11].

While our results show that the angular momentum change of all stars due to a GMC–cluster encounter is predictable analytically, we aim to explore the rotational properties of the bound subset of stars that remain in the cluster post-encounter. We consider a parameter more observationally obtainable than angular momentum, and thus define the rotational velocity at some epoch $v_{\rm rot}$ to be the mean of the tangential velocity of all bound stars about the cluster’s principal axis of rotation ($v_{\phi}$), defined by the angular momentum direction $\hat{\mathbf{L}}$ at that epoch. This follows from the approach used to detect angular momentum in M67 by [16].

In Figure 5, we show the change in the magnitude of the rotational velocity $v_{\rm rot}$ between the first and last snapshot for all simulations with $\geq 50$ bound members at simulation end. Again, we plot the median and 25th–75th percentile spread for a given simulation set. We determine boundedness in the same manner as in § 3.1. At both the initial and final epoch, we determine the cluster density center using $k=32$ nearest neighbors from all bound stars. Positions are measured relative to this density center, and velocities are measured relative to the center of mass velocity of the bound subset of stars. In the left panel of Figure 5, we include all stars energetically bound to the system when determining $v_{\rm rot}$, while in the right panel, we only consider energetically bound stars within $15 ~\mathrm{pc}$ of the cluster density center. We also show, as an orange line and band, the same difference in magnitude of initial and final rotational velocity for simulations with no GMC–cluster encounter. The vast majority of clusters undergoing GMC encounters do not gain appreciable spin in their bound subset when compared with clusters subject only to two-body relaxation, implying such encounters are not a suitable mechanism to “spin-up" clusters on a large scale. No single simulation experiences a gain in rotational velocity of $> 0.1 ~\mathrm{km ~s^{-1}}$, indicating that even the most optimal encounters do not appreciably cause an increase in cluster rotation after the stars unbound by the encounter have had sufficient time to leave the system.

Apparent Spin-Up of a Cluster

While there is no considerable “spin-up" of the cluster at simulation end, the rotational properties of the cluster a short time after the encounter are still of interest. In the top left panel of Figure 6, we instead compute the rotational velocity $5 ~\mathrm{Myr}$ before ($v_{\rm rot,-5}$) and after ($v_{\rm rot,+5}$) the time of encounter $T_{\rm enc}$ for each simulated cluster. The orange line and band represent the median and 25th–75th percentile spread of the control simulations, with rotational velocities measured $5 ~\mathrm{Myr}$ before and after the time at which a GMC would have reached closest approach. Throughout Figure 6, we define the bound subset of stars differently than in Figure 3 and Figure 5. Here, we find the density center amongst all stars, bound or unbound energetically, using $k=32$ nearest neighbors, and define a tidal radius of $r_t = 20 ~\mathrm{pc}$ [chosen as a comparative value with the Jacobi radius of M67 as reported by 17] relative to that density center; all stars within $r_t$ of the density center are considered radially bound. Positions are measured relative to this density center, and velocities are defined relative to the center of mass velocity of this radially bound subset in order to compute $\hat{\mathbf{L}}$ and $v_{\rm rot}$ at each given epoch. We consider only simulations with $\geq 50$ radially bound members. This approach is akin to how observers define bound and unbound members of a star cluster [e.g., 17] as finding the relative energy of individual stars requires accurate 6D phase space information, which is often very difficult to obtain. We see, in the top left panel of Figure 6, that when looking closely before and after the encounter and defining boundedness purely by a radial criterion, some simulations with high relative energy gains (i.e., smaller impact parameters) can produce appreciable rotation signatures when compared with the control simulations where no GMC is present.

The other three panels in Figure 6 highlight one example simulation of a cluster that has significant rotational velocity within $20 ~\mathrm{pc}$ of the density center $5 ~\mathrm{Myr}$ post-encounter. This particular simulation featured an encounter with impact parameter $b = 5 ~\mathrm{pc}$ and GMC mass $M_{\rm GMC} = 10^5 ~\mathrm{M_{\odot}}$. The top right panel shows the velocity about the rotation axis $v_{\phi}$ for each star, plotted against the 3D distance $r$ to the cluster density center. We show $v_{\phi}$ at $5 ~\mathrm{Myr}$ before (gold points) and after (navy points) the encounter, and display a rolling average as a solid line over $50$ points for each of the two epochs. Shortly after the GMC–cluster encounter, stars in the cluster outskirts beyond $10 ~\mathrm{pc}$ are characterized by a net positive value of $v_{\phi}$; the rolling mean of $v_{\phi}$ post-encounter reaches $\approx 0.5 ~\mathrm{km ~s^{-1}}$, within detectability limits for nearby OCs. The bottom two panels show the cluster in position and velocity space in the $xy$-plane at $5 ~\mathrm{Myr}$ after the encounter. This particular cluster still appears quasi-spherical in two dimensions in both position and velocity space despite the encounter, though the cluster develops a more elongated morphology in the $z$ direction due to the GMC passage. A cluster like this one, observed on the sky shortly after an energetic encounter, may appear as a spherical cluster with a detectable rotation signature, though much of this signature is produced by stars that have gained enough energy post-encounter to escape the system after a few crossing times.

Discussion

The results of our simulations confirm the theoretical framework outlined to describe the amount of angular momentum donated from a passing GMC onto a spherical star cluster in the impulsive regime. However, our results show that this angular momentum does not appreciably “spin-up" the cluster after stars energetically unbound from the system have been given ample time to leave it, as they carry away a significant fraction of the encounter-induced angular momentum. Naively, this makes sense, as the energy induced by such an encounter is proportional to $b^{-4}$ [Equation 15 in 11] while the angular momentum gain from Equation (5) is proportional to $b^{-3}$. There are two competing levers at a given relative velocity $V$ and GMC mass $M_{\rm GMC}$; reducing the impact parameter $b$ increases the energy gain, and thus the mass loss from the cluster, more quickly than it increases the angular momentum gain. In principle, many weak encounters could spin-up a cluster via a random walk process; we show below that this is infeasible given our simulation results.

From Figure 5, we see that the difference in rotational velocity for a single encounter has a spread of order $\Delta v_{\rm rot} = 0.01 ~\mathrm{km ~s^{-1}}$. Considering a random walk of encounters with the encounter parameters tested in our simulations, it would take $N \approx 10^4$ encounters to spin a cluster up from zero rotational velocity from $v_{\rm rot} = 1 ~\mathrm{km ~s^{-1}}$, taken from $v_{\rm rot} = \sqrt{N} \Delta v_{\rm rot}$. We can compare this to the expected number of interactions a cluster orbiting in the Solar neighborhood with Galactocentric radius $R_{\rm gal} \approx 8 ~\mathrm{kpc}$ will undergo with GMCs of mass $M_{\rm GMC} > 10^4 \mathrm{M_{\odot}}$, $b < 100 ~\mathrm{pc}$, and $V = 20 ~\mathrm{km ~s^{-1}}$. We estimate the number density of GMCs with masses $> 10^4 ~\mathrm{M_{\odot}}$ using the fitted surface density of GMCs as a function of Galactocentric radius given by [18], i.e., $\Sigma = 83 \exp(-R_{\rm gal} / 2 ~\mathrm{Kpc}) \mathrm{~M_{\odot} ~pc^{-2}}$, and convert this to a number density by taking $n = \Sigma / 2 \langle M_{\rm GMC}\rangle H$, where $\langle M_{\rm GMC}\rangle$ is the average mass of GMCs more massive than $10^4 ~\mathrm{M_{\odot}}$, and $H$ is the characteristic scale height of the disk. We assume $\langle M_{\rm GMC}\rangle = 8.4 \times 10^5 ~\mathrm{M_{\odot}}$ [18] and $H = 50 ~\mathrm{pc}$ [64, 81]. This corresponds to a number density of $n \approx 2 \times 10^{-8} ~\mathrm{pc^{-3}}$. The encounter rate is then computed as $\Gamma = n \sigma V = n \pi b_{\rm max}^2 V \approx 1.2 \times 10^{-2} ~\mathrm{Myr^{-1}}$, where $b_{\rm max} = 100 ~\mathrm{pc}$ and $V = 20 ~\mathrm{km ~s^{-1}}$. Over $1 ~\mathrm{Gyr}$, this encounter rate corresponds to $< 20$ such encounters. Even when ignoring the two-body relaxation processes that would “spin-down" the cluster after each successive encounter, we conclude that a long-term “spin-up" of OCs in the Solar neighborhood via a random walk of GMC–cluster encounters is not feasible; this matches the population statistics of rotating OCs, where rotators older than $1 ~\mathrm{Gyr}$ are uncommon [e.g., 46, 49, 63]. Note that this order of magnitude calculation applies only to spherical clusters; clusters with significant rotation may be flattened [e.g., 77, 34], for which $\Delta v_{\rm rot}$ may be significantly larger than $\sim 0.01 ~\mathrm{km ~s^{-1}}$, though testing this effect is outside of the scope of this work.

In this work, we make no attempt to explain the population of rotating OCs, as the population is mostly consistent with the “spin-down" of primordial rotation via coupling of two-body relaxation effects with the Galactic tides [63]. We do, however, seek to provide insight into the few anomalously rotating dynamically old OC detections that do exist, such as M67. [16] report that M67 has a rotational velocity of $0.81 \pm 0.07 ~\mathrm{km ~s^{-1}}$. This value is well above the predicted value of rotation that a fully relaxed cluster should sit at due to the Galactic tides of $\approx 0.1 ~\mathrm{km ~s^{-1}}$, following arguments by [39, 19]. This rotational velocity is also far too large to have come from primordial origin alone. [84] find that the peak of the rotation curve for a cluster at some epoch relative to its initial rotation is strongly related to the mass loss history of the cluster. Had M67 lost roughly $90\%$ of its initial mass [as predicted by the evolutionary models of 22], its initial $V_{\rm peak}$ would have needed to be $> 100 ~\mathrm{km ~s^{-1}}$ to correspond to its current $V_{\rm peak} \approx 1 ~\mathrm{km ~s^{-1}}$, which is non-physical. Other external dynamical effects that could affect the cluster’s angular momentum, like disk shocking or a binary cluster merger, both seem implausible, as disk shocks are efficient in isotropizing orbits [83], and any merger event would have needed to be within the first $\sim 100 ~\mathrm{Myr}$ of M67’s lifetime [82].

While our simulations show that inducing rotation in clusters via the integrated effect of GMC flyby encounters is infeasible, a cluster that recently underwent a sufficiently energetic interaction can appear as if it has net rotation within a predefined tidal radius. We can define a post-encounter timescale of $\Delta t = 5 ~\mathrm{Myr}$, and compare this with the timescale for a disruptive GMC–cluster encounter $t_{\rm dis}^{\rm single}$, using eq. 28 in [11]. We estimate $t_{\rm dis}^{\rm single} \approx 900 ~\mathrm{Myr}$ for a cluster of mass $500 ~\mathrm{M_{\odot}}$ and half-mass radius $3 ~\mathrm{pc}$, using a GMC mass density of $\rho = n \langle M_{\rm GMC} \rangle \approx 1.5 \times 10^{-2} ~\mathrm{M_{\odot} ~pc^{-3}}$, where all other values in eq. 28 are described by [11]. The fraction of open clusters that have undergone such a destructive interaction within the past $5 ~\mathrm{Myr}$ is then $\Delta t / t_{\rm dis}^{\rm single} \approx 0.006$. Comparing this fraction with the $\approx 3500$ observed OCs within the catalog from [17], we find that $\approx 20$ of these observed clusters have had a disruptive GMC encounter within the last $5 ~\mathrm{Myr}$; a small, but non-negligible amount. We estimate that roughly $7$, $8$, and $84$ such recently disrupted clusters exist within the OC catalogs produced by [75], [76], and [78], respectively. [17] also derive an estimated number of OCs within the Milky Way as a whole of $1.3 \times 10^5$; scaling $\Delta t / t_{\rm dis}^{\rm single}$ with this value gives $\approx 750$ very recently disrupted OCs in the Milky Way. We emphasize that these estimates are order of magnitude only, as the GMC density depends strongly on Galactocentric radius [e.g., 18, 66] and position relative to the Galactic spiral arms [e.g., 80, 65].

While M67 sits $\approx 400 ~\mathrm{pc}$ above the disk [67, 68], making a very recent GMC–cluster encounter unlikely, there is tentative evidence for some dynamical interaction. M67 likely has tidal tails [62, 69, 16], and there may be misalignment with the principle component of M67’s tidal tails and its orbital velocity vector [see Figure 2 in 16]. [70] show that old clusters on circular orbits should have tidal tails well aligned with their orbital velocity vector, and [79] show that GMC–cluster encounters can produce scatter in the relative inclination angle of a cluster’s tidal tails. We make no definitive claim of such an interaction, and only posit its plausibility as an explanation for detections of anomalously rotating, dynamically old OCs, similar to M67.

The contents of Gaia DR4, scheduled for release in December 2026, will provide extraordinary insight into the population of rotating OCs. Characterizing bulk rotation at a significant level (i.e., $> 1 \sigma$) in OCs is a challenging endeavor; one needs an accurate membership catalog, precise astrometric/radial velocity information, and enough member stars to provide a statistically significant sample. Only $\approx 26\%$ of the rotating OCs reported by [63] are found to rotate at the $1\sigma$ level, for example. The precisions of parallaxes, proper motions, and radial velocities in Gaia DR4 will improve by a factor of $\sqrt{2}$, $2 \sqrt{2}$, and $\sqrt{2}$, respectively, when compared with Gaia DR3 [61]. Moreover, the number of sources with RVS radial velocities will increase from $33$ million to over $100$ million [71], allowing for more clusters to be studied in three velocity dimensions [e.g., 43, 72]. These aforementioned data quality improvements will help to resolve differences between existing catalogs of rotating OCs [e.g., 53, 46, 49, 63], confirm or deny rotation in candidates with relatively large uncertainties, eliminate non-rotators with misattributed rotating signatures (e.g., apparent rotation due to the non-isotropic velocities of extra-tidal members), discover new rotation candidates, and shed light on the internal kinematics and dynamical histories of old rotation candidates like M67 and Ruprecht 147.

We note that the model and simulations presented in this work do not come without limitations. Our model for the approximate change in angular momentum due to a GMC–cluster encounter applies only for initially spherical clusters, and the internal energy gain of the cluster $\Delta E$ shows larger variations when compared with the analytic expectation from [11] at larger impact parameters. However, these regimes (small $M_{\rm GMC}$, large $b$) coincide with encounters that should produce more negligible changes to the cluster’s angular momentum, while both our model and simulations perform well at higher $M_{\rm GMC}$ and lower $b$. We make several simplifying assumptions in our simulations that we do not expect to alter the qualitative aspects of our results; we ignore primordial binaries, stellar evolution, Galactic tides, and the size of the GMC. Nevertheless, we do not expect any of these assumptions to change our conclusions in this work. We also emphasize that while this work is tailored towards encounters between GMCs and OCs, our model is directly applicable to encounters between dark matter subhaloes and GCs.

Conclusion

Inspired particularly by the large rotational velocity of the dynamically old OC M67, along with the pending release of Gaia DR4, we explore in this work the evolution of angular momentum, and thus rotation, of a simulated OC that underwent a nearby encounter with a point mass GMC. We derive an analytic expectation for the change in angular momentum induced by an encounter between a GMC and a star cluster, verify this expectation numerically using the $N$-body codes REBOUND and REBOUNDx, and explore the post-encounter rotational properties of the remaining subset of bound stars. Our main findings can be summarized as follows:

• The change in angular momentum of all stars in a cluster of mass $M_c = 1000 ~\mathrm{M_{\odot}}$ and $r_h = 3 ~\mathrm{pc}$ due to an impulsive GMC–cluster encounter, whether or not the stars remain bound to the cluster, is well described by Equation (5) in the tested regime of encounter parameters, i.e., impact parameter at infinity $b = 5, 10, 20, 40, 60, 80, 100 ~\mathrm{pc}$, relative velocity at infinity $V = 20 ~\mathrm{km ~s^{-1}}$, and GMC mass $M_{\rm GMC} = 10^4, 10^5 ~\mathrm{M_{\odot}}$.

• When allowing ample time for stars unbound by the encounter to escape the system, the rotational velocity of the bound subset of stars is unchanged when compared with its original pre-encounter value. Any changes to the rotational velocity of the bound subset post-encounter is consistent with the control simulations, where stars ejected from the system via two-body encounters may carry away angular momentum.

• Because these encounters produce only very small changes to the rotational velocity in the bound subset of stars post-encounters, we rule out the possibility of a long-term “spin-up" of a spherical cluster via weak perturbations from passing GMCs for clusters in the Solar neighborhood, as the necessary number of encounters to produce such a “spin-up" is orders of magnitude larger than the expected number of GMC–cluster encounters given the local number density of GMCs.

• When determining boundedness purely via a radial criterion, as opposed to an energy one, a cluster that has undergone a very recent encounter can appear as if it has significant rotation, as stars that have not yet had time to leave the system on the cluster outskirts can carry net-rotation.

• A small, but non-negligible number of OCs ($\approx 0.6\%$) have had a very recent disruptive encounter with a GMC. We estimate that roughly $20$ of these exist within the OC catalog produced by [17], and roughly $750$ of these very recently disrupted systems exist within the Milky Way at large. Such systems, in the throes of their disruption, may briefly be confused as rotating bound clusters.

• The rotation reported in M67 by [16] is truly anomalous; a primordial origin, tidally induced rotation, disk-shocking, and a binary cluster merger, all are implausible. We offer a recent perturbative encounter with a GMC as a potential explanation for M67’s apparent rotation, but make no definitive claims. Gaia DR4 will help to confirm or deny the existence of rotation within M67’s tidal boundary.

The rotational properties of open star clusters are still relatively understudied. Only very recently have efforts to simulate [e.g., 20, 21, 63] and identify [e.g., 53, 48, 46, 47, 49, 63] rotating OCs been undertaken. To our knowledge, we present the first detailed study on an origin of cluster rotation that is not primordial. At the advent of the release of Gaia DR4, studies like this one are necessary to provide insights into the underlying dynamical properties and histories of the soon-to-be-growing population of rotating open star clusters.

Acknowledgment

We thank Enrico Vesperini, Anna Lisa Varri, Carl Rodriguez, Ilya Mandel, Fred Rasio, and Brian Cook for useful discussions.

This research was supported by the Institute for Fundamental Theory scholarship through the Department of Physics at the University of Florida. The authors acknowledge UFIT Research Computing for providing computation resources and support that have contributed to the research results reported in this publication (https://it.ufl.edu/rc).

This work has made use of data from the European Space Agency (ESA) mission Gaia(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

K.T. used ChatGPT (OpenAI) to assist in developing the simulation code, and Claude (Anthropic) to assist with intermediate steps in the analytic derivations, plotting code, and editing of the manuscript text. All AI-assisted content was reviewed and verified by K.T., and the authors take full responsibility for the work.

Derivation of the Change in Angular Momentum

Consider a star with mass $m_i$ in a star cluster whose center of mass lies at the origin, undergoing an encounter with a GMC of mass $M_{\rm GMC}$. We work in the frame of the cluster center of mass, and write the star’s position as $\mathbf{r}_i(t) = \mathbf{r}_i + \delta \mathbf{r}_i(t)$, where $\mathbf{r}_i$ is its position before the encounter and $\delta \mathbf{r}_i(t)$ is its displacement relative to the cluster center of mass due purely to the GMC passage. In the impulse approximation, $\mathbf{r}_i$ may be treated as a constant. The acceleration of the star may then be written as $\mathbf{a}_i (t) = \ddot{\delta \mathbf{r}}_i (t)$, the explicit form of which we give in Equation (20). The torque this star experiences is $\boldsymbol{\tau}_i (t) = \mathbf{r}_i (t) \times m_i \ddot{\delta \mathbf{r}}_i (t)$, such that its change in angular momentum is $\Delta \mathbf{L}_i = m_i\int_{-\infty}^{\infty} \mathbf{r}_i (t) \times \ddot{\delta \mathbf{r}}_i (t) \, dt$. Expanding this integral gives

The first integral in Equation (7) is simply the final velocity kick due to the GMC passage, which may be written as

where $\lim_{t\rightarrow -\infty} \dot{\delta\mathbf{r}}_i (t) = 0$. The second integral may be evaluated by parts as

where the second integrand in Equation (9) is zero. At the lower limit, both $\delta\mathbf{r}_i (t)$ and $\dot{\delta\mathbf{r}}_i (t)$ vanish. The upper limit may be written by noting that, at late times, the star drifts freely along its final velocity kick, such that $\delta \mathbf{r}_i(t) \rightarrow \Delta \mathbf{r}_i + \Delta \mathbf{v}_i\, t$, where we define the displacement of the star due to the encounter as

The upper boundary of the first term in Equation (9) may then be written as

since $\lim_{t \rightarrow \infty} \dot{\delta \mathbf{r}}_i (t) = \Delta \mathbf{v}_i$. Equation (7) then simplifies to

To define the angular momentum change of the entire star cluster, we simply sum over all $i$ stars, such that

where, since $\mathbf{a}_i (t)$ is measured relative to the accelerating cluster center of mass, each star also experiences a torque from the uniform acceleration of the center of mass; this vanishes in the sum since $\sum_i m_i \mathbf{r}_i = 0$ and, to the order considered here, $\sum_i m_i\, \delta \mathbf{r}_i = 0$. For simplicity, we write $\Delta \mathbf{L} = \Delta \mathbf{L}^{(v)} + \Delta \mathbf{L}^{(r)}$, where $\Delta \mathbf{L}^{(v)} \equiv \sum_i m_i \mathbf{r}_i \times \Delta \mathbf{v}_i$ and $\Delta \mathbf{L}^{(r)} \equiv \sum_i m_i \Delta \mathbf{r}_i \times \Delta \mathbf{v}_i$, and derive each term individually.

To derive $\Delta \mathbf{L}^{(v)}$, we note that [12] derive the final velocity kick of some mass $M$, which we set equal to $M_{\rm GMC}$, due to a hyperbolic encounter with a star of mass $m_i$, positioned at $\mathbf{r}_i = (x_i, y_i, z_i)$, in the center of mass frame, as

where $V$ is the relative velocity at infinity, $b_i$ is the impact parameter with respect to the $i$-th, $b_{90} \equiv G M_{\rm GMC} / V^2$, and $\hat{\mathbf{V}}$ and $\hat{\mathbf{b}}_i$ are the directions of the relative velocity and impact parameter at infinity, respectively. Note that $b_{90}$ in its exact form features the sum of the two masses involved in the encounter, but we consider the GMC to be sufficiently more massive than any individual stellar mass $m_i$ or the cluster mass $M_c$, such that $b_{90} = G M_{\rm GMC} / V^2$ is a reasonable approximation. In this work, we let $\mathbf{V} = (0, 0, V)$ and $\mathbf{b}_i = (-x_i, b-y_i, 0)$. We approximate the GMC’s trajectory as a straight line, unperturbed by the cluster’s gravity. Since we ignore any motion of the star due to the cluster potential, $M_{\rm GMC} \Delta \mathbf{v}_{M_{\rm GMC}} + m_i \Delta \tilde{\mathbf{v}}_i = 0$. Approximating $m_i + M_{\rm GMC} \approx M_{\rm GMC}$, we may write the final velocity kick of the star $\Delta \tilde{\mathbf{v}}_i$ in Cartesian coordinates as

Defining the center of mass to be a particle of mass $M_c$, centered at the origin, we can similarly write the velocity kick of the center of mass to be

We may write the final velocity kick of the $i$-th star relative to the cluster center of mass, i.e., $\Delta \mathbf{v}_i \equiv \Delta \tilde{\mathbf{v}}_i -\Delta \tilde{\mathbf{v}}_{\rm com}$, as

Assuming the cluster follows a spherically symmetric Plummer density profile, all but one term in $\Delta \mathbf{L}^{(v)}$ vanishes under summation; the $x$-component of $\Delta \mathbf{L}^{(v)}$ is $\Delta L_x^{(v)} = \sum_i m_i (y_i \Delta v_{z,i} - z_i \Delta v_{y,i})$, where the first term may be non-zero since $y_i \Delta v_{z,i}$ is even in $x_i$ and $z_i$ and has no definite parity in $y_i$ due to the GMC’s passage at $y=b$. $\Delta \mathbf{L}^{(v)}$ can be written as

where the center of mass term vanishes since $\sum_i m_i y_i = 0$. For a large number of particles $N$, we may treat the sum in Equation (18) as an integral, such that

where $b_{\rm eff} = \sqrt{x^2 + (b - y)^2}$, $\rho(r)$ is the Plummer density profile $\rho(r) = (3 M_c / 4\pi a^3 ) \left[ 1 + r^2 / a^2 \right]^{-5/2}$, and $a$ is the Plummer radius, related to the half-mass radius by $r_h = 1.305 a$. Note that in the weak-encounter limit ($b_{\rm eff} \gg b_{90}$), Equation (19) is $\mathcal{O}(b_{90}^2/b_{\rm eff}^2)$, since $[1 + (b_{\rm eff}/b_{90})^2]^{-1} \simeq b_{90}^2 / b_{\rm eff}^2$.

Equation (17) gives only the velocity kick of the $i$-th star long after the encounter ($t \rightarrow \infty$), so to obtain $\Delta \mathbf{r}_i$ (and thus, $\Delta \mathbf{L}^{(r)}$), we instead follow the star’s response to the GMC over the duration of the encounter. The acceleration of the $i$-th star relative to the cluster center of mass, introduced above, is

where $\mathbf{X}(t) = (0, b, Vt)$ is the position of the GMC relative to the origin at time $t$, such that the closest approach of the GMC to the cluster center is at $t=0$. While integrating Equation (20), we hold $\mathbf{r}_i$ fixed, neglecting the effect of the star’s displacement due to the GMC passage on the acceleration it experiences. Writing the star’s position during the encounter as $\mathbf{r}_i + \delta \mathbf{r}_i(t)$ and Taylor expanding Equation (20) about $\mathbf{r}_i$, the relative acceleration of the $i$-th star may be written as

Since $\mathbf{a}_i (t) \propto M_{\rm GMC}$ and $\delta \mathbf{r}_i (t)$ is itself $\mathcal{O}(b_{90}/b_{\rm eff})$, the second term is $\mathcal{O}(b_{90}^2/b_{\rm eff}^2)$. Holding $\mathbf{r}_i$ fixed retains only the first term, such that the neglected terms contribute to $\Delta \mathbf{L}^{(r)}$ only at $\mathcal{O}(b_{90}^3/b_{\rm eff}^3)$, negligible compared with the $\mathcal{O}(b_{90}^2/b_{\rm eff}^2)$ terms retained in $\Delta \mathbf{L}$. We note that the fixed-$\mathbf{r}_i$ treatment is only used to obtain $\Delta \mathbf{r}_i$. For $\Delta \mathbf{L}^{(v)}$ we use the exact final two-body velocity kick (Equation (17)), since $\Delta v_{z,i}$, which alone contributes to Equation (19), first appears at $\mathcal{O}(b_{90}^2/b_{\rm eff}^2)$. In $\Delta \mathbf{L}^{(r)}$, by contrast, only the $\mathcal{O}(b_{90}/b_{\rm eff})$ part of $\Delta \mathbf{v}_i$ contributes since $\Delta \mathbf{r}_i$ is also $\mathcal{O}(b_{90}/b_{\rm eff})$, for which Equation (17) and Equation (22) agree.

Holding $\mathbf{r}_i$ fixed, we may write the separation between the GMC and $i$-th star as $\mathbf{X}(t) - \mathbf{r}_i = (-x_i, b-y_i, Vt - z_i)$ with magnitude $R_{i}(t) \equiv \left[ b_i^2 + (Vt - z_i)^2 \right]^{1/2}$, and similarly $R(t) \equiv (b^2 + V^2 t^2)^{1/2}$. Using the substitution $u = Vt'-z_i$, integrating Equation (20) with $\mathbf{r}_i$ fixed from $t' = - \infty$ to $t' = t$ yields the $i$-th star’s velocity kick at time $t$, $\Delta \mathbf{v}_i(t) \equiv \dot{\delta \mathbf{r}}_i(t)$, given by

where the components follow from the anti-derivatives of $(b_i^2 + u^2)^{-3/2}$ and $u\,(b_i^2 + u^2)^{-3/2}$, giving the perpendicular ($x, y$) and parallel ($z$) components of Equation (22), respectively. The center of mass terms follow identically with $z_i \rightarrow 0$ and $b_i \rightarrow b$. Note that expanding to first order in $b_{90}/b_i$ and $b_{90}/b$ in Equation (17) is equivalent to Equation (22) under the limit of $t \rightarrow \infty$.

Under our definition of $\Delta \mathbf{r}_i$ in Equation (10), we integrate Equation (22) from $t' = -\infty$ to $t' = t$, subtract the free drift $\Delta \mathbf{v}_i\, t$, and take $t \rightarrow \infty$, giving

where we again use the substitution $u = Vt' - z_i$ and utilize the anti-derivatives of $1 + u\,(b_i^2 + u^2)^{-1/2}$ for the $x$ and $y$ components of $\Delta \mathbf{r}_i$, and $(b_i^2 + u^2)^{-1/2}$ for the $z$ component.

The spherical symmetry of the system again causes all but one term in $\Delta \mathbf{L}^{(r)}$ to cancel due to odd parity in either $x_i$ or $z_i$. Specifically, the second term in $\Delta L_x^{(r)} = \sum_i m_i (\Delta y_i \Delta v_{z,i} - \Delta z_i \Delta v_{y,i})$ is non-zero, as this term is even in $x_i$ and $z_i$ and has no definite parity in $y_i$. We may write $\Delta \mathbf{L}^{(r)}$ as

where $\Delta v_{y,i} (x_i, y_i)$ is the $y$ component of Equation (17). Approximating Equation (24) as an integral over the density profile, and combining with Equation (19), we arrive at Equation (3), and write the change in internal angular momentum of the star cluster in the $x$-direction as