{'fraction_non_alphanumeric': 0.031076916879821153, 'fraction_numerical': 0.018569615919115427, 'mean_word_length': 4.688165883835491, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The impact of binaries on the evolution of star clusters from turbulent molecular clouds

Most of massive stars form in binary or higher-order systems in clumpy, sub-structured clusters. In the very first phases of their life, these stars are expected to interact with the surrounding environment, before being released to the field when the cluster is tidally disrupted by the host galaxy. We present a set of N-body simulations to describe the evolution of young stellar clusters and their binary content in the first phases of their life. To do this, we have developed a method that generates realistic initial conditions for binary stars in star clusters from hydrodynamical simulations. We considered different evolutionary cases to quantify the impact of binary and stellar evolution. Also, we compared their evolution to that of King and fractal models with different length scales. Our results indicate that the global expansion of the cluster from hydrodynamical simulations is initially balanced by the sub-clump motion and accelerates when a monolithic shape is reached, as in a post-core collapse evolution. Compared to the spherical initial conditions, the ratio of the 50% to 10% Lagrangian radius shows a very distinctive trend, explained by the formation of a hot core of massive stars triggered by the high initial degree of mass segregation. As for its binary population, each cluster shows a self-regulating behaviour by creating interacting binaries with binding energies of the order of its energy scales. Also, in absence of original binaries, the dynamically formed binaries present a mass dependent binary fraction, that mimics the trend of the observed one.


INTRODUCTION
Most stars form as members of clusters or associations, that present a clumpy spatial distribution and may also contain sub-structures (Larson 1995). Understanding the early evolution of these star-forming complexes is of fundamental importance for the comprehension of the properties of young (< 100 Myr) and open clusters (Portegies , where the presence of sub-structures and fractality is observed (e.g., Cartwright & Whitworth 2004;Sánchez & Alfaro 2009;Parker & Meyer 2012;Kuhn et al. 2019). Also, these systems are characterized by complex internal kinematics, such as sub-clump relative motions and mergers, cluster expansion, gas dispersal (Kuhn et al. 2019;Cantat-Gaudin et al. 2019) and rotation (Hénault-Brunet et al. 2012). In particular, gas dispersal due to stellar winds and supernova explosions drives the cluster out of dynamical equilibrium, leading to an expansion phase, where most stars become unbound and disperse into the field (Hills 1980;Goodwin & Bastian 2006;Baumgardt & Kroupa 2007;Pfalzner 2009). Some of these natal properties might even survive the successive evolution of the stellar system and leave an imprint on the observed properties of older, relaxed stellar clusters (e.g., they may contribute to the signatures of rotation visible in some globular clusters, van Leeuwen et al. 2000;Pancino et al. 2007;Bianchini et al. 2013;Kamann et al. 2018;Bianchini et al. 2018).
In the first phases of their life, the dynamical evolution of young stellar clusters is deeply influenced by their stellar and binary content, and vice versa. In particular, a large fraction of the most massive stars is part of binary and higher order systems (Moe & Di Stefano 2017) that can actively exchange energy and angular momentum with the host environment, thanks to the very high density (ρ ∼ 10 4 M pc −3 ) of the cluster core. On the one hand, original binary stars (i.e., stars that form as members of a binary system) contain a large reservoir of internal energy, that can be transferred to other stars in the host star cluster, through three-and multi-body encounters (e.g., Heggie 1975;Hut 1983), preventing or reversing the gravothermal collapse of the core of the cluster (Chatterjee et al. 2013;Fujii & Portegies Zwart 2014). On the other hand, the global evolution of the cluster affects the properties of the binary population: for example, core collapse leads to the formation of new binary systems and to their dynamical hardening (Spitzer & Hart 1971). On top of this, binary stars are also affected by mass transfer, common envelope, supernova kicks, tides and other evolutionary processes (e.g., Hut 1981;Webbink 1984;Portegies Zwart & Verbunt 1996;Hurley et al. 2002). All these processes are crucial for the ejection of stars from their host star cluster (e.g., runaway stars, Fujii & Portegies Zwart 2011), and for the formation of intermediate-mass black holes (e.g., Ebisuzaki et al. 2001;Portegies Zwart et al. 2004;Giersz et al. 2015;Mapelli 2016). Finally, the interplay between dynamical interactions and binary evolution (Banerjee et al. 2010;Ziosi et al. 2014;Banerjee 2017;Fujii et al. 2017;Di Carlo et al. 2020b;Kumamoto et al. 2019;Antonini & Gieles 2020;Trani et al. 2021) can explain the properties of the binary compact objects observed through gravitational wave detection by LIGO and Virgo (Abbott et al. 2020a,b).
Direct N −body simulations are usually adopted to integrate the collisional dynamics of gas-free star clusters, where length-scales of different orders of magnitude, from binary separations of some solar radii to several parsecs, need to be included. However, studies of this type often lack realistic initial conditions. For example, state-of-the-art direct N −body simulations of star clusters include realistic stellar mass functions and stellar evolution, but most of them start from spherical idealized models, such as Plummer (1911) or King (1966) models. In some recent work, fractal initial conditions were adopted to mimic the initial clumpiness of star clusters (e.g., Goodwin & Whitworth 2004;Schmeja & Klessen 2006;Allison et al. 2010;Küpper et al. 2011;Parker et al. 2014;Di Carlo et al. 2019;Daffern-Powell & Parker 2020). Few studies tried to re-simulate with a direct N −body code the initial conditions obtained from hydrodynamical simulations of star cluster formation (Moeckel & Bate 2010;Moeckel et al. 2012;Parker & Dale 2013;Fujii & Portegies Zwart 2015), but most of them do not include stellar evolution or realistic stellar mass functions or original binary populations. A recent attempt to couple magneto-hydrodynamics and direct N −body star cluster formation simulations, also considering the presence of original binaries, was proposed by Cournoyer-Cloutier et al. (2020), who developed a binary generation algorithm consistent with observations of mass dependent binary fraction and distributions of orbital periods, mass ratios and eccentricities. They found that binary systems formed dynamically do not have the same properties as the original ones, and that the presence of an initial population of binaries affects the properties of dynamically-formed binaries. An adequate modelling of the original binary population is thus necessary for a realistic description of dynamical interactions in the early stages of star clusters' evolution.
Recently, Ballone et al. (2020) and Ballone et al. (2021) proposed a new approach to connect hydrodynamics and stellar dynamics that can be used to provide more realistic initial conditions for direct N −body simulations. This approach includes a number of the ingredients necessary to self-consistently study this problem: realistic phase-space distributions of stars, drawn from sink particle distributions of collapsing molecular clouds, and a realistic stellar mass function, which is fundamental to assess the impact of stellar evolution. This method is based on the assumption that the gas, in which the newly formed star cluster is embedded, is almost instantaneously expelled by feedback (radiation, winds and, most of all, supernova explosions) from the young most massive stars (e.g., Vázquez-Semadeni et al. 2010;Dale et al. 2014;Pfalzner et al. 2014;Gavagnin et al. 2017;Chevance et al. 2020a,b;Pang et al. 2020). From that moment on, the evolution of the newly born stellar system is mainly driven by gravitational dynamics. A necessary step towards a more realistic description is the insertion of binary stars in the original stellar population.
The aim of this paper is to offer a realistic, self-consistent description of the complex interplay between binaries and their host cluster in the first phases of a cluster's life after gas expulsion, by considering the effects of dynamics, stellar and binary evolution simultaneously. To do this, we insert original binaries in the joining/splitting method introduced in Ballone et al. (2021), to generate realistic initial conditions for N −body simulations starting from hydrodynamical simulations. Also, we study the evolution of the phase-space distribution of star clusters generated by hydrodynamical simulations and we compare it to other, more idealized, initial configurations.
This paper is organized as follows. In Sect. 2, we introduce our binary generation algorithm. Section 3 describes the initial conditions of the N −body simulations. In Sect. 4, we report the results of the simulation of a stellar cluster under different evolutionary conditions and compare it to other initial phase-space distributions. In Sect. 5, we discuss the peculiar aspects of the evolution of the stellar clusters from hydrodynamical simulations. Finally, in Sect. 6 we report our conclusions.

Binary generation algorithm
We developed a new algorithm to generate a realistic initial mass function (IMF) and a realistic population of original binaries, based on observations (Sana et al. 2012;Moe & Di Stefano 2017). This algorithm can be easily coupled to different phase-space generation codes to obtain a variety of initial conditions for N −body simulations. The method consists of the following steps. a) First, the algorithm randomly draws a population of stars from a Kroupa (2001) IMF between 0.1 M and 150 M , for an assigned value of the total mass of the population. b) The stars are paired up to each other in order to obtain a distribution of mass ratios q = m2/m1 following Sana et al. (2012): The coupling is set to generate a binary fraction f bin =  N bin /(Nsing + N bin ), where N bin is the number of binary systems and Nsing is the number of single stars, which depends on the mass of the primary star, following the observational results of Moe & Di Stefano (2017). For simplicity's sake, we do not include triple systems, but we take into account their presence when evaluating the binary fraction by labeling a certain number of single stars as third components of the existing binary systems (following Moe & Di Stefano 2017). This results in a fraction of binaries counted as triples (ftrip), and prevents from having an excessive number of binary systems among the most massive stars. By this procedure, we obtain a distribution of single stars and of binary particles.
For this work, we assume the binary fraction goes to zero in the mass range 0.1 − 0.8 M : the observations indicate that the percentage of binary stars in this mass range is low anyway (Moe & Di Stefano 2017), and including these low-mass binary stars would have dramatically increased the computational cost of the simulations. The resulting binary fraction for stars with mass m > 0.8 M is 0.4. c) The single and binary particles are assigned a phase-space distribution by coupling the aforementioned algorithm to a phase-space distribution generator. For this work, we considered two choices of the phase-space distribution generator.
In the first case, we coupled our algorithm with the joining/splitting procedure summarized in the next sections and described in detail in Ballone et al. (2021). In the second case, the phase-space distribution is created with the code McLuster (Küpper et al. 2011). d) Finally, the binary particles are split into separate stars and their orbital period (P ) and eccentricity (e) distributions are generated following Sana et al. (2012): and where, for a given orbital period, we set the upper limit of the eccentricity distribution emax(P ) according to eq. (3) of Moe & Di Stefano (2017): The orbital properties of the binaries are then converted into positions and velocities by considering an isotropic distribution for the orbital planes. Figure 1 shows an example of the binary populations generated by means of this algorithm. These initial conditions can be used to study the evolution of binary stars at the early and later stages of their host stellar cluster's life with a great variety of initial configurations. In addition, the generation of initial conditions through this algorithm has negligible computational cost. Finally, the described procedure is also suited to generate initial conditions for population synthesis studies.

