© The Authors 2025

1 Introduction

Star clusters with large stellar densities ≳100 M_⊙ pc⁻³ create very active dynamical environments where binary and single stars interact on short timescales of order a few megayears (e.g., Leigh & Sills 2011; Leigh & Geller 2012). Such star clusters evolve due to intracluster dynamics (Reipurth et al. 2010) and eject stars into the Galactic field via a combination of two-body relaxation and the disintegration of higher-order multiples (usually binaries) interacting (often chaotically) due to single-binary and binary-binary interactions. This binding and subsequent dissociation of temporary groupings of three or more stars occurs in all star clusters, including nuclear star clusters, globular clusters (GCs), and open clusters (OCs). Analogously, it is also possible for stable triple systems to become unstable and eject a binary star and a single star. We expect these events to be relatively rare.

Specifically, only 2–4% of all low- and intermediate-mass triples become unstable, and of those, 42–45% lead to an unbound escaper, while the rest lead to a merger event (Toonen et al. 2022). For especially young clusters, stars can also be ejected due to supernova (SN) explosions in binaries (Perets & Subr 2012), which act to destroy the SN progenitor and launch its binary companion at the instantaneous orbital velocity.

Fast stars ejected from a cluster due to either of the above mechanisms (i.e., excluding two-body relaxation and hence tidal stripping) are known as runaway stars (RSs). For the purposes of this paper, we call any object ejected from a cluster due to either the disintegration of triples (formed primarily from single-binary interactions) or the SN-induced ejection scenario an RS object. Typical velocities of RSs escaping from OCs must be sufficiently high to identify them as RSs, but also low enough that they remain sufficiently close to their host cluster to have a non-negligible probability of being associated with it. Observationally, RSs are typically defined as stars with velocities above 30 km s⁻¹ (e.g., Blaauw 1961; Renzo et al. 2019) relative to the host cluster center-of-mass velocity, which is a factor of a few greater than the typical GC escape velocity. But they can have lower velocities as well (i.e., corresponding to ejection velocities roughly equivalent to the escape velocity), in which case they are typically referred to as walkaway stars (Schoettler et al. 2020).

Fast-moving stars have been the subject of extensive studies using observations and simulations. There have been many studies on dynamically ejected massive stars from young star clusters using numerical integrations (Perets & Subr 2012; Fujii & Portegies Zwart 2014; Oh et al. 2015; Andersson et al. 2020; Dall’Amico et al. 2021). In most cases, the studies were performed with few-body scattering experiments (e.g., Oh et al. 2015; Pflamm-Altenburg & Kroupa 2006; Gvaramadze & Gualandris 2010; Gualandris et al. 2004; Ryu et al. 2017c; Ryu et al. 2017a; Ryu et al. 2017b) and focused on young star clusters. Interestingly, however, the contribution of RSs coming from old star clusters relative to the total population of observed fast-moving stars in our Galaxy has been argued to be non-negligible, mostly due to the ejection of singles during three-body interactions (e.g., Weatherford et al. 2023; Grondin et al. 2023; Grondin et al. 2024b). Observationally, studies of RSs were initially limited to bright massive stars. However, this changed with the launch of the Gaia mission (Gaia Collaboration 2023), which has provided precise measurements of positions and velocities for more than 1.5 billion stars in the Milky Way. As a result, numerous new RSs have been discovered (e.g., Carretero-Castrillo et al. 2023; Liao et al. 2023; Li et al. 2023; Igoshev et al. 2022), including ones at the low-mass end, highlighting the fact that they exist across all stellar masses.

The first scenario for producing RSs involves an SN in a binary system (Blaauw 1961). In a single degenerate system (i.e., a compact object plus a non-degenerate companion), the primary star will explode, with the consequent mass loss causing the companion to be launched at the instantaneous orbital velocity (e.g., Tauris 2014; Renzo et al. 2019; Leigh et al. 2020). In a double degenerate system (i.e., those containing two white dwarfs), a dynamically driven double-degenerate double-detonation (D⁶) scenario could occur (Shen et al. 2018b,a). Here, unstable mass transfer leads to the coalescence of the binary, yielding a Type Ia SN. However in the D⁶ scenario, instead of being fully destroyed during the SN, the companion WD survives and becomes a hypervelocity star ejected at a velocity greater than 1000 km s⁻¹. Recently, there have been large efforts to identify hypervelocity stars produced from the D⁶ scenario (Shen et al. 2018a; El-Badry et al. 2023). Associating these stars with stellar clusters is challenging, however, as the extreme speeds would cause these stars to travel far from their birth clusters.

The second mechanism for RS production thought to operate in old star clusters is the decay of gravitationally bound systems of three or four stars. There are two ways to create an RS star via a three-(or four-) body disintegration. In the first scenario, a hierarchical triple (or quadruple) is born stable but later becomes unstable due to stellar or binary evolution or due to even internal secular dynamical processes (e.g., Lidov-Kozai cycles) (e.g., Leigh et al. 2020; Toonen 2021; Hamers 2021). In the second scenario, three (or four) bodies meet and become gravitationally bound temporarily via a single-binary (or binary-binary) interaction. In both scenarios, the triple disintegrates, launching one RS and leaving behind a runaway binary (RB) with a recoil velocity decided by linear momentum conservation (Valtonen & Karttunen 2006; Leigh & Wegsman 2018). Hence, given an observed RS and using the available observational data, it should be possible to compute exactly where the RB is if it formed from the disintegration of a triple (e.g., Reipurth et al. 2010). It is then possible to search for the RB in large stellar catalogs. Throughout this work, we focus on triple disintegrations formed via single-binary interactions, which dominate over other types of interactions in dense clusters (Leigh & Sills 2011).

In this paper, we propose a novel observational methodology for the identification of RS-RB pairs formed from the decay of three-body systems due to single-binary interactions in a star cluster. As we show, it is these RSs that should be identifiable as being associated with a progenitor cluster with the highest confidence relative to other methods that are only able to work with a single RS (e.g., the binary SN-based RS formation mechanism).

We leverage the extensive Gaia database to identify the most probable RS-RB pairs. We focus in this paper on the old OC M67 (NGC 2682) as a case study, given its close distance to the Sun, moderate age, and dynamically active environment. Located at about 889 pc from the Sun (Cantat-Gaudin et al. 2020), the old OC M67 is made up of stars of solar age and composition with an isochronal age of ∼4.0 Gyr (VandenBerg & Stetson 2004; Balaguer-Nunez et al. 2007; Viani & Basu 2017). Proper-motion and radial-velocity surveys of M67 (Yadav et al. 2008; Pasquini et al. 2011; Geller et al. 2015) have confirmed kinematic cluster memberships and estimated the population size to be on the order of 1000. Studies have shown that close stellar encounters involving binary stars should be frequent (i.e., ≳ 1 Myr⁻¹; see the next section for a more detailed calculation) in M67 (Leonard 1996; Hurley et al. 2005; Leigh & Sills 2011) such that we expect to be able to find RS-RB pairs in the field of view surrounding the cluster (e.g., Leigh et al. 2016).

A method for finding RS-RB pairs and identifying their progenitor cluster has several important applications. First, in principle, detected RS-RB pairs can be used to say something about what must be going on dynamically inside the cluster. For example, given enough RS-RB pairs, particularly when applying the method to additional clusters, if the current orbital separation of the RBs can be measured, then energy conservation can be applied to constrain the initial binary orbital separation. In turn, it could be possible to place constraints on the underlying distribution of binary orbital separations, for example, by addressing whether or not the observed RS-RB properties are anomalous or expected for a given assumed distribution of binary orbital separations.

Runaway binaries also carry essential information about binary evolution. Three-body interactions can harden or soften the binary orbits; therefore, the properties of an ejected binary (i.e., orbital period, separation, and eccentricity) may be different from its state when it was in the cluster (Leigh et al. 2024). In the case of compact-object binaries, interactions that lead to compact separations may provide information about the formation channels of transient phenomena, such as Type Ia SN in the case of a white dwarf, or gravitational wave progenitors involving neutron stars or BHs. While many of those events are commonly associated with binary evolution pathways, such as the common envelope channel for compact binary formation (Paczynski 1976; Grondin et al. 2024a), the identification of dynamically formed short-period binaries would yield complementary evidence that could be used to disentangle which dynamical or evolutionary processes provide the main contribution to such compact systems. These results point to a broader perspective of the RBs as key tools for probing cluster dynamics, binary evolution, and the interaction between them. That is, we can potentially identify events commonly associated with pure binary evolution but that may have some clear contribution from dynamics, thereby facilitating our understanding of how typical exotica may have formed and to what extent their origins are owed to dynamics.

