We simulate mergers between star clusters embedded within their natal giant molecular cloud. We extract initial conditions from cloudscale simulations of cluster formation and introduce different prescriptions for primordial binaries. We find that simulations that do not include primordial binaries result in a larger fraction of unbound stars than simulations that include a prescription for binaries based on observations. We also find a preferred direction of motion for stars that become unbound during the merger. Subcluster mergers within realistic gas environments promote binary disruption while mergers between idealized gasrich spherical clusters do not produce the same disruption. Binary systems with smaller semimajor axes are disrupted in simulations of subcluster mergers within their natal environment compared to simulations that do not include the realistic gas environment. We conclude that binary disruption and the production of an anisotropic distribution of unbound stars are the natural consequences of subcluster mergers during star cluster assembly.
['Young star clusters', 'Star forming regions', 'Star clusters', 'Stellar dynamics', 'Binary stars', '1833', '1565', '1567', '1596', '154', 'Astrophysics - Astrophysics of Galaxies']
Binary Disruption and Ejected Stars from Hierarchical Star Cluster Assembly
{'Binary Disruption and Ejected Stars from Hierarchical Star Cluster Assembly': 'Claude Cournoyer-Cloutier , 1, ∗ Jeremy Karam , 1, ∗ Alison Sills , 1 Simon Portegies Zwart , 2 and Maite J. C. Wilhelm 2  \n1 Department of Physics and Astronomy, McMaster University, 1280 Main Street West, Hamilton, ON, L8S 4M1, Canada 2 Sterrewacht Leiden, Leiden University, Niels Bohrweg 2, 2333 Leiden, The Netherlands', 'ABSTRACT': 'We simulate mergers between star clusters embedded within their natal giant molecular cloud. We extract initial conditions from cloud-scale simulations of cluster formation and introduce different prescriptions for primordial binaries. We find that simulations that do not include primordial binaries result in a larger fraction of unbound stars than simulations which include a prescription for binaries based on observations. We also find a preferred direction of motion for stars that become unbound during the merger. Sub-cluster mergers within realistic gas environments promote binary disruption while mergers between idealized, gas-rich spherical clusters do not produce the same disruption. Binary systems with smaller semi-major axes are disrupted in simulations of sub-cluster mergers within their natal environment compared to simulations that do not include the realistic gas environment. We conclude that binary disruption and the production of an anisotropic distribution of unbound stars are the natural consequences of sub-cluster mergers during star cluster assembly.  \nKeywords: Young star clusters (1833) - Star forming regions (1565) - Star clusters (1567) - Stellar dynamics (1596) - Binary stars (154)', '1. INTRODUCTION': "Star clusters form hierarchically within giant molecular clouds (GMCs, e.g. Lada & Lada 2003; Portegies Zwart et al. 2010). Simulations (e.g. Dobbs et al. 2022; Fujii et al. 2012, 2022; Grudi'c et al. 2018; Howard et al. 2018; Rieder et al. 2022) have shown that massive star clusters assemble through the merger of smaller subclusters. Observations also indicate that Orion (Fujii et al. 2022), Westerlund 2 (Sabbi et al. 2012; Zeidler et al. 2021) and R136 (Fahrion & De Marchi 2024; Fujii et al. 2012), among others, are assembling hierarchically. Simulations further suggest that sub-cluster mergers result in clusters that better match the density profiles (Fujii et al. 2012) of observed young massive clusters than isolated cluster formation, as well as providing a natural explanation for the light element abundance variations (Howard et al. 2019; Lah'en et al. 2024) observed in globular clusters. Sub-cluster mergers are also  \nCorresponding author: Claude Cournoyer-Cloutier  \ncournoyc@mcmaster.ca  \n- ∗ Claude Cournoyer-Cloutier and Jeremy Karam are co-first authors  \ndynamically rich processes, that can enhance the formation of new stars from compression of the surrounding gas (Fujii et al. 2022) and impart dynamical signatures in the stellar component of the clusters involved in the merger (e.g. Fujii et al. 2022; Karam & Sills 2024). Studying mergers between embedded star clusters therefore allows us to gain physical insights about the process of massive cluster formation.  \nMost stars and protostars in star-forming regions are part of a binary (or higher-order) system, with a multiplicity fraction increasing strongly with stellar mass (see Moe & Di Stefano 2017; Offner et al. 2023, for recent reviews). The interplay between stellar multiplicity and stellar clustering is non-trivial. Observations reveal that - for solar-mass stars, at least - stellar multiplicity in young clusters depends on cluster density: higher-density environments (e.g. Orion, Duchˆene et al. 2018) have fewer wide binaries than low-density environments (e.g. Taurus, Kraus et al. 2011). Recent simulations of embedded (Cournoyer-Cloutier et al. 2021) and gas-free (Torniamenti et al. 2021) cluster assembly further suggest that changes to populations of binaries take place while the cluster is forming but that binary properties are then stable following gas expul-  \nion. Similar trends with environment are observed in older clusters. Deacon & Kraus (2020) compare binaries in moving groups and open clusters, and find that stars formed in looser associations have more wide companions than stars formed in cluster-forming regions. For massive, dense globular clusters, stellar multiplicity is anti-correlated with cluster luminosity (and implicitly cluster mass, Milone et al. 2016), but not with central density. Taken together, those observations suggest that the differences in the populations of binaries are a result of a cluster's assembly process, rather than only its present-day properties.  \nAccounting for the high multiplicity fraction of O and B stars, which are short-lived and rare, is simultaneously both more complicated and more important. Massive OB stars regulate the subsequent star formation in their natal clouds via feedback: they produce winds, radiation, and supernova explosions that heat the nearby gas and drive it away from the central, starforming regions, thus preventing the formation of other stars within the cluster. Massive stars routinely escape from their natal cluster; for example, recent observations show OB-stars moving away from the young massive clusters M16 (Stoop et al. 2023), M17 (Stoop et al. 2024), and NGC 3603 (Kalari et al. 2019). The removal of massive stars from their birth cluster can have consequences on galactic scales, due to the enrichment and supernova feedback from stars into the interstellar medium (Andersson et al. 2020).  \nRunaway stars are often classified as stars with radial velocities v r ≳ 30 km/s. They can be produced through few-body dynamical interactions in a cluster, involving at least one binary (e.g. Poveda et al. 1967; Hoogerwerf et al. 2000; Fujii & Portegies Zwart 2011; Gvaramadze & Gualandris 2011). Through the interaction, the binary becomes more tightly bound and the energy lost is transferred to another star(s) as kinetic energy; if the amount of energy is large enough, the star's velocity may exceed the cluster's escape velocity and the star may become unbound from the cluster. Two populations of OB runaways are observed in 30 Doradus (Sana et al. 2022), attributed to binary disruption following a supernova (originally proposed by Blaauw 1961) and to few-body encounters within a star cluster. Simulations have shown that merger between sub-clusters may also enhance the production of runaway stars (Fujii et al. 2022; Polak et al. 2024). These simulations however did not include primordial binaries, and so could not compare the relative effects of sub-cluster mergers and few-body encounters on the production of runaway stars during cluster assembly.  \nAny understanding of star cluster formation therefore requires an understanding of how a population of massive binaries evolves during cluster formation. Massive OB binaries are formed within embedded clusters, and the population of massive binaries in clusters is less dynamically processed than that of massive stars in the field, which are often runaways. There is however growing evidence supporting significant hardening of massive binaries during the process of cluster formation (e.g. Sana et al. 2017; Ram'ırez-Tannus et al. 2021; Bordier et al. 2022). It is clear that there is some highly nontrivial interplay between cluster mergers, binaries, and runaway stars. Some previous studies have started exploring this interplay (e.g. Fujii et al. 2012), using idealized initial conditions.  \nIn this work, we utilize initial conditions drawn from a cluster formation simulation (Wilhelm 2024, Wilhelm et al. in prep.) 1 that includes a sub-cluster merger resulting in the production of runaway stars. The simulation only contains two sub-clusters, allowing us to cleanly isolate the merger. It was performed with Torch (Wall et al. 2019, 2020) and includes star formation and stellar feedback, but does not however have any population of primordial binaries 2 . We introduce primordial binaries to those initial conditions, and run a suite of simulations with different prescriptions for binaries. We run our simulations using only stellar dynamics and hydrodynamics. Those computational methods that are simpler and less expensive than those used in the Torch framework, allowing us to explore a larger range of parameters. In this paper, we study how sub-cluster mergers within giant molecular clouds influence a population of binary stars, and study the production of unbound stars during sub-cluster mergers. In Section 2, we describe our simulation methods and initial conditions. We present the results in Section 3 and discuss their implications for hierarchical cluster formation in Section 4. We summarize our key conclusions in Section 5.", '2. METHODS': 'We use initial conditions extracted from a simulation of star formation within a collapsing cloud with an initial mass of 10 4 M ⊙ , which naturally gives rises to a merger between two sub-clusters of stars with stellar masses of 1684.3 and 239.5 M ⊙ in 2603 and 444 stars. The merger takes place 2.13 Myr after the start of the simulation, and 0.66 Myr after the formation of the first star. The  \n)  \n(  \nFigure 1. Initial conditions from the full Torch run, zoomed in to show the merging clusters. The full domain width is 17.5 pc, and the central 7 pc are shown. The surface density shown is calculated using the resolution of the gas cells, 0.136 pc.  \n<!-- image -->  \nsnapshot we use as our initial conditions is shown in Figure 1. In Section 2.1, we list the physics included in the full simulation. In Sections 2.3 and 2.2, we describe the numerical methods used in our simulations, while the initial conditions are outlined in Section 2.4.', '2.1. Torch simulation': 'The simulation from which we draw our initial conditions was done using Torch (Wall et al. 2019, 2020) and is presented in detail in Wilhelm (2024). In the text below, we describe the relevant physics implemented in the Torch simulation, and highlight how they lead to realistic initial conditions for the simulations presented in this work.  \nThe Torch simulation uses magneto-hydrodynamics (MHD) coupled with stellar dynamics (including a treatment of dynamically-formed binaries), star formation, stellar evolution and feedback, and a sub-grid model for protoplanetary disks coupled through the Amuse framework (Portegies Zwart et al. 2009; Pelupessy et al. 2013; Portegies Zwart et al. 2013; Portegies Zwart & McMillan 2019). MHD is handled by the adaptive mesh refinement code Flash (Fryxell et al. 2000; Dubey et al. 2014), with a maximum spatial resolution of 0.136 pc. In the Torch simulation, star formation takes place in bound regions with high density gas ( ⩾ 3 . 819 x 10 -21 g cm -3 ) and converging flows. This results in the formation of subclusters with shapes consistent with observed embedded clusters (Cournoyer-Cloutier et al. 2023). The stars are sampled from a Kroupa (2001) initial mass function (IMF) from 0.08 M ⊙ to 100 M ⊙ . Stars more massive  \nthan 7 M ⊙ provide feedback to the simulation in the form of momentum-driven winds, ionizing radiation and radiation pressure (Wall et al. 2020), resulting in a nonsmooth distribution of gas as seen in real star-forming regions.', '2.2. N-body dynamics': 'We use the stellar masses, positions and velocities from the full Torch run directly. We handle stellar dynamics with Petar (Wang et al. 2020a), which is optimized to handle large numbers of binaries and few-body encounters. Petar handles long-range interactions with a Barnes-Hut tree (Barnes & Hut 1986, as implemented by Iwasawa et al. 2016), short-range interactions with a fourth-order Hermite integrator (Makino & Aarseth 1992), and stable binaries and close encounters with the slow-down algorithmic regularization method (Wang et al. 2020b, sdar ). The longest timestep dt soft of the simulation, used for long-range interactions, must be chosen in conjunction with the changeover radii r in and r out between short- and long-range interactions. We set dt soft to 1/100 th of the shortest orbital period we could have for a binary in the long-range interaction regime. We set r out = 1 . 2120 x 10 -2 pc, r in = 1 . 2120 x 10 -3 pc, and dt soft = 27 . 948 years for our simulations. We also set r bin = 100 au for the radius under which we use slow-down algorithmic regularization.', '2.3. Hydrodynamics': 'We use the smoothed particles hydrodynamics (SPH) code Gadget-2 (Springel 2005) with a particle mass of 0.01 M ⊙ leading to 6.86 × 10 5 particles in total. We convert the gas from Flash to Gadget-2 using the method described in Karam & Sills (2024). The gravitational force from the gas on the stars and vice-versa is calculated with BHT ree (based on Barnes & Hut 1986, as implemented by Jun Makino) and the Bridge (Fujii et al. 2007) scheme. We use a bridge timestep of 2 5 dt soft = 894 . 336 years for our simulations.  \nOur simulations do not include stellar feedback, stellar evolution, star formation or magnetic fields, in contrast with the Torch simulation from which we draw our initial conditions. It is however important to note that all stars in our simulations are younger than 3 Myr by the end of the runs, and that none of the stars would therefore have exploded in a supernova.', '2.4. Initial Conditions': "We summarize the different sets of initial conditions in Table 1 and in the text below. For our fiducial initial conditions (M0), we take the stars from the Torch simulation run at 2.13 Myr. The distribution of the  \nTable 1. Overview of our simulations. Columns: the name of each simulation, whether the simulation has a primordial binary prescription, the prescription used for primordial binaries, and the initial conditions for the gas and cluster shape (see text for more details). Different numbers in the run names correspond to different random seeds.  \nstars is shown in Figure 1. The fiducial initial conditions have 3047 stars for a total mass of 1923.8 M ⊙ . The stars are split into two clusters using a simple positional argument. For the runs in which we introduce binaries, each system (i.e. single star or binary) is placed so that its centre of mass is at the position of the star from M0 that has the mass closest to the system mass. This preserves the shape of the clusters, which is independent of the presence of binaries for embedded clusters (see Cournoyer-Cloutier et al. 2023). We also include two other runs without primordial binaries, M1 and M2 , where we randomly changed the mass of each star slightly from its mass in the Torch run but preserved each cluster's total mass.  \nWe use three models for binaries, based on the observations presented in Moe & Di Stefano (2017, for stars above ∼ 1 M ⊙ ) and Winters et al. (2019, for M-dwarfs). In the all binaries models (AB1, AB2, AB3), we sample the full distribution of companions in the Galactic field using the technique presented in CournoyerCloutier et al. (2021). In the inner binaries models (IB1, IB2, IB3), we sample a distribution of inner companions; for intermediate- and high-mass stars, with a large triple fractions, this shifts the distribution to shorter periods. For each model, we vary the random seed to get the three different samplings. In the wide binaries models (WB1, WB2, WB3), we impose a minimum period cut-  \nfollowing Ram'ırez-Tannus et al. (2021),  \nP min = 10 (5 ± 1) exp ( -t/ 0 . 19 +0 . 06 -0 . 04 ) +1 . 40 (1)  \nwhere P min is the minimum period in days and t is the cluster age in Myr. They derive this minimum period from the observed velocity dispersion of OB stars in young star clusters, which follows  \nσ 1D = (13 . 3 ± 1 . 1) t +(2 . 0 ± 0 . 1) (2)  \nWe pick three possible values for the age of the cluster: 0.66 Myr (for the formation time of the first star), 0.21 Myr (for the average age of the simulated stars at the start of our zoom-in simulations) and 0 Myr (since the mean formation time of the stars in the full simulation is slightly larger than 2.13 Myr). Using those ages, we modify the IB1 initial conditions, so that all systems made up of two OB stars (i.e. all systems where both stars have a mass of at least 2 M ⊙ ) have an orbital period of at least P min . We report the minimum period for an OB binary at the start of each simulation in Table 2, as P 0 OB , min . Although M0, M1 and M2 do not include any primordial binaries, they do include binaries formed through dynamics in the original Torch simulation. Those binaries have longer orbital periods than primordial binaries (see e.g Cournoyer-Cloutier et al. 2021, for a comparison between the properties of primordial and dynamically-formed binaries).  \nTo investigate the effects of the background gas on the binaries (as it was shown to be important for the merger itself by Karam & Sills 2024), we also simulate an idealized merger and an isolated cluster of the same mass, both without background gas. We simulate a merger of two Plummer spheres of gas and stars, with the same stellar and gas masses as the merging clusters (PM). The relative velocity and impact parameter are also kept fixed, as well as the total gas mass, but the gas and stars are initialized as pairs of relaxed Plummer (1911) spheres with densities consistent with young massive clusters, following Karam & Sills (2022). We also simulate an isolated Plummer sphere with the total gas and stellar mass (PB). Both simulations use exactly the same stars and binaries as the IB1 simulation. A visual comparison of the gas between the merger inherited from the full initial conditions and the merger of two Plummer spheres is shown in Figure 2.", '3.1. Binary disruption': 'The binary fraction as a function of time is shown in the upper panel of Figure 3 for simulations with primor-  \nFigure 2. Gas surface density for the merger inherited from the Torch simulation (left) and for the merger of two isolated Plummer spheres (right). The maximum surface density in the Plummer spheres simulation is about one order of magnitude higher than in the merger from the Torch simulation.  \n<!-- image -->  \ndial binaries. We calculate the binary fraction as  \nF b = B S + B (3)  \nwhere B is the number of stars with a bound companion within 10,000 au and S is the number of stars without a companion within 10,000 au. Binaries wider than 10,000  \nFigure 3. Binary fraction as defined in Equation 3 as a function of time for runs with primordial binaries. In merger runs with Torch initial conditions, about half of the total decrease in binary fraction takes place during the merger, which ends at 0.41 Myr.  \n<!-- image -->  \nau are not included in our statistics as they are very short-lived in dense stellar environments.  \nThe initial binary fraction is ∼ 27% in all simulations with binaries. In the merger simulations with initial conditions inherited from the full Torch simulation (all the IB, AB and WB runs), the final binary fraction is between ∼ 19% and ∼ 21%, and the binary fraction at the end of the merger ( ∼ 0.40 Myr) is around 23-24%. For each of these simulations, about half of the decrease in the binary fraction takes place during the merger. The binary fraction therefore decreases by the same amount over the first 0.40 Myr and the last 1.60 Myr of the simulation, which indicates more rapid changes during the merger itself. The largest decrease takes place during the first ∼ 0.1 Myr of the simulation for all simulations realized with those initial conditions. We also note that the WB runs, which have all the same stars as IB1 and all the same binary systems for primaries less massive than 6 M ⊙ , follow almost exactly the same decrease in binary fraction as IB1 during the merger. When the sampling for the binaries is changed, some scatter is introduced but the overall trend remains. The AB runs show the largest amount of scatter. Since the binaries are sampled from the parameter space of all companions rather than inner companions only, there is a larger range in the distributions of semi-major axes and bind-  \nFigure 4. Cumulative number of disrupted systems for IB1, PB and PM. The three simulations have the same initial population of binaries.  \n<!-- image -->  \nof the binaries, leading to more scatter in the fraction of disrupted systems.  \nIn the Plummer sphere merger (PM), the binary fraction decreases from ∼ 27% to ∼ 23% over the course of the simulation, with a steeper decrease in the first 0.9 Myr. This is the simulation that has the highest binary fraction at the end of the simulation. The binary fraction in the isolated Plummer sphere (PB) decreases from ∼ 27% to ∼ 22%. The binary fractions in PM and PB follow the same decrease in the first 0.9 Myr of the simulation, but diverge in the second half of the simulation.  \nWe plot the cumulative number of disrupted binaries as a function of time for all runs that use the IB prescription for primordial binaries, and present it in Figure 4. For merger simulations, most of the disruption takes place during the merger. It continues throughout the entire merger, but peaks in the early stages of the merger, during the first ∼ 0.1 Myr. We see a similar behaviour in runs with the other two binary prescriptions. More than half of the total disruptions take place during the merger, with roughly 100 binary systems being disrupted in the first 0.4 Myr. We also note that the distribution of disruptions as a function of time is almost identical for the IB runs during the merger, while they diverge at later times. Similar behaviour is observed in the AB and WB runs. The PM and PB runs also show  \nFigure 5. Semi-major axes of binary systems disrupted in IB1, PB and PM during the first 0.41 Myr of the simulation, which corresponds to the time required for the merger to complete in IB1. The three simulations have the same initial population of binaries.  \n<!-- image -->  \ndisruption at early times, although not to the extent observed in the IB, AB, and WB runs.  \nAll runs with primordial binaries and gas initial conditions inherited from a GMC-scale simulation show an excess of disruption during the merger. This excess is caused by the disruption of systems with smaller semi-major axes during the merger, as shown in Figure 5. We compare the distributions of semi-major axes for the binaries disrupted within the first 0.41 Myr of the simulation in IB1, PB and PM using a two-sample Kolomogorov-Smirnov (KS) test. We find that the semimajor axes of disrupted systems for PB and PM are consistent with being drawn from the same distributions. We are however confident at respectively 98.3% and 98.8% that the semi-major axes of the disrupted systems in IB1 are smaller than those in PB and PM.', '3.2. Unbound stars': "We present the unbound mass fraction for stars for all simulations in Figure 6. The fraction of unbound stars increases during the merger for all merger simulations (as expected from Karam & Sills 2022, 2024), while almost all stars remain bound to their host cluster in PB (the isolated Plummer sphere). A smaller fraction of stellar mass becomes unbound at early times in PM than in the other simulations with a primordial binary prescription, which behave very similarly to one  \nFigure 6. Mass fraction of unbound stars as a function of time, for all simulations.  \n<!-- image -->  \nanother. By the time the merger ends in PM , at 0.60 Myr, the mass fraction of unbound stars has reached the same value as in the other simulations with primordial binaries. The most important difference is between the runs with primordial binaries, and those without primordial binaries. M0, M1 and M2 all have a similar mass fraction of unbound stars at the end of the merger, which is about 2.5 times higher than the mass fraction of unbound stars in the runs with Torch initial conditions and primordial binaries. This suggests that the underlying physical process responsible for the production of unbound stars during sub-clusters mergers is the same that leads to the disruption of binaries. In the absence of primordial binaries, the energy from the merger is converted to kinetic energy of the stars, leading to an increase in the fraction of unbound stars. In the presence of primordial binaries, however, the energy from the merger is used to disrupt binaries and is not sufficient to also unbind stars from the cluster.  \nFor each unbound star, we calculate its threedimensional radial velocity v r relative to the cluster's center of mass. The distributions of velocities in the simulations with and without primordial binary distributions are similar. We label unbound stars with v r ≥ 30 km/s as runaways. We report the number of runaways N runaways in Table 2, along with the number of runaways more massive than 6 M ⊙ (N runaways , M ), the number of unbound stars more massive than 6 M ⊙ (N unbound , M ) and the number of runaways initially in  \nFigure 7. Projected direction of motion in the xy plane of all unbound stars in simulations with Torch initial conditions, at the end of the merger. The length of the arrow is proportional to the projected velocity and the origin of the arrow corresponds to the location of the star. Runaway stars are color-coded by simulation while stars with v r < 30 km/s are shown in light grey.  \n<!-- image -->  \nFigure 8. Mollweide projection of the direction of motion of all unbound stars in simulations with Torch initial conditions, at the end of the merger.  \n<!-- image -->  \nbinaries (N runaways ,b ). The only merger simulation that does not produce any runaway stars is AB2 , which however has two unbound stars with v r ≥ 25 km/s, just below our runaway cut-off. PB, on the other hand, does not produce any runways, and the maximum radial velocity for an unbound star is 12.8 km/s.  \nThe runaways produced during the merger for simulations with Torch initial conditions are shown in Figure 7, along with the unbound stars with v r ≤ 30 km/s. We also calculate the direction of motion of all unbound  \nstars and the runaways stars after the mergers. We show the positions and directions of motion of all runaways produced during the merger in Figure 7. There are no runaways moving towards the (-x, +y) quadrant, which is the original location of the less massive cluster. If we consider all unbound stars produced during the merger for all Torch runs, however, the distribution of directions of motions peaks in the -x direction. Comparing the distribution of angles to uniform distributions in ϕ and θ with a KS test, we are confident at > 99.9% that the direction of motion of the unbound stars is not uniform. We present the three-dimensional direction of motion of unbound stars in Figure 8. The directions of motions are clustered around an angle ( ϕ, θ ) = (60 o , 0 o ).  \nWe also report the binary fraction for stars that are unbound at the end of the simulation in Table 2. We calculate the binary fraction at the start of the simulation, at the end of the merger, and at the end of the simulation. The binary fraction for those stars decreases in all simulations, consistent with our expectation of disrupted binaries leading to unbound stars. We note that for simulations without primordial binaries, which only include wide, dynamically-formed bianries, the binary fraction of unbound stars at the end of the simulation is 0, while it ranges from 7 to 17% for all runs with primordial binaries (except PB).", '3.3. Outcomes for massive binaries': "We now consider specifically how massive OB binaries evolve throughout the simulations. We consider OB binaries to be systems with a primary mass above 6 M ⊙ (the lower mass threshold considered by Ram'ırezTannus et al. 2021) and companion mass above 2 M ⊙ (the lower mass limit for B-type main sequence stars). We are interested in two main features: can cluster mergers with background gas harden massive binaries, and can they disrupt them?  \nWe first turn our attention to the OB binary with the shortest orbital period in each simulation. The minimum orbital period for an OB binary at the beginning and at the end of each simulation is reported in Table 2. There are no significant changes to the orbital period of the OB star with the shortest period in IB1, IB2, AB1, AB2, PB or PM. In IB3, the OB binary with the shortest period at the end of the simulation originally had a an orbital period of 123 days, and was hardened to an orbital period of 32.9 days through an exchange at 1.92 Myr. In AB3, the OB binary that had initially had the shortest orbital period, at 40.2 days, got hardened to an orbital period of 34.7 days, likely through a few-body encounter.  \nIn simulations with initially wide massive binaries (WB1, WB2 and WB3), the minimum orbital period for OB binaries decreases significantly. In WB1, it decreases by a factor of > 2, going from ∼ 3100 days to ∼ 1440 days. The OB binary that has the shortest period at the end of the simulation however initially had an orbital period slightly longer than 13,000 days, and was hardened through an exchange at 1.23 Myr followed by a few-body interaction at 1.59 Myr. In WB2 and WB3, the OB binary that has the shortest orbital period at the end of the simulation also had the shortest possible period at the beginning of the simulation. Those periods were shorted by a factor of 7.4 and 25 in WB2 and WB3 respectively; we also note that the shortest period is shorter in WB3 than WB2 at the end of the simulation. This hardening is the expected behaviour from observations of OB binaries in young clusters (Ram'ırezTannus et al. 2021, orbital periods decrease during cluster formation). It however must be noted than none of the OB binaries in WB runs reaches an orbital period below 1000 days. Additionally, exchanges and interactions that significantly shorten the orbital period of OB binaries all take place well after the merger.  \nWe look for disrupted OB binaries in all simulations that include primordial binaries. In the runs with Torch initial conditions, all the disruption of OB binaries takes place after the merger. There is one early disruption in PM. We stress that energetic binaries can be disrupted early if they encounter other energetic binaries. If we calculate the hard-soft limit for the resultant clusters, we find values around 2 x 10 44 erg, with a maximum value of 3.10 x 10 44 erg. Almost all the disrupted OB binaries were more energetic than this limit: this suggests that the short-period OB binaries observed in young star clusters might not be fully representative of the primordial population, despite their high binding energies. We also note that OB binaries were disrupted in all merger runs with primordial binaries, but not in the isolated Plummer sphere.  \nOnly one massive star becomes unbound from its host cluster the during the merger: a star with mass 68 M ⊙ in M0 , the fiducial set of initial conditions. This star is the second most massive in the simulation. It forms a binary with semi-major axis ∼ 300 au with a star with mass 63 M ⊙ during the merger. It is then disrupted by an encounter with the most massive star in the simulation, which has a mass of 85 M ⊙ ; after the encounter, the 85 M ⊙ and 63 M ⊙ stars form a binary with semimajor axis ∼ 80 au. M0 , M1 and M2 produce the most massive unbound stars: by the end of the simulation, M0 has unbound stars with masses 22, 23, 50 and 68 M ⊙ , M1 has unbound stars with masses 49 and 66 M ⊙ ,  \nTable 2. Properties of unbound stars and OB binaries in the different simulations. Columns: run name, number of runaways, number of runways with mass above 6 M ⊙ , number of unbound stars with mass above 6 M ⊙ , number of runaways originally in binaries, fraction F b, i of stars unbound at the end of the simulation originally in binaries, fraction F b, m of stars unbound at the end of the simulation in binaries at the end of the merger, fraction F b, f of stars unbound at the end of the simulation in binaries at the end of the simulation, minimum initial period for an OB binary, minimum period for an OB binary at the end of the simulation.  \nand M2 has an unbound stars with mass 29 M ⊙ . They have all been identified as members of a wide binary in at least one snapshot, but those did not persist between snapshots, indicating that a few-body interaction took place. This illustrates how wide, dynamically-formed binaries can results in the ejection of massive stars or large numbers of ejected stars, due to their larger dynamically cross-sections. Indeed, the other two simulations without primordial binaries, M1 and M2 , also have the largest number of runaway stars produced during the merger, and during the full simulation. On the other hand, the mass of the most massive unbound star in any simulation with primordial binaries is ∼ 11 M ⊙ .", '4.1. Causes of binary disruption during hierarchical cluster assembly': "There are several differences between the simulations using the initial conditions from Torch and those using idealized Plummer models. We discuss the effects of the stellar and gas distributions below.  \nThe stellar density profiles for IB1, PB, and PM are very similar. All three simulations reach densities ≳ 10 4 M ⊙ pc -3 (and therefore number densities > 10 4 pc -3 ) in their central regions. If high stellar densities were the cause of the excess disruption we observe in the Torch runs, we would expect the mass or number density to be significantly higher in the Torch runs than in the idealized models, which is not the case. The gas distribution however differs strongly between the Torch runs  \n- which were extracted from a GMC simulation - and the idealized Plummer models. The central gas density is about an order of magnitude higher in the Plummer models than in the Torch models, as illustrated in Figure 2. Conversely, the gas distribution is more extended in the Torch models, and includes background gas.  \nThe more realistic gas distribution contributes in two ways to the excess disruption of binaries observed in the Torch models. First, the presence of background gas promotes the merger between the sub-clusters, as found by Karam & Sills (2024). This means that the stars from the smaller cluster are accreted unto the larger cluster much more quickly in the Torch simulations than in the Plummer spheres' merger. This rapid accretion means that binaries in either cluster are more likely to encounter other binaries and stars, and therefore have a larger chance to be disrupted. On the other hand, the Plummer spheres do not immediately merge, and the two clusters can be easily identified after the first passage. Second, the potential is shallower in the Torch runs, due to the more extended gas distribution. This allows binaries with higher binding energies (and therefore smaller semi-major axes) to be disrupted by allowing the stars to move apart more easily after an encounter. We note that the potential increases sharply during the first ∼ 0.2 Myr of the simulations, during which the disruption rate of binaries peaks in all Torch runs.", '4.2. Unbound stars in embedded star clusters': 'Recent observations have found a possible anisotropy in the distribution of runaway stars around the young massive star cluster R136 (Stoop et al. subm.). Using simulations of massive cluster formation, Polak et al. (2024) attribute this effect to the production of runaways from hierarchical cluster assembly inside a GMC. They find that an entire sub-cluster can become a system of runaway stars, due to the large tidal forces involved in the merger which completely destroys an incoming cluster. Our simulations however probe a lower cluster mass range, resulting in a smaller number of runaways, which is not sufficient to investigate whether those stars have a preferred direction. Considering the full population of unbound stars - of which runaways are only a subset - however allows us to conclude that stars ejected during sub-cluster mergers inside GMCs have a preferred direction of motion, even in a lowermass regime.', '4.3. Implications for massive cluster formation': "Our results offer several insights into the process of (massive) star cluster formation. First, the disruption of binaries during hierarchical cluster assembly within giant molecular clouds offers a natural explanation for the lower wide binary fraction found in the ONC than in other young star-forming regions (e.g. Duchˆene et al. 2018). As the ONC has been proposed to be the recent site of a merger between gas-rich sub-clusters (see e.g. Fujii et al. 2022), the disruption of the wide binaries during the merger process would be consistent with our findings. Howard et al. (2018) found that massive clusters acquire about half of their stellar mass via gas-rich mergers. Such a mode of cluster assembly could contribute to the low binary fractions observed in globular clusters (e.g. Milone et al. 2016).  \nBinary disruption during gas-rich or gas-driven subcluster mergers leading to a shift towards smaller semimajor axes could also explain the apparent discrepancy between the results of Cournoyer-Cloutier et al. (2021) and Parker & Meyer (2014) or Torniamenti et al. (2021) regarding the evolution of field-like population of binaries in star clusters. In Cournoyer-Cloutier et al. (2021), we found that a field-like distribution of semi-major axes, like the one used in the AB runs, shifts to smaller values during hierarchical star cluster formation. On the other hand, for a similar distribution of binaries, Parker & Meyer (2014) use fractal initial conditions for their stars, while Torniamenti et al. (2021) use stellar positions inherited from the gas distribution. Neither, however, included gas along with the N-body dynamics while investigating the binaries. Both found the distribution of semi-major axes statistically unchanged by  \nthe process of (gas-free) cluster assembly. Along with the results presented in this paper, this suggests that changes in the gas potential that coincide with subclusters merger play a critical role in disrupting binaries during cluster formation.  \nAnother important implication of our results concerns the hardening of massive binary systems. The density in the inner 0.3 pc of IB1 is roughly 1.2 x 10 4 M ⊙ pc -3 at the start of the simulation, and 8.6 x 10 3 M ⊙ pc -3 after 0.9 Myr, around the typical central densities for young massive clusters (Portegies Zwart et al. 2010). Despite those high densities, however, simulations that started without any close massive binaries did not succeed in forming any OB binary with an orbital period shorter than 1000 days, while the most conservative choice proposed by Ram'ırez-Tannus et al. (2021) in Equation 1, using a cluster age of 2 Myr, gives a period of 337 days, about a factor of 4 shorter than the shortest period OB binary we get in WB1. Although our simulated clusters explore a lower mass range than the clusters studied by Ram'ırez-Tannus et al. (2021), our results suggest that the orbital parameters of OB binaries can be modified on timescales shorter than the cluster formation timescale.", '5. SUMMARY': 'We have conducted simulations of stellar sub-cluster mergers with different realistic prescriptions for the initial distribution of binaries. The shapes and masses of the sub-clusters, as well as the background gas, were taken from larger-scale simulations of cluster formation within a giant molecular cloud. We find that massive binaries can be disrupted or undergo significant changes to their orbital periods over timescales shorter than the cluster-formation timescale, in all of our merger simulations. The observed distributions of OB binaries is likely not the same as their formation distribution even in very young clusters and associations.  \nWe also find that the merger process - and by extension, hierarchical cluster formation - lowers the binary fraction and disrupts a large number of binaries during the merger, in excess of what is observed for an isolated cluster or an idealized model of cluster merger. The disrupted systems also tend to have smaller semimajor axes than in idealized models. This excess disruption is attributed to the rapid merger driven by the background gas distribution, as well as the shallower potential. Sub-cluster mergers result in stars becoming unbound from their host cluster, with a larger mass fraction of unbound stars when primordial binaries are taken into account. Simulations without primordial binaries, in which the most massive stars tend to pair up  \nin a wide, dynamically-formed binaries, are capable of ejection higher-mass stars for the cluster mass and density regime we consider. We further note that stars that become unbound during the merger tend to move in the same direction.  \nWe conclude that the production of unbound stars with a preferred direction of motion is a natural consequence of sub-cluster mergers within a GMC, extending the results of Polak et al. (2024) to lower cluster masses. We further argue that viewing binary disruption as a by-product of sub-cluster mergers within a giant molecular cloud offers a natural explanation for the lower binary fraction observed in denser star-forming environments and in globular clusters, which form hierarchically through subsequent mergers.  \nCCC is grateful for the hospitality of Leiden University, where this work was started during a visit in June 2023. The authors thank the referee for comments that improved the clarity of the paper. The authors warmly thank Eric Andersson, Sabrina Appel, Mordecai-Mark Mac Low, Stephen McMillan and Brooke Polak for ongoing discussions about Torch . The authors also thank Veronika Dornan, William Harris and Marta ReinaCampos for helpful discussions. CCC and JK are grateful to Gwendolyn Eadie for very valuable insights regarding the choice of initial stellar distributions.  \nCCC is supported by a Canada Graduate Scholarship - Doctoral (CGS D) from the Natural Sciences and Engineering Research Council of Canada (NSERC). The visit to Leiden University was made possible by a Michael Smith Foreign Study Supplement (CGS MSFSS) held by CCC at the Max Planck Institute for Astrophysics in Summer 2023. JK and AS are supported by NSERC. This research was enabled in part by support provided by Compute Ontario (https://www. computeontario.ca/) and the Digital Research Alliance of Canada (alliancecan.ca) via the research allocation FT #2665: The Formation of Star Clusters in a Galactic Context. The Torch simulation used for initial conditions was run on Snellius through the Dutch National Supercomputing Center SURF grants 15220 and 2023/ENW/01498863.  \nSoftware: matplotlib (Hunter 2007), numpy (Harris et al. 2020), pynbody (Pontzen et al. 2013), yt (Turk et al. 2011)', 'REFERENCES': 'Andersson, E. P., Agertz, O., & Renaud, F. 2020, MNRAS,  \n```\n494, 3328, doi: 10.1093/mnras/staa889 Barnes, J., & Hut, P. 1986, Nature, 324, 446, doi: 10.1038/324446a0 Blaauw, A. 1961, BAN, 15, 265 Bordier, E., Frost, A. J., Sana, H., et al. 2022, A&A, 663, A26, doi: 10.1051/0004-6361/202141849 Cournoyer-Cloutier, C., Tran, A., Lewis, S., et al. 2021, MNRAS, 501, 4464, doi: 10.1093/mnras/staa3902 Cournoyer-Cloutier, C., Sills, A., Harris, W. E., et al. 2023, MNRAS, 521, 1338, doi: 10.1093/mnras/stad568 Deacon, N. R., & Kraus, A. L. 2020, MNRAS, 496, 5176, doi: 10.1093/mnras/staa1877 Dobbs, C. L., Bending, T. J. R., Pettitt, A. R., & Bate, M. R. 2022, MNRAS, 509, 954, doi: 10.1093/mnras/stab3036\n```  \nDubey, A., Antypas, K., Calder, A. C., et al. 2014, The International Journal of High Performance Computing Applications, 28, 225, doi: 10.1177/1094342013505656  \nDuchˆene, G., Lacour, S., Moraux, E., Goodwin, S., &  \nBouvier, J. 2018, MNRAS, 478, 1825,  \ndoi: 10.1093/mnras/sty1180  \nFahrion, K., & De Marchi, G. 2024, A&A, 681, A20,  \ndoi: 10.1051/0004-6361/202348097  \nFryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273, doi: 10.1086/317361  \nFujii, M., Iwasawa, M., Funato, Y., & Makino, J. 2007,  \nPASJ, 59, 1095, doi: 10.1093/pasj/59.6.1095  \nFujii, M. S., & Portegies Zwart, S. 2011, Science, 334, 1380,  \ndoi: 10.1126/science.1211927  \nFujii, M. S., Saitoh, T. R., & Portegies Zwart, S. F. 2012, ApJ, 753, 85, doi: 10.1088/0004-637X/753/1/85  \nWilhelm, M. J. C. 2024, Doctoral dissertation, Leiden University. https://hdl.handle.net/1887/3717680 Winters, J. G., Henry, T. J., Jao, W.-C., et al. 2019, AJ, 157, 216, doi: 10.3847/1538-3881/ab05dc  \nZeidler, P., Sabbi, E., Nota, A., & McLeod, A. F. 2021, AJ, 161, 140, doi: 10.3847/1538-3881/abd6ca'}