Hydrodynamical simulations
The star clusters studied in this work are obtained by applying our algorithm to the output of the hydrodynamical simulations of turbulent molecular clouds presented in Ballone et al. (2020) and Ballone et al. (2021). These hydrodynamical simulations are performed with the smoothed-particle hydrodynamics code gasoline2 (Wadsley et al. 2004(Wadsley et al. , 2017. For this work, we consider the hydrodynamical simulation initialized with a total mass of 2 × 10 4 M . The cloud has an initial uniform density of 250 cm −3 , an initial temperature of 10 K and it is in an initial marginally bound state, with a virial ratio αvir ≡ 2 T /|V | = 2, where T and V are the kinetic and potential energy, respectively. The turbulence consists of a divergence-free Gaussian random velocity field, following a Burgers (1948) power spectrum. The gas thermodynamics has been treated by adopting an adiabatic equation of state with the addition of radiative cooling (Boley 2009). Stellar feedback was not included. Star formation is implemented through a sink particle algorithm adopting the same criteria as in Bate et al. (1995).
At 3 Myr (for a discussion of this choice see Ballone et al. 2021), we instantaneously remove all the gas from the simulations, mimicking the impact of a supernova explosion. We apply the joining/splitting algorithm to the properties of the sink particles at 3 Myr, as detailed in the next sub-section. We refer to Ballone et al. (2020) and Ballone et al. (2021) for more details on the hydrodynamical simulations.
2.3 The joining/splitting algorithm Ballone et al. (2021) introduced a new algorithm to generate stellar populations from sink particles obtained through hydrodynamical simulations. This algorithm consists in either joining or splitting the sink particle masses, which are affected by non-physical effects (such as the simulation resolution and the adopted sink particle algorithm), so to obtain a new, more realistic mass function of "children" stars. In this way the obtained stellar population inherits the turbulent phase space distribution generated from hydrodynamical simulations, but features a realistic mass function. Here we summarize the main steps of the joining/splitting process.
First, a population of stars with a chosen IMF is created, for an assigned value of their total mass. The joining algorithm is used when a star is more massive than the most massive sink particle. According to the joining algorithm, we select the densest region of the sink particle distribution and merge the neighbour sinks until we obtain the mass of the star. The position and the velocity of the star are assigned as the position and the velocity of the centre of mass of the joined sinks. The joining algorithm tends to enforce mass segregation in the central regions of the simulated star clusters.
The splitting branch of the algorithm, instead, is applied if a massive sink is more massive than any left star. In this case, we subtract the mass of individual stars from the massive sink particle, until a mass smaller than 0.1 M is left. The leftover mass is reassigned to the closest sink, so to enforce local and total mass conservation. The children stars of each sink particle are then distributed around the position and velocity of their parent sink according to a virialized Plummer distribution (for this step we make use of the new plummer model module in amuse, Pelupessy et al. 2013). In Ballone et al. (2021), we considered a Plummer half-mass radius of 10 −3 pc, that allowed a good energy and virial ratio conservation for all the hydrodynamical simulations of the sample. For this work, we prefer a Plummer half-mass radius of 10 −2 pc because, for this specific star cluster, this choice allows a better conservation of the total energy and a smaller variation of the virial ratio. The process of joining/splitting is cycled until either all the sink particles or the stars are consumed.

Direct N -body simulations
For our direct N −body simulations, we made use of the direct summation N −body code nbody6++gpu (Wang et al. 2015) coupled with the population synthesis code mobse (Mapelli 2017;, an upgraded version of bse (Hurley et al. 2000(Hurley et al. , 2002. nbody6++gpu implements a 4thorder Hermite integrator, individual block timesteps (Makino & Aarseth 1992) and a Kustaanheimo-Stiefel regularization of close encounters and few-body subsystems (Stiefel et al. 1965;Mikkola & Aarseth 1993). A neighbour scheme (Nitadori & Aarseth 2012) is used to compute the force contributions at short time intervals (irregular force/time steps), while at longer time intervals (regular force/time steps) all the members in the system contribute to the force evaluation. The irregular forces are evaluated using CPUs, while the regular forces are computed on GPUs using the CUDA architecture. The force integration includes a solar neighbourhoodlike static external tidal field (Wang et al. 2016). In all our cases, we consider a star as an escaper if it reaches a distance from the centre of density greater than four times the tidal radius of the cluster. The value chosen for the removal distance avoids the presence of potential escapers in the calculation (Takahashi & Baumgardt 2012;Moyano Loyola & Hurley 2013). mobse includes up-to-date prescriptions for massive star winds , for core-collapse supernova explosions (Fryer et al. 2012; and for pair instability . nbody6++gpu and mobse are integrated, as described by Di Carlo et al. (2019Carlo et al. ( , 2020a.

INITIAL CONDITIONS FOR N -BODY SIMULATIONS
The initial conditions for the N −body simulations from hydrodynamical simulations (hereafter labeled as Hydro) are obtained by combining the binary generation algorithm described in Sect. 2.1 and the the joining/splitting procedure (Sect. 2.3). The main properties of the initial conditions for the star cluster are reported in Table 1. The system has a total mass of Mtot = 6687 M , a half-mass radius (defined as the 50% Lagrangian radius centred in the centre of density) r50 = 1.70 pc, and a core radius (defined as the 10% Lagrangian radius centred in the centre of density) r10 = 0.06 pc. After the instantaneous removal of the gas, the system is left in a super-virial state, with αvir ≡ 2 T /|V | = 1.53. In order to quantify the impact of different physical ingredients on the dynamical evolution of a stellar cluster, we take into account four different evolutionary cases: • Bin: evolution with original binary stars and without stellar evolution. • Bin+SE : evolution with original binary stars and with stellar evolution. We assumed solar metallicity (Z = 0.02, Anders & Grevesse 1989), in order to match the young star clusters of the Milky Way (Portegies ) and to maximize the difference with respect to the case without stellar evolution, because mass loss by stellar winds is extremely high at solar metallicity (e.g., Vink et al. 2001;Kudritzki 2002). • NoBin: case with no original binary stars and no stellar evolution. • NoBin+SE : case without original binary stars but with stellar evolution.
The comparison between the aforementioned four different cases allows us to have a complete view of the impact of binaries and of stellar evolution on the dynamical evolution of a cluster with a realistic phase-space distribution of stars.
For each case, we ran 10 simulations of different joining/splitting realizations in order to filter out stochastic fluctuations.

Comparison with other initial conditions
We compared the evolution of the Hydro initial conditions to that of other initial phase-space distributions, which are commonly used in studies of star cluster dynamics. In order to have a fair comparison, we set initial conditions that match the mass scale and either the central length scale (r10) or the global length scale (r50) of our Hydro clusters. All the initial conditions are generated by coupling our binary generation code to McLuster (Küpper et al. 2011) as descibed in Section 2.1. We considered three cases: • King: a King (1966) model matching the core radius of the hydrodynamical initial conditions. We generated a King model with a half-mass radius of r50 = 0.25 pc and a central concentration of W0 = 9. The chosen value for the central concentration is high, typical of clusters that are believed to have undergone core-collapse. For this reason a post-core collapse evolution may be expected for both this case and for the central regions of the hydrodynamical case. • Loose Fract: a fractal sphere, with the same total mass and half-mass radius as the Hydro case. For this case, we selected a fractal dimension D = 1.6. • Dense Fract: a fractal sphere with the same mass and fractal dimension as the previous case, but with a half-mass radius set according to the Marks & Kroupa (2012) relation: In this case, we have r50 ≈ 0.3 pc and r10 ≈ 0.07 pc. Interestingly, the core radius results very similar to that of the Hydro initial conditions.
For all these initial conditions we set the same virial ratio as the Hydro case. The physical properties for all the initial conditions are summarized in Table 1.

Initial clumpiness of the stellar cluster
The initial space distribution of the Hydro simulation is clumpy and sub-structured, as can be seen in Fig. 2. The stellar cluster mainly consists of two very dense main subclumps and some minor and irregular clusters and filaments. We first defined the two main sub-clumps by using the dbscan (Density-Based Spatial Clustering of Applications with Noise) algorithm (Ester et al. 1996) 1 . This algorithm allows to group together points in high-density regions: these are labeled as core points, and they are distinguished from points in low density areas, that are labeled as noise The result of 1 The implementation we referred to is that of the python library Scikit-learn (sklearn.cluster.DBSCAN, Pedregosa et al. 2011). dbscan requires to define two parameters, and minP ts. The parameter minP ts is the number of points within the reference distance needed for a point to be considered as a core point. Otherwise, it is labeled as noise. For our case, we set these parameters based on the half-mass radius of the cluster and on the total number of stars: minP ts = Ntot/10 and = r 50 /5. the clustering procedure is shown in the top panel of Fig. 2: the algorithm manages to identify the two main sub-clusters. The main sub-clump presents a mass of Msc ≈ 2304 M (35% of the total mass) and a half-mass radius of rSC = 0.15 pc, while the second sub-clump has a mass of Msc ≈ 1132 M (17% of the total mass) and a half-mass radius of rSC = 0.16 pc. We checked if the sub-clump masses and halfmass radii are consistent with eq. (5), that is the relation between total mass and half-mass radius found in star-forming cloud cores by Marks & Kroupa (2012). Recently, Fujii et al. (2021) found that this relation holds in N -body/SPH simulations for embedded clusters with mass up to about 10 3 M and it is preserved after gas expulsion. In this case, we consider the sub-clumps found in the sample of stellar clusters considered by Ballone et al. (2021), that extend to higher masses, between 10 3 and 10 4 M . As the lower panel of Fig.  2 shows, this sample is well consistent with eq. (5).

Global evolution
4.2.1 Early evolution (t < 1 Myr) Figure 3 shows the very first phase (t ≤ 1 Myr) of the evolution for one representative cluster. At t = 0 Myr, the centre of density is located well within the main clump, while the second main sub-clump is out of the sphere defined by the half-mass radius. At t = 0.5 Myr, the cluster structure has significantly evolved. On the one hand, at small scales, each sub-clump rapidly expands, as a consequence of the instantaneous gas removal, thus lowering its local density. On the other hand, the two main sub-clumps get closer to each other, thus balancing the small scale expansion on a larger scale. These competing mechanisms characterize the first ≈ 1 Myr of the simulation.
At t = 1 Myr the cluster has nearly a monolithic shape. The half-mass radius is slightly larger (r50 ≈ 2 pc) than at the beginning of the simulation (when r50 ≈ 1.7 pc), while the core radius has grown much faster, as can be easily seen from Fig. 3. Typically, a realization reaches a monolithic shape after 1 − 1.5 Myr (only in a limited number of cases, this condition is fulfilled at about 2 − 2.5 Myr), after a short period in which the two sub-clumps tidally interact without merging. The resultant cluster presents an elongated shape, as a consequence of the strong tidal interaction and the relative motion between the sub-clumps.
The range of merger timescales is in agreement with the results by Fujii (2015), whose simulations can simultaneously reproduce the properties of different types of young star clusters, from massive and dense ones to open clusters and looser OB associations. In this sense, when N -body simulations are exploited to study the early evolution of stellar clusters, the timescale of sub-clump mergers is strongly dependent on the initial energetic state of the molecular cloud, as can be inferred by comparing the results in Fujii & Portegies Zwart (2014) and Ballone et al. (2021), who initialized their clouds in a marginally bound state. On the observational side, this kind of mergers between sub-clumps seems to be disfavoured to explain the formation of young star clusters like NGC 3603 (Banerjee & Kroupa 2013), whose observational properties require either a monolithic formation channel or a prompt assembly in t < 1 Myr. However, the results by Sabbi et al. (2012) hint that ongoing mergers between very young clus-  ters (such as R136 and the Northeast Clump in NGC 2070) may also occur.

Cluster expansion
In order to consider both the initial clumpy evolution and the successive monolithic expansion, we evolved the clusters for 10 Myr. Figure 4 shows the expansion of the cluster, described by r50 and r10, for all the four evolutionary cases. As a consequence of the mechanism described in Section 4.2.1, the half-mass radius initially grows, reaches a peak at about 0.5 Myr, that is when the secondary sub-clump enters the sphere of the half-mass radius of the main sub-clump, and then decreases. At 1 Myr, r50 reaches a minimum and then grows monotonically. The expansion of r50 is no longer influenced by the relative sub-clump motion, which at this time have merged or are very close to each other, but is due to the small scale expansion that has now reached larger scales. In contrast, the core radius grows rapidly since the very beginning of the simulation, because the sub-clump motion has no effect at these small scales.
The impact of binary stars is evident in the second phase of the evolution of the cluster, during the monolithic expansion. In fact, clusters with original binaries expand faster after 1 Myr: at this point the large-scale interaction of the subclumps is no longer present, and the density in the central region is still high enough (of the order of 10 3 −10 4 M pc −3 ) to allow efficient interactions and energy exchange between the binary stars and their surrounding environment. In the very first phases, instead, the faster expansion due to binary stars is balanced by the global evolution of the sub-clumps.
As explained in Sect. 3.1, the central regions of the cluster are matched by a King model with W0 = 9, that is typical of stellar clusters that are thought to have undergone corecollapse. We thus compared the monolithic expansion of the cluster with that expected based on a self-similar evolution, at constant mass (Spitzer 1987): where B is a proportionality constant. If the evolution of the cluster is a post-core collapse expansion, the time increase of r50 should be roughly consistent with eq. (6). We performed a fit to the median values of r50 from 1 Myr curves by using eq. (6). The resulting best-fit curves are the dashed lines in the left-hand panel of Fig. 4. We show the curves for the cases Bin and NoBin, where the lack of stellar evolution should avoid the presence of additional effects (e.g., mass loss) and make the dynamical effect by binaries more evident. The curves of both cases seem to be consistent with a post-core collapse phase until 10 Myr.

Mass loss
As the cluster expands, stars get further away from its centre, until they are eventually removed from the cluster dynamics by the tidal field of the host galaxy. This makes the total mass of the stellar system decrease. The presence of binary stars enhances the number of escaping stars, by powering a faster expansion. Also, close interactions between binary stars and single (or other binary) stars may lead to the ejection of stars, and possibly also of binary systems. In addition, stellar evolutionary processes (e.g. stellar winds, supernova explosions) make single stars, and thus the cluster, lose mass. Figure 5 shows the variation of the total mass of the cluster. Stellar evolution gives the main contribution to mass loss in the early stages of the simulation, resulting in a steeper slope of the mass evolution. After 10 Myr, the mass loss in the cases with stellar evolution is twice as large as that in the cases without stellar evolution. The absence of primordial binaries delays the mass loss, because the cluster needs to form its binaries dynamically before they start ejecting other stars. Figure 6 shows the evolution of the total kinetic energy (E k ), the total potential energy (Eg) of the centres of mass and the total binding energy of binary systems (E b ). Binary stars produce an initial sharp increase of the kinetic energy by yielding their internal energy to the surrounding stars. This results in the fast cluster expansion seen in Fig. 4. After this initial sharp increase, the kinetic energy of the clusters with original binaries decreases at a fast rate as a consequence of the ejection or evaporation of high velocity stars. After the first ∼ 5 Myr, the kinetic energy of the star clusters with original binary systems becomes similar to that of the other clusters, and they evolve in the same way for the rest of the simulation.

Energy variation
The total binding energy of the initial binary population is much higher than the typical gravitational energy of the centres of mass. Our original binary stars are, in fact, mostly hard and a small fraction of their total internal energy is sufficient to deeply affect the evolution of the cluster. The decrease of the total binding energy springs from two factors. Firstly, some binary stars escape from the system. This causes the slow decrease of the black line in Fig. 6. Secondly, stellar and binary evolution tend to remove binary stars from the population, via mergers, supernova explosions but also direct collisions between stars. This process is important since the very first stages, because the binary fraction is very high for the most massive stars and because the initial semi-major axes from Sana et al. (2012) are skewed to small values. By comparing the Bin models and the Bin+SE models, one can infer that this second factor is the main responsible for the variation of the total binding energy. If there are no original binary systems (NoBin and NoBin+SE models), the cluster creates its own population, with binding energies of the order of the gravitational energy scale. The case without stellar evolution is characterized by a monotonic increase of the binding energy, where binaries form and the hardest ones tend to harden. In the end, the total binding energy is dominated by the binding energy of a very few binaries. In presence of stellar and binary evolution, after an initial increase, the total binding energy decreases when stellar and binary evolution processes take over.

Binary populations
In order to understand how binary populations evolve and interact with the host cluster, we must estimate how their binding energy distribution is related to the mean energy of the cluster. Figure 7 shows the distribution of binding energies for one representative simulation at four different snapshots, in presence of original binary systems and stellar evolution.
At the beginning of the simulation, binding energies are very large if compared to the mean kinetic energy. In particular, the hardest part of the distribution is about five orders of magnitude higher than the typical energy scale of the star cluster. This means that the other stars in the cluster "see" the hardest binary systems as if they were single stars: the cross section of the hardest binary systems is so small that these can hardly interact with single stars.
In absence of original binary systems with a sufficiently large cross section, the star cluster creates new binary systems, with a larger semi-major axis and, thus, a large cross section for three-body encounters. This is the reason behind the large number of binary systems created at successive snapshots, that are close to the mean kinetic energy of the cluster. Finally, the loosely bound tail of the binary distribution consists of soft binaries that are continuously created and destroyed by dynamical interactions with their neighbours. Figure 8 shows the evolution of probability density function (PDF) of the binary semi-major axes and of the mass ratios and the multiplicity fraction, defined as the sum of the fraction of binaries and the fraction of triples. We consider two representative populations, one for simulations with original binaries (the same as in Fig. 7) and one for simulations without original binaries, in presence of stellar evolution.

Orbital parameters and multiplicity fraction
In presence of original binaries, the PDFs significantly change with time, because of the creation of a large number of dynamical binaries. In particular, the distribution of semi-major axes extends to higher values, and shows a secondary peak at 10 3 AU, the typical value at which dynamical binaries form. This value corresponds to ≈ 5 × 10 −3 pc, that is the lowest distance scale (it is the typical distance of stars split into Plummer spheres). As explained above, the cluster responds to the absence of interacting binaries by creating its own. This also explains why the distributions of the dynamically formed semi-major axes and mass ratios are very similar to those that form in absence of original binaries (as shown in the lower-left panel of Fig. 8).
As for the mass ratios' (q) distribution, dynamical interactions produce a steep increase of the PDF at high values, because the new binaries are typically formed by the low mass stars in the Plummer spheres. Also, the distribution of mass ratios extends towards lower values than the initial lower limit (q = 0.1). Most of the variations in the PDFs take place in the first 1 Myr, that is when the environment is dynamically active. Since then, the binary distributions remain almost unchanged. Also, the large number of dynamicallycreated small-mass binaries increases the total multiplicity fraction from ≈ 6% to ≈ 26%. In particular, these systems populate the lowest mass bin of Figure 1, by increasing the binary fraction from 0 to 20%.
In the absence of original binaries (NoBin+SE case), dynamical interactions produce a distribution of semi-major axes that is similar to the distribution of dynamically formed binary systems in the Bin+SE case, but cannot reproduce the hardest part of the Sana et al. (2012) binary distribution. Also, dynamical mechanisms tend to create equal-mass binaries. Remarkably, the binary fraction of dynamically formed binaries in the NoBin+SE case is mass-dependent: it grows with the mass of the primary star and mimics the trend of the observed distribution (Moe & Di Stefano 2017). Hence, in the absence of original binary stars, the cluster is able to produce a mass-dependent binary fraction. However, there is not sufficient energy at small scales to reproduce the hardest part of the initial distribution of Sana et al. (2012).

Exchanges
The degree of interactions between the binary systems and their host cluster can be quantified by evaluating the number of exchanges that take place. Fig. 9 shows the variation of the incremental number of exchanges. The original binaries take part in a limited number of exchanges, most of which are in the first 2 Myr of the cluster's life, when densities allow an efficient interaction with the other stars. In the following evolution, the original binaries interact much less, as indicated by the flatness of the curve. Nonetheless, because the original binaries are very hard, the few interactions they undergo exchange a sufficient amount of energy to affect the global evolution of the cluster, as shown by the evolution of r50 (Fig. 4). Interestingly, the total number of exchanges is about two orders of magnitude higher than that of original binaries and does not depend on the presence of an initial population of binary stars. This aspect indicates that the cluster under consideration is a very active environment for binary interactions and confirms that the most interacting binaries are dynamically created by the cluster itself. However, most of these exchanges involve binaries that are loosely bound (see also Fig. 7) and thus their energy exchange is quite low with respect to that of the original binaries.

Comparison with other initial conditions
The novelty of the Hydro initial conditions can be better understood if we compare their evolution to that of other, more idealized initial conditions. To this purpose, we ran simulations with the initial conditions presented in Section 3.1. Since we want to focus on the dynamical evolution with different initial phase-space distributions, we decided to run these simulations without stellar evolution. Figure 10 shows the evolution of the medians of the distributions of r50, of r10, and of the ratio r50/r10, that measures the concentration of the system. In the initial conditions, the Hydro clusters have a much larger ratio r50/r10 than the other models. Hence, they have very dense cores and rather extended halos, because of the scale of the sub-structures. For these intrinsic differences, the evolution of the characteristic radii of the Hydro simulations is considerably different from that of the other distributions.

Cluster Expansion
In the first Myr, the Hydro case is the only one that does not show a monotonic increase of r50 because of the initial sub-cluster motion (as discussed in Section 4.2.1). All of the other initial conditions present a monotonic increase of r50 and r10, but with different slopes. The Loose Fract case, that is initialized with the same half-mass radius as the Hydro case, shows a mild expansion on both scales, due to its supervirial state. The low density of the central regions (the initial value of r10 is larger than in the Hydro case by a factor of 5, see Tab. 1) does not allow efficient star-star interactions, that would power a faster expansion. The King and Dense Fract models, that are set to match the core radius of the initial Hydro simulations, undergo a stronger expansion from the very beginning of their evolution. These two different initial conditions present a very similar behavior.
The peculiarity in the evolution of the Hydro case is evident when the evolution of the ratio r50/r10 is taken into account. All the cases except the Hydro present a monotonic slow decrease for the r50/r10 ratio, that indicates that the systems expand at a similar rate at both scales. The Hydro initial conditions, instead, show an initial steep decrease, because the growth of r50 is balanced by the sub-clump motion (see Fig. 4), while at smaller scales the cluster expands rapidly. Even when it has reached a nearly monolithic shape, the evolution of its r50/r10 ratio is very different from the other initial conditions: this ratio rapidly increases until it reaches a maximum at about 5 Myr. Such a difference may be explained in terms of the stronger mass segregation that features the Hydro simulation (we will discuss this point in Sect. 5). Figure 11 shows the evolution of the total binding energy for different initial conditions. In absence of original binaries, every initial configuration creates its own population and the resulting total binding energy is strictly connected to the initial energy scale of the system. In particular, the Hydro, King and Dense Fract final binding energies are similar to each other as they are initialized with similar core radii, whereas the total binding energy of the Loose Fract systems is about one order of magnitude lower. Most of the total binding energy is contained in a limited number of binaries (from 2 to 5) that go on hardening as the simulation proceeds. This relation between the total binding energy and the global scales of the clusters confirms that star clusters are self-regulating systems with respect to their binary populations (Goodman & Hut 1989: in absence of binaries, each system creates its own population of binaries, with binding energies of the order of its global energy scales.

DISCUSSION
The Hydro star clusters present a very distinctive evolution of the r50/r10 ratio. We studied what factor determines the growth of this ratio during the monolithic phase. In particular, we focus on the impact of the initial degree of mass segregation. In fact, a high degree of mass segregation would allow the most massive stars to rapidly form a centrally concentrated core that is dynamically separated from the rest of the cluster, the scenario usually referred to as Spitzer instability (Spitzer 1969). If this happens, the distribution of massive stars is hotter than the rest of the cluster, because they remain more concentrated and the local value of the velocity dispersion decreases with the distance from the centre. Previous N −body simulations have found evidence that, for a wide range of initial conditions, the most massive stars in a system do not move slower than the low-mass stars ( as one would expect based on the tendency of stellar systems towards energy equipartition (Trenti & van der Marel 2013;Bianchini et al. 2016). A confirmation that the most massive stars can present higher velocities has also been found in proper motion observations of the open cluster NGC 6530 . Wright & Parker (2019) showed that this aspect can be explained by the combination of Spitzer instability and a cool collapse. If the most massive stars remain more concentrated than the rest of cluster, then the core, that is mostly populated by these massive stars, is expected to expand slower than the rest of the cluster.
To quantify the impact of mass segregation on the evolution of the cluster, we selected the 30 most massive stellar particles 2 and evaluated the ratio between their velocity dispersion σmass and the velocity dispersion of all the stellar particles σ all . For these calculation, only stars inside two half-mass radii are considered, as done in Wright & Parker (2019). The evolution of the ratio between these two velocity dispersions is shown in the upper panel Fig. 12. In all the phase-space configurations except the Hydro, the velocity dispersion ratio is about one and does not change very much with time. In the Hydro case, instead, the high initial value of the velocity dispersion ratio suggests that the stellar cluster presents a strong initial mass segregation. Also, during the monolithic phase, the velocity dispersion ratio grows because, after the merger of the two main sub-clumps, their most massive stars rapidly segregate towards the centre, while the system globally expands. The segregation of the massive stars towards the centre of the potential well may be enhanced by the fact that, in each sub-clump, the stars have already formed a massive core that segregates as one single, very massive particle. In the case with binaries, the velocity dispersion value grows enough to match the observed value for NGC 6530.
The connection between the growth of the velocity disper-sion ratio and the degree of mass segregation is confirmed by the trend shown by the ratio of the half-mass radii of the 30 most massive stellar particles rmass and the overall half-mass radius r50, shown in the bottom panel of Fig. 12. The Hydro simulations show an initial strong degree of mass segregation. The initial small scale expansion makes this ratio instantly grow; but, then, it rapidly decreases because of the strong segregation at the centre of the cluster. The initial degree of mass segregation seems to be the most important factor in the growth of the velocity dispersion ratio in the Hydro case: a stronger initial degree of mass segregation triggers the rapid formation of a dense core that expands more slowly than the rest of the cluster. Also, the rapid formation of a dense core could influences the interaction rate between binaries and the host cluster. If the primordial hard binaries live in a denser environment, they are more likely to interact: this explains why the Hydro initial conditions present different expansions for the cases with and without binaries (Fig. 10).

SUMMARY AND CONCLUSIONS
We studied the early dynamical evolution (t < 10 Myr) of young stellar clusters with realistic populations of binaries and different initial phase-space distributions. The initial conditions for our N −body simulations are obtained by combining a new algorithm to generate realistic stellar and binary distributions (Sana et al. 2012;Moe & Di Stefano 2017) with the joining/splitting algorithm defined in Ballone et al. (2021), to derive initial conditions from hydrodynamical simulations.
For the hydrodynamical initial conditions (Hydro cluster), we considered different evolutionary cases by switching on and off the presence of original binary stars and stellar evolution in order to weight their contribution to the dynamical evolution. Our results show that the evolution of the cluster is characterized by two distinct evolutionary phases: first, the global expansion of the cluster is balanced by the approaching of its main sub-clumps, while at small scales the cluster expands instantaneously. After 1 Myr, the cluster has reached a nearly monolithic shape and expands as a whole, following a post-core collapse expansion. Binaries tend to speed up the expansion of the cluster in this phase, making the half-mass radius expand faster, while stellar evolution has a minor impact on the early dynamical evolution of the cluster, but has a major impact on mass loss.
We compared the evolution of the Hydro star cluster to that of star clusters with spherical distributions of stars (King, Loose Fract, Dense Fract). The main difference between the Hydro cluster and the others relies in the evolution of the r50/r10 ratio, that measures the concentration of the system. The Hydro cluster, in fact, shows a distinctive trend of r50/r10. At the beginning of the simulations, r50/r10 is much larger in the Hydro cluster than in the other models, because the Hydro cluster is an aggregate of several sub-clumps, resulting in a large total half-mass radius, but its core radius is very small, since it basically coincides with the core radius of the densest sub-clump. The r50/r10 ratio decreases very fast (< 1 Myr) in the Hydro cluster, reaching values similar to the other clusters, because of the hierarchical merger of the sub-clumps, which reduces the total half-mass radius. However, at t > 1 Myr the value of r50/r10 keeps decreasing Lower panel: evolution of the ratio between the half-mass radius of the 30 most massive star particles (rmass) and the overall half-mass radius, r 50 . Lines and colours are the same as Fig. 10.
in the spherical models, while it grows again in the Hydro case. The late growth of r50/r10 in the Hydro cluster is due to its initial high degree of mass segregation, which allows it to form a centrally concentrated core of massive stars. As this core expands more slowly than the rest of the cluster, the ratio between the velocity dispersion of the most massive stars and that of all the stars increases. In the case with binaries, it grows enough to match the observed value for NGC 6530 .
The initial binary stars we set based on observational constraints (Sana et al. 2012;Moe & Di Stefano 2017) are generally too hard to interact in an efficient way with the host environment. The stellar systems recover from the lack of interacting binaries by dynamically creating additional binaries with binding energy of the order of their kinetic energy. Also, in the absence of primordial binaries, the binary fraction of the dynamically formed binaries show a binary fraction that increases with the mass of the primary star. This behaviour spontaneously reproduces the relation between binary fraction and stellar mass found in observations (Moe & Di Stefano 2017).