We also wish to test our theoretical understanding of three-body disintegration and hence the chaotic nature of the three-body problem. That is, we theoretically know the outcome distributions (see Valtonen & Karttunen 2006; Stone & Leigh 2019, etc.), but we have never been able to test if these are consistent with data. To determine this, we can ideally generate a larger sample size by looking at other clusters. Such a study becomes most interesting if we are looking for specific outcomes, such as those producing X-ray binaries. The rates will depend sensitively on the number of compact objects in a cluster, how many of these objects are in binaries, and so on. Therefore, the number of expected pairs depends on quantities that are not always directly observable but can be constrained using observations of runaway single-binary pairs. Hence, we obtain an independent measure of the underlying frequency and properties of the compact multiple star systems in a cluster relative to, for example, radial velocity surveys.

Following this initial study, a potential future application lies in the detailed study of black hole (BH)-main sequence (MS) binaries ejected from clusters. We can investigate the frequency of BHs in GCs since the predicted ejection rates depend on how many BHs are in a given cluster. Hence, this offers an independent means of constraining the number of BHs in clusters relative to X-ray observations of accreting BH binaries or gravitational wave observations. Once ejected from a dense GC, a BH-MS binary can be studied further using radial velocity observations, whereas this is not typically the case when inside a dense crowded cluster since radial velocity surveys cannot be performed due to crowding. This approach will allow us to verify the properties of the identified binary systems and potentially uncover systems that may contain BHs.

Another application of our method for detecting RS-RB pairs is in the study of high-velocity stars, many of which are observed in the Galaxy but do not have trajectories that point back toward the Galactic center. One hypothesis to explain these stars is that they originate from three-body ejection events in GCs. By analyzing the frequency of these interactions, we can determine whether all such high-velocity stars are consistent with having come from dense GCs, as opposed to the Galactic center. In future work, we aim to search for BH-MS binaries ejected by single MS stars. Even without population statistics, we can compare the observed properties of these binaries with theoretical distributions from three-body interactions to assess their rarity. By including more GCs in our sample, we can further strengthen the connection between theory and observation and address the origin of all high-velocity stars in the Galaxy.

In Sect. 2, we present and discuss the relevant physical concepts needed to understand how and why our selection criteria are constructed and the software, Corespray, used to compute the predicted outcome distributions. In Sect. 3, we present our methodology for identifying the most probable RS-RB pairs. In Sect. 4, we present the results of applying our method to the old OC M67 using the Gaia data and identify the most probable pairs of RSs-RBs for further study (e.g., radial-velocity observations). Finally, in Sect. 5, we conclude and discuss improvements on the present analysis that can be implemented when applying our method to other clusters.

2 Theoretical expectations for three-body disintegrations

In this section we present the basic theoretical expectations for three-body disintegrations that go into deciding how our final selection criteria are constructed, in addition to the software Corespray that we use to obtain predicted outcome distributions from theory (e.g., ejection velocities, directions, etc.). Corespray provides the theoretical information we need to be able to directly compare our observational results to theory.

2.1 Causality and conservation of momentum and energy

In this section, we discuss in more detail the expected observed properties of RS-RB pairs ejected from clusters due to single-binary interactions and the underlying physical mechanisms that produce their properties. We show that for the triple disintegration scenario, a larger number of selection criteria can be applied to the data to constrain the origins of RS-RB pairs relative to scenarios where only a single RS is expected to be produced (e.g., the SN-based mechanism).

To quantify what we expect for typical velocities of RSs formed from the triple disintegration scenario, we can compute a theoretically expected distribution of escaper velocities. Let the velocity of the escaper in the center-of-mass coordinate system be v_s. Then, a typical escape velocity distribution for the single star in three-body disintegrations can be described by the formula (Valtonen & Karttunen 2006) $$f\left(v_{\mathrm{s}}\right)d{v}_{\mathrm{s}}=\frac{\left(3.5{\left|E_0\right|}^{7/2}{m}_{\mathrm{s}}M/m_{\mathrm{B}}\right)v_{\mathrm{s}}d{v}_{\mathrm{s}}}{{(\left|E_0\right|+\frac{1}{2}\left(m_{\mathrm{s}}M/m_{\mathrm{B}}\right)v_{\mathrm{s}}^2)}^{9/2}},$$(1)

which is suitable to the isotropic case. Here, the total energy of the three-body system is $$E_0=\frac{1}{2}m{\dot{r}}_{\mathrm{s}}^2+\frac{1}{2}\mathcal{M}{\dot{r}}^2-G\frac{m_{\mathrm{a}}{m}_{\mathrm{b}}}{r}-G\frac{m_{\mathrm{s}}{m}_{\mathrm{B}}}{r_{\mathrm{s}}},$$(2)

where the separation between the components of the binary system is r and r_s is the separation between the escaper and the binary center of mass. For the masses of the bodies we have: m_a and m_b are the masses of the stars in the binary system, m_B is the total binary mass (m_a + m_b), m_s is the mass of the single star escaper and M is the sum of all masses in the three-body system. We also have $\mathcal{M}=\frac{m_{\mathrm{a}}{m}_{\mathrm{b}}}{M}$ and $m=\frac{m_{\mathrm{B}}{m}_{\mathrm{s}}}{M}$, which correspond to the reduced masses of the relative motions of the binary and the third body, respectively. The peak of this distribution is obtained by differentiating Eq. (1) and setting it equal to zero: $${\left(v_{\mathrm{s}}\right)}_{\text{peak\hspace{0.17em}}}=\sqrt{\frac{2(M-m_{\mathrm{s}})}{5m_{\mathrm{s}}M}}\sqrt{\left|E_0\right|}.$$(3)

The most likely escape velocity assuming an initial binary separation of 100 AU (integrating over Öpik’s Law from twice the orbital separation corresponding to a contact state to twice the hard-soft boundary yields an average orbital separation of ∼ 168 AU) is of the order of a few to a few tens of kilometers per second. We note, however, that deriving a predicted velocity distribution requires assuming an underlying distribution of binary orbital parameters, both for the triple disintegration scenario as well as for the SN-based mechanism.

Given an assumed ejection velocity for an RS candidate coming from the triple disintegration mechanism, both timescales as well as momentum conservation-based arguments can be used to predict the exact position and velocity of the associated RB candidate. In particular, since momentum conservation should be upheld during the decay of triple systems, the ratio of the RS and RB velocities should be equal to the inverse of the ratio of their masses, or $$v_{\mathrm{R}\mathrm{B}}=\frac{m_{\mathrm{R}\mathrm{S}}}{m_{\mathrm{R}\mathrm{B}}}{v}_{\mathrm{R}\mathrm{S}},$$(4)

where v_RS/v_RB and m_RS/m_RB are, respectively, the RS-RB velocities and masses. Given a mass-luminosity ratio, the masses can be converted to observable luminosities (or magnitudes). Additionally, the candidate RS and its associated candidate RB should have equal travel times from the position where the velocity vectors intersect. We call this position the cross-point, and it should lie within the tidal radius of the host cluster. Finally, using energy conservation, it is in principle possible to also constrain the most likely initial binary orbital separation of the putative single-binary interaction, given observational constraints on the final binary orbital separation as well as all three stellar masses. We emphasize that none of the preceding constraints can be applied to RSs coming from the SN-based mechanism (or any other mechanism producing only a single RS) because they do not have an associated RB, which should in principle translate into an ability to more confidently associate an identified RS-RB pair candidate with a host (old) star cluster for the triple disintegration mechanism.

Few constraints on the orbital properties of the underlying binary population are available, particularly in old, massive clusters like GCs with crowded fields of view. We argue in this paper that by looking observationally for the products of single-binary interactions, limits can be placed on the distribution of binary orbital separations. It is possible to obtain these constraints because the mean timescale between single-binary interactions in star cluster cores (where most ejection events are thought to occur) is given by (Leigh & Sills 2011) $$\begin{array}{rl}{\tau }_{1+2}=3.4\times {10}^7{(1-f_b)}^{-1}{f}_b^{-1}{\left(\frac{1\mathrm{p}\mathrm{c}}{r_c}\right)}^3{\left(\frac{{10}^3{\mathrm{p}\mathrm{c}}^{-3}}{n_0}\right)}^2& \\ \left(\frac{v_{\mathrm{r}\mathrm{m}\mathrm{s}}}{5\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}\right)\left(\frac{m}{1{\mathrm{M}}_{\odot }}\right)\left(\frac{1\mathrm{A}\mathrm{U}}{a}\right)\mathrm{y}\mathrm{r}& \end{array}$$(5)

