Two main mechanisms have classically been proposed for the formation of runaway stars. In the binary supernova scenario BSS a massive star in a binary explodes as a supernova ejecting its companion. In the dynamical ejection scenario a star is ejected during a strong dynamical encounter between multiple stars. We propose a third mechanism for the formation of runaway stars the subcluster ejection scenario SCES where a subset of stars from an infalling subcluster is ejected out of the cluster via a tidal interaction with the contracting gravitational potential of the assembling cluster. We demonstrate the SCES in a starbystar simulation of the formation of a young massive cluster from a 10SUP6SUP MSUBSUB gas cloud using the TORCH framework. This star cluster forms hierarchically through a sequence of subcluster mergers determined by the initial turbulent spherical conditions of the gas. We find that these mergers drive the formation of runaway stars in our model. Lateforming subclusters fall into the central potential where they are tidally disrupted forming tidal tails of runaway stars that are distributed highly anisotropically. Runaways formed in the same SCES have similar ages velocities and ejection directions. Surveying observations we identify several SCES candidate groups with anisotropic ejection directions. The SCES is capable of producing runaway binaries two wide dynamical binaries in infalling subclusters were tightened through ejection. This allows for another velocity kick via subsequent via a subsequent BSS ejection. An SCESBSS ejection is a possible avenue for the creation of hypervelocity stars unbound to the Galaxy. The SCES occurs when subcluster formation is resolved. We expect nonspherical initial gas distributions to increase the number of calculated runaway stars bringing it closer to observed values. The observation of groups of runaway stars formed via the SCES can thus reveal the assembly history of their natal clusters.
['stars: kinematics and dynamics', 'ISM: clouds', 'galaxies: star clusters: general', 'Astrophysics - Astrophysics of Galaxies']
Massive star cluster formation II. Runaway stars as fossils of subcluster mergers
{'II. Runaway stars as fossils of subcluster mergers': 'Brooke Polak 1 , 2, ⋆ , Mordecai-Mark Mac Low 2 , Ralf S. Klessen 1 , 4 , Simon Portegies Zwart 5 , Eric P. Andersson 2 , Sabrina M. Appel 6 , Claude Cournoyer-Cloutier 7 , Simon C. O. Glover 1 , and Stephen L. W. McMillan 3  \n- 1 Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Heidelberg, Germany e-mail: bpolak@amnh.org\n- 2 Department of Astrophysics, American Museum of Natural History, New York, NY, USA\n- 3 Department of Physics, Drexel University, Philadelphia, PA, USA\n- 4 Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Heidelberg, Germany\n- 5 Sterrewacht Leiden, Leiden University, Leiden, the Netherlands\n- 6 Department of Physics and Astronomy, Rutgers University, Piscataway, NJ, USA\n- 7 Department of Physics and Astronomy, McMaster University, Hamilton, ON, Canada  \nReceived 17 May 2024; accepted 5 August 2024', 'ABSTRACT': 'Two main mechanisms have classically been proposed for the formation of runaway stars. In the binary supernova scenario (BSS), a massive star in a binary explodes as a supernova, ejecting its companion. In the dynamical ejection scenario, a star is ejected during a strong dynamical encounter between multiple stars. We propose a third mechanism for the formation of runaway stars: the subcluster ejection scenario (SCES), where a subset of stars from an infalling subcluster is ejected out of the cluster via a tidal interaction with the contracting gravitational potential of the assembling cluster. We demonstrate the SCES in a star-by-star simulation of the formation of a young massive cluster from a 10 6 M ⊙ gas cloud using the torch framework. This star cluster forms hierarchically through a sequence of subcluster mergers determined by the initial turbulent, spherical conditions of the gas. We find that these mergers drive the formation of runaway stars in our model. Late-forming subclusters fall into the central potential, where they are tidally disrupted, forming tidal tails of runaway stars that are distributed highly anisotropically. Runaways formed in the same SCES have similar ages, velocities, and ejection directions. Surveying observations, we identify several SCES candidate groups with anisotropic ejection directions. The SCES is capable of producing runaway binaries: two wide dynamical binaries in infalling subclusters were tightened through ejection. This allows for another velocity kick via subsequent via a subsequent BSS ejection. An SCES-BSS ejection is a possible avenue for the creation of hypervelocity stars unbound to the Galaxy. The SCES occurs when subcluster formation is resolved. We expect nonspherical initial gas distributions to increase the number of calculated runaway stars, bringing it closer to observed values. The observation of groups of runaway stars formed via the SCES can thus reveal the assembly history of their natal clusters.  \nKey words. Star clusters - star formation - ISM: clouds - globular clusters: general - stars: kinematics and dynamics', '1. Introduction': "Young stars displaced from their birthplace in the star-forming spiral arms of the Galaxy and moving away from the Galactic disk were first observed by Blaauw & Morgan (1954). Since then, many more stars moving rapidly away from their formation sites have been observed throughout the Galaxy. These stars are typically classified as runaways when their velocity relative to their associated nebula or cluster is ≥ 30 km s -1 (Gies & Bolton 1986).  \nDespite their ubiquity, the ejection mechanism remains unknown for most runaway stars. Currently, there are two popular proposed mechanisms for producing runaway stars: the binary supernova scenario (BSS; sometimes referred to as Blaauw kicks; Blaauw 1961) and the dynamical ejection scenario (DES; Poveda et al. 1967; Hoogerwerf et al. 2000; Fujii & Portegies Zwart 2011). According to the BSS, when two massive stars are  \nin a binary system and one explodes as a supernova (SN), loss of the ejecta from the system reduces the gravity on the companion, which begins moving through space at a velocity comparable to its orbital velocity. According to the DES, runaway stars are ejected during a strong dynamical encounter involving at least one binary, in which orbital binding energy is converted to kinetic energy. In some cases, a runaway system is formed moving in the opposite direction of the system that ejected it (in the center-of-mass frame of the encounter; Poveda et al. 1967; Fujii & Portegies Zwart 2011). BSS and DES runaways can be distinguished by their rotational and linear velocities: the BSS produces slow-moving, rapidly rotating stars, while the DES produces fast-moving, slowly rotating stars (Sana et al. 2022).  \nIn this work we propose a third mechanism for producing runaway stars: the subcluster ejection scenario (SCES), in which part of a subcluster is ejected from the cluster after it is tidally disrupted by the contracting potential of the assembling cluster. Star clusters form hierarchically, with giant molecular clouds (GMCs) fragmenting into dense clumps that form subclusters  \nof stars. These subclusters merge, forming a single central cluster. The combined feedback eventually blows away the unused gas, leaving a gas-free star cluster (see, e.g., Rahner et al. 2017; Grudi'c et al. 2018; Rahner et al. 2019; Wall et al. 2019, 2020; Cournoyer-Cloutier et al. 2021; Lewis et al. 2023; Wilhelm et al. 2023; Cournoyer-Cloutier et al. 2023). Lucas et al. (2018) shows through N -body simulations that merging groups of stars can produce a tidal tail of unbound stars. The SCES describes when this mechanism occurs between young subclusters during initial cluster assembly. To thoroughly investigate the origin of SCES runaway stars, the entire complex dynamical history of young star clusters must be modeled consistently. This necessitates modeling the entire star cluster formation process from the birth of stars within subclusters.  \nWepresent a star-by-star simulation of cluster formation performed with the torch 1 framework (Wall et al. 2019, 2020). torch follows gas dynamics, the N -body dynamics of stars, subgrid star formation via sink particles, stellar evolution, and stellar feedback in the form of winds, radiation, and SNe. With this simulation, we can determine the origin of runaway stars formed self-consistently in the cluster environment.  \nIn this paper we analyze the runaway stars formed in the M6 star cluster presented in Polak et al. (2023, hereafter Paper I), which forms from a molecular cloud with an initial mass of M cloud = 10 6 M ⊙ and an initial radius of R cloud = 11 . 7 pc. We briefly describe our model and the initial conditions in Sect. 2. In Sect. 3 we present the properties of the runaway stars formed in M6, as well as our evidence that the SCES mechanism produces runaway stars. We discuss the implications of our results in Sect. 4 and present some observational candidates for runaways formed via the SCES. We conclude in Sect. 5.", '2. Methods': "Weused torch (Wall et al. 2019, 2020) to model the formation of a star cluster from a turbulent spherical cloud of gas. The torch star cluster formation framework is built with the Astrophysical MUltpurpose Software Environment (AMUSE; Portegies Zwart et al. 2009, 2013; Pelupessy et al. 2013; Portegies Zwart & McMillan 2018) framework, linking several physics codes addressing hydrodynamics, stellar evolution, stellar dynamics, star formation via sink particles, and stellar feedback in the form of winds, radiation, and SNe. In this section we provide a brief overview of the numerical methods used for these simulations. A more detailed description of the applied methods is provided in Paper I.  \nThe adaptive mesh refinement (AMR) magnetohydrodynamics code flash (Fryxell et al. 2000) handles the gas dynamics. The gas dynamics are evolved with an HLLD Riemann solver (Miyoshi & Kusano 2005) and third-order piecewise parabolic method reconstruction (Colella & Woodward 1984). The selfgravity of the gas is evolved with a multi-grid Poisson solver (Ricker 2008), and the gravitational interaction between the stars and gas is modeled using a leapfrog gravity bridge scheme (Fujii et al. 2007; Wall et al. 2019).  \ntorch uses a sub-grid star formation model via sink particles in flash . Sink particles form when the local gas density exceeds a threshold value ρ sink based on the Jeans density at the maximum refinement level, and when several other criteria are met (see Federrath et al. 2010). When a sink particle forms, the gas  \nTable 1. Simulation parameters.  \nNotes. Rows: initial cloud mass, central density, average density, surface density, peak initial vertical magnetic field (see Eq. 1), initial central Jean's length, initial free-fall time, cloud radius, half-width of box, virial parameter, average velocity and velocity dispersion of the gas within R cloud, minimum cell width, maximum cell width, sink radius, sink threshold density, approximate initial sink mass, minimum feedback star mass, agglomeration mass of low-mass stars, mass sampling range of Kroupa IMF.  \nwithin the sink accretion radius r sink and above the threshold density is added to the sink's mass. For each sink, the Kroupa initial mass function (IMF; Kroupa 2002) is randomly sampled to create a list of stellar masses for the sink to form sequentially. From this list, the sink forms star particles until it runs out of mass. The sink continues to form stars whenever it accretes enough mass to form the next star on the list. When stars form, they are placed randomly in a uniform spherical distribution within the sink accretion radius and given a velocity equal to the sink's velocity plus an additional velocity component with the direction sampled isotropically and the magnitude sampled from a Gaussian distribution with a standard deviation equal to the local sound speed. The local sound speed is taken to be the average sound speed of gas with temperature T ≤ 100 K in cells ≤ 2 r sink from the sink. If no gas fulfills this criterion within the region, the sink is blocked from forming stars. Sink momentum is set by adding the current sink momentum to the momentum of the gas it accretes.  \nStellar feedback in the form of stellar winds, radiation, and SNe, are also implemented as additional physics in flash (Wall et al. 2020). Radiative transport is modeled in Flash using the ray-tracing routine Fervent (Baczynski et al. 2015). Stars are evolved from the zero-age main sequence until their death via S e B a (Portegies Zwart & Verbunt 1996), which informs the stellar feedback properties in flash . petar (Wang et al. 2020) evolves the N-body dynamics of the star particles.  \nTo reduce the computational expense caused by the large number of stars formed in a 10 6 M ⊙ cloud, we applied three modifications to the standard torch model:  \n- 1. Star particles forming below M agg = 4 M ⊙ are agglomerated until their summed mass is above M agg and then are assigned to a single super-star particle. This reduces the number of stars from ∼ 10 6 to ∼ 10 5 to make the N-body calculations feasible.\n- 2. The mass-loading of stellar winds is artificially raised to limit the temperature of wind-blown bubbles to T w = 3 × 10 5 K thereby alleviating the timestep constraint set by the Courant condition. Although this causes wind bubbles to be smaller and cooler, the primary e ff ect of winds is to clear out dense gas, allowing ionizing radiation to escape and form H ii regions. The clearing of dense gas still occurs with the massloaded winds. At the densities in our simulation, hot stellar winds rapidly cool regardless of mass-loading.\n- 3. Ray-tracing on the AMR grid quickly becomes expensive for a large number of sources. To reduce the cost of the raytracing calculations, we only include feedback (radiation, winds, and SNe) from stars with M feedback ≥ 20 M ⊙ . Most of the mechanical wind energy and ionizing radiation in a star cluster comes from massive stars, so with this limit < 20% of the total feedback energy is lost (see Paper I). SNe of stars below this mass limit take ≥ 10 Myr , which is much longer than the scope of these simulations.  \nA detailed discussion of the e ff ect of these modifications is provided in Paper I.  \nWe used the M6 simulation introduced in Paper I of initial mass M cloud = 10 6 M ⊙ , with properties listed in Table 1. The initial density profile is a Gaussian (Bate et al. 1995; Goodwin et al. 2004) with ρ edge /ρ c = 1 / 3. We used outflow boundary conditions that allow inflow from ghost zones. The turbulence was set by imposing a Kolmogorov (1941) velocity spectrum on the gas. The initial magnetic field, B = Bz ˆ z , is uniform in z and decreases radially with the mid-plane density ρ ( x , y, z = 0) in the x -y plane:  \nBz ( x , y ) = B 0 , z exp h -( x 2 + y 2 ) ln(3) / R 2 cloud i , (1)  \nwhere B 0 , z = 18 . 5 µ G.", '3. Results': 'In our analysis, we made an initial selection for runaway stars by filtering for unbound stars, 2 with E total , i = Ek , i + U ⋆, i j + U gas > 0 , where Ek , i is the kinetic energy, U ⋆, i j is the gravitational potential of other star particles, and U gas is the gravitational potential of the gas. Our analysis only covers the early runaways ejected during the hierarchical assembly of the cluster. The free-fall time of the M6 cloud is t ff = 0 . 67 Myr, and the final simulation time is 1 . 36 t ff = 0 . 91 Myr. The cluster is fully assembled by the end of the simulation, and the SCES only occurs during active subcluster formation. This run time is therefore su ffi cient for analysis of the SCES. At the final time, the fraction of runaway stars is 0 . 51%. We cannot yet determine the final runaway star fraction due to the short run time of the simulation. Simulating the total runaway fraction requires integration times of ≳ 2 Myr, which is outside the scope of this paper introducing a production mechanism for early-forming runaway stars in young clusters.  \nFigure 1 shows a histogram of the runaway star count binned by the angular direction of the runaway velocities in the centerof-mass frame. There are two distinct groupings of runaway stars  \nFig. 1. Ejection velocity angular directions of M6 runaway stars with θ u and ϕ u binned in 9 · and 18 · angle bins, respectively. There are two distinctly peaked angular regions, which are labeled as group 1 and group 2. The rectangles indicate the angular bins of runaways we selected for trajectory analysis.  \n<!-- image -->  \nejected in the same directions, and the overall directional distribution of runaways is highly anisotropic. This is caused by the tidal disruption of two infalling subclusters by the cluster center of mass, which then form two groups of runaway stars as tidal tails. White boxes in Fig. 1 indicate the selection regions used to investigate the origin of these two groups of runaways. These regions were chosen by eye with the intent of performing a broad initial analysis.  \nThe properties of the runaway stars are depicted in Fig. 2, which shows a Mollweide spherical projection of the runaway stars colored according to radial velocity from the center-ofmass, mass, and age. The group 2 population is moving faster than group 1. The average velocity and standard deviation for group 1 and 2 are ¯ v 1 = 38 . 6 km s -1 , σ 1 = 8 . 2 km s -1 , ¯ v2 = 87 . 1 km s -1 , σ 2 = 20 . 7 km s -1 . Both groups are composed mostly of low-mass stars with a few high-mass stars scattered within. The stars in a given group were formed at similar ages.  \nFigure 3 shows the age and velocity of the runaways in a scatter plot colored according to stellar mass, with the selected angle bins in Fig. 1 marked by symbol type. The runaways are distinctly clustered in both radial velocity and stellar age. Group 1 is distinguished from the rest of the runaways in the upper left; stars in group 1 are older and have lower velocities ( v r ≈ 30 km s -1 ). Group 2, however, consists of two populations of runaways seen as an age gap in the diamond points in Fig. 3. This indicates stars from two subclusters or two generations of stars within one subcluster were ejected in the same direction. The mass distributions of both groups are the same, though, with both mostly containing low mass stars and a few stars above 10 M ⊙ . Figure 4 shows the Kroupa IMF power law with α = -2 . 3 and the histograms of the stellar masses in each runaway group as well as the entire cluster. The overpopulation of low mass bins ( < 8 M ⊙ ) is a result of the mass agglomeration of stars below < 4 M ⊙ . This suggests that each group consists of a random sampling of the IMF, which makes sense for a subcluster. Both runaway groups roughly resemble the IMF indicating that the SCES mechanism is unbiased toward stellar mass. We propose that this is because the stars in the subcluster closest to  \nFig. 2. Mollweide map projection of the ejection directions of the runaway stars in M6. The colors correspond to the radial velocity from the cluster center of mass (top), the star mass (middle), and the star age (bottom). The highest values are plotted over lower values for visibility.  \n<!-- image -->  \nTable 2. Runaway group filters.Notes. Selection ranges for runaway stars in groups 1 and 2. Rows: ϕ and θ directions of velocity, radial velocity from center of mass, and stellar age.  \nthe potential are tidally captured while the outer stars are ejected, and the spatial distribution of stellar masses within a subcluster is random by construction. Most of the runaway stars are agglomerated, which is expected from a sampling of the broader stellar population. This is acceptable as there is no apparent mass preference or cuto ff for ejection via the SCES.  \nWe applied a second filter to select the runaways belonging to the two groups using stellar age and radial velocity indicated in Fig. 3. We selected the younger of the two populations in group 2 to focus on a specific subcluster. The three most important characteristics for identifying groups of runaways are ejection direction, velocity, and stellar age. The complete list of  \nFig. 3. Scatter plot of the runaway stars, with the radial velocity and age colored by mass. Squares are group 1, diamonds are group 2, and dots are the rest of the runaways. These groups are the runaways selected by ejection direction only. The two groups are distinctly separate in both age and velocity. Group 2 contains two subgroups ejected in the same direction but at slightly di ff erent times. For the purpose of demonstrating the SCES mechanism, we focus on the later forming subgroup in group 2. The red rectangles indicate the additional filtering of the runaway groups by age and velocity.  \n<!-- image -->  \nFig. 4. Particle mass distributions of the entire star cluster (gray) and the runaway groups along with the Kroupa IMF power law. Note that at the low-mass end ( ≤ 8M ⊙ ), the particle masses are skewed higher than the sampled IMF due to the agglomeration of stars below < 4M ⊙ .  \n<!-- image -->  \nfilters for groups 1 and 2 is listed in Table 2. We find that selecting groups based on direction, age, and radial velocity provides a su ffi cient selection of runaways from a particular subcluster merger.  \nUsing these age, velocity, and direction filters, we traced the stars in the two groups back through the evolution and formation of the star cluster to determine their ejection mechanism. Figure 5 is a still from a movie (available in the online version of  \nFig. 5. Movie of the formation and ejection of runaway star group 1. Left: Stars as they are forming in the cluster. Right: Gravitational potential of both the stars and gas. The trajectories of the runaway stars are shown as lines. Once each runaway star forms, it appears with a star marker. This plot shows the same variables as shown in Figs. 6 and 7.  \n<!-- image -->  \nthe paper) showing the formation and ejection of runaway group 1 as well as the formation of the star cluster. In addition to the angular selection outlined in Table 2, we also imposed an age restriction to select stars that formed at the same time. In this way, we selected only one of the two age groups in group 2.  \nFigures 5-8 reveal tidal interactions of the subcluster with the main cluster as the physical mechanism for ejecting these two groups, explaining their shared ejection direction, age, and velocity as well as their uniformly sampled IMF mass composition. Group 1 and 2 stars form from distinct subclusters. The subcluster merges with the assembled central massive cluster becoming tidally disrupted. Some stars are tidally captured while others are slingshot around the central potential and escape as tidal tails.  \nThe timing of the contraction and expansion of the central cluster potential is crucial for the SCES. In a static potential, the stars falling in could not escape the cluster due to conservation of energy. In the SCES, a subcluster approaches the central potential as the potential contracts due to the assembly of other subclusters and gas infall. After the tidal interaction, as the SCES stars are leaving the central cluster, the potential well expands as star formation from the gas proceeds and the cluster virializes. The timing of the subcluster formation and cluster potential evolution determines whether the SCES stars become unbound.  \nGroup 2 is more populous than group 1. This is not because subcluster 2 formed more stars, but because the potential driving the tidal force on group 2 is stronger than when group 1 fell in, resulting in more stars being ejected. From Fig. 6 we see many stars that do not become runaways forming in the vicinity of group 1 stars, whereas the majority of stars forming near group 2 stars in Fig. 8 become runaways.  \nWe find that SCES runaways from later interactions with a deeper potential well result in a higher fraction of stars ejected within a subcluster. The tidal interaction is highlighted further by the figures showing trajectories with the evolving potential. At the onset of infall for group 1, the potential is significantly shallower than when group 2 forms. The peak of the potential has just formed and is still migrating when group 1 forms, but it is at its final destination when group 2 forms. However, the potential reaches its peak by the time of slingshot (point of maximum stellar velocity) for both groups. The greater integrated acceleration results in more stars in subcluster 2 being ejected than in subcluster 1. This correlation allows us to probe the assembly history of a star cluster from the ejection directions, ages, and velocities of runaway stars.', '4.1. Persistence of runaway groups through projection effects': 'In our simulated cluster, the anisotropy of runaway stars is directly linked to the subcluster merger history. To apply this to observations, we must determine whether the anisotropy can be clearly seen from any given observation angle, and the degree to which the groups of runaways appear to be distinct. We tested this by rotating the cluster about a randomly generated unit vector by a randomly generated angle and then recalculating the escape directions. Figure 9 shows a histogram of 1000 random projections (top) and histograms of the average, maximum, and minimum values of the projections as well as the cumulative sum of runaway fractions per bin (bottom). The bins are centered such that the central bin is the most prevalent escape angle for each projection. The histograms are strongly peaked, and in almost every projection (each row in the top panel) a fainter secondary peak in the escape direction is still visible. The secondary peak indicates a second SCES runaway group. The strongly peaked histograms, paired with the steep cumulative sum distribution, tells us that if there is a significant subcluster merger history in a cluster resulting in a population of SCES runaways, observers will see an anisotropic distribution of runaway stars regardless of viewing angle.', '4.2. Fossils of cluster assembly history': "Groups of runaways ejected by the SCES are distinctly grouped kinematically, temporally, and directionally. We can identify clusters ejected by this mechanism by finding corresponding grouping in age-velocity space and ϕ u -θ u space. With wellresolved velocities and ages of runaway stars from Gaia, this technique can be applied to find runaways originating from the same subcluster. If paired with other groups of SCES runaways, one can unwind the dynamical assembly history of the natal cluster. Furthermore, in the case of multiple runaway groups, the relative size and velocity of each group can indicate which subcluster merger happened first. The groups with more stars at higher velocities likely were ejected later. Since SCES groups can only be ejected while the central potential is contracting, the kinetic ages of multiple SCES groups can constrain the timescale of their natal cluster's assembly.  \nObservations typically only look for runaway O and B-type stars because of their short lifetime. If an O or B star is far from a star forming region, it must be a runaway whereas low-mass stars could just be field stars. Furthermore, there is a bias toward massive stars for ejection via BSS and DES : Massive stars have a multiplicity fraction of > 90% (Sana et al. 2014), and because of mass segregation massive stars are clustered in the stellar core and are more likely to have close dynamical encounters (Oh & Kroupa 2016). Therefore, the presence of many low-mass runaway stars is a unique feature of the SCES. There are 3 O stars ( ≥ 15 M ⊙ ) in group 1 and 7 O stars in group 2, so these SCES runaway groups would be detectable using just the O stars. In the case of anisotropic runaways, we predict these high-mass runaways to be accompanied by many more low-mass stars with similar velocities and ages. Depending on how many massive stars are in a runaway grouping, one can estimate how many low-mass stars accompany them using the IMF.  \nFig. 6. Trajectories of runaway star group 1 originating from a subcluster merging into the central cluster, with color indicating the stellar velocity (light = v max, dark = v min). The dots are star particles at the given time, and the star markers indicate the runaway stars in the group as they form. The total number of stars in the cluster is listed in each panel.  \n<!-- image -->", '4.3. Observational examples': 'The most promising example of the SCES is a group of runaways ejected to the north of R136 in the 30 Doradus star-forming region (Stoop et al. 2024a). Of the 18 runaways ejected < 1 Myr ago, 16 are ejected in the same direction. Additionally, there is an older group of runaways ejected more isotropically, which suggests that they were ejected by DES or BSS. There is a distinct separation in age-velocity space between the SCES and BSS / DES groups, just as we see in our model. This also con-  \nfirms detectability of SCES runaways using age and velocity. Stoop et al. (2024a) suggest that the anisotropic runaways were ejected by an interaction with another cluster. Our results further suggest that these ejected stars are the tidal tails produced in the tidal disruption of a late-forming, infalling subcluster. Observations show that there is an ongoing merger between two subclusters within R136 (Sana et al. 2012). This confirms that R136 formed via hierarchical assembly and supports the case for the SCES as the ejection mechanism of its northern runaway star group.  \nFig. 7. Trajectories of runaway star group 1 originating from a subcluster merging into the central cluster, with color indicating the stellar velocity (light = v max, dark = v min). The isosurfaces show four values of the gravitational potential of both the stars and gas, | U g | = 10 47 , 48 , 49 , 50 erg, with dark to light from lowest to highest. As the runaway stars form, they are plotted with star markers. The total number of stars in the cluster is listed in each panel.  \n<!-- image -->  \nAnother example of possible SCES ejection is seen in the runaway OB stars ejected from the young massive cluster NGC 6618, which seem to have a preferential direction. In Fig. 6 of Stoop et al. (2024b), the directions of runaways are plotted showing 7 of 13 stars ejected in a < 90 · region of the sky. Furthermore, in Fig. 8 of Stoop et al. (2024b), there is a distinct grouping of stars in velocity and stellar age. These stars are ejected from a cluster ≤ 1 Myr old (Hanson et al. 1997; Ho ff meister et al. 2008;  \nPovich et al. 2009; Ramírez-Tannus et al. 2017), which indicates a dynamical ejection, as no stars have exploded as SNe yet. We argue that this group of runaway stars was ejected via subcluster merger rather than individual dynamical interactions.  \nFigure 12 of Drew et al. (2021) shows the ejection directions of runaway OB stars in NGC 360 and Westerlund 2, colored according to the time of their ejection. The runaway populations from both clusters display some anisotropy, with clusters of run-  \nFig. 8. Trajectories of runaway star group 2 originating from a subcluster merging into the central cluster, with color indicating the stellar velocity (light = v max, dark = v min). The dots in the top row (analogous to Fig. 6) represent star particles, and the star markers indicate the runaway stars in the group as they form. The bottom row (analogous to Fig. 7) shows the same runaway star trajectories, with the potential of gas and stars plotted as isosurfaces. The isosurfaces show four values of the gravitational potential of both the stars and gas, | U g | = 10 47 , 48 , 49 , 50 erg , with dark to light from lowest to highest.  \n<!-- image -->  \naways moving in roughly the same direction. While it is possible these have all been ejected by the DES, we argue the anisotropy indicates a subcluster merger as their origin.', '4.4. Runaway binary stars': 'Two dynamical binaries 3 were ejected via the SCES in our model, both in group 1. We note that our model did not include any primordial binaries, a topic we will explore in future work. The time evolution of their velocities and orbital separations from the system center of mass is shown in Fig. 10 along with the velocity of all SCES runaways. Both systems formed as extremely wide dynamical binaries. After passing through the central potential, they became significantly harder as orbital energy is lost during the dynamical encounter(s) in the dense core. This suggests that tight primordial binaries ejected via the SCES are unlikely to be ionized. On the contrary, we show that ejection through the SCES can further harden binaries. This needs to be confirmed by analyzing a similar model with primordial binaries.  \nThe hardening e ff ect occurred in both a circular and a slightly eccentric binary. One binary transitioned into a circular ( e = 0 . 089) orbit after passing through the central potential, gaining orbital angular momentum through dynamical encounters in the  \ndense core. The other was irregular, settling into a steady eccentric ( e = 0 . 276) orbit roughly ≈ 0 . 1 Myr after peak acceleration.  \nWith the expectation that binaries are preserved and hardened through this mechanism, the SCES has the potential to produce runaway binaries with a broad range of properties. This is because runaways produced by the SCES are a mostly uniform sampling of the stars formed in a subcluster.  \nObservations of runaway binaries are fairly rare. The binary fraction of observed OB runaway stars is ≈ 8% (Gies & Bolton 1986; Mason et al. 1998). Simulations of equal mass binarybinary collisions also found that ≃ 10% of O type runaways were binaries (Leonard & Duncan 1990). The ejection mechanism producing the observed runaway binaries is not always certain.  \nThe ejection mechanism is unknown for two confirmed runaway O-star binaries HD 14633 and HD 15137 4 (Gies & Bolton 1986). They both are in a tight ( P ≈ 15 , 30), eccentric ( e ≈ 0 . 7 , 0 . 48) orbits with low mass companions (1-3 M ⊙ ) (McSwain et al. 2007). The low-mass companions could be neutron stars, but neither pulsars nor X-ray emissions are detected, so ejection via BSS cannot be demonstrated 5 . This leaves the case of a low-  \n<!-- image -->  \nFig. 9. Histograms of runaway star directions from 1,000 di ff erent projections. Top: 2D histogram showing the 1,000 projections, sorted vertically and centered horizontally by peak number of runaway stars, N run, and wrapped around the x-axis. For every projection there is a strong peak. Looking at each row individually, a fainter peak can be seen to the right or left of the peak value. Bottom: Average value of the projections (black histogram) and the range from maximum to minimum (gray area), with values given by the left axis. The dashed line shows the cumulative sum of the runaway stars per bin, with values given by the right axis.  \n<!-- image -->  \ntar companion. Most primordial O star binaries with lowmass companions 6 have wide orbits ( a > 100 AU) (Sana et al. 2014; Moe & Di Stefano 2017). Given our results, it is possible that HD 14633 and HD 15137 originated as wide O-B pairs that were hardened and ejected via the SCES.  \nThe SCES ejects binary systems while they are young, and they most likely survive. This means binaries ejected via the SCES can potentially undergo two velocity kicks, with a second kick from the BSS following the SCES 7 . Furthermore, because the binaries are hardened during the SCES, the velocity kick from the BSS will be even higher. This could be a channel for producing hypervelocity stars (HVSs). HVSs are stars unbound from the Galaxy, which requires v ≳ 400 km s -1 in the Galactic rest frame depending on location and direc-  \nFig. 10. Radial velocity, v r , relative to the center of mass over time and the orbital properties of the SCES runaway singles and binaries. Top: Time evolution of all the runaways in groups 1 and 2. Middle and bottom: Velocity (solid) and semimajor axes, ai (dotted) of each member in the two runaway binary stars formed in our simulation. Both are from group 1. The dark lines are the primary stars, and the light lines are the companion stars. Orbital properties for each system are listed in the corresponding panel, calculated for the stable orbits after t = 0 . 85 Myr, which is indicated by the vertical black line.  \n<!-- image -->  \non. Tauris (2015) found the BSS capable of producing HVSs. For HVSs with masses M ⋆ = 0 . 9 , 3 . 5 , and 10 M ⊙ they found kick velocities in the Galactic rest frame up to v max g r f = 1 , 280 , 770 , and 550 km s -1 , respectively. These maximum velocities are only possible via the BSS for particularly favorable binary parameters, such as closer orbits.  \nBinaries ejected via the SCES become more favorable HVS progenitors via subsequent BSS: they are closer in orbit, further from the galactic disk, and already locally unbound. Therefore, SCES binaries are more likely to produce HVSs via a subsequent BSS than if a binary undergoes the BSS while still bound to its parent cluster. This two-step mechanism is a likely channel for producing HVSs with trajectories not pointing toward the supermassive black hole at the Galactic center 8 (Hills 1988).  \nThe runaway binaries produced in our model could not become HVSs, as their orbital velocities are too low. Rather, we argue that a close primordial binary could become further hardened by ejection through the SCES and thus acquire a fast  \nenough orbital velocity to become a HVS after the SN of its companion.', '4.5. Importance of initial conditions': 'Our initial conditions are simplified: a spherical cloud of almost uniform density with isotropic turbulence. In reality, GMCs forming star clusters are filamentary (see Heyer & Dame 2015; Klessen & Glover 2016; Hacar et al. 2023). Subclusters forming from spherical clouds are closer to the center of mass and more evenly distributed than those forming from filamentary clouds. Accounting for filamentary initial conditions could dramatically increase the e ff ectiveness of the SCES.  \nThe two subclusters that formed furthest from the center of the cluster in M6 were ejected. Subclusters that form late or far enough away that they approach the assembling cluster after most of the other subclusters have already merged will be most e ff ectively accelerated in the SCES. The subcluster must approach during the period when the central potential is still contracting. With more realistic filamentary initial conditions, there will be more subclusters formed as the dense gas will be more distributed. We therefore predict a much higher fraction of SCES runaways from star cluster models with more realistic initial conditions due to the increase in the number and infall distance of subcluster merger events.  \nIn future work, we plan to import GMCs formed in largescale galaxy simulations into torch to assess the extent to which realistic initial GMCs a ff ect the fraction of runaway stars formed via the SCES. This will give us a better idea of how many observed runaways can be attributed to the SCES versus the BSS or the DES. This issue highlights the necessity of using realistic initial conditions, particularly for realistic subcluster dynamics in star cluster formation models.', '4.6. High peak stellar density': 'The stellar encounter rate per unit volume γ in a star cluster is given by γ ∝ ρ 2 /σ where ρ is the stellar density and σ is the velocity dispersion (Verbunt & Hut 1987). Due to the high initial density and low virial parameter of the initial cloud for our M6 model, the star cluster collapses to a high stellar density, resulting in a high stellar encounter rate. We compare the half-mass density of our model to young ( < 10 Myr) Milky Way and Local Group star clusters in Fig. 11. We also indicate the time periods of ejection (peak velocity) for the two runaway groups.  \nThe stellar density of our model when the runaway groups are ejected is higher than the observed values. However, we note that this is the density just after core collapse. The cluster will relax to a lower stellar density as it continues to evolve, and, indeed, the stellar density has already begun to decrease at a steady rate following the ejection of group 1. Considering their ages (the youngest is 3 Myr old), the observed clusters used to make Fig. 11 are most likely already further post-collapse than our model. Also, aside from R136, none are as massive as our M6 model, which is more characteristic of clusters in starburst galaxies.  \nRegardless, the high stellar density in our model compared to observations implies that the interaction rate could be inflated. The fraction of stars ejected via the SCES in clusters that do not reach as high stellar densities during core collapse would likely be lower than in our model. Additional models with more selfconsistent initial conditions are needed to constrain the e ffi ciency  \nFig. 11. Half-mass stellar density of our model over time. Regions of observed ρ hm values in young clusters ( < 10 Myr) are indicated by horizontal gray shaded regions. The dashed border lines show extragalactic clusters, and the solid border lines show Milky Way clusters. The vertical shaded regions show the time windows when runaways are dynamically ejected at peak velocity in our M6 model. The observational data were taken from Portegies Zwart et al. (2010).  \n<!-- image -->  \nand frequency of SCES runaways in clusters of varying mass and density.', '5. Conclusions': "Using a star-by-star simulation of star cluster formation from a gas cloud, we have identified the SCES as a new channel for the origin of runaway stars. This scenario occurs when a subcluster forms late, after the rest of the cluster is mostly assembled, and the subcluster then falls into the contracting central potential, becoming tidally disrupted and ejecting the majority of the stars in the subcluster as runaways in tidal tails.  \nWe believe the reason this phenomenon had not been identified earlier is the oversimplification of initial conditions in many star cluster formation simulations. Using a smooth spherical cloud instead of a cloud more realistically structured by supersonic turbulence results in less energetic subcluster mergers since star formation is more centrally concentrated in a spherical cloud. This emphasizes the need for linking scales and using realistic clouds formed in galaxy formation simulations as initial conditions for models of cluster formation.  \nTo determine the fraction of runaways formed by the SCES, BSS, and DES, we must import a self-consistently formed GMC to model the hierarchical assembly correctly. We also need to include primordial binaries, as they are essential for ejection via the DES or BSS. Furthermore, the simulations must be run for long enough to account for long-term dynamics and the timing of SNe. The M6 simulation in this study was only run for 1 . 35 t ff , equivalent to less than 1 Myr, and therefore no SNe occurred during this time.  \nThe runaways formed by the SCES are highly correlated with respect to their velocities, ages, and ejection directions. Finding groupings of runaway stars with similar values of these  \nthree properties is a straightforward method for detecting runaways formed from the same subcluster merger. If anisotropy is observed in a group of runaway stars, this indicates a distinct subcluster merger history. Conversely, if the runaways are isotropic, this indicates a mild cluster assembly history and ejection primarily via the DES or BSS channels. A multimodal spatial or velocity distribution of runaway stars could indicate multiple distinct subcluster mergers, as seen in our model. Depending on the resolution of observations, these groups of runaways can be used as fossils to trace the assembly history of a star cluster.  \nAn important caveat for confirming an anisotropic distribution of runaway stars concerns the number of detectable runaways: if a cluster ejects a small number of runaway stars isotropically, even fewer will be massive and detectable. This could falsely imply anisotropy in runaway star ejection directions. A su ffi ciently large sample of massive runaway stars is needed to exclude this possibility and confirm the anisotropy of the ejection direction. On the other hand, if grouped runaways in the small sampling have the same velocities and kinetic ages, it is likely that they were ejected by the SCES.  \nWe surveyed observational work and found several anisotropic populations of runaway stars. Stoop et al. (2024a) show that the group of runaways moving north of R136 has properties consistent with ejection via the SCES. We also find runaway groups in other clusters that could be produced by the SCES, though more analysis needs to be done to confirm this. Regardless, observations of runaway stars must be looked at through an additional lens of possible SCES origins.  \nThe SCES is capable of producing runaway binaries. Two wide binaries in our model were ejected, and through ejection, their orbits hardened significantly. This suggests that primordial binaries can survive ejection and end up with harder orbits. The SCES can thus explain observations of runaway binaries. This scenario also allows for a subsequent BSS ejection, which will be more energetic than a birthplace BSS as the binary progenitor is already unbound and in a hardened orbit. A two-step SCESBSS ejection is therefore a potential production mechanism for HVSs unbound from the Galaxy.  \nAcknowledgements. We thank the referee Mark Gieles for his insightful comments, one of which led to Sect. 4.6. We thank Yuri Levin for pointing out runaway binary stars as a potential signature of SCES. B.P. was partly supported by a fellowship from the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). M.M.M.L., B.P., and E.P.A were partly supported by NSF grants AST18-15461 and AST23-07950. E.P.A. also received support from NASA Astrophysical Theory grant 80NSSC24K0935. C.C.-C. is supported by a Canada Graduate Scholarship - Doctoral (CGS D) from the Natural Sciences and Engineering Research Council of Canada (NSERC). This work used Stampede 2 at TACC through allocation PHY220160 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants 21-38259, 21-38286, 21-38307, 21-37603, and 2138296. The code development that facilitated this study was done on Snellius through the Dutch National Supercomputing Center SURF grants 15220 and 2023 / ENW / 01498863. S.M.A. acknowledges the support of NSF grant AST2009679. R.S.K. and S.C.O.G. acknowledge financial support from the European Research Council via the ERC Synergy Grant 'ECOGAL' (project ID 855130), from the Heidelberg Cluster of Excellence (EXC 2181 - 390900948) 'STRUCTURES', funded by the German Excellence Strategy, and from the German Ministry for Economic A ff airs and Climate Action in project 'MAINN' (funding ID 50OO2206). The team in Heidelberg also thanks The Länd and the German Science Foundation (DFG) for computing resources provided in bwHPC supported by grant INST 35 / 1134-1 FUGG and for data storage at SDS@hd supported by grant INST 35 / 1314-1 FUGG. The authors acknowledge Interstellar Institute's program 'II6' and the Paris-Saclay University's Institut Pascal for hosting discussions that nourished the development of the ideas behind this work.", 'References': "Baczynski, C., Glover, S. C. O., & Klessen, R. S. 2015, Monthly Notices of the  \nRoyal Astronomical Society, 454, 380 Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362 Blaauw, A. 1961, Bull. Astron. Inst. Netherlands, 15, 265 Blaauw, A. & Morgan, W. W. 1954, ApJ, 119, 625 Brown, W. R., Geller, M. J., & Kenyon, S. J. 2014, ApJ, 787, 89 Colella, P. & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174 Cournoyer-Cloutier, C., Sills, A., Harris, W. E., et al. 2023, MNRAS, 521, 1338 Cournoyer-Cloutier, C., Tran, A., Lewis, S., et al. 2021, MNRAS, 501, 4464 Drew, J. E., Monguió, M., & Wright, N. J. 2021, MNRAS, 508, 4952 Federrath, C., Banerjee, R., Clark, P. C., & Klessen, R. S. 2010, The Astrophys-  \n- ical Journal, 713, 269\n- Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJs, 131, 273\n- Fujii, M., Iwasawa, M., Funato, Y., & Makino, J. 2007, PASJ, 59, 1095\n- Fujii, M. S. & Portegies Zwart, S. 2011, Science, 334, 1380\n- Gies, D. R. & Bolton, C. T. 1986, ApJS, 61, 419  \nGoodwin, S. P., Whitworth, A. P., & Ward-Thompson, D. 2004, A&A, 414, 633  \nGrudi'c, M. Y., Guszejnov, D., Hopkins, P. F., et al. 2018, MNRAS, 481, 688  \nHacar, A., Clark, S. E., Heitsch, F., et al. 2023, in Astronomical Society of the  \nPacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inut-  \nsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 153 Hanson, M. M., Howarth, I. D., & Conti, P. S. 1997, ApJ, 489, 698 Heyer, M. & Dame, T. M. 2015, ARA&A, 53, 583 Hills, J. G. 1988, Nature, 331, 687 Ho ff meister, V. H., Chini, R., Scheyda, C. M., et al. 2008, ApJ, 686, 310 Hoogerwerf, R., de Bruijne, J. H. J., & de Zeeuw, P. T. 2000, ApJ, 544, L133 Klessen, R. S. & Glover, S. C. O. 2016, in Saas-Fee Advanced Course,  \nVol. 43, Saas-Fee Advanced Course, ed. Y. Revaz, P. Jablonka, R. Teyssier, & L. Mayer, 85  \nKolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301 Kroupa, P. 2002, Science, 295, 82 Leonard, P. J. T. & Duncan, M. J. 1990, AJ, 99, 608 Lewis, S. C., McMillan, S. L. W., Low, M.-M. M., et al. 2023, ApJ, 944, 211 Lucas, W. E., Rybak, M., Bonnell, I. A., & Gieles, M. 2018, MNRAS, 474, 3582 Mason, B. D., Gies, D. R., Hartkopf, W. I., et al. 1998, AJ, 115, 821 McSwain, M. V., Ransom, S. M., Boyajian, T. S., Grundstrom, E. D., & Roberts,  \n- M. S. E. 2007, ApJ, 660, 740\n- Miyoshi, T. & Kusano, K. 2005, Journal of Computational Physics, 208, 315 Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15  \nOh, S. & Kroupa, P. 2016, A&A, 590, A107  \nPelupessy, F. I., van Elteren, A., de Vries, N., et al. 2013, A&A, 557, A84  \nPflamm-Altenburg, J. & Kroupa, P. 2010, MNRAS, 404, 1564  \n- Polak, B., Mac Low, M.-M., Klessen, R. S., et al. 2023, arXiv e-prints, arXiv:2312.06509\n- Portegies Zwart, S. & McMillan, S. 2018, Astrophysical Recipes, 2514-3433 (IOP Publishing)\n- Portegies Zwart, S., McMillan, S., Harfst, S., et al. 2009, New Astronomy, 14, 369\n- Portegies Zwart, S., McMillan, S. L. W., van Elteren, E., Pelupessy, I., & de Vries, N. 2013, Computer Physics Communications, 184, 456  \nPortegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431  \nPortegies Zwart, S. F. & Verbunt, F. 1996, A&A, 309, 179  \n- Poveda, A., Ruiz, J., & Allen, C. 1967, Boletin de los Observatorios Tonantzintla y Tacubaya, 4, 86\n- Povich, M. S., Churchwell, E., Bieging, J. H., et al. 2009, ApJ, 696, 1278\n- Rahner, D., Pellegrini, E. W., Glover, S. C. O., & Klessen, R. S. 2017, MNRAS, 470, 4453\n- Rahner, D., Pellegrini, E. W., Glover, S. C. O., & Klessen, R. S. 2019, MNRAS, 483, 2547\n- Ramírez-Tannus, M. C., Kaper, L., de Koter, A., et al. 2017, A&A, 604, A78 Ricker, P. M. 2008, ApJS, 176, 293\n- Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444\n- Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15\n- Sana, H., Ramírez-Agudelo, O. H., Hénault-Brunet, V., et al. 2022, A&A, 668, L5\n- Stoop, M., de Koter, A., Kaper, L., et al. 2024a, Nature, subm.\n- Stoop, M., Derkink, A., Kaper, L., et al. 2024b, A&A, 681, A21 Tauris, T. M. 2015, MNRAS, 448, L6\n- Verbunt, F. & Hut, P. 1987, in The Origin and Evolution of Neutron Stars, ed. D. J. Helfand & J. H. Huang, Vol. 125, 187\n- Wall, J. E., Mac Low, M.-M., McMillan, S. L. W., et al. 2020, ApJ, 904, 192 Wall, J. E., McMillan, S. L. W., Mac Low, M.-M., Klessen, R. S., & Portegies Zwart, S. 2019, ApJ, 887, 62\n- Wang, L., Iwasawa, M., Nitadori, K., & Makino, J. 2020, MNRAS, 497, 536 Wilhelm, M. J. C., Portegies Zwart, S., Cournoyer-Cloutier, C., et al. 2023, MNRAS, 520, 5331"}