where f_b is the binary fraction, r_c is the core radius in parsecs, n₀ is the number density in the core in pc⁻³, v_rms is the root-mean-square velocity in kilometers per second, m is the average mass (in M_⊙) and a is the mean binary semi-major axis (in AU). We note that Eq. (5) suggests that the rate of single-binary interactions scales linearly with the binary semi-major axis, which is the most difficult quantity in this equation to measure observationally. It follows that, by constraining the frequency of RSs formed from single-binary interactions, limits can be placed on the underlying binary orbital separation distribution. Based purely on Eq. (5), it is specifically the mean semi-major axis that is constrained.

2.2 Corespray

To obtain more robust predictions, including outcome distributions for disintegrating triples, we use Corespray¹ to generate a distribution of RS-RB pair escaper velocities, as a function of our assumptions for the properties of the underlying binary population. Corespray is a PYTHON-based three-body particle spray code that simulates extra-tidal stars ejected from a star cluster’s core due to single-binary interactions (Grondin et al. 2023). All orbits in corespray are defined and integrated using galpy² a Python-package for galactic dynamics (Bovy 2015). We use KingPotential and MWPotential2014 from galpy to encompass influences from both the cluster and the galaxy. Corespray uses the formalism from Valtonen & Karttunen (2006) to generate the outcome distributions from three-body disintegrations. The code allows us to generate a large sample of ejected RS-RB pairs, from which we obtain the expected distribution of ejection velocities.

We simulate with Corespray a sample of 5000 pairs ejected randomly over a period of 10 Myr. The cluster M67 is modeled with a King profile for the gravitational potential (King 1962) and positioned in the Galaxy such that its center of mass and center-of-mass motion are in agreement with observational constraints (see Sect. 3.1.2 for the exact values). We assume a range of stellar masses between 0.1 M_⊙−1.4 M_⊙ adopting a Salpeter mass function (i.e., with a power-law of α = −1.35; Salpeter 1955). For the minimum and maximum orbital separations, or a_min and a_max respectively, we assume limits between twice the orbital separation corresponding to a contact state and twice the hard-soft boundary (assuming two solar mass stars). Using the data from Corespray we obtain the expected velocity distribution for the ejected pairs, as shown in Fig. 1.

Next, we consider how many RSs we expect to currently observe around M67 due to triple disintegration coming from single-binary interactions in the cluster core. Specifically, consider a circular field of view of a radius of 150 pc at the distance of M67. We chose this value for the radius initially to be a balance between computational expense (which increases with the number of candidate pairs) and the probability of finding an ejected RS-RB pair in the field of view. Moving at a tangential velocity of 30 km s⁻¹ (where Corespray predicts an average velocity of 29.4 km s⁻¹ for RSs and 6.44 km s⁻¹ for RBs; see Fig. 1), a star will cross this field-of-view in roughly 10 Myr, as shown from our Corespray simulations in the bottom right panel of Fig. 2. But the time between single-binary interactions in M67 is 55 Myr. These timescales are not directly correlated; however, comparing them reveals that, theoretically, the detection of one or more pairs is indeed expected. This comparison highlights the consistency between the theoretical framework and the anticipated observational outcomes, reinforcing the plausibility of such detections. To calculate the mean time between interactions, we used the same parameters as adopted in Sect. 3.1 of Leigh & Geller (2013). We adopted Öpik’s Law (Öpik 1924) for the initial orbital parameter distribution assumed in Corespray, and we integrated over this from twice the separation corresponding to a contact state to twice the hard-soft boundary yields an average orbital separation of ∼ 168 AU. Assuming all of these interactions resulted in an ejection velocity for the single star of ∼30 km s⁻¹ by a binary with a semi-major axis of 100 AU (i.e., close to the mean value of 168 AU for our assumed initial semi-major axis distribution), we therefore expected to observe on the order of two RS-RB pairs originating from the cluster in this field of view at the present time due solely to single-binary interactions. An observation of more or less than this amount could imply that some of our assumptions for our timescale calculation are incorrect, such as the mean binary orbital separation and/or the binary fraction. Based on the above calculation for the expected number of pairs in our chosen field of view, we expect few RS-RB pairs to come from individual clusters, but this limitation on statistical analyses can be somewhat mitigated by including other clusters in the sample.

To summarize this section, given the additional constraints available when searching for RS candidates coming from the triple disintegration scenario, we intuitively expect to be able to more confidently associate them with their host cluster relative to runaway-producing ejection mechanisms that produce only a single RS. Hence, the inclusion of the RB in the analysis increases the number of selection criteria we have to work with.

Fig. 1 Histogram showing the predicted velocity distributions of escapers ejected in the Corespray simulation. The x-axis shows the ejection velocity in kilometers per second.

Fig. 2 Four different snapshots in time of a simulation using Corespray. We show the spatial distribution of the ejected pairs in the plane of the sky at 0.1, 1, 10, and 100 Myr.

3 Methods

In this section, we develop a methodology to identify RS candidates from clusters coming from the disintegration of three-body systems formed from single-binary interactions, which we apply to the old OC M67 as a case study. To this end, we present our methodology for identifying not only candidate RSs, but also candidate RBs that are consistent with having interacted with them in M67. The methodology results in a (short) list of RS-RB pair candidates and their corresponding observed properties.

3.1 Data acquisition

Here we introduce the data used to search for RS-RB pairs in the vicinity of M67, as well as our preliminary cuts in distance and position on the sky (but defer more in-depth calculations to Appendix A).

3.1.1 Gaia DR3 data

Ideally, we must build a 6D phase space, requiring right ascension, declination, 2D proper motions, radial velocity and parallax (i.e., distance). In our case, we are able to construct a nearly complete 5D phase space using the Gaia DR3 catalog.

Excluding radial velocities quantifies how constraining the data are without the third dimension in velocity, which can be exceedingly difficult to obtain for many, especially more distant, clusters and especially for those with high densities and hence crowded fields of view, such as GCs. It is also important to note that the Gaia parallaxes are somewhat limited in providing a precise distance for distant clusters like M67. The median parallax uncertainties are 0.02–0.03 mas for G < 15, 0.07 mas at G = 17, 0.5 mas at G = 20, and 1.3 mas at G = 21 mag³. Hence, we adopted a course-grained cut on distance for our candidates, as described in the subsequent section.

This data release represents a major advance with respect to Gaia DR2 and Gaia EDR3 because of the unprecedented quantity, quality, and variety of astrophysical data. Even though Gaia DR3 includes the sixth data release from the Radial Velocity Experiment (RAVE) survey (Steinmetz et al. 2020), we do not include these data in our analysis. From the total number of sources in Gaia only ∼1.8% have a radial velocity value³.

3.1.2 Preliminary sample selection

In order to identify our initial sample of candidate RSs, we extract from the Gaia data (Gaia Collaboration 2023) a complete 5D kinematic data set (i.e., proper motions, parallax and position on the sky) that covers our chosen field of view centered on the motion of the center of M67 (as reported by Jacobson et al. 2011) with a projected radius of 150 pc (i.e., forming a circular field of view in 2D). This limit and hence our first cut in the data is decided, as calculated in Sect. 2, based on the computed time between single-binary interactions weighed against the added computational expense of increasing our field of view and hence the total number of candidate pairs. In the end, we expect order unity RS-RB pairs to be observed within our chosen field of view at any given time.

First, we apply a cut in distance. For this, we adopt a cylindrical volume with a diameter and a length of 300 pc (in the direction of the cluster), centered on M67. We select all stars that are inside this 3D cylinder. We do not incorporate the uncertainties in parallax for this preliminary cut, which can be substantial in some cases.

Within this cylinder, we identify 203 190 objects not factoring in the uncertainties in parallax. Including objects with parallax uncertainties that make them consistent with lying within the domain of our cylinder to within 1σ (one standard deviation) increases this number by about a factor of two. We do not include these additional objects in our initial sample since most of them are dim with large uncertainties in much of our 5D parameter space. It is the dimmest stars that tend to have the largest uncertainties, making it more likely that they will satisfy one of our criteria at the 1σ level. It is therefore arguably a better choice to omit these additional sources since we do not expect it to be possible in many cases to securely identify them as disrupted triples. We note that the change in sample size we obtain upon performing a similar cut in the 2D position on the sky (i.e., including all sources that lie within 1σ of falling inside our circular field of view) is negligible given the small uncertainties. Hence, we do not include these additional objects in our initial sample either.

For our second cut in the data, we do not consider objects located within the cluster in projection on the plane of the sky. Using for the cluster center α = 8^h51^m23^s.3, β = +11^∘49^′02^′′ (J2000) (Jacobson et al. 2011), we identify those RSs located outside the cluster tidal radius. In previous studies, the tidal radius was found to be either 14 ± 2.5 pc (Keenan et al. 1973), 16.8 pc (Davenport & Sandquist 2010) or 17.98 pc (Kharchenko et al. 2013). For simplicity and to choose a balance between these values, we adopt a tidal radius of 17 pc. For all of our calculations, we assume the following additional parameters for M67. We assume a total mass of the cluster M_cl = 2100 ± 600 M_⊙ (Geller et al. 2015), a core radius of r_c = 1 pc (Bonatto & Bica 2003), a half-light radius r_h = 2.6 pc (Balaguer Núñez 2006) and a central 1D velocity dispersion of 0.59 ± 0.07 km s⁻¹ (Geller et al. 2015).

Finally, we remove from our sample all objects with 2D velocities smaller than the escape velocity from M67 relative to the cluster center of mass. For this, we adopt an escape velocity of 2.6 km s⁻¹, calculated using the method adopted in Georgiev et al. (2009) by assuming a correction factor of f_c,t(1.2) = 0.09154 (corresponding to the velocity once the object reaches the tidal radius). We also consider only sources for which tracing back the velocity vector passes within the tidal radius of M67. This cut reduces our sample size from 203 190 to 15 346. As we show in Sect. 5, our chosen criteria are in some ways conservative since the inclusion of the third dimension in velocity (i.e., radial velocities) could yield velocity vectors that are inconsistent with having originated from M67.

3.2 Final sample selection

In this section, we select from our preliminary sample to identify the most likely candidate RS-RB pairs that were ejected from M67 due to single-binary interactions. We begin, as described in Sect. 2, with a total of 15346 objects, yielding in total 117864981 unique combinations.

3.2.1 Traceback position and time

First, we identify every pair of objects whose 2D velocity vectors trace back to intersect within the tidal radius of M67. This is done as follows. Using Eq. (A.1), we first transform to a coordinate system in velocity space that is centered on the motion of M67’s center, as reported by Jacobson et al. (2011). We then calculate in this frame the 2D proper-motion vector for every object in both magnitude and direction using the Gaia proper motions in both RA and Dec (van Leeuwen 2009; Gaia Collaboration 2018a; Getman et al. 2019; Kuhn et al. 2019, see Eqs. (A.2, A.3, and A.4)), along with the corresponding uncertainties.

Using the proper motions and their uncertainties in RA and Dec, we calculate three vectors for each object. The “primary” vectors are defined by the actual parameters in position and proper motion reported in the catalog for every pair. For illustrative purposes, we show two additional vectors defined by the uncertainties in position and proper motion as shown in Fig. 3, calculated using Eqs. (A.10), (A.12), and (A.11), yielding maximum and minimum distances in this parameter space that separate these points of intersection from that of the primary. With this, we then calculate three points of intersection, called cross-points, by applying the equations in Appendix A.3. We keep in our sample all pairs with at least one of these three cross-points lying within the tidal radius of M67 and discard the rest.

Next, all pairs of candidate objects with cross-points located within the domain of M67 must coincide in not only space but also time. For a chosen cross-point that is located within the M67 tidal radius, we calculated the time expected for both sources to get to that point tracing back in time. Hence, we calculated for both objects in every pair the traceback time (TBT) to the location of the cross-point, defined as the time taken to travel from the cross-point to the presently observed object position. This is done using the measured values in both distance traveled and velocity, assuming the presently observed velocity is constant and has not changed over the extent of its travel across the field of view. We then calculate an uncertainty by propagating the provided uncertainties in both position and velocity. Hence, the time it will take object a (and it is the same for object b) to reach the cross-point is given by $$TB{T}_{\mathrm{a}}=\frac{|{\mathit{r}}_{\mathrm{a}}-{\mathit{r}}_{\mathrm{c}\mathrm{p}}|}{\left|{\mathit{v}}_{\mathrm{a}}\right|}=\sqrt{\frac{{(r_{\mathrm{a}}^{\mathrm{x}}-r_{\mathrm{c}\mathrm{p}}^{\mathrm{x}})}^2+{(r_{\mathrm{a}}^{\mathrm{y}}-r_{\mathrm{c}\mathrm{p}}^{\mathrm{y}})}^2}{{\left(v_{\mathrm{a}}^{\mathrm{x}}\right)}^2+{\left(v_{\mathrm{a}}^{\mathrm{y}}\right)}^2}},$$(6) $$TB{T}_{\mathrm{b}}=\frac{|{\mathit{r}}_{\mathrm{b}}-{\mathit{r}}_{\mathrm{c}\mathrm{p}}|}{\left|{\mathit{v}}_{\mathrm{b}}\right|}=\sqrt{\frac{{(r_{\mathrm{b}}^{\mathrm{x}}-r_{\mathrm{c}\mathrm{p}}^{\mathrm{x}})}^2+{(r_{\mathrm{b}}^{\mathrm{y}}-r_{\mathrm{c}\mathrm{p}}^{\mathrm{y}})}^2}{{\left(v_{\mathrm{b}}^{\mathrm{x}}\right)}^2+{\left(v_{\mathrm{b}}^{\mathrm{y}}\right)}^2}},$$(7)

where r_a is the distance separating the object’s current position from the center of mass of M67 and v_a is the velocity relative to the center-of-mass motion of M67. Similarly, r_cp is the vector from the center of mass of M67 to the cross-point with $r_{\mathrm{c}\mathrm{p}}^{\mathrm{x}}$ and $r_{\mathrm{c}\mathrm{p}}^{\mathrm{y}}$ being the corresponding values for each coordinate. If the difference in TBT values for a given pair of objects is consistent with zero to within 1σ, we keep the pair in our sample and discard it otherwise.

Fig. 3 Diagram illustrating the selection method for the angle of intersection. We have the two limiting vectors described by the error in position and proper motion (pma ± Δ pma and pmb ± Δ pmb) for the RS and RB, which yield a minimum (lower) and maximum (upper) possible angle (using the 1σ uncertainties in both directions; see Eqs. (A.10), (A.11), and (A.12)). If a given pair has an angle of intersection consistent with 180∘ to within these upper and lower bounds, we keep it in our sample and discard it otherwise. We note that the RB is not resolved and the cross point is at the center of light.

3.2.2 Angle of intersection between velocity vectors

We further required that in the plane of the sky, the angle between the velocity vectors (to the hypothetical cross-point) of a given pair must be consistent with 180^∘, as illustrated in Fig. 3. To evaluate this consistency, we used the two limiting vectors described above when evaluating the possible range of cross-points for each pair, which yield a minimum (lower) and maximum (upper) possible angle (using the 1σ uncertainties in both directions; see Fig. 3). We calculate the angle between the intersection of each pair of vectors using Eqs. (A.10), (A.11), and (A.12). If a given pair has an angle of intersection consistent with 180^∘ to within these upper and lower bounds, we keep it in our sample and discard it otherwise.

3.2.3 Velocity ratio

We also required that both objects have velocities consistent with linear momentum conservation. This means that when assuming the slower moving object is a binary at least twice as massive as its single star counterpart (and that it was the lowest-mass star that was ejected from the three-body interaction as a single; Aarseth & Heggie 1976), the ratio of the object velocities should exceed two. We applied this cut to our data using the provided proper motions in RA and Dec by calculating a total magnitude for each object velocity in 2D (relative to the center-of-mass motion of M67), and then use these values to compute a ratio. We propagate the 1σ uncertainties in proper motions to compute an uncertainty for each velocity ratio. We then keep in our sample all pairs that have a velocity ratio greater than two by 1σ or more and discard them otherwise.

3.2.4 Position in the color-magnitude diagram of the progenitor cluster

Finally, we assess whether or not both objects in each pair have a location in the color-magnitude diagram (CMD) that is consistent with having originated from M67 using the Gaia DR3 G_BP and G_RP magnitudes. Specifically, we required that one object (the candidate single star) lie within 1σ of an isochrone for M67 in both color and brightness simultaneously (but we note that our results do not change significantly if we instead use both magnitudes instead of a color), while the other object (the candidate binary) lies within 1σ (in both color and brightness) of the area in color-brightness-space defined as the region between the isochrone and an identical track shifted upward by 0.75 mag (corresponding to the equal-mass binary sequence).

Additionally, it is the slower object that should be the brighter one from linear momentum conservation (assuming the brighter object is more massive) since the binary is most likely to include the most massive star as a dynamical outcome. To evaluate this, we use the PARSEC isochrone fit provided in Childs et al. (2024). We consider that the pair is a valid candidate if the fastest object is the dimmest and it simultaneously lies within 1σ of the isochrone in both magnitude and color. At the same time, the slower object is the brightest and is consistent with lying between the isochrone and the equal-mass binary sequence to within 1σ in both magnitude and color. If both objects in a given pair satisfy these criteria then we keep it in our sample and discard it otherwise.

To summarize, the criteria, we implement are as follows, independent of the order in which they appear:

• Traceback position: identifies candidate pairs whose velocity vectors, rewound back in time, intersect within the cluster. Ensures a spatial coincidence.

• Traceback time: checks if these candidates share the same travel time within errors to reach that intersection point. This approach ensures temporal coincidence.

• Angle of intersection: angle between their velocity vectors, considering uncertainties, is consistent with 180 degrees in the plane of the sky.

• Velocity ratio: only consider pairs with a velocity ratio exceeding two, with the brighter object moving more slowly. The velocity ratio criterion is because it is the least massive object that is the most likely to be ejected as a single, and it should be moving faster which follows causally from linear momentum conservation given our assumptions.

• Position in the CMD: examination of locations in a CMD and comparison with a model for the cluster. One object (the single star) should match the isochrone, while the other (the binary) should be brighter and lie within the region separating the isochrone from the equal-mass binary sequence.

4 Results

In this section, we first describe the results of our analysis, focusing on the fraction or number of pairs that simultaneously satisfy all of the criteria described in the previous section. Our initial preliminary sample size, described in Section 3.1.2, is 117 742 185.

In the end, we present the only pair that simultaneously satisfies all of our criteria, which is the pair with the smallest calculated uncertainties. The uncertainties in angle, TBT, and position in the CMD quickly become large when looking at the other eight candidates that pass all of our criteria (see Table 1).

For each criterion, however, we evaluated to what extent it reduces our sample size and, wherever possible, compare to the results of our Corespray simulations to assess how each part of the criteria and our method performs when tested with simulated data.

Fig. 4 Histograms showing the distribution of the difference in TBTs (first and second panel) as well as their ratio (third and fourth panel) for every pair in our observational data and a 2D velocity projection in corespray. The TBTs are provided in years.

4.1 Traceback position and time distributions

One of the cuts applied to the data ensures that the velocity vectors of a given pair of objects rewound in time (typically ≲1 Myr), intersect in 2D space within the tidal radius of M67. This cut alone reduces the sample size by 71%. We also required that these traceback vectors intersect not only in space but also in time, which represents a 90% reduction. If the TBTs are consistent within the 1σ uncertainties, then both objects in a given pair are taken as consistent with having left the cluster at the same time. The mean uncertainty on the TBTs is ± 0.201 Myr. These cuts combined reduce our initial sample size by 97.46%.

The distributions in both the ratio and differences and the ratios of TBTs for each pair are shown in Fig. 4. In Fig. 4, one can see that most TBTs are equal to within a factor of a few. We note that we do not compare to the Corespray results here since the difference in TBTs is always zero for these simulations.

The reason that these cuts only reduce our sample size by roughly an order of magnitude is that most objects have velocities that follow an overall Galactic flow and hence tend to point in similar directions with comparable magnitudes. Additionally, all objects lie a roughly similar distance away from the center of mass of M67 (i.e., to within a factor of a few given our field of view and the limiting tidal radius of M67). Given these limitations, it is not so surprising that a relatively large fraction of pairs are consistent with having originated from M67 while also having similar TBTs.

4.2 Distribution of angles of intersection between velocity vectors

We further required that a given pair have an angle at the crosspoint between the two velocity vectors that are consistent with 180^∘, as predicted from linear momentum conservation. To this end, we computed the angle of intersection for every pair of objects in our sample. The resulting distribution is shown in Fig. 5.

As is clear, the least probable angle observed is 180^∘, making this cut especially constraining. The typical uncertainty in this angle is highly dependent on how dim the two objects in a given pair are but is around a degree for pairs having both members brighter than G_bp < 20. Hence, most sources do not fall within 1σ of 180^∘. Specifically, this cut reduces our initial sample size by 99.78%.

For comparison, we also show the predicted distribution of intersection angles we obtain from Corespray in Fig. 5. Corespray predicts that the distribution should be much more clustered around 180^∘, given that these pairs are the result of a three-body disintegration. Hence, Corespray predicts that the more likely angle to be obtained is 180^∘.

4.3 Velocity ratio distribution

In order to keep a given pair in our sample, we required that the velocity ratio of a given pair exceed two by at least 1σ. The distribution of velocity ratios for every pair in our sample is shown in Fig. 6, where the velocities are computed with respect to the cluster center-of-mass motion. The vertical dashed line indicates the critical velocity ratio of two. As is clear, this cut is not especially constraining. It reduces our initial sample size by 69%.

Curiously, most of the pairs in the Corespray results show a significantly broader velocity ratio than is observed, extending to much higher ratios. Most of the pairs in the Corespray simulation tend to have velocity ratios >2, as theory predicts (given that Corespray only simulates the products of three-body disintegration). The surplus of velocity ratios >2 relative to the observed distribution could be due to three reasons. First, the effects of both the cluster and Galactic potential can significantly alter the ratio of velocities post-ejection, and Corespray adopts potentials that are not identical to the true underlying potentials. Second, in order to properly compare theory to observations, the Corepray simulations would need to include foreground contaminants, which tend to have comparable velocities. Third, Corespray provides 3D velocities which we compare to our 2D velocities projected onto the plane of the sky, and the third dimension can increase the velocity ratio. We come back to the velocity ratio criterion in Section 5.2, but for now, we note that using our proposed velocity cut by choosing only sources with velocity ratios greater than 2 is likely not a useful criterion to incorporate in future work without a more sophisticated calculation that accounts for sources of acceleration/deceleration (e.g., cluster and/or Galactic potentials). However, removing this cut from our criteria for identifying the most probable pairs does not change our results significantly.

Finally, since we have assumed that the least massive object is always ejected from a three-body interaction, we also required that for each pair, the slower moving object be the brighter one. If we apply this simple cut after the one made previously in this section, we further reduce our sample size by one order of magnitude from the initial sample. Hence, including a cut that combines the velocity ratio with the relative brightness of each object is much more constraining than a simple velocity cut alone.

Fig. 5 Histogram showing the distribution of intersection angles between the primary vector of every pair. The x-axis shows the angle in degrees. The top panel shows the projected angle for every pair from our Corespray simulations. The bottom panel shows the distribution with observational data. The simulation shows that most of the pairs tend to have angles near 180∘. In contrast, the pairs from the observational data set have a very low number of pairs near 180∘. This difference between the simulation (only pairs resulting from three-body disintegration) and observational (random sources in a field of view) distributions shows the angle to be one of the most important filters for the method.

4.4 Position in the color-magnitude diagram

Both objects should lie within 1σ (simultaneously in both color and brightness) of the area in color–brightness-space defined by the region between the isochrone and an identical track shifted upward by 0.75 mag (corresponding to the equal-mass binary sequence). Finally, we required that for a given pair, the fastest source (the candidate single) must lie within 1σ (in both color and brightness) of an isochrone for the progenitor cluster. Further, considering the fact that the fastest star must lie on the isochrone to within 1σ confidence, we eliminate approximately 98.54% of the total sample.

Fig. 6 Histogram showing the distribution of velocity ratios. The top panel represents every pair from a simulation using Corespray, and the bottom panel shows all the pairs in our observational data. The vertical dashed line shows a critical ratio of two. Theory predicts that it is very likely that RS-RBs have velocity ratios greater than two, as can be seen in the simulated data set. On the other hand, the observational data show that more than half of the pairs have velocity ratios of less than two.

4.5 Number of candidates remaining

After applying all of the above criteria, we are left with only one pair that satisfies them all simultaneously. This one pair is shown in Fig. 7, and summarized below. The small dots show high-probability cluster members obtained from Childs et al. (2024), whereas the two large circles correspond to the two components of our one candidate pair. The lines show an isochrone fit and the equal-mass binary sequence (obtained by shifting the isochrone up in brightness by 0.75 mag). The color coding indicated in the provided legend shows which of the two objects is the fastest, confirming that it is indeed the brighter object that is moving slower.

Therefore, after applying all of the criteria available to us from momentum conservation and causality, our total sample size is reduced by eight orders of magnitude, from ∼ 10^∘ to 1.

This resulting pair has the following properties. For the RS, the properties are ID_a: 599619181604349440, v_a = 66.9 km s⁻¹, g_bp,a = 18.36, and g_rp,a = 16.57. And for the RB, they are ID_b: 604856670884640384, v_b = 27.6 km s⁻¹, g_bp,b = 16.95, and g_rp,b = 15.46. Between these two sources, we have the following: v_ratio = 2.43, TBT_diff = 0.005 ± 0.007 [Myr], and ${\theta }_{a,b}^{\circ }=179.6\pm 0.5$. The RS candidate is consistent with lying on an isochrone in the CMD 7, while the RB is consistent with lying within the binary sequence just above the MS.

A few additional properties of the RS-RB pair can also be inferred. Given where the RS falls on the isochrone, we can infer that it has a mass of roughly 0.67 M_⊙. For the RB, we can think of the following limit. We assume that the vertical displacement of the binary from the isochrone could be explained if the primary component has a mass of 0.73 M_⊙ (i.e., as inferred from the vertical line connecting the RB’s CMD location to the isochrone, a magnitude for the primary of 17.5). We then took the magnitude required for the secondary in order to end up with a total magnitude equal to that of the binary, or M_bin ∼ 16.95. Using this approach, we compute a magnitude for the secondary of 17.95, corresponding to a mass of 0.7 M_⊙ and hence a mass ratio of q ∼ m_b/m_a ∼ 0.96, whereas before we assumed m_a and m_b are respectively the masses of the primary and secondary of the binary. Given these component masses, and the fact that the most likely semi-major axis is around 100 AU (see Section 2.2), the period of the binary would be 836.2 years. This corresponds to a low orbital velocity where detection with radial velocity is not feasible. However, at a distance of 1 kpc, 100 AU corresponds to 0.1 arcsec in the plane of the sky, which is resolvable from the ground.

Fig. 7 Color-magnitude diagram of the old OC M67 plotted using the Gaia data. Only stars indicated as high-probability cluster members in Childs et al. (2024) have been selected and are plotted as small blue dots. The solid red lines show the best-fitting isochrone (lower line) as well as the equal-mass binary sequence shifted 0.75 mag above it (upper line). The larger dots indicate the two components of our one candidate pair that passes all of our criteria. The red and black dots correspond to, respectively, the faster and slower moving objects. The cluster members do not include blue stragglers and other binary evolution products because they are far from a standard isochrone fit.

5 Discussion

This section examines potential limitations of our methodology and suggests possible enhancements for future research. Addressing these areas could strengthen the approach and improve results.

5.1 False positives

How many false positives do we expect from our analysis? Said another way, if we choose a random position on the sky for the cluster center-or-mass, how many pairs in a field of view comparable to the one assumed herein do we expect to simultaneously satisfy all of our selection criteria? To address this, we choose a position on the sky shifted in either RA or Dec by 150 pc and then repeat our analysis for all pairs in this field of view. This generates four different tests corresponding to the four different RA/Dec combinations (i.e., a positive/negative shift in RA and/or a positive/negative shift in Dec). In total, this generates an initial sample of 468 586 715 pairs to which we apply our methodology, which is four times greater than the initial sample size we obtained. Out of all these pairs, zero objects satisfy all of our criteria simultaneously (mostly due to the CMD location requirements).

Based on this analysis, we expect the frequency of false positive detection to be lower than 2 × 10⁻⁹ and hence lower than the detection frequency found in this paper (i.e., 10⁻⁸). The most likely false positives should in general be associated with dimmer objects, due to the larger uncertainties. Despite having roughly four times more pairs compared to the original test, none passed all of our selection criteria.

5.2 The effects of the galactic potential

Our method neglects the effects of the cluster and Galactic potentials. This is because the expectation is that any changes imparted to the velocity vectors of escapers from M67 should be small over our chosen field of view. This assumption turns out to be a good approximation when both objects in a given pair have high velocities and are located close to the cluster such that their travel times to their currently observed locations are short. However, this assumption can fail substantially when the object velocities are low and their positions are far from the host cluster. This is illustrated in Fig. 8, which is meant as an independent calculation to what is shown from Corespray in Figs. 6 and 5. Ultimately, Fig. 8 confirms what was found in Figs. 6 and 5, which suggests that velocity ratios significantly larger than 2 are common theoretically (see the top panel of Fig. 6) but much less common observationally (see the bottom panel of Fig. 6).

Hence, a simple cut in only the velocity ratio need not be implemented in future work applying the method presented in this paper to other star clusters. For now, we again note that removing this cut from our final sample selection does not affect our results.

Figure 8 shows the results of simulating 500 ejection events from M67. For each, the trajectories of the ejected objects are integrated into the Galactic potential using AMUSE (Astrophysical Multipurpose Software Environment, Pelupessy et al. 2013), a component library for performing astrophysical simulations involving different physical domains and scales, and galpy: a Python Library for Galactic Dynamics (Bovy 2015). We work over a 50 pc distance (i.e., a distance representative of that for our identified candidate pair) from the center of mass of the progenitor cluster (for the single star only), and the final angle between their velocity vectors as well as the change in magnitude of the velocities are calculated. The initial velocity of the binary is set equal to the cluster escape velocity, and the initial single star speed is then decided assuming linear momentum conservation. Hence, we adopted the minimum possible ejection velocity, which maximizes the effects of the Galactic potential in changing the object velocities over our chosen field of view. We note that we work with 3D velocities, but these results remain unchanged using 2D velocities.

As indicated by the color coding in the bottom panel of Fig. 8, most pairs experience the most significant change in velocity when ejected parallel to the orbital motion of the cluster and the smallest change in velocity when ejected perpendicular to it. We note as well that objects tend to be decelerated in the former case and accelerated in the latter. As is clear, for those objects leaving the cluster with small escape speeds, it is indeed possible for the final angle of intersection between their velocity vectors to deviate significantly from 180^∘. This becomes even more so if a larger field of view is considered, and the orbits are integrated out to 150 pc or more. If we instead adopt for the initial binary velocity ten times the cluster escape speed, then Fig. 8 changes dramatically, and 100% of the pairs have a final intersection angle within 5^∘ of 180^∘ since the typical uncertainty on this angle are a distribution around 4.4^∘. We note that, for the one candidate pair identified in this paper, the observed velocities are of order a factor of ten larger than the cluster escape speed, and both objects are close to M67 in projection on the sky.

This general behavior is also found in our Corespray simulations. As shown in Fig. 9, including all sources with intersection angles >175^∘ (i.e., roughly 5σ away from 180^∘) reduces our sample size by less than a factor of two, whereas only including those sources with angles of intersection >179^∘ reduces our sample by more than two orders of magnitude. Importantly, of those remaining 31 simulated pairs, only four have velocity vectors that are consistent with pointing orthogonally away from the direction of the orbital motion of the cluster to within 5^∘. The final sub-sample from Corespray then yields four pairs that satisfy all of our criteria while also having trajectories that should minimize the effects of the Galactic potential by traveling orthogonally to the orbital motion of the cluster (by ±5^∘).

Thus, Corespray predicts four pairs out of an initial sample of 5000 that meet the same criteria as our single best candidate, for which the observed parameters are summarized in Sect. 4.5. This potentially introduces another factor into our analysis, which suggests that many ejection events are likely needed before one will adhere to all of our selection criteria. How to calculate a correction factor coming from our theoretical analysis is non-trivial, and should be one of the focuses of future work. But, the theoretical predictions from Corespray confirm that some candidate pairs adhering to our strictest selection criteria should exist, and typically be moving tangentially away from the cluster center-of-mass motion in the Galactic potential since this minimizes possible changes to the RS-RB velocity vectors postejection from the cluster. The importance of taking into account the Galactic field is highlighted by these results and ultimately summarized in Figs. 8 and 9: the Galactic potential can erase any signature of a three-body decay significantly and quickly, and it becomes the dominant effect when within a few tidal radii from the host cluster.

A calculation using conservation of energy suggests that a single star escaping the cluster from the central core will only have its velocity changed by less than 1 km s⁻¹. Hence, we do not expect that taking the cluster potential into account will change our results significantly, suggesting as well that the Galactic potential is the more important effect, at least for such low-mass OCs as M67 and candidates with low velocities that lie far from the cluster tidal radius in projection.

To summarize, the results of this preliminary analysis suggest that the effects of the Galactic potentials should be considered in future work when searching for ejected RS-RB pairs. The effect is minimized for faster objects located closer to the progenitor cluster, as is the case for our one resulting pair and confirmed by the Corespray and AMUSE simulations.

Fig. 8 Results of integrating the trajectories of escapers from M67 in the Galactic potential until the least massive object reaches 50 pc from the cluster center. Each inset is plotted with respect to the final angle between the 3D velocity vectors of the escapers, θvf. (The initial angle between the velocity vectors is 180∘.) The top panel shows a histogram of the distribution of θvf for all pairs. The second panel shows the total change in velocity magnitude for each object in each pair. The third panel shows the final distances from the cluster center. In the bottom panel, the y-axis shows the initial direction of ejection of the lighter star, where θGC = 0∘ corresponds to ejection directly toward the Galactic center, and θGC = 90∘ corresponds to ejection in the direction of motion of the cluster.

5.3 The sixth dimension

Unfortunately, a RAVE radial velocity is not provided for our candidate. In order to further. In order to further improve our method, radial-velocity follow-up observations can be used. This cannot only constrain whether or not a given pair is consistent with having originated from the progenitor cluster along the line-of-sight, but also potentially the binary nature of the RB candidate and even its orbital properties. This makes additional predictions based on momentum conservation and causality that can be tested using radial velocity measurements (to confirm/reject that it is consistent with having originated from M67 and confirm/reject that it is a binary).

Fig. 9 Resulting velocity ratio, angle of intersection, and TBT ratio distributions from Corespray. From top to bottom, we isolate pairs with angles greater than 179∘, 175∘, and 170∘ (but less than 180∘).

6 Summary

In this paper, we have presented a systematic methodology to identify RSs ejected from star clusters due to single-binary interactions. This was done by exploiting the fact that we expect such disintegrations to produce a pair of objects (i.e., a single star and a binary) whose relative velocity vectors, positions on the sky, and positions in the progenitor cluster CMD are consistent with linear momentum conservation and causality. In other words, both objects are consistent with having originated from within the cluster at the same time.

Using the old OC M67 as a benchmark test case, we compared our theoretical expectations to a 5D kinematic data set obtained from the Gaia DR3 catalog via the provided parallaxes, proper motions, and positions on the sky. Our analysis yields, out of an initial sample size of roughly 10⁸ candidate pairs, only one single-binary pair that simultaneously satisfies all of our selection criteria. We have presented the observed distributions for our entire sample to illustrate which of our criteria are the most constraining, and we compared these distributions to theoretical predictions, coming mostly from Corespray. We find that a small fraction (<1%) of our initial sample of candidate pairs has an angle of intersection between their velocity vectors that is consistent with 180^∘ to within 1σ. Additionally, on the order of 10% of our initial sample simultaneously has a ratio of velocities that exceed two, with the brighter of the two objects in the pair traveling the slowest. Both of these criteria are expected from linear momentum conservation and tend to be more constraining than simple causality arguments related to the relative travel times and distances of both objects in a given pair being commensurate. With that said, the velocity ratio alone is not very constraining since velocity ratios greater than two are common (see the top panel in Fig. 6).

While the results obtained in this study for M67 are promising, it is important to bear in mind that our candidate pair is still preliminary. It should be further constrained using, for example, radial velocities and abundances as well as a more rigorous false positive test. Also, it is important to validate our method by applying it to other star clusters. This would not only further test the effectiveness of our method but would also potentially provide valuable insight into the dynamical evolution of other clusters. Furthermore, by studying multiple clusters, we could potentially identify trends in the properties of RS-RB pairs and test current theoretical predictions for the internal dynamical evolution of the progenitor clusters. That is, our method can reveal differences in runaway demographics between different clusters, which would give valuable information on the dynamical evolution of different types of clusters. In future work, we hope to apply our method to a large sample of clusters in order to generate the required statistics to address these questions.

Acknowledgements

AEHR acknowledges financial support from Millenium Nucleus NCN19_058 (TITANs). NWCL gratefully acknowledges the generous support of a Fondecyt General grant 1230082, as well as support from Millenium Nucleus NCN2023_002 (TITANs) and funding via the BASAL Centro de Excelencia en Astrofisica y Tecnologias Afines (CATA) grant PFB-06/2007. NWCL also thanks support from ANID BASAL project ACE210002 and ANID BASAL projects ACE210002 and FB210003. AMG acknowledges financial support from NSF AAG grant no. AST-2107738. SMG acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) and is partially funded through a NSERC Postgraduate Scholarship – Doctoral. RDM gratefully acknowledges financial support by the Fulbright U.S. Scholar Program, which is sponsored by the U.S. Department of State and the Chile-American Fulbright Commission. Its contents are solely the responsibility of the author and do not necessarily represent the official views of the sponsors.

Data availability

The data underlying this article will be shared on reasonable request to the corresponding author.

Appendix A Reference frame

Appendix A.1 Positions

We begin with the positions for all the sources in the catalog, which we first consider in a two dimensional frame. Using μ_α* ≡ μ_α cos δ, the coordinate transformations for the positions frame are $$\begin{array}{rl}x& =\cos \delta \sin (\alpha -{\alpha }_0)\\ y& =\sin \delta \cos {\delta }_0-\cos \delta \sin {\delta }_0\cos (\alpha -{\alpha }_0)\end{array}$$(A.1)

The positions of the sources are denoted by α (RA) and δ (Dec). The center of the cluster is identified by its RA and Dec, taking values α₀ = 8^h51^m23^s.3, δ₀ = +11^∘49^′02” (J 2000) (Jacobson et al. 2011)

Appendix A.2 Proper motions

Our goal is to determine two components of stellar velocities based on the proper motions μ_α and μ_δ.

The apparent motions of stars within stellar systems can be influenced by perspective effects and the motion of the system’s center of mass, as noted in van Leeuwen (2009). In particular, the motion of a cluster in the radial direction can produce an effect known as “perspective expansion,” which has been observed in Gaia measurements of GCs, many of which have radial velocities (RVs) reaching hundreds of kilometers per second Gaia Colaboration et al. (2018b). A first-order approximation of perspective expansion can be derived from Eq. 13 in van Leeuwen (2009). The equations for the additional shift in proper motion of a star are $$\begin{array}{rl}\mathrm{\Delta }{\mu }_{\alpha \ast ,\text{\hspace{0.17em}per\hspace{0.17em}}}& \approx \mathrm{\Delta }{\alpha }_{\mathrm{i}}({\mu }_{\delta ,0}\sin {\delta }_0-\frac{v_{\mathrm{r}}{\varpi }_0}{\kappa }\cos {\delta }_0),\\ \mathrm{\Delta }{\mu }_{\delta ,\text{\hspace{0.17em}per\hspace{0.17em}}}& \approx -\mathrm{\Delta }{\alpha }_i{\mu }_{\alpha \ast ,0}\sin {\delta }_0-\mathrm{\Delta }{\delta }_{\mathrm{i}}\frac{v_{\mathrm{r}}{\varpi }_0}{\kappa }.\end{array}$$(A.2)

Here, Δα_i and Δδ_i are the differences in right ascension and declination between the system center of mass and an individual star. The proper motions in RA and Dec are, respectively, μ_δ,0 = −10.97 mas yr⁻¹ and μ_δ,0 = −2.9396 mas/yr⁻¹. We assume v_r = 33.64 ± 0.03 km s⁻¹ (Geller et al. 2015) for the radial velocity of the center-of-mass motion of M67 and ϖ₀ = 860 pc for the cluster’s parallax (Cantat-Gaudin et al. 2018). The conversion factor that converts milliarcseconds per year to kilometers per second at a distance of 1 kpc is κ = 4.74. The first term in each equation accounts for motion within a spherical coordinate system, while the second term accounts for the perceived expansion or contraction of an object as its distance increases or decreases.

These contributions can be subtracted from the observed proper motions (relative to the system center of mass) to isolate the influence of internal kinematics (Brown et al. 1997). We calculated the components of the corrected velocities parallel to lines of constant right ascension and declination: $$\begin{array}{rl}{\mu }_{\alpha ,\text{rest}}& \approx -\kappa \left(\frac{\mathrm{\Delta }{\mu }_{\alpha \ast ,\text{obs}}-\mathrm{\Delta }{\mu }_{\alpha \ast ,\text{per}}}{{\varpi }_0}\right),\\ {\mu }_{\delta ,\text{rest}}& \approx \kappa \left(\frac{\mathrm{\Delta }{\mu }_{\delta \ast ,\text{obs}}-\mathrm{\Delta }{\mu }_{\delta \ast ,\text{per}}}{{\varpi }_0}\right)\end{array}$$(A.3)

with Δμ_α*,obs = μ_α* − μ_α*,0 and Δμ_δ*,obs = μ_δ − μ_δ,0 are the differences in proper motion in right ascension and declination, respectively, between an individual star and the system center of mass. We have as a result the μ_α,rest and μ_δ,rest as the proper motions in the rest of frame calculated for RA and DEC. We then transformed the velocities to a 2D Cartesian coordinate system, v_x and v_y, using the orthographic projection (Gaia Colaboration et al. 2018b): $$\begin{array}{rl}& v_{\mathrm{x}}={\mu }_{\alpha ,\text{\hspace{0.17em}rest\hspace{0.17em}}}\cos (\alpha -{\alpha }_0)-{\mu }_{\delta ,\text{\hspace{0.17em}rest\hspace{0.17em}}}\sin \delta \sin (\alpha -{\alpha }_0)\\ & v_{\mathrm{y}}={\mu }_{\alpha ,\text{\hspace{0.17em}rest\hspace{0.17em}}}\sin {\delta }_0\sin (\alpha -{\alpha }_0)+\\ & {\mu }_{\delta ,\text{\hspace{0.17em}rest\hspace{0.17em}}}(\cos \delta \cos {\delta }_0+\sin \delta \sin {\delta }_0\cos (\alpha -{\alpha }_0)).\end{array}$$(A.4)

Appendix A.3 Cross-point and angle of the stars

We used the following equation to find the cross-point. With X equal to all the possible points on the line connecting each velocity vector to the center-of-mass motion of the cluster.: $$\mathit{X}=\lambda \hat{n}+{\mathit{X}}_0,$$(A.5)

The factor λ determines the cross-point. The unit vector $\hat{n}$ that we adopt satisfies $\mathit{n}=\frac{\mathit{v}}{|\mathit{v}|}$. The vector equations of the line for sources a and b are: $${\mathit{X}}_{\mathrm{a}}={\lambda }_{\mathrm{a}}{\hat{n}}_{\mathrm{a}}+{\mathit{r}}_{\mathrm{a}}={\lambda }_{\mathrm{a}}\frac{{\mathit{v}}_{\mathrm{a}}}{\left|{\mathit{v}}_{\mathrm{a}}\right|}+{\mathit{r}}_{\mathrm{a}}$$(A.6) $${\mathit{X}}_{\mathrm{b}}={\lambda }_{\mathrm{b}}{\hat{n}}_{\mathrm{b}}+{\mathit{r}}_{\mathrm{b}}={\lambda }_{\mathrm{b}}\frac{{\mathit{v}}_{\mathrm{b}}}{\left|{\mathit{v}}_{\mathrm{b}}\right|}+{\mathit{r}}_{\mathrm{b}},$$(A.7)

where X_a need to be equal to X_b. We needed to find the λ_a and λ_b values that satisfy the following condition: $${\lambda }_{\mathrm{a}}=\frac{n_{\mathrm{a}}^{\mathrm{y}}(r_{\mathrm{a}}^{\mathrm{x}}-r_{\mathrm{b}}^{\mathrm{x}})-n_{\mathrm{b}}^{\mathrm{x}}(r_{\mathrm{s}}^{\mathrm{y}}-r_{\mathrm{b}}^{\mathrm{y}})}{n_{\mathrm{a}}^{\mathrm{y}}{n}_{\mathrm{b}}^{\mathrm{x}}-n_{\mathrm{a}}^{\mathrm{x}}{n}_{\mathrm{b}}^{\mathrm{y}}}$$(A.8) $${\lambda }_{\mathrm{b}}=\frac{n_{\mathrm{b}}^{\mathrm{y}}(r_{\mathrm{a}}^{\mathrm{x}}-r_{\mathrm{b}}^{\mathrm{x}})-n_{\mathrm{a}}^{\mathrm{x}}(r_{\mathrm{a}}^{\mathrm{y}}-r_{\mathrm{b}}^{\mathrm{y}})}{n_{\mathrm{a}}^{\mathrm{y}}{n}_{\mathrm{b}}^{\mathrm{x}}-n_{\mathrm{a}}^{\mathrm{x}}{n}_{\mathrm{b}}^{\mathrm{y}}}.$$(A.9)

To calculate the angle between the RSs-RBs, we used the following equation: $${\theta }_{\mathrm{a},\mathrm{b}}=\arccos \frac{(v_{\mathrm{x}}^{\mathrm{a}}{v}_{\mathrm{x}}^{\mathrm{b}}+v_{\mathrm{y}}^{\mathrm{a}}{v}_{\mathrm{y}}^{\mathrm{b}})}{v^{\mathrm{a}}{v}^{\mathrm{b}}}.$$(A.10)

By propagating the errors, we obtained a value corresponding to 1σ of θ_1σ. Then the upper and lower angles were obtained using $${\theta }_{\mathrm{a},\mathrm{b}}^{+}=\arccos \frac{(v_{\mathrm{x}}^{\mathrm{a}}{v}_{\mathrm{x}}^{\mathrm{b}}+v_{\mathrm{y}}^{\mathrm{a}}{v}_{\mathrm{y}}^{\mathrm{b}})}{v^{\mathrm{a}}{v}^{\mathrm{b}}}+{\theta }_{1\sigma }$$(A.11) $${\theta }_{\mathrm{a},\mathrm{b}}^{-}=\arccos \frac{(v_{\mathrm{x}}^{\mathrm{a}}{v}_{\mathrm{x}}^{\mathrm{b}}+v_{\mathrm{y}}^{\mathrm{a}}{v}_{\mathrm{y}}^{\mathrm{b}})}{v^{\mathrm{a}}{v}^{\mathrm{b}}}-{\theta }_{1\sigma }.$$(A.12)

The numerator of the fraction in the above equations, $(v_{\mathrm{x}}^{\mathrm{a}}{v}_{\mathrm{x}}^{\mathrm{b}}+v_{\mathrm{y}}^{\mathrm{a}}{v}_{\mathrm{y}}^{\mathrm{b}})$, represents the dot product of the velocity vectors of objects a and b. The denominator, v^a v^b, represents the product of the magnitudes of the velocity vectors of objects a and b. Taking the fraction $\frac{v_x^{\mathrm{a}}{v}_x^{\mathrm{b}}+v_y^{\mathrm{a}}{v}_y^{\mathrm{b}}}{v^{\mathrm{a}}{v}^{\mathrm{b}}}$, we obtain the cosine of the angle between the two velocity vectors.

References

All Tables

All Figures

	Fig. 1 Histogram showing the predicted velocity distributions of escapers ejected in the Corespray simulation. The x-axis shows the ejection velocity in kilometers per second.
In the text

	Fig. 2 Four different snapshots in time of a simulation using Corespray. We show the spatial distribution of the ejected pairs in the plane of the sky at 0.1, 1, 10, and 100 Myr.
In the text

Fig. 3 Diagram illustrating the selection method for the angle of intersection. We have the two limiting vectors described by the error in position and proper motion (pma ± Δ pma and pmb ± Δ pmb) for the RS and RB, which yield a minimum (lower) and maximum (upper) possible angle (using the 1σ uncertainties in both directions; see Eqs. (A.10), (A.11), and (A.12)). If a given pair has an angle of intersection consistent with 180∘ to within these upper and lower bounds, we keep it in our sample and discard it otherwise. We note that the RB is not resolved and the cross point is at the center of light.

In the text

	Fig. 4 Histograms showing the distribution of the difference in TBTs (first and second panel) as well as their ratio (third and fourth panel) for every pair in our observational data and a 2D velocity projection in corespray. The TBTs are provided in years.
In the text

Fig. 5 Histogram showing the distribution of intersection angles between the primary vector of every pair. The x-axis shows the angle in degrees. The top panel shows the projected angle for every pair from our Corespray simulations. The bottom panel shows the distribution with observational data. The simulation shows that most of the pairs tend to have angles near 180∘. In contrast, the pairs from the observational data set have a very low number of pairs near 180∘. This difference between the simulation (only pairs resulting from three-body disintegration) and observational (random sources in a field of view) distributions shows the angle to be one of the most important filters for the method.

In the text

Fig. 6 Histogram showing the distribution of velocity ratios. The top panel represents every pair from a simulation using Corespray, and the bottom panel shows all the pairs in our observational data. The vertical dashed line shows a critical ratio of two. Theory predicts that it is very likely that RS-RBs have velocity ratios greater than two, as can be seen in the simulated data set. On the other hand, the observational data show that more than half of the pairs have velocity ratios of less than two.

In the text

Fig. 7 Color-magnitude diagram of the old OC M67 plotted using the Gaia data. Only stars indicated as high-probability cluster members in Childs et al. (2024) have been selected and are plotted as small blue dots. The solid red lines show the best-fitting isochrone (lower line) as well as the equal-mass binary sequence shifted 0.75 mag above it (upper line). The larger dots indicate the two components of our one candidate pair that passes all of our criteria. The red and black dots correspond to, respectively, the faster and slower moving objects. The cluster members do not include blue stragglers and other binary evolution products because they are far from a standard isochrone fit.

In the text

Fig. 8 Results of integrating the trajectories of escapers from M67 in the Galactic potential until the least massive object reaches 50 pc from the cluster center. Each inset is plotted with respect to the final angle between the 3D velocity vectors of the escapers, θvf. (The initial angle between the velocity vectors is 180∘.) The top panel shows a histogram of the distribution of θvf for all pairs. The second panel shows the total change in velocity magnitude for each object in each pair. The third panel shows the final distances from the cluster center. In the bottom panel, the y-axis shows the initial direction of ejection of the lighter star, where θGC = 0∘ corresponds to ejection directly toward the Galactic center, and θGC = 90∘ corresponds to ejection in the direction of motion of the cluster.

In the text

	Fig. 9 Resulting velocity ratio, angle of intersection, and TBT ratio distributions from Corespray. From top to bottom, we isolate pairs with angles greater than 179∘, 175∘, and 170∘ (but less than 180∘).
In the text