# The Physics of Star Cluster Formation and Evolution

Martin Krause, Stella Offner, Corinne Charbonnel, Mark Gieles, Ralf Klessen, Enrique Vázquez-Semadeni, Javier Ballesteros-Paredes, Philipp Girichidis, J. Diederik Kruijssen, Jacob Ward, et al.

► To cite this version:

Martin Krause, Stella Offner, Corinne Charbonnel, Mark Gieles, Ralf Klessen, et al. The Physics of Star Cluster Formation and Evolution. Space Science Reviews, 2020, 216 (4), pp.64. 10.1007/s11214-020-00689-4. hal-04894664

HAL Id: hal-04894664

https://hal.science/hal-04894664v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Physics of Star Cluster Formation and Evolution

Martin G. H. Krause · Stella S. R. Offner ·

Corinne Charbonnel · Mark Gieles ·

Ralf S. Klessen · Enrique Vázquez-Semadeni ·

Javier Ballesteros-Paredes Philipp Girichidis ·

J. M. Diederik Kruijssen · Jacob L. Ward ·

Hans Zinnecker

Received: 31 Jan 2020 / Accepted: date

Martin G. H. Krause

Centre for Astrophysics Research, School of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane, Hatfield, Hertfordshire AL10 9AB, UK

E-mail: M.G.H.Krause@herts.ac.uk

Stella S. R. Offner

Department of Astronomy, The University of Texas, Austin TX, 78712, U.S.A.

Corinne Charbonnel

Department of Astronomy, University of Geneva, Chemin de Pegase 51, 1290 Versoix, Switzerland; IRAP, CNRS & Univ. of Toulouse, 14, av.E.Belin, 31400 Toulouse, France

Mark Gieles

Institut de Ciències del Cosmos (ICCUB-IEEC), Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain; ICREA, Pg. Lluis Companys 23, 08010 Barcelona, Spain

Ralf S. Klessen

Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany

Enrique Vázquez-Semadeni

Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Campus Morelia, Apdo. Postal 3-72, Morelia 58089, México

Javier Ballesteros-Paredes

Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Campus Morelia, Apdo. Postal 3-72, Morelia 58089, México

Philipp Girichidis

Leibniz-Institut für Astrophysik (AIP), An der Sternwarte 16, 14482 Potsdam, Germany

J. M. Diederik Kruijssen

Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstraße 12-14, 69120 Heidelberg, Germany

Jacob L. Ward

Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstraße 12-14, 69120 Heidelberg, Germany

Hans Zinnecker

Nucleo de Astroquimica y Astrofísica, Universidad Autonoma de Chile, Avda Pedro de Valdivia 425, Providencia, Santiago de Chile, Chile
**Abstract** Star clusters form in dense, hierarchically collapsing gas clouds. Bulk kinetic energy is transformed to turbulence with stars forming from cores fed by filaments. In the most compact regions, stellar feedback is least effective in removing the gas and stars may form very efficiently. These are also the regions where, in high-mass clusters, ejecta from some kind of high-mass stars are effectively captured during the formation phase of some of the low mass stars and effectively channeled into the latter to form multiple populations. Star formation epochs in star clusters are generally set by gas flows that determine the abundance of gas in the cluster. We argue that there is likely only one star formation epoch after which clusters remain essentially clear of gas by cluster winds. Collisional dynamics is important in this phase leading to core collapse, expansion and eventual dispersion of every cluster. We review recent developments in the field with a focus on theoretical work.

**Keywords** galaxies: star clusters: general · ISM: kinematics and dynamics · open clusters and associations: general · stars: formation

# 1 Star clusters: more than a collection of stars

Star clusters have caught human attention since ancient times, as evidenced for example by depictions of the Pleiades on cave walls and the Nebra Disk (Rappenglück, 2001; Mozel, 2003). They continue to be a fascinating topic today, thanks to new and puzzling observations challenging the theoretical models.

Spitzer has traced the dense gas in a number of nearby young clusters (Fig 1) and shown its connection to young stellar objects (e.g., Gutermuth et al., 2011). Thanks to GAIA (Gaia Collaboration et al., 2018), we now know the kinematics of many clusters on a star-by-star basis (e.g., Ward & Kruijssen, 2018; Karnath et al., 2019; Kuhn et al., 2019). Chemistry is traced by spectroscopic and photometric surveys (Bastian & Lardo, 2018; Gratton et al., 2019); cluster winds have been detected in spaceborne X-ray observations (Kavanagh et al., 2011) and young super star clusters show evidence of MASER emission (Gorski et al., 2019).

These observations place strong constraints on theoretical modelling. The latter has been typically attempted from different angles with a view on explaining a particular subset of observations. A simulation that includes gas dynamics, stellar dynamics and chemistry to sufficient accuracy and from cloud collapse to cluster dispersal remains beyond reach for the foreseeable future. Approaches that focus on each aspect separately, or combine some aspects making some approximations therefore have to form the basis of our understanding of stellar clusters.

This review aims to provide an overview of the different theoretical approaches, puts them in context with each other, and aims to paint a comprehensive and coherent picture of the physics of star cluster formation and evolution. We are not aware of past projects with such an ambition, but previous reviews that have significant overlap with the present one include Mac Low & Klessen (2004); Zinnecker & Yorke (2007); Portegies Zwart et al. (2010); Gratton et al. (2012); Renzini (2013); Kruijssen (2014); Krumholz (2014); Longmore et al. (2014); Charbonnel (2016); Klessen & Glover (2016); Bastian & Lardo (2018); Gratton et al. (2019) and Krumholz et al. (2019). After defining star clusters in §2, we first review the onset of star formation in molecular clouds (§3) and the formation of stars in clusters
(§4). In this formation period, physical processes and multiple scales are coupled. From the end of the star formation epoch, gas dynamics (§5), stellar dynamics (§6) and nucleosynthesis (§7) evolve independently. Each section includes, however, links to the other fields. In particular, the chemistry in the predominantly old, multiple population clusters, discussed in §7, refers back to the formation epoch, where all the different processes are coupled. We conclude with a summary and outlook in §8.

## 2 What is a star cluster?

We adopt the ontological definition that a star cluster is a gravitationally bound group of stars inside a closed tidal surface if this volume is

1. not dark matter-dominated and

2. contains at least 12 stars.

The first condition distinguishes star clusters from galaxies. The second one from multiple star systems. This definition essentially follows Krumholz et al. (2019, though we do not distinguish here between different overdensities required in different environments). Groups of stars that are not gravitationally bound are called associations (Blaauw, 1964; Gieles & Portegies Zwart, 2011; Gouliermis, 2018, and Adamo et al. 2020, in prep.). For the Milky Way, bound star clusters have been subdivided into open clusters in the disc and globular clusters associated with the bulge and halo. Open clusters are generally young ($\lesssim 1$ Gyr) and have low mass ($\lesssim 10^5 M_\odot$) while globular clusters are generally old ($> 1$ Gyr) and massive ($\gtrsim 10^4 M_\odot$), quite typically survivors from the early Universe, representing the relics of star formation at high redshift. In fact, the oldest globular clusters in the Milky Way have a likely age $> 13$ Gyr and provide an important constraint for the age of the Universe (Krauss & Chaboyer, 2003; O'Malley et al., 2017). The distinction between open and globular clusters happens to correspond closely to a fundamental distinction in photometric properties and chemical abundance patterns: Open clusters are mostly single population clusters with a single main sequence in the colour-magnitude diagram, while almost all globular clusters have multiple main sequences and strong star to star variations in light-element abundances, i.e., multiple stellar populations. A more useful classification of star clusters is therefore between single and multiple population clusters (Carretta et al., 2010; Bastian & Lardo, 2018).

## 3 The onset of star formation in molecular clouds

Star clusters form from molecular clouds, which are the densest regions in the interstellar medium, and consist mostly of molecular hydrogen and several other molecules, which are used as tracers for observing these regions and their substructure. Molecular clouds range in mass from $\sim 10^3$ to $\sim 10^7 M_\odot$, and have extremely complex hierarchical (or fractal) morphologies (Elmegreen & Falgarone, 1996), with the densest regions embedded in larger, lower-density ones, and so on (e.g., Blitz & Williams, 1999). It has been suggested that the internal structure
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Three-color Spitzer images (3.6 (blue), 5.8 (green), and 24 µm (red)) of young, nearby ($d < 1000$ pc) clusters. Cluster source catalogs are contained in the Spitzer Extended Solar Neighbourhood Archive (SESNA, Gutermuth et al. in prep). Top left: Young (4-5 Myr), massive cluster Cep OB3b, which is part of the Cep OB3 molecular cloud complex (Gutermuth et al., 2011; Allen et al., 2012). Top right: IC 348 cluster (2-3 Myr), which is forming in a sub-region of the Perseus molecular cloud (Gutermuth et al., 2009). Bottom left: Extended field containing exposed cluster L988-e (courtesy of R. Gutermuth). Bottom right: Small, dense cluster IC 5146, where protostars are forming around a bright PAH emission bubble (Gutermuth et al., 2009).</figcaption></figure>

and dynamics of molecular clouds is instrumental in determining the early structure and kinematics of star clusters (e.g., Klessen et al., 2000; Klessen & Burkert, 2000, 2001; Offner et al., 2009; Kruijssen et al., 2012; Girichidis et al., 2012b; Vázquez-Semadeni et al., 2017, hereafter VS17).

## 3.1 The Gravoturbulent (GT) scenario

Molecular clouds are known to have internal supersonic non-thermal motions (Wilson et al., 1970), which follow a relation between the observed linewidth and the spatial scale (Larson's relation Larson, 1981; Hennebelle & Falgarone, 2012), although with substantial scatter (Ballesteros-Paredes et al., 2011; Miville-Deschênes et al., 2017). These supersonic motions were originally interpreted as large-scale radial motions, likely to originate from global collapse (Liszt et al., 1974; Goldreich & Kwan, 1974). However, this interpretation was soon rejected because, as it was argued, it would lead to excessively large star formation rates (SFRs) and should produce systematic velocity differences (i.e., red or blue line shifts) between
emission lines produced by HII regions located at the centers of the clouds and absorption lines produced at the outskirts of the clouds. Since such shifts were not observed, the supersonic motions were then interpreted as *small-scale* supersonic turbulence that produces a turbulent pressure capable of supporting the clouds against their self-gravity (Zuckerman & Palmer, 1974; Zuckerman & Evans, 1974). The requirement for the motions to be confined to small scales was necessary in order to avoid the generation of the unobserved line shifts and to produce an isotropic pressure that could support the clouds.

Since then, the prevailing paradigm for molecular clouds is that they are supported by some agent against their self-gravity, typically turbulence and/or magnetic fields. These are invoked in part to explain the observed star-formation efficiencies per free-fall time, $\epsilon_{\text{eff}}$ – the fraction of gas mass converted to stars over a free fall time – of ~1% for most giant molecular clouds (GMCs) (Zuckerman & Evans, 1974; Krumholz et al., 2019). Since turbulence is known to dissipate rapidly, typically in a crossing time, it was first proposed that the motions consisted of Alfvénic turbulence, because Alfvén waves were thought to be less dissipative than shocks (e.g., Shu et al., 1987). However, subsequent numerical simulations of MHD turbulence showed that it dissipates as rapidly as hydrodynamic turbulence (Mac Low et al., 1998; Stone et al., 1998; Padoan & Nordlund, 1999), implying that constant driving of the turbulence must be present to maintain it. In this *gravoturbulent* (GT) scenario (e.g., Klessen et al., 2000; Vázquez-Semadeni et al., 2003; Mac Low & Klessen, 2004), the clouds are supported globally by the pressure of the (continuously driven) supersonic, small-scale, isotropic turbulence, while locally, shocks are produced that in turn generate density fluctuations (sheets, filaments, and clumps), which may locally become Jeans unstable and collapse. Moreover, if magnetised turbulence provides support such that clouds are neither dispersing nor globally collapsing, it is assumed that the clouds are in approximate virial equilibrium between turbulence and self-gravity (Krumholz et al., 2006; Goldbaum et al., 2011). This assumption is consistent with observations (e.g., Larson, 1981; Heyer et al., 2009). In particular, the Larson (1981) linewidth-size relation observed in molecular clouds is interpreted as the manifestation of the energy spectrum, $E(k) \propto k^{-2}$, corresponding to strongly compressible, highly supersonic turbulence.

## 3.2 The Global Hierarchical Collapse (GHC) scenario

On the other hand, there is evidence that the process of formation of the molecular clouds is important for their subsequent dynamical evolution. The clouds seem to form by accreting tenuous ($n \sim 10 \text{ cm}^{-3}$) atomic gas, which often appears gravitationally bound to the molecular gas it surrounds (Fukui et al., 2009). Moreover, molecular clouds exhibit a hierarchical structure, so that their internal dynamics are governed by very similar processes. On smaller scales, star-forming cores accrete material from the scale of their parent clumps (i.e., the cores are said to be *clump-fed*; Liu et al., 2015; Yuan et al., 2018), and longitudinal, multi-parsec scale flows are routinely observed along filamentary clouds, which feed the main cores (or *hubs*) within the filaments (e.g., Myers, 2009b; Schneider et al., 2010; Kirk et al., 2013; Peretto et al., 2014; Wyrowski et al., 2016; Hacar et al., 2017; Chen et al., 2019b). Additionally, numerical simulations of the formation and evolution
of cold, dense atomic clouds from large-scale compressions in the warm, diffuse gas also suggest that the clouds engage into global, hierarchical collapse (GHC; Vázquez-Semadeni et al., 2019) soon after they reach their thermal Jeans mass (Vázquez-Semadeni et al., 2007, 2009; Heitsch et al., 2008). In what follows, we focus on this scenario, as it provides a direct link between the processes occurring in the gas during the collapse and the structural properties of the resulting stellar cluster(s).

### 3.2.1 Onset of large-scale gravitational contraction and turbulence generation

The clouds are expected to rapidly reach and exceed their thermal Jeans mass because the Jeans mass in the dense, cold gas is ~ $10^4$ times smaller than in the diffuse, warm gas (Gómez & Vázquez-Semadeni, 2014) and simulations indicate that the clouds actively accrete from their diffuse environment (Ballesteros-Paredes et al., 1999; Hartmann et al., 2001; Vázquez-Semadeni et al., 2006; Heitsch & Hartmann, 2008; Banerjee et al., 2009; Heiner et al., 2015; Wareing et al., 2019). This accretion implies that the clouds generally grow in mass, allowing them to become magnetically supercritical (i.e., unsupported by the magnetic field), gravitationally unstable, and molecular at roughly the same column density (~ $10^{21}$ cm$^{-2}$) for solar-neighbourhood pressures and metallicities (Hartmann et al., 2001; Heitsch et al., 2009; Vázquez-Semadeni et al., 2011; Heiner et al., 2015).

Simulations of the self-consistent formation and evolution of clouds by converging streams of diffuse gas (e.g., Heitsch et al., 2005, 2006; Audit & Hennebelle, 2005, 2010; Vázquez-Semadeni et al., 2006, 2007; Hennebelle et al., 2008; Banerjee et al., 2009) show that the very formation process of the cloud causes the generation of moderately supersonic (with respect to the sound speed in the cold, dense gas) turbulence by the combined action of various instabilities, such as the non-linear thin-shell instability (Vishniac, 1994), thermal instability (Field, 1965) and Kelvin-Helmholtz instability (see Heitsch et al., 2006, and Klessen & Hennebelle, 2010a, for further discussions). Similar effects have also been shown for shells of expanding bubbles (Krause et al., 2013a). The turbulence generates nonlinear density fluctuations in which the free-fall time $\tau_{ff} = \sqrt{3\pi/(32G\rho)}$ is significantly shorter than the average in the cloud.

The energy in the turbulent motions generated by the instabilities, which are only moderately supersonic with respect to the cold gas, and subsonic with respect to the warm gas, quickly becomes overwhelmed by the gravitational energy of the whole cloud (actually, a cloud complex), which then begins to undergo global gravitational contraction. In the GHC scenario, thus, the apparent near-virial state of molecular clouds and their substructures is not due to turbulent support, but rather to the infall motions driven by the self-gravity (Ballesteros-Paredes et al., 2011). It should be noted, however, that the infall is highly chaotic and so a truly random (turbulent) component is in fact maintained by the collapse (Klessen & Hennebelle, 2010b; Robertson & Goldreich, 2012; Murray & Chang, 2015; Li, 2018), although it is apparently not able to significantly delay the collapse, possibly due to the rapid dissipation. More experiments are needed to clarify exactly how much turbulence is generated by the collapse, especially in the context of the formation of the first stars, where energy loss via radiative cooling is suppressed due to the low metallicity of the gas (Sur et al., 2012; Latif et al., 2013; Schober
et al., 2012; Bovino et al., 2013; Federrath et al., 2014; Schober et al., 2015; Klessen, 2019).

The presence of turbulent density fluctuations with nonlinear amplitudes, together with the generally amorphous and flattened or filamentary shape of the clouds, has the important implication that realistic collapse is far from homogeneous (uniform spherical configurations, all material in the sphere reaching the center at the same time). It is well known that already in non-uniform spherical configurations (“cores”) with centrally-peaked radial density profiles, the central, densest parts terminate their collapse (i.e., reach protostellar densities) earlier than the outer parts, and then the rest of the material, which was initially at lower densities, continues to accrete onto the previously collapsed material (e.g., Larson, 1969; Penston, 1969; Shu, 1977; Hunter, 1977; Whitworth & Summers, 1985; Foster & Chevalier, 1993; Mohammadpour & Stahler, 2013; Keto et al., 2015; Naranjo-Romero et al., 2015). In a turbulent system, the nonlinear density fluctuations have free-fall times significantly shorter than that of the whole cloud, and so they can collapse faster, as soon as they become locally gravitationally unstable (compare Vázquez-Semadeni et al., 2019).

Under this regime, the cloud evolves towards containing a large number of thermal Jeans masses, in agreement with the observation that molecular clouds typically have masses $M_c$ upwards of $10^3 M_\odot$ (e.g., Mac Low & Klessen, 2004, and references therein). Thus, the cloud becomes a system of *collapses within collapses*, with an ever-larger hierarchy of collapsing scales, each one accreting from the next larger scale (Vázquez-Semadeni et al., 2019). This is a mass cascade, in some senses similar to the turbulent energy cascade (Field et al., 2008). This is also essentially Hoyle's fragmentation (Hoyle, 1953), except with nonlinear density fluctuations and non-spherical geometry of the clumps (Vázquez-Semadeni et al., 2019). Also, it can be considered as an extension of the *competitive accretion* scenario (Bonnell et al., 2001; Bonnell & Bate, 2006), with the accretion extending to cloud scales (~10 parsecs or more), and with the added ingredient that a whole hierarchy of chaotic, gravitational contraction flows is present.

### 3.2.2 Filament formation and filamentary accretion

At sufficiently advanced stages of a cloud's evolution, when its mass $M_c$ is much larger than the Jeans mass, it must behave essentially as a pressureless collapse, because precisely the meaning of $M_c \gg M_J$ is that the gravitational energy over-whelms the internal energy of the cloud. But it is known that pressureless collapse amplifies anisotropies, so that a triaxial ellipsoid contracts first along its shortest dimension to form a sheet, and then an elliptical sheet contracts again along its shortest dimension to form a filament (Lin et al., 1965). Therefore, it is expected that multi-Jeans mass molecular clouds should evolve to develop filaments, which are actually akin to “rivers” funnelling the mass from large to small scales (Gómez & Vázquez-Semadeni, 2014). This is consistent with the observation that dense molecular cloud cores appear as “hubs” at the intersection of filaments (e.g., Myers, 2009a), with the filaments feeding material to the hubs (e.g., Schneider et al., 2010; Sugitani et al., 2011; Kirk et al., 2013; Peretto et al., 2014; Chen et al., 2019b).

Since the majority (~90%) of pre- and protostellar cores in molecular clouds are located either in filaments or in hubs (Könyves et al., 2015), it follows that
star formation is initiated already in the flows feeding the hubs. This mechanism was referred to as “conveyor belt cluster formation” by Longmore et al. (2014) and in Krumholz et al. (2019), in opposition to “monolithic cluster formation”, in which the gas first collapses and subsequently forms stars in a centrally-concentrated cluster. The conveyor-belt mechanism is also observed in simulations of self-consistent cloud formation and evolution, in which the filaments form spontaneously by anisotropic gravitational contraction (Gómez & Vázquez-Semadeni, 2014; Vázquez-Semadeni et al., 2019).

### 3.2.3 Acceleration of star formation and delayed formation of massive stars

Another expected consequence of the global collapse and continued accretion onto the star-forming hubs observed in the simulations is an acceleration of the star formation before massive stars form. This increase in the star formation rate (SFR) is routinely observed in simulations of cloud evolution (e.g., Vázquez-Semadeni et al., 2010, 2017; Hartmann et al., 2012; Colín et al., 2013; Lee et al., 2015; Li et al., 2018), and predicted by models of clouds dominated by gravity (e.g., Zamora-Avilés et al., 2012; Zamora-Avilés & Vázquez-Semadeni, 2014; Caldwell & Chang, 2018). Observational evidence of the acceleration is provided by, for example, a) the age histograms of young embedded clusters, which systematically show a maximum at either the smallest ages, or at a certain, relatively recent age, together with a tail of older stars, of ages up to several Myr (e.g., Ballesteros-Paredes et al., 1999; Palla & Stahler, 2000, 2002; Huff & Stahler, 2006; Da Rio et al., 2010); b) a superlinear ($\sim t^2$) temporal dependence of the total number of stars formed at time $t$ in several young clusters (Caldwell & Chang, 2018).

In the GHC scenario, the increase of the SFR during the early stages of star-forming regions is due to the growth in mass, density, and size of the regions due to accretion from their parent structures. The increase in density implies an increase of the SFR because a larger fraction of the mass is at densities high enough that their free-fall time is much shorter than that of the mean density of the parent structure (Zamora-Avilés & Vázquez-Semadeni, 2014). But, additionally, the larger mass of the more evolved regions provides a larger mass reservoir, allowing for the formation of more massive stars. So, the star-forming regions evolve towards forming more massive stars, meaning that the formation of massive stars is delayed with respect to that of the first low mass stars, by several Myr in moderate-mass regions, according to the simulations (Vázquez-Semadeni et al., 2009, 2017, top left panel of Fig. 2). Note that low-mass stars always form, but the maximum mass of the stars that can form is capped by the instantaneous mass of the hub where they form, and increases with time as long as the hub’s mass increases. Eventually, however, the stellar feedback begins to erode the hub in the simulations, decreasing its density and mass, and also eroding the filamentary accretion flow, decreasing the maximum stellar mass that can form. A model for the development of the high-mass slope of the IMF based on the same principle, of the mass of the most massive star being bounded by the mass of the hub in which it forms, has been developed by Oey (2011). This delayed formation of massive stars also implies that the age range of the massive stars is smaller (and they are younger) than that of the low-mass stars, which begin to form since the onset of the star formation activity in the region. Equivalently, the mass range of the younger stars
<figure><img src="image_3.png" /></figure>

**Fig. 2** *Top left:* Normalised cumulative stellar mass histograms of star forming regions in a numerical simulation of GHC (Vázquez-Semadeni et al., 2017) at various times. As time proceeds, a larger fraction of the stars are seen to be massive, until feedback begins to disrupt the cloud. *Top right:* Mass versus age of the members of a stellar group in the same simulation, at time $t = 22.4$ Myr, corresponding to roughly 3.5 Myr after the onset of star formation in the region. *Bottom left:* Mass versus distance from center of mass of the group members at $t = 23.7$ Myr in the simulation (≈ 4.7 Myr after the onset of star formation. The more massive stars are seen to be located near the center of mass. *Bottom right:* Groups constituting the cluster in the simulation at time $t = 30.0$ Myr (11 Myr after the onset of SF), as identified by a friends-of-friends algorithm. Stars belonging to the same group have the same colours. Figures from Vázquez-Semadeni et al. (2017).

is larger and extends to higher masses than that of the older stars (top right panel of Fig. 2).

### 3.2.4 Mass and age radial gradients

Star formation occurring in the filaments generally involves lower-mass cores, because they are themselves part of the flow falling onto the main hubs, which are the main accreting centers. That is, the stars formed in the filaments do so in secondary gravitational potential wells, while the hubs are the primary wells. Therefore, the more extended secondary star formation in the filaments generally produces lower-mass stars, and thus tends to produce a primordial mass segregation in the cluster (bottom left panel of Fig. 2), independent of any N-body processes that may occur afterwards (§6.5).

Because the secondary star formation in the filaments, involves lower-mass stars that have been forming for a longer time, the median age of the more distant
stars tends to be larger than that of those nearer the hub. In the simulation from Fig. 2 a median-age gradient of ~1 Myr pc⁻¹ is found (Getman et al., 2018), consistent with the gradient observed in the clusters in the MYStIX (Feigelson et al., 2013) and SFiNCs (Getman et al., 2017) star-forming region catalogs.

The hierarchical and filamentary structure of the collapse flow is imprinted on the structure of the cluster itself, which therefore adopts a self-similar, fractal-like spatial distribution, and retains traces of a filamentary morphology, as seen in the bottom right panel of Fig. 2, which shows the groups constituting the cluster in the simulation at time *t* = 30.0 Myr, as identified by a friends-of-friends algorithm. The groups are seen to be strung along a long filamentary structure, and also, when the linking parameter of the algorithm is varied, the number of identified groups varies, indicating a hierarchical structure (see Vázquez-Semadeni et al., 2017, for further details).

The above structural properties of the nascent cluster are blurred to some degree when the massive stars lose their mass, the gas is cleared by the stellar feedback and the stars have time for *n*-body interactions (see below). Therefore, these “primordial” structural features are expected to be more prominent in younger clusters.

## 4 The formation of stars in clusters

The formation and evolution of star clusters is a multi-scale process that depends in detail on the formation of the constituent individual stars. In turn, forming stars influence the accretion and dynamics of their neighbours through stellar feedback, including winds, radiation and supernovae.

Stars accreting from a common gas clump or protostellar core may compete with one another for fuel, while accretion disk properties depend on the local ionising flux, which is set by the distribution of nearby massive stars. Ionising radiation, winds and jets from protostars interact with their own accretion streams, as well as those of other stars. Finally, there is the puzzling observation of multiple populations in clusters (§7), where massive star ejecta may affect the accretion flows of lower mass stars. Understanding how an individual star in a cluster grows by accretion of gas and why it reaches a particular mass is therefore inseparable from the larger cluster context.

During the earliest stages of accretion, stars are hidden from view, and direct probes of accretion, like stellar spectral lines, which are exploited to study accretion in T-Tauri stars, are unavailable. Instead, indirect evidence of accretion, such as protostellar outflows and luminosities must be used to reconstruct the magnitude and history of accretion. In this section, we first discuss the observational signatures of accretion and their implications. We then summarise a variety of theoretical models for protostellar accretion. We sub-divide these into three categories: models based on the properties of cores and/or filaments that host the protostar, models that focus on the protostar-disk relationship and models that depend on feedback and the larger protostellar environment. This division is mainly for convenience, since in practice, accretion is determined by a variety of nonlinear processes that span a broad range of times and physical scales. Finally, we discuss models for how and why accretion ultimately ceases, which is critical for understanding the accretion histories of individual stars as well as global
properties such as the star formation efficiency, star formation rate and lifetime of molecular clouds.

## 4.1 Observational Signatures of Accretion: Protostellar Luminosities, Outflows and Spectral Lines

Protostellar outflows are a direct byproduct of the accretion process (see Chapter Processes for more details). If a fixed fraction of accreting material is flung outwards in an outflow, then in principle by measuring the outflow mass flux it is possible to reverse engineer the accretion rate and history. Observations of protostellar outflows suggest protostellar accretion rates of $10^{-4} - 10^{-9} M_{\odot} yr^{-1}$, where younger or more massive protostars have higher inferred accretion rates (Bally, 2016).

Outflow morphology also gives important insights in the accretion process. Outflows and jets (highly-collimated flows, usually observed in optical emission), frequently exhibit regularly spaced clumps along the outflow axis, “bullets” (Bally, 2016; Zhang et al., 2016). The spacing of the bullets indicates that accretion is variable on timescales of hundreds to thousands of years (Bachiller et al., 1991; Lee et al., 2009; Arce et al., 2013). In the most extreme events, the accretion rate, and hence the source luminosity, rises by several orders of magnitude over a period of years in a brief “episodic” accretion burst (Audard et al., 2014).

Outflows also exhibit precession or changes in direction, providing a window into the angular momentum of accreting material and the impact of binarity on the accretion process (Shepherd et al., 2000; Hirano et al., 2010; Lee et al., 2017a). A number of outflows appear to have two components: a highly-collimated component, likely launched close to the protostar, and a wider-angle, slower component that likely arises from the accretion disk (Hirano et al., 2010; Arce et al., 2013).

A variety of uncertainties underpin the connection between outflows and accretion (Dunham et al., 2014b). Much of the outflowing material is entrained core material (Offner & Chaban, 2017), so the outflow is not a direct measure of accreting gas. The typical outflow dynamical time, as measured by the outflow extent and gas velocity, $t_{dyn} \sim L_{out}/(2v_{out}) \sim 10^3$ yr, is shorter than the expected protostellar lifetime, and thus provides only a narrow window into the total accretion history (Bally, 2016). This is probably related to accretion physics (§4.3) rather than protostellar dynamics in clusters: Since the velocity dispersions of dense gas and young stars in star-forming environments are typically of the order $0.1-1$ km s$^{-1}$ (e.g., Kirk et al., 2010; Foster et al., 2015), a protostar would require of the order of $10^5$ years to move out of a filament or core of 0.1 pc thickness.

Protostellar luminosities provide another constraint on accretion (see Chapter processes). At early times and for low masses, i.e., before the intrinsic luminosity of the protostar becomes significant, the luminosity is directly proportional to the accretion rate. By assuming reasonable properties for the protostellar mass and radii, it is possible to set limits on the accretion rate. However, protostellar evolution remains uncertain in part because it is itself sensitive to the accretion history (Palla & Stahler, 1991; Baraffe et al., 2009; Hosokawa et al., 2011). Observations of clusters of protostars show orders of magnitude scatter, such that on average the luminosity is weakly dependent, at best, on the protostellar class (Dunham et al., 2014a; Fischer et al., 2017). The difficulty of mapping classes to evolutionary
stage further confuses accretion trends over time (Robitaille et al., 2006; Dunham et al., 2010; Offner et al., 2012).

Time-domain studies of protostellar luminosities are more informative and support the highly variable nature of accretion suggested by outflow observations. Changes in luminosity are observed on timescales of days to decades spanning changes from as little as a few percent to several orders of magnitude in brightness (Rebull et al., 2014; Audard et al., 2014). Low magnitude, shorter timescale variations, which are quite common, are likely caused by stellar activity or disk oc- cultations, while more extreme and rarer luminosity changes can only be explained by accretion fluctuations (Hillenbrand & Findeisen, 2015).

Early observations of protostars noted that they were on average about 10 times dimmer than simple accretion models and timescale arguments would suggest (Kenyon et al., 1990; Kenyon & Hartmann, 1995; Evans et al., 2009). This became known as the “protostellar luminosity problem”. A variety of theoretical solutions have been proposed that resolve this problem, including episodic and slow accretion (§4.3).

Spectroscopic measurements, which become possible once young stellar objects are more than ≈ 2 × 10⁵ yr old, indicate that accretion declines steeply at late times (Hartmann et al., 2016). However, the accretion rate depends on both age and mass, which are difficult to disentangle due to measurement and model uncertainties. Observations of Balmer continuum, photometry and emission lines suggest $\dot{M} \propto M_*^\alpha$, where $\alpha = 1.5 - 3.1$ and $\dot{M} \propto t^\beta$, where $\beta = -1.6 - -1.2$ (Hartmann et al., 2016).

## 4.2 Core-Regulated Accretion Models

Core-regulated models assume that collapsing gas (compare §3) efficiently proceeds from ~ 0.1 pc to au scales such that the infall rate is equal to the protostellar accretion rate. Under this assumption, the details of accretion depend only upon the properties of the local gas reservoir.

In the simplest model, collapse is regulated by the interplay of thermal pressure and self-gravity. The accretion rate due to the collapse of a uniform isothermal sphere of gas is known eponymously as the Larson-Penston solution (Larson, 1969; Penston, 1969):

$$ \dot{m} = 46.9 \frac{c_s^3}{G} = 7.4 \times 10^{-5} \left( \frac{T}{10 \text{ K}} \right)^{3/2} M_{\odot} \text{yr}^{-1}, \quad (1) $$

where $c_s$ is the sound speed. If the gas is isothermal and centrally condensed the infall solution is self-similar and can be written (Shu, 1977):

$$ \dot{m} = 0.975 \frac{c_s^3}{G} = 1.5 \times 10^{-6} \left( \frac{T}{10 \text{ K}} \right)^{3/2} M_{\odot} \text{yr}^{-1}. \quad (2) $$

Figure 3a shows the isothermal sphere accretion rate and a variety of other analytical predictions as defined below. In these models, the accretion is by nature time-invariant and independent of stellar mass.

These idealised solutions, however, gloss over a great deal of important physics. Cores are observed to be magnetised and turbulent (Crutcher, 2012; Kirk et al.,
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Stellar accretion rates as a function of time. Left (a): different analytic model predictions for protostellar accretion. Right (b): accretion of individual protostars in an MHD simulation of protostars forming in a turbulent giant molecular cloud (figure adapted from Li et al. (2018) and reproduced with permission). Turbulence, dynamics and protostellar outflows together significantly modulate the accretion histories.</figcaption></figure>

2017). They also exhibit velocity gradients indicative of rotation (Chen et al., 2019a,c). Moreover, cores are not spheres but are frequently asymmetric, where elongation is likely dictated by the greater filamentary environment that hosts them (Pineda et al., 2010; Arzoumanian et al., 2011). A variety of theoretical models have attempted to address these complications (Terebey et al., 1984; Fatuzzo et al., 2004; Adams & Shu, 2007). The simplest way of modifying equation (2) is to treat turbulence and magnetic fields as contributions to the total pressure support against gravity. For example, replace the thermal sound speed with $c_{\text{eff}} = c_s(1 + 2\alpha + \beta)^{1/2}$, where $\alpha = P_B/P_{\text{th}}$ and $\beta = P_{\text{turb}}/P_{\text{th}}$ are the ratio of magnetic and turbulent pressure to thermal pressure, respectively (Stahler et al., 1980). However, both turbulence and magnetic fields are intrinsically anisotropic, which suggests this approach over-simplifies their true impact on the accretion rate and over-estimates their contributions to pressure support. Also, these models implicitly assume that star formation can be represented by discrete collapsing regions and thus, arguably, are applicable only for isolated, low-mass star formation (cf. §3).

The “turbulent core model” developed by McKee & Tan (2003) treats cores as high-column density centrally condensed objects and includes turbulence as an effective pressure. This model predicts that $\dot{M} \propto M^{1/2} M_f^{1/4}$, where M is the instantaneous stellar mass and $M_f$ is the final mass of the star at the end of accretion. This naturally implies that high-mass stars have higher accretion rates than low-mass stars, form faster and that their accretion *increases* in time. Hydrodynamic calculations of high-mass star formation, which adopt high-column density, high-mass cores as initial conditions exhibit these trends (Krumholz et al., 2012; Rosen et al., 2016).

Hydrodynamic simulations of forming star clusters paint a very dynamical picture, particularly in clusters with high-stellar densities. In the “competitive accretion” model protostars begin as small seeds formed by local collapse, which compete with one another for the available gas (Bonnell et al., 2001). Birth lo-
cation and dynamical interactions determine the protostellar locations within the gravitational potential well and thus their rate of gas accretion. More massive stars naturally form in the center of the cluster and are best positioned to rapidly accrete gas (Fig. 2, Bonnell et al., 2001; Bonnell & Bate, 2006). In this scenario accretion continues until the gas runs out, which occurs on ~ a global free-fall time. As a result, all stars have the same formation time, which is set by the cluster environment. Analytically, this corresponds to accretion rates of $\dot{M} \propto M^{2/3}$ for stars in gas-dominated potentials (Bonnell et al., 2001). At late times when the stellar mass exceeds the gas mass, accretion limits to $\dot{M} \propto M^2$ (compare Ballesteros-Paredes et al., 2015; Kuznetsova et al., 2017, 2018). Numerical simulations following the formation of massive stars from a 250 pc scale interstellar medium region suggest that massive stars form over a longer time period via converging, filamentary gas flows (Padoan et al., 2019). Their cores are less massive than the final stellar mass at any given time, i.e., massive stars do not form from progenitor massive turbulent cores. Inflow rather than competition drives the accretion behavior.

Stellar accretion is highly variable but does not increase with stellar mass as predicted by both the turbulent core and competitive accretion models. All core-regulated models fall somewhere in the continuum between constant accretion rate and constant accretion time, between the highly dynamical and isolated star formation paradigms.

## 4.3 Disk-Regulated Accretion Models

Observations suggest that mass does not pass smoothly from the outer envelope to the protostar but instead accretes in a more variable process as mediated by an accretion disk. Thus, the semi-analytic models outlined above only describe the *time-average* accretion behavior. Variability in models that do not explicitly include disk physics arises purely from variation in the environment and evolution of the host gas reservoir (Padoan et al., 2014; Li et al., 2018; Padoan et al., 2019).

The formation of a disk is a direct consequence of angular momentum in the star formation process. In the absence of angular momentum or in the limit of perfectly efficient angular momentum transport accretion disks would not exist. However, observations tell us disks are common (Tychoniec et al., 2018; Andrews et al., 2018). They act as a repository for high-angular momentum gas and effectively sort low-angular momentum material, which moves inwards towards the protostar, and high-angular momentum material, which moves outwards. The two dominant processes for angular momentum transport in disks, viscous torques due to turbulence (Balbus & Hawley, 1994) and gravitational instability (GI) (Toomre, 1964; Laughlin & Bodenheimer, 1994), both produce variability in the accretion flow.

Viscous torques require the activation of the magnetorotational instability (MRI), which depends on the local ionization fraction (Balbus & Hawley, 1994). If the gas is not sufficiently ionized then the magnetic field is poorly coupled, reducing the efficacy or shutting off the MRI entirely (Blaes & Balbus, 1994). Thus, MRI-regulated accretion disks may undergo periods with little or no accretion during which material builds up in the disk, followed by periods when the gas is thermally ionized initiating a burst of accretion (Zhu et al., 2009b). During
these bursts accretion may be elevated by several orders of magnitude, similar to
observed FU Ori bursts (Audard et al., 2014).

The requisite ionization for the MRI may be provided by the parent star and
its environment including FUV radiation, x-rays, and cosmic rays (Umebayashi &
Nakano, 1981; Semenov et al., 2004; Glassgold et al., 2007; Perez-Becker & Chiang,
2011). Consequently, disk surface layers are generally strongly ionized such that
accretion continues in a layered fashion, where gas accretes in the surface layers
while the disk mid-plane remains predominantly neutral and is an MRI “dead
zone” (Gammie, 1996). High periods of accretion may in turn increase the x-ray
and cosmic-ray ionization towards the disk mid-plane, prompting accretion deeper
in the disk and boosting the magnitude of the accretion burst (Offner et al., 2019).

Gravitational torques, which are the dominant transport mechanism in the outer disk, may also prompt large accretion variations (Kratter & Lodato, 2016). If mass builds up in the inner disk, the disk may undergo GI and form small clumps. If these clumps migrate inwards they produce burst events as they accrete onto the star (Vorobyov & Basu, 2005, 2006). Mild GI, in the form of spiral arms, to severe GI, which causes catastrophic disk fragmentation, produce accretion variability from factors of a few to orders of magnitude (Audard et al., 2014).

Dynamical interactions between stars or close binary companions can also pro-
duce accretion variability (Adams & Lin, 1993). Close passage gravitationally per-
turbs the disk, prompting instability and elevated accretion (Bonnell, 1994). Fi-
nally, variation of the angular momentum of the infalling gas on larger scales may
also create luminosity variations, either through direct accretion (Padoan et al.,
2014) or by affecting disk properties (Lee et al., 2017b).

The frequency and magnitude of disk-mediated accretion bursts depend both
on disk microphysics and the larger disk environment (Kratter et al., 2010). Curr-
rent observations show a heterogeneous distribution of large, small, smooth and
structured disks. Likely both, GI and MRI, play a role in disk evolution (Armitage
et al., 2001; Zhu et al., 2009a). The corresponding scatter in protostellar luminosi-
ties provides one solution for the protostellar luminosity problem (Kenyon et al.,
1990; Offner & McKee, 2011; Dunham & Vorobyov, 2012; Padoan et al., 2014).

4.4 Feedback-Regulated Accretion Models

Stellar feedback, in the form of protostellar outflows, winds and radiation, also
shapes the accretion process, either by reducing the mass reservoir available for
accretion (as found in observational and numerical work, see e.g. Dale et al., 2015;
Ginsburg et al., 2016) or by dispersing bound gas and halting accretion alto-
gether. The earliest semi-analytic model for feedback-regulated accretion weighed
the competition between accretion and outflow feedback (Norman & Silk, 1980).
Feedback-regulated models are often formulated more generally in terms of a distri-
bution of stopping times or probabilities, which has the advantage that the model
can be agnostic about the particular mechanism halting accretion. For example,
several more recent models assume accretion durations follow the probability dis-
tribution, f(t) = 1/τe⁻t/τ, where τ is the mean accretion time (of the order of
10⁵ yr). Such models can reproduce the stellar IMF and match the observed pro-
tostellar luminosity distributions without appealing to overly long accretion times
or significant periods of episodic accretion (Basu & Jones, 2004; Myers, 2009b, 2012).

A variety of hydrodynamic simulations of accreting protostars including protostellar outflows have been carried out, which demonstrate that outflows can indeed efficiently expel 30-60% of the dense core material and reduce overall star formation efficiencies by ~30% (Hansen et al., 2012; Machida & Hosokawa, 2013; Offner & Arce, 2014; Federrath, 2015; Offner & Chaban, 2017; Tanaka et al., 2017). Simulations of isolated dense cores including protostellar outflows find that the main phase of accretion continues for 0.3-0.5 Myr, depending on the degree of turbulence and magnetic field strength (Machida & Hosokawa, 2013; Offner & Arce, 2014; Offner & Chaban, 2017). The accretion rate of a protostar accreting within a turbulent, magnetised dense core can be described in terms of the current protostellar mass, $m$, and its final mass, $m_f$: $\dot{m} = m_0 \left(\frac{m}{m_f}\right)^{1/2} m_f^{3/4} \left[1 - \left(\frac{m}{m_f}\right)^{1/2}\right]^2$,
where $m_0 \propto \Sigma_c^{3/4}$ is a constant coefficient related to the surface density of the core, $\Sigma_c$, and both $m$ and $m_f$ are in solar masses (Offner & Chaban, 2017). This is effectively the predicted turbulent core model accretion rate (McKee & Tan, 2003), tapered by a multiplicative factor. While the final masses are influenced by the core magnetic field and turbulence, the accretion history can be analytically described independently of the gas physical properties. Simulations of the impact of outflows on accretion within forming star clusters find wide variation in the accretion histories as shown for example in Figure 3b with some accretion rates steadily declining over time to $\dot{m} = 10^{-8} M_\odot \text{yr}^{-1}$ and others declining and then rising again due to protostellar dynamics (Li et al., 2018).

Feedback, turbulence and gravitational interactions may all play important roles in setting the accretion histories of individual stars. These same processes also drive the global evolution of the molecular cloud, gas dispersal (§5) and star cluster dynamics (§6). Thus, it is not possible to separate the formation of individual stars from the structure and evolution of the larger star cluster, i.e., whether it is a strongly bound cluster or a quickly dispersing association.

However, once the initial formation of the stars is completed, stars and gas effectively de-couple. In the following two sections we therefore first review studies that focus on the evolution of the gas, and then ones that treat the dynamics of the stars. The limitations of these approaches will become obvious when the transition from the star formation epoch will be considered and contact to models that specifically target the transition phase will be made.

# 5 Feedback and gas dynamics

After the initial star formation process, clusters become exposed, i.e., no dense gas is found in clusters from this stage onward. The process of a cluster becoming exposed may be driven by collective stellar feedback and may influence the dynamics of the stars. Later, star cluster winds can convey feedback energy to larger scales. Cooling flows have been discussed in the context of secondary star formation episodes, although age spreads in clusters are small, such that secondary star formation is likely restricted to associations.

Star clusters have a closed tidal surface and usually contain a focal point, the minimum of the gravitational potential. It is therefore generally expected that a
global pattern for the gas dynamics will form, which may in principle be inflow, outflow, or hydrostatic equilibrium. Contrary to galaxies, (even approximate) hydrostatic equilibrium is probably not relevant for star clusters.

## 5.1 Impossibility of hydrostatic equilibrium

To see this, we present a simple argument and show that starting from a situation close to hydrostatic equilibrium, stellar feedback would alter the gas properties quickly. Either cooling would take over leading to inflow, or heating, leading to outflow. Let us start with the hydrostatic equilibrium condition <fnref n="1" />:

$$ \frac{d\Phi}{dr} = -\frac{1}{\rho}\frac{dp}{dr} \qquad (3) $$

Approximating gradients by the absolute change out to the half-mass radius, we can write eq. (3) as<fnref n="2" />: $G(M/2)/r_h = p/\rho = k_B T/(\mu m_p)$. Radiative cooling will reduce the gas pressure. To maintain hydrostatic equilibrium, the cooling time<fnref n="3" /> $t_c = k_B T/(n\Lambda)$ therefore must at least exceed the crossing time<fnref n="4" /> $t_x = 2r_h/\sigma$ (Krause et al., 2019). Using also the definition

$$ \sigma^2 = GM/(\eta r_h) \qquad (4) $$

with $\eta = 7.5$ for a Plummer (1911) model, we arrive at the constraint

$$ \begin{aligned} n < & \frac{G^{3/2} \mu m_p M^{3/2}}{\eta^{1/2} \Lambda r_h^{5/2}} \\ = & 65 \, \text{cm}^{-3} \left( \frac{\Lambda}{10^{-27} \text{erg cm}^3 \text{s}^{-1}} \right)^{-1} \left( \frac{M}{10^5 M_{\odot}} \right)^{3/2} \left( \frac{r_h}{3 \text{pc}} \right)^{-5/2} . \end{aligned} \qquad (5) $$

For the relevant densities, $\Lambda$ is of the order of $10^{-27} \text{ erg cm}^3 \text{s}^{-1}$ (Bialy & Sternberg, 2019), which we have used for the scaling in eq. (5).

The immediate effect of stellar feedback is to add mass and energy to the intracluster gas. Hydrostatic equilibrium may only be maintained, if the energy input matches the energy loss via gas cooling. The particle density in the cluster increases at a rate (e.g., Mathews & Brighenti, 2003):

$$ \dot{n} = \frac{\alpha M/2}{4\pi r_h^3/3} = 56 \, \text{cm}^{-3} \text{Myr}^{-1} \left( \frac{\alpha}{10^{-16} \, \text{s}^{-1}} \right) \left( \frac{M}{10^5 M_{\odot}} \right) \left( \frac{r_h}{3 \, \text{pc}} \right)^{-3} . \qquad (6) $$

Here, we have scaled the mass loss factor $\alpha = \dot{M}/M$ to $10^{-16} \, \text{s}^{-1}$, a value that would be expected for a very young (≈ Myr) stellar population (Leitherer et al., 1999; Gaibler et al., 2005; Krause et al., 2013b).

Therefore, within a short timescale compared to the timescale of a cluster's evolution, stellar feedback would increase the gas density beyond the cooling limit. Hydrostatic equilibrium could then only be maintained, if the energy input was

<footnote n="1">¹ $\Phi$: gravitational potential, r: radius, p: gas density, p: pressure</footnote>

<footnote n="2">² M: cluster mass, μ: mean molecular weight</footnote>

<footnote n="3">³ k_B: Boltzmann constant, T: temperature, n: particle density, Λ: cooling function.</footnote>

<footnote n="4">⁴ r_h: half-mass radius, σ: line-of sight velocity dispersion.</footnote>
spatially fine-tuned and arranged to increase in time as required for the increasing
cooling rates. Since cooling rates are determined by atomic physics and energy
input by stellar physics, this will not be the case.

The late time evolution of α can be approximated as $\alpha = 4.7 \times 10^{-20} s^{-1}$ ($t/13 \text{ Gyr}$)$^{-1.3}$ (Mathews & Brighenti, 2003). Therefore, if at late times the cluster was for some reason in a state of hydrostatic equilibrium, it would take longer for the stellar feedback to increase the gas density beyond the stability limit. However, the relevant timescales also grow, such that the cluster would always become unstable on a timescale that is shorter than its age. The analysis depends only weakly on cluster radius. For smaller masses, the argument becomes stronger. Therefore, we can conclude that stellar feedback generally inhibits hydrostatic equilibrium in star clusters at all times.

## 5.2 Initial gas clearance

How much gas is left over from the star formation process and how violently this is
removed from the cluster has wide-ranging implications for star formation. From
abundances and ages of clusters and associations, Lada & Lada (2003) concluded
that stars generally form in dense clusters and get dispersed due to violent gas
expulsion and the associated change in gravitational potential (infant mortality).
Subsequent work has superseded this initial picture, showing that the statistics
depend crucially on the surface density threshold for the definition of star clusters
(Bressert et al., 2010), as well as the initial gas density at which stars are forming
(Kruijssen, 2012). Many recent studies show that star formation proceeds at a
variety of densities and spatial scales (e.g., Bastian et al., 2007; Sun et al., 2018;
Rodríguez et al., 2019) and detailed analysis of OB associations shows that they
did not evolve from significantly smaller structures (Wright et al., 2014; Ward &
Kruijssen, 2018; Ward et al., 2019). A good example is the Gaia study of the closest
OB association, Sco-Cen OB2 (Wright & Mamajek, 2018), for which alternative
formation scenarios based on multi-wavelength observations have been suggested
(Krause et al., 2018). While the formation of bound clusters and their dispersal
may be less common than once thought (compare also Kruijssen, 2012; Krumholz
et al., 2019), it is still interesting to ask what fate the gas experiences and what
roles it can play in any given cluster.

### 5.2.1 Gas expulsion

We use the term *gas expulsion* to refer to a special kind of gas removal, where a
significant mass of gas (≥ 50% of the total mass) is removed quickly (compared to
the crossing time for stars, i.e. impulsively) from the cluster such that some or all
stars are left unbound and escape (Hills, 1980). The only situation where this can
happen is at the end of the initial formation of the star cluster from the primordial
gas cloud, hence the frequent use of the term primordial gas expulsion.

Assuming that the gas retains the same spatial profile as the young stellar
cluster, the effect of primordial gas expulsion has been studied extensively in pure
N-body simulations, where the stars are represented by a large number of gravi-
tationally interacting bodies and the gas by a smooth potential that is varied in
time (e.g., Portegies Zwart et al., 2010; Banerjee & Kroupa, 2017, for reviews).
<figure><img src="image_5.png" /></figure>

**Fig. 4** Two time frames of a numerical simulation of the evolution of a molecular cloud with stellar feedback. Upper panels are maps of the column density. Lower panels are the corresponding *x*-profiles of the gravitational potential at the *y* position of each sink particle. As the expansion of the HII region proceeds, the gravitational potential is flipped-up, and thus, the stars are pulled out toward the edges. Figure from Zamora-Avilés et al. (2019).

Baumgardt & Kroupa (2007) show in a large parameter study that most clusters are completely destroyed or lose a substantial number of stars. Those that survive have expanded by a typical factor of 3-4. More recent N-body simulations vary, e.g., the kinematic state at gas expulsion or the level of substructure (Smith et al., 2013; Farias et al., 2015) and find that cluster dispersal becomes more difficult in more realistic scenarios (Farias et al., 2018).

Hydrodynamic simulations by Geen et al. (2018) and Zamora-Avilés et al. (2019, Fig. 4) have shown that the dispersal of the parental molecular cloud could have a “gravitational feedback” effect on the newborn stellar cluster: feedback from the newborn massive stars expels the gas from the collapse centre. Since neither the parental clouds, nor the formed shells are distributed symmetrically around the HII region, net forces can even accelerate the stars towards the edges of the cavity and may produce a “Hubble flow-like” ($v \propto r$) expansion.

### 5.2.2 Observed kinematics in young star clusters

Several candidates for stellar groups undergoing expansion or dispersal related to gas expulsion have been found with Gaia Data Release 2, for clusters with
masses up to $10^4 M_\odot$: Kuhn et al. (2019) study the kinematics of 28 young stellar groups with typically 100 stars with proper motion measurements each. For 75% of their objects, they find a positive offset of the generally Gaussian distributions of the cluster-centric radial velocities, i.e. an expansion of the system. Some of their groups are likely unbound and may have formed as associations (compare also Bravi et al., 2018; Wright et al., 2019), while some could have undergone an expansion phase and are settling in virial equilibrium. There is an interesting variety in expansion states also for compact systems. Kuhn et al. (2019) find the partially embedded Orion Nebula cluster (ONC) ($r_h = 0.9$ pc) to be only slowly expanding ($v_{\text{out}} = 0.43 \pm 0.20$ km s$^{-1}$), but the likewise partially embedded cluster Cep B ($r_h = 1.4$ pc) to be expanding at more than twice this rate and to clearly show a Hubble-law-like behavior (compare §5.2.1) as expected to develop in dispersing star clusters. Karnath et al. (2019) provide detailed kinematics for two expanding sub-clusters of Cep OB3b, arguing that the clusters might lose 25-65% of their stars, before re-settling in virial equilibrium.

Young ($\approx 10$ Myr), exposed, massive ($\gtrsim 10^4 M_\odot$) star clusters also appear frequently to have velocity dispersions above the expectation for virial equilibrium, given the mass expected for the observed luminosity and age (e.g., Bastian & Goodwin, 2006; Goodwin & Bastian, 2006; Gieles et al., 2010a; Portegies Zwart et al., 2010). This has been discussed as evidence for dissolution after gas expulsion (Goodwin & Bastian, 2006). However, N-body simulations show that many of these clusters would have re-virialised by the time of observation (Baumgardt & Kroupa, 2007; Gieles et al., 2010a; Portegies Zwart et al., 2010). An interpretation in terms of a large contribution from binaries to the velocity dispersion (compare, e.g., Leigh et al., 2015; Oh et al., 2015) seems more plausible (Gieles et al., 2010a; Cottaar et al., 2012; Hénault-Brunet et al., 2012).

### 5.2.3 Gas expulsion in massive star clusters

It can be shown that there exist a critical compactness $M/r_h$ above which gas expulsion with associated dispersal of stars can no longer work in a star cluster even if the gas dominates the gravitational potential at the time when massive star feedback becomes effective: while the gravitational binding energy $E_b$ is proportional to$^5$ $(1 - \epsilon_{\text{SF}})M^2/r_h$, the cumulative feedback energy by winds and supernovae at any given cluster age is only linear in the mass: $E_f \propto \epsilon_{\text{SF}}M$. Therefore, gravity must eventually win.

If we demand that for successful gas expulsion to happen, the provided feedback energy must exceed a critical energy proportional to the binding energy, i.e.,

$$E_f > a^{-1}E_b, \qquad (7)$$

with a constant $a^{-1}$ that will depend on the details of the feedback physics, then we can derive a critical star formation efficiency, defined here as the ratio of stellar mass in the cluster to its total mass during the embedded phase, for gas expulsion to succeed:

$$\epsilon > \epsilon_{\text{crit}}(C_5) = aC_5 \left( -\frac{1}{2} + \sqrt{\frac{1}{4} + \frac{1}{aC_5}} \right), \qquad (8)$$

<footnote>$^5$ $\epsilon_{\text{SF}}$: stellar mass M over total mass (stars + gas) of an embedded cluster</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Gas dynamical constraint on the star formation efficiency for successful gas expulsion. If the star formation efficiency, i.e. the ratio between stellar and total mass in the embedded cluster, is less than the given value, the given type of feedback will not be able to expel the gas on the crossing timescale of the cluster. N-body simulations find an upper limit of the star formation efficiency of 50%, if gas expulsion is to significantly affect the stars (non-shaded region). This results in an upper limit on the compactness of a star cluster, $C_5$, above which gas expulsion cannot lead to significant expansion, loss of stars or dispersal. Curves for three assumptions on the type of feedback responsible for the gas expulsion are shown: Stellar winds at metallicity [Fe/H] = -1.5 (black), supernovae (green) and a burst of hypernovae (red), as an extreme upper limit for gas expulsion via stellar feedback. Based on Fig. 4 in Krause et al. (2016).</figcaption></figure>

where we have defined the compactness index as

$$C_5 = \frac{M/r_h}{10^5 M_{\odot} \text{pc}^{-1}} = \left( \frac{\sigma}{7.5 \text{ km s}^{-1}} \right)^2 \quad (9)$$

and used eq. (4) in the final equality above.

The function $\epsilon_{crit}(C_5)$ tends towards zero for small $C_5$ ($\sigma^2$) and towards one for very high cluster compactness. Krause et al. (2016) have shown that a thin-shell superbubble model reproduces this equation (Fig. 5).

Thin-shell superbubble models (Krause et al., 2012, 2016) compute the kinematics of the supershell assuming some prescription for the energy input and spherical symmetry. They can, however, take 3D effects into account by evaluating a criterion for the shell's acceleration. The shell will be destroyed by the Rayleigh-Taylor instability as soon as modes comparable to the size of the shell become unstable. The hot, pressurised bubble interior then escapes through holes in the shell, and the dense shell gas falls back. The more stars there are compared to the amount of gas, the stronger the feedback, and the easier to push out the gas without making the shell unstable.

Successful gas expulsion therefore requires the star formation efficiency to be above a certain limit. There is, however, also an upper limit ($\approx 50\%$), if one wants gas expulsion to affect the stellar population. Both constraints together imply that only clusters with a compactness below a critical value can suffer expansion or even dispersal due to gas expulsion. For both, solar metallicity winds and supernovae, Krause et al. (2016) show that the critical compactness index is $C_5 \approx 1$ ($\sigma = 7.5 \text{ km s}^{-1}$, also compare Fig. 5).

The thin shell models effectively correspond to the assumption of maximum efficiency for stellar feedback: the hot gas is assumed to be always more central than the cold gas, thus maximising the outward push on the cold gas. The spherical shell prevents hot gas from escaping, thus all of it can be used to act on the
cold gas. Finally, Krause et al. (2012, 2016) assume 80% of the released feedback energy to be radiated away, thus 20% to be available for gas dynamics. This is likely a generous assumption, given the high efficiency of mixing and associated radiative losses seen in recent 3D superbubble simulations with time-dependent driving (Krause et al., 2013a; Vasiliev et al., 2017; Gentry et al., 2019).

Stellar winds become less efficient at low metallicities $Z$, their energy output scaling with $Z^{0.7}$ (Maeder & Meynet, 2012). If stellar winds (augmented by photoionisation and radiation pressure effects further away from the massive stars, compare below) dominate feedback in young clusters, rather than supernovae, for which there is some evidence from the timescales observed for massive clusters to become exposed (Hollyhead et al., 2015; Sokal et al., 2016; Kruijssen et al., 2019; Chevance et al., 2019), the critical compactness index becomes smaller at low metallicities, $C_5 = 0.3$ at $[\text{Fe/H}] = -1.5$. It is also possible to increase it by extreme assumptions on stellar feedback. If the most massive stars in a cluster exploded as hypernovae, all releasing ten times the conventional supernova energy output of $10^{51}$ erg (e.g., Mazzali et al., 2014; Lü et al., 2018), this would increase the critical compactness index to $C_5 \approx 30$. A more comprehensive analytic treatment by Matzner & Jumper (2015) that takes into account accretion and various feedback processes separately find the threshold at $3 \text{ km s}^{-1}$ ($C_5 = 0.2$).

### 5.2.4 Slow gas clearance

Crocker et al. (2018) consider the effect of the radiation pressure taking into account re-radiated infrared radiation due to the presence of dust. They argue that indirect radiation pressure on dust would first expand the gas gently, and that direct radiation pressure would later, but still before the first supernova, expel the gas on the dynamical timescale. For favourable assumptions, they find a maximum stellar surface density of $10^4 M_{\odot} \text{pc}^{-2}$ at which up to $\epsilon_{\text{SF}} = 50\%$ of the mass in a cluster can be stars without them forcing the remaining gas mass out of the cluster. Taking a typical cluster radius of 1 pc converts this result to a compactness index $C_5 = 0.3$. Hence, radiation pressure is expected to be somewhat less effective at expelling gas than stellar winds, but of comparable order of magnitude (compare also Reissl et al., 2018).

Rahner et al. (2017) use a self-gravitating thin-shell model to predict gas removal in clusters with $0.05 < C_5 < 100$. They find comparable contributions from radiation pressure, winds and supernovae, with radiation pressure dominating at the high-mass end (compare also Kim et al., 2016). Rahner et al. (2019) use an updated treatment of the hot gas pressure, similar to Krause et al. (2012, 2016). They also investigate the effect of the power-law slope of the initial gas distribution outside the core radius. For $0.025 < C_5 < 2.5$ they find minimum star formation efficiencies for gas removal in the percent range, with the exception of their most concentrated clouds with their steepest gas density power-law index of -2, where it can reach 50%. As they consider only gas removal on any timescale, and not the specific condition for gas expulsion on the dynamical timescale, their critical star formation efficiencies are somewhat lower than the ones of Krause et al. (2012, 2016) despite them including radiation pressure into their calculations.
### 5.2.5 Gas clearance in 3D hydrodynamics simulations

Multi-dimensional simulations of star cluster formation that take into account the actual formation of the stars and follow feedback from individual massive stars typically find a reduced effect of feedback compared to the much more idealised works above. Dale et al. (2015) study the evolution of a turbulent molecular cloud with photoionisation and conservatively implemented stellar wind feedback using smooth particle hydrodynamics. They report a variety of conditions for star formation, including tenuous and very dense regions, with the overall number of expelled stars remaining low. Gavagnin et al. (2017) conducted a similar study using adaptive mesh refinement hydrodynamics together with photoionisation from individual stars in an initially subvirial cloud. They report runs with different feedback strength. The fraction of unbound stars depends only weakly on the feedback strength, and ejections are mainly due to gravitational star-star interactions. Surprisingly, their star cluster without feedback disperses at the end of the simulation, whereas the cluster with the strongest feedback forms a subvirial system, despite 80% of the gas being ejected. This is, because the feedback efficiently slows the overall collapse, such that the stellar density remains lower, and less dynamical interactions between stars take place.

The accuracy with which the strength of feedback is predicted by these models may be subject to further improvement. The different feedback processes (accretion radiation, protostellar jets and outflows, photoionisation, radiation pressure, stellar winds and supernovae) require very different computational methods. The simulations discussed above all include photoionisation. Dale et al. (2015) exclusively use the momentum from stellar winds. This is an underestimate, because the energy in the winds will not be entirely radiated away, but produce some additional momentum. That the simulations generally underestimate feedback is underlined by the fact that many runs do not terminate star formation (e.g. Dale, 2017) within an observationally required time frame of 3 Myr (Chevance et al., 2019). This is particularly relevant, given that the timescale of gas loss strongly affects any expansion or dispersal (Smith et al., 2013). The virial state is expected to have a strong influence on the fraction of bound stars (Farias et al., 2015). Hence, simulations with subvirial clouds, only, (Gavagnin et al., 2017) cannot provide the full picture.

More recently, Li et al. (2019) simulated star cluster formation from turbulent clouds in different kinematic states with a moving mesh hydrodynamics code. In each run, they form a variety of stellar structures, hierarchically merging into bigger ones. They apply feedback via mass and momentum deposition around each star, which is varied within a factor of 20. The latter range reflects the still existing uncertainty on the feedback strength. Gas expulsion with associated dispersal of stars seems to occur in some of their simulations with the highest level of feedback. Most of their simulations do, however, not show a strong unbinding of stars due to bound structures being generally subvirial prior to the gas expulsion treatment.

### 5.2.6 Gas clearance summary

Analytical and semi-analytical models combined with N-body simulations tend to overestimate the effects of stellar feedback and predict strong gas expulsion effects with cluster expansion and dispersal for a virial velocity dispersion $\sigma <$
3–7 km s⁻¹. They firmly exclude strong effects of gas expulsion on the cluster stars for σ > 7 km s⁻¹ unless one assumes non-standard mechanisms. Multi-dimensional simulations tend to underestimate feedback and usually see little effects of gas expulsion and a small amount of unbound stars. However, tuning up the feedback strength, such effects have also been reported (Zamora-Avilés et al., 2019; Li et al., 2019). *Gaia* stellar kinematics observations suggest that gas expulsion may be responsible for expansion and possibly dispersal of some stellar groups, while compact clusters appear unaffected, implying slow gas outflow or exhaustion of gas turned into newly formed stars. All reported cases where expansion or dispersal may take place have σ < 3 km s⁻¹, consistent with the theoretical constraints.

These results are in good agreement with direct observations of gas and stars in young massive clusters and progenitor clouds: there are no clouds compact enough, so that a young massive cluster could form with the same structure as that cloud (Longmore et al., 2014; Walker et al., 2015, 2016). The implication is that star formation has to proceed as the cloud collapses (compare §3), which has been termed the ‘conveyor belt’ model of cluster formation (Longmore et al., 2014). As gas clouds fall in, they can already be forming stars, with a star formation efficiency that peaks within the regions of the highest densities. This leads to local gas exhaustion in these regions and limits the effects of gas expulsion on the virial state of the resulting cluster, producing high bound fractions and thus compact cluster formation (Kruijssen, 2012). Low central gas fractions in an embedded cluster are indeed reported by Ginsburg et al. (2016). Gas exhaustion has also been seen in cluster formation simulations (Girichidis et al., 2012a; Kruijssen et al., 2012; Dale et al., 2015).

5.3 Steady-state cluster winds

Star cluster winds form, where the gas in the cluster is heated faster than it can
cool. This is usually expected for the epoch just after the gas has been cleared from
the cluster, either because most of it has been accreted on to the stars (exhaustion)
or, because the feedback processes removed it from the cluster.

If there is a significant number of massive stars in the cluster, all driving winds and exploding as supernova in more or less regular intervals, one can assume that the energy is efficiently thermalised in the local interactions, and the mass input from the various wind sources to be smoothed out over the size of the cluster. Chevalier & Clegg (1985) developed a classical steady-state wind model that applies to this situation. Important assumptions in the model are:

1. Spherical symmetry.

2. The region of interest comprises a large number of stars (sources of mass and energy), such that individual sources interact locally and we can describe the gas physics using smooth mass and energy input functions *q(r)* and *Q(r)*, respectively.

3. Top-hat flat source profile, i.e.<fnref n="6" />, $q(r) = \dot{M}/V$, $Q(r) = \dot{E}/V$ for $r < R$, and $Q(r) = q(r) = 0$, otherwise.

4. Gravity is negligible.

<footnote n="6">⁶ $\dot{M}$: total mass loss rate; $\dot{E}$ total energy release rate; $V = 4\pi r^3/3$.</footnote>
Assumption 2 above restricts the theory effectively to massive star clusters (and galaxies, of course). For the wind phase, this is because of the strong dependence of the stellar wind strength on the stellar mass. For example, Krause et al. (2013a) show in 3D hydrodynamics simulations that for a group that harbours star of 25, 30 and $60M_{\odot}$ (typical for a 1000$M_{\odot}$ cluster using the initial mass function from Kroupa et al., 2013), the $60M_{\odot}$ star completely dominates the gas dynamics as long as it exists. Interacting stellar winds and supernovae will heat the cluster to typically, $10^7$ K, which corresponds to a sound speed of $\approx 500$ km s$^{-1}$. For a typical cluster diameter of, say, 10 pc, the dynamical timescale is then 20,000 yr. If we require one supernova per dynamical timescale, we need roughly 1500 massive stars, which we expect for a cluster with $\approx 10^5 M_{\odot}$. Steady-state winds in clusters are therefore frequently referred to as super star cluster winds. Cantó et al. (2000) show using 3D hydrodynamics simulations that in a cluster with 30 massive stars with similar properties, the 1D case with smooth source functions is approximately recovered.

Given these conditions, the 1D hydrodynamics equations can be solved analytically (see also Zhang et al., 2014). Pressure, density and outward velocity are given by

$$ \begin{pmatrix} p \\ \rho \\ u \end{pmatrix} = \begin{pmatrix} p_* \dot{M}^{1/2} \dot{E}^{1/2} R^{-2} \\ \rho_* \dot{M}^{3/2} \dot{E}^{-1/2} R^{-2} \\ u_* \dot{M}^{-1/2} \dot{E}^{1/2} \end{pmatrix}, \qquad (10) $$

where the functions containing the radial dependencies are given by$^7$:

$$ u_*^2 = \frac{2M^2}{M^2 + \frac{2}{\gamma-1}} \qquad (11) $$

$$ \rho_* = \frac{r_*^a}{4\pi u_*} \qquad (12) $$

$$ p_* = \frac{2\rho_*}{\gamma \left( M^2 + \frac{2}{\gamma-1} \right)} \qquad (13) $$

with $r_* = r/R$, $a = 1$ (-2) for $r_* < 1$ ($r_* > 1$), and the implicit definition of the Mach number $M$:

$$ r_* = \begin{cases} \left( \frac{\gamma-1+2M^{-2}}{\gamma+1} \right)^{\frac{\gamma+1}{2+10\gamma}} \left( \frac{3\gamma+M^{-2}}{1+3\gamma} \right)^{-\frac{3\gamma+1}{5\gamma+1}} & r_* < 1 \\ \left( \frac{\gamma-1+2M^{-2}}{\gamma+1} \right)^{\frac{\gamma+1}{4\gamma-4}} M^{\frac{1}{\gamma-1}} & r_* > 1 \end{cases} \qquad (14) $$

The solution is shown graphically in Fig. 6.

The Chevalier & Clegg (1985) solution is characterised by a slow, hot and subsonic flow inside the star cluster. The flow turns supersonic at the boundary of the source region and then continues to accelerate towards an asymptotic value of $1.414\sqrt{\dot{E}/\dot{M}}$.

As an example we give here parameters for the Arches cluster, one of the most massive, young (2-3 Myr, Lohr et al., 2018) star clusters in the Milky Way. Clark et al. (2019) estimate $\gtrsim 50$ stars with masses $\gtrsim 60M_{\odot}$. Using the initial mass function from Kroupa et al. (2013), this translates to a total mass of $5 \times 10^4 M_{\odot}$

<footnote>$^7$ $\gamma$: adiabatic index, 5/3 for the usual monatomic ideal gas</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Steady-state wind solution by Chevalier & Clegg (1985). Shown are the dimensionless quantities given in eqs. (11-13). Section 5.3 for details.</figcaption></figure>

(consistent with the kinematic measurement, Clarkson et al., 2012), and a total number of massive stars > $8M_\odot$ of ≈ 1000 (compare also Figer et al., 1999, 2002, who further give a cluster radius of 0.2 pc). Population synthesis of stellar mass loss and energy output then yields (Voss et al., 2009): $\dot{M} \approx 10^{-3} M_\odot \text{ yr}^{-1}$ and $\dot{E} \approx 2 \times 10^{39} \text{ erg s}^{-1}$ (similar estimates can be found in Stevens & Hartwell, 2003). Within the cluster this yields particle densities and temperatures of

$$n_0 = 300 \text{ cm}^{-3} \left( \frac{\dot{E}}{2 \times 10^{39} \text{ erg s}^{-1}} \right)^{-1/2} \left( \frac{\dot{M}}{10^{-3} M_\odot \text{ yr}^{-1}} \right)^{3/2} \left( \frac{\dot{R}}{0.2 \text{ pc}} \right)^{-2} \quad (15)$$

$$T_0 = 9 \times 10^7 \text{ K} \left( \frac{\dot{E}}{2 \times 10^{39} \text{ erg s}^{-1}} \right) \left( \frac{\dot{M}}{10^{-3} M_\odot \text{ yr}^{-1}} \right)^{-1} \quad (16)$$

Arches and similar young, massive and compact clusters are therefore expected to be faint diffuse X-ray emitters (e.g., Cantó et al., 2000; Añorve-Zeferino et al., 2009), which has been confirmed by X-ray observations for the Arches cluster, Westerlund 1 and possibly also the Quintuplet cluster (Yusef-Zadeh et al., 2002; Wang et al., 2006; Kavanagh et al., 2011). Generally the temperature is somewhat lower than predicted by the Chevalier & Clegg (1985) model (Stevens & Hartwell, 2003). This may be related to unaccounted for effects of non-equilibrium ionisation (Ji et al., 2006) or limitations of our understanding of mass loading and thermalisation efficiency of the winds. Also, the spatial distribution of the massive stars plays a role. X-ray emission is also expected from the interaction of the cluster wind with the surrounding gas. The superbubble is particularly bright in soft, ≈ 1 keV, X-rays, whenever individual supernova shock waves interact with the shell (Krause et al., 2014). Cluster winds can, however, be identified by their harder spectra and co-location with the optical star cluster (Silich et al., 2005).

The basic wind solution has been modified and enhanced by many authors for example to include cooling (Silich et al., 2004) and more sophisticated shapes for the source functions (Silich et al., 2011; Palouš et al., 2013; Zhang et al., 2014). Wünsch et al. (2011) present detailed models for the entire phase when massive stars (> $8M_\odot$) are present in the cluster. They find that clusters are always in a steady outflow regime similar to the Chevalier-Clegg model, unless the energy input is significantly overestimated (factor ≳ 20) by current population synthesis models or the wind loads a significant amount of gas that was leftover from the star formation event. If those conditions applied, part of the massive star ejecta
would cool, be compressed by the remaining hot gas into UV-shielding filaments and form stars in an extended or second star formation episode (Palouš et al., 2014). The total cold gas dropout from the wind can reach 1-6% of the total stellar mass of a cluster (Wünsch et al., 2017). The population of stars formed would be (moderately due to the mass loading) enriched in He-burning products ejected from Wolf-Rayet stars and supernovae (Wünsch et al., 2011). The latter two predictions disfavour this mechanism as explanation for the frequently observed chemically distinct populations in globular clusters (compare §7).

For a normal stellar population, type Ia supernovae appear from about 100 Myr after a star formation event at a rate characterised by the delay time distribution (Heringer et al., 2019):

$$ \text{DTD}(t) = 7 \times 10^{-13} M_{\odot}^{-1} \text{yr}^{-1} \left( \frac{t}{\text{Gyr}} \right)^{-1.34} . \qquad (17) $$

For a $10^6 M_{\odot}$ cluster at an age of 100 Myr, this yields about 15 events per Myr. This is near the limit where one might consider the energy injection continuous and apply cluster wind models. D’Ercole et al. (2008) show in a 1D hydrodynamic simulation with individual SN Ia that even one event can turn the cluster into an outflow state. A SN Ia rate comparable to the one in the field would leave star clusters in a continuous outflow state.

SN Ia occur in binary systems (e.g., Diehl et al., 2014), which may be more frequent in massive star clusters (Leigh et al., 2015). Direct searches for type Ia SNe in massive star clusters have, however, so far only produced upper limits (Washabaugh & Bregman, 2013). Dynamical effects should lead to a high net destruction rate for binaries in clusters, so that, at least at late times, they may have actually fewer SN Ia than the field (Cheng et al., 2018; Belloni et al., 2019).

## 5.4 Cooling flows?

After the end of the type II supernova phase (≈ 30–40 Myr after star formation) and before SN Ia start to occur (≈ 100 Myr, e.g., Liu & Stancliffe, 2018) a star cluster has little internal energy production and may in principle have a cooling flow. D’Ercole et al. (2008) show that mass loss and energy injection are dominated by AGB stars:

$$ \alpha = \dot{M}/M = 3 \times 10^{-17} \text{s}^{-1} \qquad (18) $$

$$ \dot{E}/M = \alpha v_w^2 / 2 = 10^{29} \text{ erg s}^{-1} \text{ M}_{\odot}^{-1} \qquad (19) $$

The cooling flow is robust for these parameters as they also show that more than ten times higher energy input would be required to turn the cooling flow into a wind.

D’Ercole et al. (2010, 2012, 2016) and others have argued that the gas mass flowing to the centre may initiate secondary star formation. Challenges in explaining the chemically distinct multiple populations in globular clusters (compare §7) in this way include the small mass of enriched material available and fine-tuning required for the dilution of the ejecta.

Conroy & Spergel (2011) have conjectured that the gas in the cooling flow would actually not be able to form stars for several 100 Myr, because the UV flux
of the remaining intermediate-mass stars would keep the gas photo-dissociated and too warm for star formation. The accumulating gas would only form stars when the UV luminosity has declined enough to allow the formation of molecular hydrogen. However, gas cooling can also take place very efficiently in atomic gas via $C^+$ (Glover & Clark, 2012). Also, it is unclear if type Ia SNe would be delayed sufficiently for the model to work (Lyman et al., 2018).

Dense gas or late star formation as postulated in the above cooling flow models is generally not observed in star clusters (e.g., Cabrera-Ziri et al., 2015; Longmore, 2015; Bastian & Lardo, 2018). This calls into question our understanding of the gas dynamics in star clusters with ages between the type II and type Ia supernova phases. One possibility is that the cooling flow gas accretes on to the dark remnants, i.e., the stellar mass black holes and neutron stars (Krause et al., 2013b; Roupas & Kazanas, 2019). The energy released in jets, winds and radiation could then drive a cluster wind. D'Ercole et al. (2008) derive a critical energy input of $6 \times 10^{37}$ erg s$^{-1}$ for their $10^7 M_{\odot}$ cluster. Even the typical luminosity of one X-ray binary (few $10^{38}$ erg s$^{-1}$, Jordán et al., 2004) would be sufficient to accomplish this. Pulsar winds can keep star clusters in an outflow state (Naiman et al., 2020).

# 6 Collisional dynamics and long-term evolution

After the gas cloud a star cluster formed from has been partially transformed into stars and dispersed, its fate is governed by gravity (i.e. collisional dynamics and tidal perturbations) and mass loss of the stars due to stellar evolution. Here we discuss the various physical processes separately, but it is important to keep in mind that most processes act simultaneously and an important area of research is understanding the interplay between them, which is often non-linear.

## 6.1 Stellar evolution

Stellar evolution leads to a decrease of the total cluster mass, at a rate that is slow compared to the orbital frequencies of the stars, such that the cluster can approximately maintain its virial equilibrium. The removal of mass leads to a reduction of the binding energy and an increase of the cluster radius. If the stellar mass loss happens throughout the cluster with no preferred location then the cluster radius is inversely proportional to the mass (Hills, 1980). For mass segregated clusters most stellar mass loss occurs in the cluster centre where the binding energy is larger, resulting in a faster expansion.

## 6.2 Tidal shocks

During the first 0.1-1 Gyr, cluster dissolution is likely dominated by tidal 'shocks', i.e. impulsive tidal perturbations from Galactic substructure, such as transient spiral arms and molecular gas clouds (e.g. Gieles et al., 2006; Elmegreen & Hunter, 2010; Kruijssen et al., 2011). These perturbations boost the energy of stars in the cluster, some of which will exceed the escape energy and will therefore become unbound (Spitzer, 1958). The rate of shock-driven mass loss scales inversely with
the mass volume density of the cluster, and is proportional to the surface density of
the individual clouds and the ISM density in the host galaxy disc (Spitzer, 1958).

When integrated over the lifetime of a cluster, this mass loss mechanism could dominate the total mass loss budget (Elmegreen, 2010; Kruijssen, 2015), even in environments of relatively low gas density such as the solar neighbourhood (e.g. Spitzer, 1958; Gieles et al., 2006; Lamers & Gieles, 2006), where a single encounter with a GMC ($\gtrsim 10^5 M_{\odot}$) can completely disrupt a modest open cluster ($\sim 10^3 M_{\odot}$, Wielen, 1985; Terlevich, 1987). By scaling $N$-body models of individual encounters, it was found that in gas-rich environments like galaxy discs, GMCs dominate the disruption of clusters (Gieles et al., 2006; Webb et al., 2019), decimating the initial globular cluster population to the survivors that remain at the present day (Elmegreen, 2010; Kruijssen, 2015).

Tidal shocks do not only drive considerable mass loss, they also dominate the structural evolution of stellar clusters: after an initial phase of expansion due to the escape of unbound stars (Webb et al., 2019), the remaining cluster of bound stars may shrink due to energy conservation (centrally concentrated clusters shrink, while low-concentration clusters expand, see Gieles & Renaud, 2016). When ignoring other effects, a density increase makes tidal shocks self-limiting (Gnedin & Ostriker, 1999). However, a higher density makes two-body relaxation more important which tends to reduce the cluster density, thereby counteracting the shock-induced density increase. Under the assumption that statistical equilibrium is reached, eventually the ratio of the shock dissolution time-scale and the relaxation time-scale will become constant, resulting in a shallow mass-radius relation ($r_h \propto M^{1/9}$, Gieles & Renaud, 2016). The normalisation of the predicted mass-radius relation depends on the environment, such that clusters are smaller at higher ISM densities (which is likely already the case at formation, see Choksi & Kruijssen 2019), slowing down their shock-driven disruption. However, even for correspondingly more compact clusters, Kruijssen (2015) predict that the total shock-driven mass loss dominates over relaxation-driven mass loss when considering the dynamical evolution of globular clusters over a Hubble time.

To obtain a complete understanding of the interplay between shocks and relaxation, a comprehensive parameter study of $N$-body simulations including both processes is required, which again highlights that this is an important area for future research. A complementary approach to controlled $N$-body experiments would be to use direct $N$-body simulations with realistic particles numbers ($N \gtrsim 10^6$) and evolve them in the time-dependent tidal field extracted from models of galaxy formation at the epoch of GC formation.

## 6.3 Two-body relaxation

The importance of collisional dynamics in the evolution of star clusters depends on the evolutionary stage of the cluster and the timescale that is considered. Given sufficient time, all clusters will dissolve due to collisional effects, even the clusters that are not in a Galactic tidal field (Baumgardt et al., 2002).

Globular clusters are the archetypical collisional systems, that survived the initial phase of tidal shock-driven disruption, meaning that orbital energy diffusion via gravitational interactions – so-called two-body relaxation, or collisional dynamics – plays an important role in their evolution. This is because the velocities of
stars are relatively low (~10 km/s) and stellar densities are high (~10$^{4-6}$ pc$^{-3}$), making two-body encounters frequent and long-lasting. Another way of saying this is that the relaxation timescale is short (few Gyr) compared to their ages (10-12 Gyr). Two-body relaxation is also relevant during the formation phase of clusters, contrary to some propositions made in the literature (e.g. Fall & Zhang, 2001; Krumholz et al., 2014). To explain this, we start by painting a broad-brush picture of the classical theory of relaxation that was developed for (old) globular clusters.

The consequences of two-body relaxation are reasonably well understood for the idealised case of a single-mass cluster, which is often regarded to be a reasonable approximation for globular clusters, because their stars are confined to a narrow range of masses. A single-mass cluster, with stars initially in hydrostatic equilibrium, without primordial binaries, develops a radial energy flow as a result of energy diffusion such that energy flows through the half-mass radius ($r_h$) at a rate $\sim |E|/\tau_{rh}$ (ignoring constants of order unity). Here $E \sim -GM^2/r_h$ is the total energy of the cluster, with $G$ the gravitational constant and $M$ the total cluster mass. The timescale $\tau_{rh}$ is the half-mass relaxation time, which we shall define below (equation 22). This energy flow originates from the core, where stars lose kinetic energy to the stars outside the core via two body interactions. As a result, the velocity dispersion of the stars in the core reduces – i.e. they ‘cool’ – and the core radius contracts, while the stars outside the core heat up. The stars in the core now experience a higher binding energy and due to the virial theorem, the stars now move faster. This somewhat paradoxical result is a direct consequence of the negative heat capacity of self-gravitating systems. The time evolution of the cluster structure can be solved in various ways. Lynden-Bell & Eggleton (1980) used a set of equations that are strikingly similar to the stellar structure equations (so-called gaseous models, or continuum models) and found that this process of core contraction continues until the core has an infinite density and zero mass: the core has collapsed. In reality this mathematical endpoint is never reached, and a binary star forms when the number of stars in the core has reduced to $\sim 10$ (Goodman, 1987). To sustain the two-body relaxation process after core collapse, the freshly formed binary absorbs the negative energy that continues to flow into the core, by hardening in interactions with other stars. From that moment onwards, the evolution of the cluster is approximately self-similar, with the rate of energy absorption of the central binaries determined by the global energy flow through the cluster (Hénon, 1961)

$$ \dot{E}_{\text{bin}} \propto -\frac{|E_{\text{ext}}|}{\tau_{\text{rh}}} \qquad (20) $$

Here $E_{\text{ext}}$ is the external energy of the cluster, which is the total energy excluding the negative energy locked up in binaries (i.e. $E_{\text{ext}} \sim -GM^2/r_h$). Because of energy conservation, the external energy must increase at a rate

$$ \dot{E}_{\text{ext}} = -\dot{E}_{\text{bin}} \qquad (21) $$

The insight that there must exists a balance between the rate of energy production in the core and the energy flow through the cluster came from Michel Hénon (Hénon, 1975) and is a fundamental building block in the theory of cluster evolution. It allowed him to derive two models for the post-collapse evolution: in
<figure><img src="image_8.png" /><figcaption>Fig. 7 The half-mass density of Milky Way globular clusters as a function of their mass (x-axis) and Galactocentric radius (colour coding). ‘Isochrones’ from the cluster evolution model of Gieles et al. (2011) are shown. This diagram is the equivalent of the Hertzsprung-Russell diagram of stars. Collisional dynamics is not yet important for objects that are well below the model lines.</figcaption></figure>

the absence of a Galactic tidal field, the energy increase leads to an expansion of the cluster at an approximately constant mass (Hénon, 1965), while tidally limited clusters lose mass over the tidal boundary (sometimes referred to as ‘evaporation’) at a constant rate, while maintaining a constant density (Hénon, 1961). These two solutions describe the two extreme ends of the life cycle of tidally limited star clusters with high initial density. By smoothly ‘stitching’ the two models the relaxation driven evolution of star clusters can be described from the moment they emerge from a gas-rich environment (i.e. once tidal shocks no longer dominate the instantaneous disruption rate) to their eventual dissolution (Gieles et al., 2011). Thanks to its analytic nature, this simple model for the relaxation driven evolution of star clusters readily provides expressions for $M(t)$ and $\rho_h(t)$ at different Galactocentric radii, which can be used to construct evolutionary ‘tracks’ and ‘isochrones’ of globular cluster radius (or density) and mass as a function of location in the Galaxy, which – despite their first order nature – provide a satisfactory match with the observed mass-density distribution of globular clusters (see Fig. 7), supporting Hénon’s suggestion that collisional dynamics is important for almost all globular clusters and shapes these relations.

## 6.4 Two-body relaxation in young clusters

We now turn to the relevance of relaxation for young clusters. To get an idea for the physical time it takes for relaxation to become important, we write the expression for $\tau_{rh}$ from Spitzer & Hart (1971) as

$$ \tau_{rh} \approx 18 \text{ Myr} \psi^{-1} \frac{M}{10^4 M_{\odot}} \left( \frac{\rho_h}{10^4 M_{\odot}/\text{pc}^3} \right)^{-1/2}, \quad (22) $$
where $\rho_h = 3M/(8\pi r_h^3)$ is the average mass density within $r_h$. To derive this, we assumed an average mass of stars of 0.5 $M_\odot$ and assumed that the slight dependence of the Coulomb logarithm on the number of stars can be neglected ($\ln \Lambda = 8$). The term $\psi \ge 1$ is a factor that depends on the mass spectrum within $r_h$, and often the assumption of single-mass clusters is made (i.e. $\psi = 1$). It takes about 16 initial $\tau_{rh}$ for a single-mass cluster to reach core collapse (e.g. Cohn, 1980) which even for the relatively low-mass and high-density scaling adopted in equation (22) corresponds to a relatively long time of $\sim 300$ Myr, i.e. well after star formation ceased, suggesting that relaxation plays no role in the early evolution of clusters.

However, there are several important differences between young clusters and their older counterparts that are important for this discussion and make collisional dynamics important at young ages. Most importantly, young clusters have a well populated (initial) stellar mass function between $\sim 0.1 M_\odot$ and $\sim 100 M_\odot$. The presence of the high-mass stars significantly speeds up the relaxation process, because the energy transfer from high-mass stars to low-mass stars is more efficient than between stars of the same mass$^8$. With numerical N-body experiments, Portegies Zwart & McMillan (2002) find that core collapse happens after $0.2\tau_{rh}$ (with $\tau_{rh}$ defined with $\psi = 1$), i.e. 2 orders of magnitude faster than for single-mass clusters (Gieles et al., 2010b), which corresponds to $\sim 3$ Myr for the fiducial case in equation (22), and earlier for lower-mass/denser clusters. The definition of core collapse is often taken to be the moment of dynamical formation of the first hard binary, which means that in the first 3 Myr the most massive stars are migrating to the centre of the cluster (if they did not form there), making collisional dynamics important already before core collapse.

There are several reasons to assume that binary activity starts even earlier. Firstly, clusters may form hierarchically from mergers of lower-mass and denser sub-clusters, which can each dynamically form a binary. To quantify the importance of this, lets assume that the cluster formed from two clumps with each half the mass of the cluster. If they merge with negligible orbital energy, we find from conservation of energy that the clumps have half the radius of the final cluster, i.e. 4 times shorter relaxation time, such that the clumps undergo core collapse within a Myr. Considering an additional step in the hierarchy would reduce $\tau_{rh}$ of the first sub-clumps to form to $\sim 10^5$ yr, i.e. well within the timescale of cluster formation itself. In addition, the massive stars may form in the centre of the sub-clumps/cluster (see Section 6.5), resulting in a cluster forming in a core collapsed state, setting off binary activity immediately. Finally, several models of the formation of star clusters (Longmore et al., 2014; Vázquez-Semadeni et al., 2017; Gieles et al., 2018; Krumholz & McKee, 2019) and massive stars (Padoan et al., 2019) suggest that gas inflow from larger scale is important (compare §3). Accretion of low-angular momentum gas on proto-stars leads to an efficient contraction of the parent cluster (Bonnell et al., 1998), driving the cluster into core collapse by reducing $\tau_{rh}$ as the cluster gains mass (Moeckel & Clarke, 2011).

<footnote n="8">⁸ The single-mass approximation, therefore, also breaks down for globular clusters with a stellar-mass black hole population (Breen & Heggie, 2013; Giersz et al., 2019; Kremer et al., 2019; Wang, 2020; Antonini & Gieles, 2020).</footnote>
## 6.5 Mass segregation in young clusters

Observational support for the importance of collisional dynamics at young ages comes from the (low) densities of young-stellar objects in the solar neighbourhood, which is consistent with the density distribution of a population of dynamically expanding (i.e. post-collapse) low-mass star clusters (few 100 stars, Gieles et al., 2012).

Collisional dynamics will equalise “temperature” between stars and hence the kinetic energy per star. More massive stars will therefore acquire lower velocities to have the same kinetic energy as the lower mass stars. The high-mass stars will therefore be found deeper in the potential well, i.e., closer to the centre of the cluster. This is called mass segregation. Dynamical mass segregation is a consequence of two-body relaxation, whereas primordial mass segregation describes a situation where the massive stars already form in the centre.

The observational assessment of mass segregation is mixed. Young star clusters show signs of mass segregation (Hillenbrand & Hartmann, 1998; de Grijs et al., 2002; Littlefair et al., 2003; Stolte et al., 2005, 2006; Kim et al., 2006; Harayama et al., 2008; Espinoza et al., 2009; Bontemps et al., 2010; Gennaro et al., 2011). However, observations of pre-main sequence stars in star-forming regions do not indicate mass segregation (Parker et al., 2011, 2012; Gennaro et al., 2017; Parker & Alves de Oliveira, 2017; Dib et al., 2018). Plunkett et al. (2018) find the prestellar cores in Serpens South to be mass-segregated, whereas the pre-main sequence stars are not. The complexity of defining the details of clusters and subclusters in combination with measurement differences complicate the problem of the distribution of masses. Kirk & Myers (2011) observe mass segregation in small groups in Taurus. Investigating the stellar distribution in total, Parker et al. (2011) find the most massive stars to be inversely mass segregated.

Parker et al. (2015) point out that the degree of mass segregation depends on how it is measured in simulations and compared to observations.

It can conveniently be computed using the method of Allison et al. (2009a,b), which is based on the minimum spanning tree (Gower & Ross, 1969). Kirk et al. (2014) showed that small clusters in hydrodynamic simulations exhibit primordial mass segregation with distributions consistent with nearby, young embedded clusters. Parker et al. (2015) found that the degree of mass segregation is reduced if the clusters form under the influence of feedback from massive stars. Girichidis et al. (2012b) report that the degree of mass segregation depends on the initial density configuration, but that no inverse mass segregation occurred. Clusters undergoing competitive accretion are expected to be primordially mass-segregated (Bonnell et al., 2001, cf. Bate 2009). In all cases the time scales are consistent with dynamical relaxation times, so all clusters had enough time to dynamically mass segregate.

## 6.6 Dynamical feedback

We conclude this section by discussing the feedback from collisional dynamics on star formation. The dynamically formed binaries consist of massive stars, which soon after formation start ejecting other (massive) stars in dynamical interactions
(Poveda et al., 1967, see also Gavagnin et al. 2017), possibly explaining the origin of the O-stars that are found with high velocities (≥ 30 km s⁻¹), far from star forming regions (Blaauw, 1961). In addition, a large fraction if not all massive stars are expected to form in binaries and higher order multiples (e.g., Sana et al., 2012), which have larger gravitational cross section than single stars making binary-binary interactions an additional channel for ejecting massive stars from an ongoing cluster formation site (Leonard & Duncan, 1990). The removal of massive stars reduces the mechanical and radiation feedback from massive stars on the cluster and the more distributed feedback in the low(er) density ISM has consequences for galaxy formation (Ceverino & Klypin, 2009). Finally, the high central density of massive stars in the centre affect the ionisation level (and thus accretion rate) and survival of discs around smaller mass stars

In conclusion, collisional dynamics is likely important from the very beginning of cluster evolution and it may have played a role in the origin of the multiple populations in GCs (see §7.3.3), with tidal perturbations being an additional important process in the early evolution. It is unclear, if observed signs of mass segregation are of dynamical or primordial origin. Mass segregation and ejection of massive stars modify ionisation levels in accretion discs and hence accretion rates on to stars and the strength of feedback from star clusters.

# 7 Nucleosynthesis

For a long time, nucleosynthesis (or, in other words, internal chemical evolution) has been ignored in star cluster modelling, based on both theoretical and observational arguments. Galaxies have a deep potential well and are hence expected to retain even some of the ejecta that massive stars shed at high velocity. This has recently been confirmed by measurements of Doppler kinematics of the radioactive decay line of unstable $^{26}$Al, which traces high-mass star ejecta (Kretschmer et al., 2013; Krause et al., 2015; Rodgers-Lee et al., 2019). The high observed velocities suggest that a large fraction of the ejecta is blowing away from their birth places at high speeds. The scale height of the order of kpc is in agreement with expectations from fountain-flow super-bubbling disc models where ejecta diffuse into the hot halo and return in part on a Gyr timescale (Pleintinger et al., 2019; Rodgers-Lee et al., 2019).

The need for a sufficiently deep potential well to retain the gas despite the energetic feedback from the massive stars and eventually recycle it internally to make new stars is supported by the fact that open clusters present no (within measurement uncertainties) spread in Fe-peak, α, and s-process elements (hereafter heavy metals). These specific species actually vary only in the most massive globular clusters (hereafter GCs), with the extreme case being Omega Cen which is thought to be the remnant of a dwarf galaxy nucleus (Butler et al., 1978; Zinnecker et al., 1988; Marino et al., 2011). Such rare objects (recently called Type II GCs; see e.g. Milone et al. 2017 and Marino et al. 2018) possibly make the link between star clusters (open clusters and Type I GCs) and chemically evolved dwarf galaxies.

Interestingly, large variations in carbon, nitrogen, oxygen, sodium, magnesium, and aluminium (C, N, O, Na, Mg, Al, hereafter light elements) were discovered in GCs already in the 1970's (among bright red giants, Osborn 1971; see Kraft
1979, and references therein), but they were initially attributed to internal deep
mixing processes occurring along the evolution of the stars themselves (Sweigart
& Mengel e.g. 1979; Denisenkov & Denisenkova e.g. 1990; Langer et al. e.g. 1993;
Cavallo et al. e.g. 1996; Denissenkov & Weiss e.g. 1996; Da Costa e.g. 1997; Weiss
et al. e.g. 2000; but see e.g. Peterson 1980 and Brown & Wallerstein 1992 who
already advocated for a primordial origin). Hence, the “classical paradigm” became
established, presenting individual GCs as the archetype of a single, coeval, and
chemically homogeneous stellar population, i.e., a system that did not undergo
any internal chemical evolution.

7.1 Light element abundance variations

The surprise came in the 2000's from studies with 8-10m class telescopes which
opened a spectroscopic window on less evolved stars down to the main sequence
turnoff in the case of the closest GCs (Gratton et al., 2001; Thévenin et al., 2001).
All the Galactic and extra-galactic GCs that were scrutinised this way were shown
to host multiple stellar populations (hereafter MSP) located all along the color-
magnitude diagram and exhibiting similar variations in light elements (for reviews
see Gratton et al., 2012, 2019; Charbonnel, 2016; Bastian & Lardo, 2018). Exten-
sive surveys established that while every Galactic and extra-galactic GC contains
main sequence and red giant stars with field-like abundances (the so-called first
population, 1P), 30 to 90 % of their hosts (the second population, 2P) actually line
up along an O-Na anticorrelation (e.g. Carretta et al., 2009a,b, 2010; Lind et al.,
2009; Gratton et al., 2012; Wang et al., 2016), a C-N anticorrelation (e.g. Norris
& Freeman, 1979; Norris & Pilachowski, 1985), and a Mg-Al anticorrelation with
a possible linkage towards Si depletion in the case of the most massive and/or
most metal-poor GCs (e.g. Norris et al., 1981; Ivans et al., 1999; Yong et al.,
2003; Carretta, 2015; Mészáros et al., 2015, 2019; Pancino et al., 2017; Mucciarelli
et al., 2018; Masseron et al., 2019). This definitively argued for a primordial origin
of these anomalies, calling for GC self-enrichment by the hot hydrogen-burning
yields of short-lived massive stars before or during the formation of the low-mass
stars that we observe today (Prantzos et al., 2007, 2017). The light element abun-
dance variations within GCs reflect in the photometric properties of their MSPs.
When one uses specific combinations of optical and ultraviolet HST filters, the
MSPs spread out along broadened or multiple sequences of the color-magnitude
diagram (CMD) and of the so-called chromosome map (e.g. Bedin et al., 2004;
Piotto et al., 2007, 2012, 2015; Marino et al., 2008, 2019; Han et al., 2009; Milone
et al., 2015a, 2017; Bellini et al., 2017; Zennaro et al., 2019, compare Adamo et
al. 2020, in prep.).

7.2 Multiple sequences in the colour-magnitude diagram

Color variations and/or separations in the CMD are used to infer He abundance
variations among MSPs (typically between 0.003 and 0.19 in mass fraction, with He
enrichment increasing with the present-day mass; Norris e.g. 2004; Piotto et al. e.g.
2005; King et al. e.g. 2012; Sbordone et al. e.g. 2011; Milone e.g. 2015; Nardiello
et al. e.g. 2015; Milone et al. e.g. 2018; Lagioia et al. e.g. 2019). Importantly,
the photometric approach revealed the presence of MSPs similar to GC ones in extragalactic massive star clusters with ages down to ~2 Gyr (Larsen et al., 2014; Bastian, 2016; Dalessandro et al., 2016; Niederhofer et al., 2017a,b; Martocchia et al., 2018, 2019; Gilligan et al., 2019; Nardiello et al., 2019), with the extent of the MSP increasing with the age of the clusters (Martocchia, 2019). It thus seems that the formation of MSPs was not restricted to old GCs, but that it continued to occur in sufficiently massive star clusters down to a redshift of ~0.17. This provides a very strong link between YMSC and ancient GCs, and suggests a common formation and evolution path.

The exact shape and the extension of these features vary from GC to GC, and depend primarily on their mass, compactness, metallicity, and age (e.g. Carretta et al., 2009b; Krause et al., 2016; Pancino et al., 2017; Masseron et al., 2019; Cabrera-Ziri, 2019; Martocchia, 2019), and these abundance patterns have never been found in open clusters (Pancino et al., 2010; Bragaglia et al., 2012) with the possible exception of NGC 6791 (Geisler et al. 2012, but see Bragaglia et al. 2014).

## 7.3 Self-enrichment scenarios

Different self-enrichment scenarios were proposed to explain the abundance patterns described above.

### 7.3.1 AGB model

In the 'AGB model' (e.g. Ventura et al., 2001; D'Ercole et al., 2012; Ventura et al., 2013; D'Antona et al., 2016) it is assumed that a second generation of stars forms from material that is polluted by AGB winds from a first generation. The model starts from the point that AGB winds are slow enough so that they may be unable to escape the potential well of a massive star cluster. It is conjectured that after the type II supernova phase, AGB winds would be the only energy source for the intracluster gas, which would be insufficient to overcome radiative losses. Consequently, a cooling flow forms, and stars would form in the centre of the cluster. This scenario has been criticised in multiple ways: AGB nucleosynthesis builds an O-Na correlation instead of the observed anti-correlation (Forestini & Charbonnel, 1997), and it releases He-burning products, thus predicting total C+N+O variations that are only observed in a few GCs (e.g. Decressin et al., 2009; Yong et al., 2015). The AGB stars need to be massive enough to undergo hot-bottom burning to reach the required temperatures (~6.5 $M_{\odot}$, Ventura et al., 2001), which implies that not enough mass is available to pollute (in some cases) >80% of the stars. This is commonly referred to as the mass budget problem (e.g. Prantzos & Charbonnel, 2006; Schaerer & Charbonnel, 2011; Gieles & Charbonnel, 2019). Some of the ideas that have been put forward to overcome the mass budget problem, such as the loss of >90% of the first generation stars over the tidal boundary (D'Ercole et al., 2008), do not address the GC mass dependence of the fraction of polluted stars, i.e. the requirement of more polluted material per unit of GC mass in more massive GCs. We refer to this as the specific mass budget problem. Whether star clusters have a cooling flow phase at all is questionable (§5.4). Also, the gas expulsion paradigm, which had been identified as the mechanism to expel the surplus 1P stars, may not be applicable to many low-mass
clusters (§5.2.1,§5.2.2), and would certainly require some very unusual assumptions for it to work in high-mass clusters (§5.2.3, Bastian & Lardo 2015). Finally, in dwarf galaxies the amount of field stars with the same metallicity is comparable to the total mass in GCs, putting an upper limit of $\lesssim 50\%$ on the amount of mass that GCs could have lost (Larsen et al., 2012). This is referred to as the external mass budget problem. Recent models combining the variety of dynamical mass loss mechanisms discussed in Section 6 predict that globular clusters were only 2–4 times more massive at birth (Kruijssen, 2015; Reina-Campos et al., 2018).

### 7.3.2 Fast-rotating massive star model

Ordinary massive stars ($\sim 20 - 100 M_{\odot}$) have the correct central temperature to create the O-Na anti-correlation, but they are not convective and the material does therefore not reach the surface. To overcome this, Decressin et al. (2007) assume fast-rotating massive stars (FRMS) which undergo rotationally induced mixing and possibly lose a large amount of material through a mechanical wind via an outflowing (decretion) disk. Krause et al. (2013b) laid out a detailed scenario, showing that 2P stars may form in the decretion discs from material spreading through the disc from the surface of the star, and accreting pristine gas during a somewhat extended embedded phase of $\approx 10$ Myr. The embedded phase would be longer in massive star clusters, because stellar feedback would not be strong enough to remove the gas. However, the FRMS scenario also faces the mass budget problem (Prantzos & Charbonnel, 2006; Krause et al., 2013b). In addition, FRMS reach high-enough temperatures to activate the MgAl chain only when the He fraction has increased significantly, predicting a larger He spread in GCs with a Mg-Al anti-correlation than observed (e.g. Chantereau et al., 2016; Nardiello et al., 2015; Milone et al., 2015b).

### 7.3.3 Supermassive star model

Finally, supermassive stars (SMS) have also been suggested as polluters. SMS models with masses between $\sim 2 \times 10^3 M_{\odot}$ and $\sim 2 \times 10^4 M_{\odot}$ reach the required central temperature to activate the MgAl chain ($\sim 72-78$ MK) already at the very beginning of the evolution on the main sequence, when the He abundance is close to pristine (Denissenkov & Hartwick, 2014; Prantzos et al., 2017). Consequently, in this early evolutionary phase the H-burning products of SMSs show remarkable agreement with the various observed abundance anti-correlations and Mg isotopic ratios (Denissenkov & Hartwick, 2014). A formation scenario has been proposed by Gieles et al. (2018). Similar to the picture described in detail in §3, a proto-GC would form at the confluence of gas filaments in gas-rich environments. At the highest density peak, inflow motions are converted to turbulence, and, where gravity dominates locally, the gas fragments into proto-stars. Inflow over at least these two hierarchies leads to cancellation of angular momentum by the time the gas reaches the proto-stars. Accretion of this low-angular momentum gas then reduces the specific angular momentum (i.e. the angular momentum per unit mass) and the stellar density increases as $\rho \propto M^{10}$ (Bonnell et al., 1998). An SMS would then form by runaway collisions. The SMS is assumed to be fully convective and the nucleosynthesis yields are efficiently brought to the surface, streaming off via the usual radiatively driven wind. The wind is initially fast $\gtrsim 1000$ km s$^{-1}$, but
brakes down quickly by interaction with dense gas in the still embedded cluster. The ejecta would then mix into the star-forming dense gas that is accreting on the proto-stars in the cluster or collapse locally to form stars independently.

If SMSs form via stellar collisions, then it may be possible to keep the He abundance low and also produce an order of magnitude more polluted material thereby addressing the mass budget problem (Gieles et al., 2018). This model also provides a pathway to solve the specific mass budget problem because the dynamical models predict a super-linear scaling between the amount of polluted material released via the SMS wind and GC mass. As of today, SMS thus appear to be the most appealing candidate, provided that these stars really exist in nature and are fully convective (cf. Haemmerlé et al., 2018) so they can release the material at the very beginning of their evolution to avoid overproduction of He.

It may come as a surprise that SMS would be difficult to find observationally. Martins et al. (2020) predicted the detectability of cool SMS in proto-GCs at high redshift through deep imaging with JWST NIRCAM camera. One problem at low redshift however is that clusters that would be massive enough are not found in the Milky Way. R136 in the Large Magellanic cloud is, at 50 kpc distance, the closest example of a young massive cluster that may just be massive and compact enough to show some signs of massive star formation via collisions. Indeed some very massive stars ($> 160M_\odot$) have been observed (Crowther et al., 2010). Since SMS are expected to occur in embedded star clusters, absorption would likely be an issue.

An interesting alternative might be MASER emission. GHz MASERs are reliable tracers of massive star formation (Ellingsen et al., 2018; Billington et al., 2019). SMS might hence be expected to show particularly bright MASER emission. Active Galactic Nuclei are an interesting analogue, as they also come with molecular tori of significant optical depth and strong central UV source. Strong MASER emission is frequently seen in heavily-absorbed AGN (Castangia et al., 2019). A good example is the archetypical nearby AGN in NGC 1068, where the outer accretion disc is spatially and kinematically resolved, and the mass of the central object has thus been measured (Murayama & Taniguchi, 1997; Gallimore et al., 2001).

Similar to the case of AGN, forming super star clusters have also been found to be associated with 22 GHz H₂O MASER emission (Gorski et al., 2019). In particular, Gorski et al. (2019) find a nuclear kilomaser in NGC 253, also associated with a forming super star cluster. The spectrum has more than one component and a total width of ≈ 170 km s⁻¹. Whether this relates to SMSs in their centres, and if rotation curves of any SMS disc can be obtained, remains to be seen.

Many other models have been proposed in the literature (de Mink et al., 2009; Bastian et al., 2013; Elmegreen, 2017; Kim & Lee, 2018; Szécsi & Wünsch, 2019). A recent example is the model by Zinnecker et al. (2020, submitted). The model suggests that first, only convective, still accreting high-mass stars form (Hosokawa & Omukai, 2009), with slow, heavily cooling winds (Vink, 2018) producing a chemically anomalous population of predominantly low-mass stars. In this model, an SMS is not needed, yet it also solves the mass budget problem and the He overproduction problem. More research into such models is required.
8 Synthesis: physical processes in complex systems of stars and gas

It is well established that a complex cycle of matter and energy takes place within galaxies. The non-linear flow patterns in the multi-scale multi-phase interstellar medium are the central engine of galaxy evolution, they determine where and at what rate stellar clusters and associations form in our Milky Way. They build up in regions of the interstellar medium that become unstable under their own weight. These are the star forming clumps and cores in the deep interior of molecular clouds.

Altogether, star cluster formation can be seen as a three-phase process. First, supersonic turbulence creates a highly transient and inhomogeneous molecular cloud structure that is characterised by large density contrasts. Some of the high-density fluctuations are transient. Others exceed the critical mass for gravitational contraction, they begin to collapse and eventually decouple from the complex cloud environment. Second, the collapse of these unstable cores leads to the formation of individual stars in clusters and associations. In this phase, a nascent protostar grows in mass via accretion from the infalling envelope until the available gas reservoir is exhausted or stellar feedback effects become important and remove the parental cocoon. Third, stellar feedback becomes so efficient that all the remaining gas is cleared. This reveals the young cluster to optical telescopes, and its subsequent secular evolution is then largely dominated by gravitational processes rather than by the complex competition between gravity and many other physical agents.

We begin our discussion with a definition of what we actually mean when talking about star clusters in Section 2. Then we present evidence from analytic studies and numerical models that indicate that the proto-cluster gas is heavily influenced by the initial conditions and the dynamical properties of the parental cloud in Section 3.

As gas contracts to form stars, the density increases by more than 25 orders of magnitude and the temperature rises by a factor of a million. The process comes to an end when nuclear burning processes set in and provide stability: a star is born. We discuss the different models and suggestions to describe this evolutionary phase in Section 4.

Star formation is controlled by the intricate interplay between the self-gravity of the ISM and various opposing agents, such as supersonic turbulence, magnetic fields, radiation and gas pressure. The evolution is modified by the thermodynamic response of the gas, which is determined by the balance between heating and cooling, which in turn depends on the chemical composition of the material and the detailed interaction of gas and dust with the interstellar radiation field. Altogether, stellar feedback provides enough energy and momentum to remove the parental gas from the cluster. It may also be responsible for clusters being in a global outflow state for the rest of their life. The various physical agents contributing to this process are discussed in Section 5.

Once gas is removed, the subsequent dynamical evolution of a star cluster becomes relatively simple. It is solely governed by the mutual gravitational attraction of the stars in the cluster, modified only by tidal forces exerted from the larger-scale galactic environment, which are weak except near the galactic center or when clusters pass nearby overdensities (such as giant molecular clouds, GMCs, spiral arms, or other clusters), and by mass loss due to the internal evolution of
the constituent stars. Large self-gravitating systems such as star clusters exhibit complex dynamical behavior which we discuss in Section 6.

The chemical composition of stars can provide important constraints on the origin of stellar clusters and help us to distinguish between different physical scenarios. We therefore introduce in Section 7 key aspects of stellar nucleosynthesis and discuss their relation to cluster formation and evolution. Specifically, we speculate about the physical reasons for the observed O-Na anti-correlation observed in globular clusters.

These different perspectives emphasise the interdependence of the different processes: How long gas remains in a state of turbulence before accreting onto a star (§3), and how accretion discs are connected to the upper hierarchies of the gas structure (§4) is crucial to understand how, and what kind of massive stars can pollute the gas out of which the low-mass stars form in massive star clusters and why this is not happening in lower-mass clusters (§7). Stellar feedback determines the abundance of gas in the cluster at all times after formation (§5) with implications for star formation (§4,§7) and the dynamics of the stars (§6).

In summary, star clusters originate from a large reservoir of diffuse gas and dust that permeates the Galaxy, the interstellar medium. The process is governed by the complex interplay of often competing physical agents such as gravity, turbulence, magnetic fields, and radiation across the entire electromagnetic wavelength range. The system is organised in a hierarchy of scales, that link the global dynamical evolution of the galactic gas, to dense, star-forming clouds, and eventually to the newly born star clusters in their interior. Stellar feedback creates highly non-linear feedback loops that strongly influence the dynamical evolution across the entire cascade of scales.

We provide an overview of the most important physical agents involved in the formation and early evolution of star clusters. We argue that stellar birth is a highly dynamical event, and that it couples a wide range of scales in space, time, and energy across the overall hierarchical structure of the galaxy. Star clusters are key constituents of the complex galactic ecosystem, where large parts evolve far from equilibrium and which exhibits highly non-linear dynamical behavior. Progress in this field rests on a comprehensive understanding of the underlying physics and chemistry. Due to the stochastic nature of problem, any theory of star formation is necessarily based on a statistical approach combined with an inventory of the different Galactic environments and knowledge of their possible variations across all scales. This is the direction in which future research efforts in this field will go.

**Acknowledgements** This review emerged from a workshop at the International Space Science Institute in Bern, Switzerland in May 2019. The authors thank the staff of ISSI for their generous hospitality and creating the inspiring atmosphere that initiated this project. We thank the referee for a constructive report which has improved the manuscript.

## Conflict of interest

J.M.D.K. gratefully acknowledges funding from the German Research Foundation (DFG) in the form of an Emmy Noether Research Group (grant number KR4801/1-1) and a DFG Sachbeihilfe Grant (grant number KR4801/2-1), from the European Research Council (ERC) under the European Union's Horizon 2020
research and innovation programme via the ERC Starting Grant MUSTANG (grant agreement number 714907), and from Sonderforschungsbereich SFB 881 "The Milky Way System" (subproject B2) of the DFG. S.S.R.O. acknowledges funding from NSF Career grant AST-1650486. P.G. acknowledges funding from the European Research Council under ERC-CoG grant CRAGSMAN-646955.

R.S.K. acknowledges support from the Deutsche Forschungsgemeinschaft via the SFB 881 The Milky Way System (subprojects B1, B2, and B8) as well as funding from the Heidelberg Cluster of Excellence STRUCTURES in the framework of Germany's Excellence Strategy (grant EXC-2181/1 - 390900948).

JBP acknowledges UNAM-DGAPA-PAPIIT support through grant number IN-111-219.

MG acknowledges support from the ERC (CLUSTERS, StG-335936).

E.V.-S. acknowledges financial support from CONACYT grant 255295.

## References

* Añorve-Zeferino G. A., Tenorio-Tagle G., Silich S., 2009, MNRAS, 394, 1284

* Adams F. C., Lin D. N. C., 1993, in Levy E. H., Lunine J. I., eds, Protostars and Planets III. p. 721

* Adams F. C., Shu F. H., 2007, ApJ, 671, 497

* Allen T. S., et al., 2012, ApJ, 750, 125

* Allison R. J., Goodwin S. P., Parker R. J., Portegies Zwart S. F., de Grijs R., Kouwenhoven M. B. N., 2009a, MNRAS, 395, 1449

* Allison R. J., Goodwin S. P., Parker R. J., de Grijs R., Portegies Zwart S. F., Kouwenhoven M. B. N., 2009b, ApJ, 700, L99

* Andrews S. M., et al., 2018, ApJ, 869, L41

* Antonini F., Gieles M., 2020, arXiv:1906.11855

* Arce H. G., Mardones D., Corder S. A., Garay G., Noriega-Crespo A., Raga A. C., 2013, ApJ, 774, 39

* Armitage P. J., Livio M., Pringle J. E., 2001, MNRAS, 324, 705

* Arzoumanian D., et al., 2011, A&A, 529, L6

* Audard M., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 387 (arXiv:1401.3368), doi:$10.2458/azu'uapress'9780816531240-ch017$

* Audit E., Hennebelle P., 2005, A&A, 433, 1

* Audit E., Hennebelle P., 2010, A&A, 511, A76

* Bachiller R., Martin-Pintado J., Fuente A., 1991, A&A, 243, L21

* Balbus S. A., Hawley J. F., 1994, MNRAS, 266, 769

* Ballesteros-Paredes J., Hartmann L., Vázquez-Semadeni E., 1999, ApJ, 527, 285

* Ballesteros-Paredes J., Hartmann L. W., Vázquez-Semadeni E., Heitsch F., Zamora-Avilés M. A., 2011, MNRAS, 411, 65

* Ballesteros-Paredes J., Hartmann L. W., Pérez-Goytia N., Kuznetsova A., 2015, MNRAS, 452, 566

* Bally J., 2016, Annual Review of Astronomy and Astrophysics, 54, 491

* Banerjee S., Kroupa P., 2017, A&A, 597, A28

* Banerjee R., Vázquez-Semadeni E., Hennebelle P., Klessen R. S., 2009, MNRAS, 398, 1082

* Baraffe I., Chabrier G., Gallardo J., 2009, ApJ, 702, L27
Bastian N., 2016, in EAS Publications Series. pp 5–37 (arXiv:1606.09468), doi:10.1051/eas/1680002

Bastian N., Goodwin S. P., 2006, MNRAS, 369, L9

Bastian N., Lardo C., 2015, MNRAS, 453, 357

Bastian N., Lardo C., 2018, ARA&A, 56, 83

Bastian N., Ercolano B., Gieles M., Rosolowsky E., Scheepmaker R. A., Gutermuth R., Efremov Y., 2007, MNRAS, 379, 1302

Bastian N., Lamers H. J. G. L. M., de Mink S. E., Longmore S. N., Goodwin S. P., Gieles M., 2013, MNRAS, 436, 2398

Basu S., Jones C. E., 2004, MNRAS, 347, L47

Bate M. R., 2009, MNRAS, 392, 590

Baumgardt H., Kroupa P., 2007, MNRAS, 380, 1589

Baumgardt H., Hut P., Heggie D. C., 2002, MNRAS, 336, 1069

Bedin L. R., Piotto G., Anderson J., Cassisi S., King I. R., Momany Y., Carraro G., 2004, ApJ, 605, L125

Bellini A., Milone A. P., Anderson J., Marino A. F., Piotto G., van der Marel R. P., Bedin L. R., King I. R., 2017, ApJ, 844, 164

Belloni D., Giersz M., Rivera Sandoval L. E., Askar A., Ciecielåg P., 2019, MNRAS, 483, 315

Bialy S., Sternberg A., 2019, ApJ, 881, 160

Billington S. J., et al., 2019, MNRAS, 490, 2779

Blaauw A., 1961, Bulletin of the Astronomical Institutes of the Netherlands, 15, 265

Blaauw A., 1964, ARA&A, 2, 213

Blaes O. M., Balbus S. A., 1994, ApJ, 421, 163

Blitz L., Williams J. P., 1999, in Lada C. J., Kylafis N. D., eds, NATO Advanced Science Institutes (ASI) Series C Vol. 540, NATO Advanced Science Institutes (ASI) Series C. p. 3

Bonnell I. A., 1994, MNRAS, 269, 837

Bonnell I. A., Bate M. R., 2006, MNRAS, 370, 488

Bonnell I. A., Bate M. R., Zinnecker H., 1998, MNRAS, 298, 93

Bonnell I. A., Bate M. R., Clarke C. J., Pringle J. E., 2001, MNRAS, 323, 785

Bontemps S., Motte F., Csengeri T., Schneider N., 2010, A&A, 524, A18+

Bovino S., Schleicher D. R. G., Schober J., 2013, New Journal of Physics, 15, 013055

Bragaglia A., Gratton R. G., Carretta E., D'Orazi V., Sneden C., Lucatello S., 2012, A&A, 548, A122

Bragaglia A., Sneden C., Carretta E., Gratton R. G., Lucatello S., Bernath P. F., Brooke J. S. A., Ram R. S., 2014, ApJ, 796, 68

Bravi L., et al., 2018, A&A, 615, A37

Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779

Bressert E., et al., 2010, MNRAS, 409, L54

Brown J. A., Wallerstein G., 1992, AJ, 104, 1818

Butler D., Dickens R. J., Epps E., 1978, ApJ, 225, 148

Cabrera-Ziri I., 2019, in A Synoptic View of the Magellanic Clouds: VMC, Gaia and Beyond. p. 14, doi:10.5281/zenodo.3472465

Cabrera-Ziri I., et al., 2015, MNRAS, 448, 2224

Caldwell S., Chang P., 2018, MNRAS, 474, 4818

Cantó J., Raga A. C., Rodríguez L. F., 2000, ApJ, 536, 896
Carretta E., 2015, ApJ, 810, 148

Carretta E., et al., 2009a, A&A, 505, 117

Carretta E., Bragaglia A., Gratton R., Lucatello S., 2009b, A&A, 505, 139

Carretta E., Bragaglia A., Gratton R. G., Recio-Blanco A., Lucatello S., D'Orazi V., Cassisi S., 2010, A&A, 516, A55

Castangia P., Surcis G., Tarchi A., Caccianiga A., Severgnini P., Della Ceca R., 2019, A&A, 629, A25

Cavallo R. M., Sweigart A. V., Bell R. A., 1996, ApJ, 464, L79

Ceverino D., Klypin A., 2009, ApJ, 695, 292

Chantereau W., Charbonnel C., Meynet G., 2016, A&A, 592, A111

Charbonnel C., 2016, in EAS Publications Series. pp 177–226 (arXiv:1611.08855), doi:10.1051/eas/1680006

Chen C.-Y., et al., 2019a, MNRAS, p. 2252

Chen H.-R. V., et al., 2019b, ApJ, 875, 24

Chen H. H.-H., et al., 2019c, ApJ, 886, 119

Cheng Z., Li Z., Xu X., Li X., 2018, ApJ, 858, 33

Chevalier R. A., Clegg A. W., 1985, Nature, 317, 44

Chevance M., et al., 2019, MNRAS in press, p. arXiv:1911.03479

Choksi N., Kruijssen J. M. D., 2019, MNRAS submitted, p. arXiv:1912.05560

Clark J. S., Lohr M. E., Patrick L. R., Najarro F., 2019, A&A, 623, A84

Clarkson W. I., Ghez A. M., Morris M. R., Lu J. R., Stolte A., McCrady N., Do T., Yelda S., 2012, ApJ, 751, 132

Cohn H., 1980, ApJ, 242, 765

Colín P., Vázquez-Semadeni E., Gómez G. C., 2013, MNRAS, 435, 1701

Conroy C., Spergel D. N., 2011, ApJ, 726, 36

Cottaar M., Meyer M. R., Andersen M., Espinoza P., 2012, A&A, 539, A5

Crocker R. M., Krumholz M. R., Thompson T. A., Baumgardt H., Mackey D., 2018, MNRAS, 481, 4895

Crowther P. A., Schnurr O., Hirschi R., Yusof N., Parker R. J., Goodwin S. P., Kassim H. A., 2010, MNRAS, 408, 731

Crutcher R. M., 2012, ARA&A, 50, 29

D'Antona F., Vesperini E., D'Ercole A., Ventura P., Milone A. P., Marino A. F., Tailo M., 2016, MNRAS, 458, 2122

D'Ercole A., Vesperini E., D'Antona F., McMillan S. L. W., Recchi S., 2008, MNRAS, 391, 825

D'Ercole A., D'Antona F., Ventura P., Vesperini E., McMillan S. L. W., 2010, MNRAS, 407, 854

D'Ercole A., D'Antona F., Carini R., Vesperini E., Ventura P., 2012, MNRAS, 423, 1521

D'Ercole A., D'Antona F., Vesperini E., 2016, MNRAS, 461, 4088

Da Costa G. S., 1997, in Bedding T. R., Booth A. J., Davis J., eds, IAU Symposium Vol. 189, IAU Symposium. pp 193–202

Da Rio N., Robberto M., Soderblom D. R., Panagia N., Hillenbrand L. A., Palla F., Stassun K. G., 2010, ApJ, 722, 1092

Dale J. E., 2017, MNRAS, 467, 1067

Dale J. E., Ercolano B., Bonnell I. A., 2015, MNRAS, 451, 987

Dalessandro E., Lapenna E., Mucciarelli A., Origlia L., Ferraro F. R., Lanzoni B., 2016, ApJ, 829, 77

Decressin T., Meynet G., Charbonnel C., Prantzos N., Ekström S., 2007, A&A,
464, 1029

Decressin T., Charbonnel C., Siess L., Palacios A., Meynet G., Georgy C., 2009, A&A, 505, 727

Denisenkov P. A., Denisenkova S. N., 1990, Soviet Astronomy Letters, 16, 275

Denissenkov P. A., Hartwick F. D. A., 2014, MNRAS, 437, L21

Denissenkov P. A., Weiss A., 1996, A&A, 308, 773

Dib S., Schmeja S., Parker R. J., 2018, MNRAS, 473, 849

Diehl R., et al., 2014, Science, 345, 1162

Dunham M. M., Vorobyov E. I., 2012, ApJ, 747, 52

Dunham M. M., Evans Neal J. I., Terebey S., Dullemond C. P., Young C. H., 2010, ApJ, 710, 470

Dunham M. M., et al., 2014a, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 195 (arXiv:1401.1809), doi:$10.2458/azu'uapress'9780816531240-ch009$

Dunham M. M., Arce H. G., Mardones D., Lee J.-E., Matthews B. C., Stutz A. M., Williams J. P., 2014b, ApJ, 783, 29

Ellingsen S. P., Voronkov M. A., Breen S. L., Caswell J. L., Sobolev A. M., 2018, MNRAS, 480, 4851

Elmegreen B. G., 2010, ApJ, 712, L184

Elmegreen B. G., 2017, ApJ, 836, 80

Elmegreen B. G., Falgarone E., 1996, ApJ, 471, 816

Elmegreen B. G., Hunter D. A., 2010, ApJ, 712, 604

Espinoza P., Selman F. J., Melnick J., 2009, A&A, 501, 563

Evans Neal J. I., et al., 2009, ApJS, 181, 321

Fall S. M., Zhang Q., 2001, ApJ, 561, 751

Farias J. P., Smith R., Fellhauer M., Goodwin S., Candlish G. N., Blaña M., Dominguez R., 2015, MNRAS, 450, 2451

Farias J. P., Fellhauer M., Smith R., Domínguez R., Dabringhausen J., 2018, MNRAS, 476, 5341

Fatuzzo M., Adams F. C., Myers P. C., 2004, ApJ, 615, 813

Federrath C., 2015, MNRAS, 450, 4035

Federrath C., Schober J., Bovino S., Schleicher D. R. G., 2014, ApJ, 797, L19

Feigelson E. D., et al., 2013, ApJS, 209, 26

Field G. B., 1965, ApJ, 142, 531

Field G. B., Blackman E. G., Keto E. R., 2008, MNRAS, 385, 181

Figer D. F., Kim S. S., Morris M., Serabyn E., Rich R. M., McLean I. S., 1999, ApJ, 525, 750

Figer D. F., et al., 2002, ApJ, 581, 258

Fischer W. J., et al., 2017, ApJ, 840, 69

Forestini M., Charbonnel C., 1997, A&A Supplement series, 123

Foster P. N., Chevalier R. A., 1993, ApJ, 416, 303

Foster J. B., et al., 2015, ApJ, 799, 136

Fukui Y., et al., 2009, ApJ, 705, 144

Gaia Collaboration et al., 2018, A&A, 616, A1

Gaibler V., Camenzind M., Krause M., 2005, Proceedings of the 1st Arizona/Heidelberg Symposium "The High Redshift Frontier"

Gallimore J. F., Henkel C., Baum S. A., Glass I. S., Claussen M. J., Prieto M. A., Von Kap-herr A., 2001, ApJ, 556, 694

Gammie C. F., 1996, ApJ, 457, 355
Gavagnin E., Bleuler A., Rosdahl J., Teyssier R., 2017, MNRAS, 472, 4155

Geen S., Watson S. K., Rosdahl J., Bieri R., Klessen R. S., Hennebelle P., 2018, MNRAS, 481, 2548

Geisler D., Villanova S., Carraro G., Pilachowski C., Cummings J., Johnson C. I., Bresolin F., 2012, ApJ, 756, L40

Gennaro M., Brandner W., Stolte A., Henning T., 2011, MNRAS, 412, 2469

Gennaro M., Goodwin S. P., Parker R. J., Allison R. J., Brandner W., 2017, MNRAS, 472, 1760

Gentry E. S., Krumholz M. R., Madau P., Lupi A., 2019, MNRAS, 483, 3647

Getman K. V., Broos P. S., Kuhn M. A., Feigelson E. D., Richert A. J. W., Ota Y., Bate M. R., Garmire G. P., 2017, ApJS, 229, 28

Getman K. V., Feigelson E. D., Kuhn M. A., Bate M. R., Broos P. S., Garmire G. P., 2018, MNRAS, 476, 1213

Gieles M., Charbonnel C., 2019, arXiv e-prints, p. arXiv:1908.02075

Gieles M., Portegies Zwart S. F., 2011, MNRAS, 410, L6

Gieles M., Renaud F., 2016, MNRAS, 463, L103

Gieles M., Portegies Zwart S. F., Baumgardt H., Athanassoula E., Lamers H. J. G. L. M., Sipior M., Leenaarts J., 2006, MNRAS, 371, 793

Gieles M., Sana H., Portegies Zwart S. F., 2010a, MNRAS, 402, 1750

Gieles M., Baumgardt H., Heggie D. C., Lamers H. J. G. L. M., 2010b, MNRAS, 408, L16

Gieles M., Heggie D. C., Zhao H., 2011, MNRAS, 413, 2509

Gieles M., Moeckel N., Clarke C. J., 2012, MNRAS, 426, L11

Gieles M., et al., 2018, MNRAS, 478, 2461

Giersz M., Askar A., Wang L., Hypki A., Leveque A., Spurzem R., 2019, MNRAS, 487, 2412

Gilligan C. K., et al., 2019, MNRAS, 486, 5581

Ginsburg A., et al., 2016, A&A, 595, A27

Girichidis P., Federrath C., Banerjee R., Klessen R. S., 2012a, MNRAS, 420, 613

Girichidis P., Federrath C., Allison R., Banerjee R., Klessen R. S., 2012b, MNRAS, 420, 3264

Glassgold A. E., Najita J. R., Igea J., 2007, ApJ, 656, 515

Glover S. C. O., Clark P. C., 2012, MNRAS, 421, 9

Gnedin O. Y., Ostriker J. P., 1999, ApJ, 513, 626

Goldbaum N. J., Krumholz M. R., Matzner C. D., McKee C. F., 2011, ApJ, 738, 101

Goldreich P., Kwan J., 1974, ApJ, 189, 441

Gómez G. C., Vázquez-Semadeni E., 2014, ApJ, 791, 124

Goodman J., 1987, ApJ, 313, 576

Goodwin S. P., Bastian N., 2006, MNRAS, 373, 752

Gorski M. D., Ott J., Rand R., Meier D. S., Momjian E., Schinnerer E., Ellingsen S. P., 2019, MNRAS, 483, 5434

Gouliermis D. A., 2018, PASP, 130, 072001

Gower J. C., Ross G. J. S., 1969, Appl. Stat., 18, 54

Gratton R. G., et al., 2001, A&A, 369, 87

Gratton R. G., Carretta E., Bragaglia A., 2012, A&A Review, 20, 50

Gratton R., Bragaglia A., Carretta E., D'Orazi V., Lucatello S., Sollima A., 2019, A&A Review, 27, 8

Guttermuth R. A., Megeath S. T., Myers P. C., Allen L. E., Pipher J. L., Fazio
G. G., 2009, ApJS, 184, 18

Guttermuth R. A., Pipher J. L., Megeath S. T., Myers P. C., Allen L. E., Allen
T. S., 2011, ApJ, 739, 84

Hacar A., Alves J., Tafalla M., Goicoechea J. R., 2017, A&A, 602, L2

Haemmerlé L., Woods T. E., Klessen R. S., Heger A., Whalen D. J., 2018, MNRAS, 474, 2757

Han S.-I., Lee Y.-W., Joo S.-J., Sohn S. T., Yoon S.-J., Kim H.-S., Lee J.-W., 2009, ApJ, 707, L190

Hansen C. E., Klein R. I., McKee C. F., Fisher R. T., 2012, The Astrophysical Journal, 747, 22

Harayama Y., Eisenhauer F., Martins F., 2008, ApJ, 675, 1319

Hartmann L., Ballesteros-Paredes J., Bergin E. A., 2001, ApJ, 562, 852

Hartmann L., Ballesteros-Paredes J., Heitsch F., 2012, MNRAS, 420, 1457

Hartmann L., Herczeg G., Calvet N., 2016, ARA&A, 54, 135

Heiner J. S., Vázquez-Semadeni E., Ballesteros-Paredes J., 2015, MNRAS, 452, 1353

Heitsch F., Hartmann L., 2008, ApJ, 689, 290

Heitsch F., Burkert A., Hartmann L. W., Slyz A. D., Devriendt J. E. G., 2005, ApJ, 633, L113

Heitsch F., Slyz A. D., Devriendt J. E. G., Hartmann L. W., Burkert A., 2006, ApJ, 648, 1052

Heitsch F., Hartmann L. W., Slyz A. D., Devriendt J. E. G., Burkert A., 2008, ApJ, 674, 316

Heitsch F., Ballesteros-Paredes J., Hartmann L., 2009, ApJ, 704, 1735

Hénault-Brunet V., et al., 2012, A&A, 546, A73

Hennebelle P., Falgarone E., 2012, A&A Review, 20, 55

Hennebelle P., Banerjee R., Vázquez-Semadeni E., Klessen R. S., Audit E., 2008, A&A, 486, L43

Hénon M., 1961, Annales d'Astrophysique, 24, 369; English translation: arXiv:1103.3499

Hénon M., 1965, Ann. Astrophys, 28, 62; English translation: ArXiv:1103.3498

Hénon M., 1975, in A. Hayli ed., Proc. IAU Symp. 69, Dynamics of the Solar Systems. Reidel, Dordrecht, p. 133.

Heringer E., Pritchet C., van Kerkwijk M. H., 2019, ApJ, 882, 52

Heyer M., Krawczyk C., Duval J., Jackson J. M., 2009, ApJ, 699, 1092

Hillenbrand L. A., Findeisen K. P., 2015, ApJ, 808, 68

Hillenbrand L. A., Hartmann L. W., 1998, ApJ, 492, 540

Hills J. G., 1980, ApJ, 235, 986

Hirano N., Ho P. P. T., Liu S.-Y., Shang H., Lee C.-F., Bourke T. L., 2010, ApJ, 717, 58

Hollyhead K., Bastian N., Adamo A., Silva-Villa E., Dale J., Ryon J. E., Gazak
Z., 2015, MNRAS, 449, 1106

Hosokawa T., Omukai K., 2009, ApJ, 691, 823

Hosokawa T., Offner S. S. R., Krumholz M. R., 2011, ApJ, 738, 140

Hoyle F., 1953, ApJ, 118, 513

Huff E. M., Stahler S. W., 2006, ApJ, 644, 355

Hunter C., 1977, ApJ, 218, 834

Ivans I. I., Sneden C., Kraft R. P., Suntzeff N. B., Smith V. V., Langer G. E., Fulbright J. P., 1999, AJ, 118, 1273
Ji L., Wang Q. D., Kwan J., 2006, MNRAS, 372, 497

Jordán A., et al., 2004, ApJ, 613, 279

Karnath N., Prchlik J. J., Gutermuth R. A., Allen T. S., Megeath S. T., Pipher
J. L., Wolk S., Jeffries R. D., 2019, ApJ, 871, 46

Kavanagh P. J., Norci L., Meurs E. J. A., 2011, New Astronomy, 16, 461

Kenyon S. J., Hartmann L., 1995, ApJS, 101, 117

Kenyon S. J., Hartmann L. W., Strom K. M., Strom S. E., 1990, AJ, 99, 869

Keto E., Caselli P., Rawlings J., 2015, MNRAS, 446, 3731

Kim J. J., Lee Y.-W., 2018, ApJ, 869, 35

Kim S. S., Figer D. F., Kudritzki R. P., Najarro F., 2006, ApJ, 653, L113

Kim J.-G., Kim W.-T., Ostriker E. C., 2016, ApJ, 819, 137

King I. R., et al., 2012, AJ, 144, 5

Kirk H., Myers P. C., 2011, ApJ, 727, 64

Kirk H., Pineda J. E., Johnstone D., Goodman A., 2010, ApJ, 723, 457

Kirk H., Myers P. C., Bourke T. L., Gutermuth R. A., Hedden A., Wilson G. W., 2013, ApJ, 766, 115

Kirk H., Offner S. S. R., Redmond K. J., 2014, MNRAS, 439, 1765

Kirk H., et al., 2017, ApJ, 846, 144

Klessen R., 2019, Formation of the first stars. pp 67-97,
doi:10.1142/9789813227958'0004

Klessen R. S., Burkert A., 2000, ApJS, 128, 287

Klessen R. S., Burkert A., 2001, ApJ, 549, 386

Klessen R. S., Glover S. C. O., 2016, Saas-Fee Advanced Course, 43, 85

Klessen R. S., Hennebelle P., 2010a, A&A, 520, A17

Klessen R. S., Hennebelle P., 2010b, A&A, 520, A17

Klessen R. S., Heitsch F., Mac Low M.-M., 2000, ApJ, 535, 887

Könyves V., et al., 2015, A&A, 584, A91

Kraft R. P., 1979, ARA&A, 17, 309

Kratter K., Lodato G., 2016, Annual Review of Astronomy and Astrophysics, 54,
271

Kratter K. M., Matzner C. D., Krumholz M. R., Klein R. I., 2010, ApJ, 708, 1585

Krause M., Charbonnel C., Decressin T., Meynet G., Prantzos N., Diehl R., 2012, A&A, 546, L5

Krause M., Fierlinger K., Diehl R., Burkert A., Voss R., Ziegler U., 2013a, A&A,
550, A49

Krause M., Charbonnel C., Decressin T., Meynet G., Prantzos N., 2013b, A&A,
552, A121

Krause M., Diehl R., Böhringer H., Freyberg M., Lubos D., 2014, A&A, 566, A94

Krause M. G. H., et al., 2015, A&A, 578, A113

Krause M. G. H., Charbonnel C., Bastian N., Diehl R., 2016, A&A, 587, A53

Krause M. G. H., et al., 2018, A&A, 619, A120

Krause M. G. H., Hardcastle M. J., Shabala S. S., 2019, A&A, 627, A113

Krauss L. M., Chaboyer B., 2003, Science, 299, 65

Kremer K., et al., 2019, arXiv:1911.00018,

Kretschmer K., Diehl R., Krause M., Burkert A., Fierlinger K., Gerhard O., Greiner J., Wang W., 2013, A&A, 559, A99

Kroupa P., Weidner C., Pflamm-Altenburg J., Thies I., Dabringhausen J., Marks M., Maschberger T., 2013, The Stellar and Sub-Stellar Initial Mass Function of Simple and Composite Populations. p. 115, doi:$10.1007/978-94-007-5612-0'4$
Kruijssen J. M. D., 2012, MNRAS, 426, 3008

Kruijssen J. M. D., 2014, Classical and Quantum Gravity, 31, 244006

Kruijssen J. M. D., 2015, MNRAS, 454, 1658

Kruijssen J. M. D., Pelupessy F. I., Lamers H. J. G. L. M., Portegies Zwart S. F., Icke V., 2011, MNRAS, 414, 1339

Kruijssen J. M. D., Maschberger T., Moeckel N., Clarke C. J., Bastian N., Bonnell I. A., 2012, MNRAS, 419, 841

Kruijssen J. M. D., et al., 2019, Nature, 569, 519

Krumholz M. R., 2014, Physics Reports, 539, 49

Krumholz M. R., McKee C. F., 2019, arXiv:1909.01565

Krumholz M. R., Matzner C. D., McKee C. F., 2006, ApJ, 653, 361

Krumholz M. R., Klein R. I., McKee C. F., 2012, ApJ, 754, 71

Krumholz M. R., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 243 (arXiv:1401.2473), doi:10.2458/azu'uapress'9780816531240-ch011

Krumholz M. R., McKee C. F., Bland-Hawthorn J., 2019, ARA&A, 57, 227

Kuhn M. A., Hillenbrand L. A., Sills A., Feigelson E. D., Getman K. V., 2019, ApJ, 870, 32

Kuznetsova A., Hartmann L., Burkert A., 2017, ApJ, 836, 190

Kuznetsova A., Hartmann L., Heitsch F., Ballesteros-Paredes J., 2018, ApJ, 868, 50

Lada C. J., Lada E. A., 2003, ARA&A, 41, 57

Lagioia E. P., Milone A. P., Marino A. F., Dotter A., 2019, ApJ, 871, 140

Lamers H. J. G. L. M., Gieles M., 2006, A&A, 455, L17

Langer G. E., Hoffman R., Sneden C., 1993, PASP, 105, 301

Larsen S. S., Strader J., Brodie J. P., 2012, A&A, 544, L14

Larsen S. S., Brodie J. P., Grundahl F., Strader J., 2014, ApJ, 797, 15

Larson R. B., 1969, MNRAS, 145, 271

Larson R. B., 1981, MNRAS, 194, 809

Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J., 2013, ApJ, 772, L3

Laughlin G., Bodenheimer P., 1994, ApJ, 436, 335

Lee C.-F., Hirano N., Palau A., Ho P. T. P., Bourke T. L., Zhang Q., Shang H., 2009, The Astrophysical Journal, 699, 1584

Lee E. J., Chang P., Murray N., 2015, ApJ, 800, 49

Lee C.-F., Ho P. T. P., Li Z.-Y., Hirano N., Zhang Q., Shang H., 2017a, Nature Astronomy, 1, 0152

Lee J. W. Y., Hull C. L. H., Offner S. S. R., 2017b, ApJ, 834, 201

Leigh N. W. C., Giersz M., Marks M., Webb J. J., Hypki A., Heinke C. O., Kroupa P., Sills A., 2015, MNRAS, 446, 226

Leitherer C., et al., 1999, ApJS, 123, 3

Leonard P. J. T., Duncan M. J., 1990, AJ, 99, 608

Li G.-X., 2018, MNRAS, 477, 4951

Li P. S., Klein R. I., McKee C. F., 2018, MNRAS, 473, 4220

Li H., Vogelsberger M., Marinacci F., Gnedin O. Y., 2019, MNRAS, 487, 364

Lin C. C., Mestel L., Shu F. H., 1965, ApJ, 142, 1431

Lind K., Primas F., Charbonnel C., Grundahl F., Asplund M., 2009, A&A, 503, 545

Liszt H. S., Wilson R. W., Penzias A. A., Jefferts K. B., Wannier P. G., Solomon P. M., 1974, ApJ, 190, 557
Littlefair S. P., Naylor T., Jeffries R. D., Devey C. R., Vine S., 2003, MNRAS, 345, 1205

Liu Z.-W., Stancliffe R. J., 2018, MNRAS, 475, 5257

Liu H. B., Galván-Madrid R., Jiménez-Serra I., Román-Zúñiga C., Zhang Q., Li Z., Chen H.-R., 2015, ApJ, 804, 37

Lohr M. E., Clark J. S., Najarro F., Patrick L. R., Crowther P. A., Evans C. J., 2018, A&A, 617, A66

Longmore S. N., 2015, MNRAS, 448, L62

Longmore S. N., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. The University of Arizona press, pp 291–314 (arXiv:1401.4175)

Lü H.-J., Lan L., Zhang B., Liang E.-W., Kann D. A., Du S.-S., Shen J., 2018, ApJ, 862, 130

Lyman J. D., et al., 2018, MNRAS, 473, 1359

Lynden-Bell D., Eggleton P. P., 1980, MNRAS, 191, 483

Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125

Mac Low M.-M., Klessen R. S., Burkert A., Smith M. D., 1998, Phys. Rev. Lett., 80, 2754

Machida M. N., Hosokawa T., 2013, MNRAS, 431, 1719

Maeder A., Meynet G., 2012, Reviews of Modern Physics, 84, 25

Marino A. F., Villanova S., Piotto G., Milone A. P., Momany Y., Bedin L. R., Medling A. M., 2008, A&A, 490, 625

Marino A. F., et al., 2011, ApJ, 731, 64

Marino A. F., et al., 2018, ApJ, 859, 81

Marino A. F., et al., 2019, MNRAS, 487, 3815

Martins F., Schaerer D., Haemmerlé L., Charbonnel C., 2020, A&A, 633, A9

Martocchia S., 2019, in A Synoptic View of the Magellanic Clouds: VMC, Gaia and Beyond. p. 11 (arXiv:1907.12440), doi:10.5281/zenodo.3472459

Martocchia S., et al., 2018, MNRAS, 477, 4696

Martocchia S., et al., 2019, MNRAS, 487, 5324

Masseron T., et al., 2019, A&A, 622, A191

Mathews W. G., Brighenti F., 2003, ARA&A, 41, 191

Matzner C. D., Jumper P. H., 2015, ApJ, 815, 68

Mazzali P. A., McFadyen A. I., Woosley S. E., Pian E., Tanaka M., 2014, MNRAS, 443, 67

McKee C. F., Tan J. C., 2003, ApJ, 585, 850

Mészáros S., et al., 2015, AJ, 149, 153

Mészáros S., et al., 2019, MNRAS, p. 3134

Milone A. P., 2015, MNRAS, 446, 1672

Milone A. P., et al., 2015a, MNRAS, 447, 927

Milone A. P., et al., 2015b, ApJ, 808, 51

Milone A. P., et al., 2017, MNRAS, 464, 3636

Milone A. P., et al., 2018, MNRAS, 481, 5098

Miville-Deschênes M.-A., Murray N., Lee E. J., 2017, ApJ, 834, 57

Moeckel N., Clarke C. J., 2011, MNRAS, 410, 2799

Mohammadpour M., Stahler S. W., 2013, MNRAS, 433, 3389

Mozel P., 2003, Journal of the Royal Astronomical Society of Canada, 97, 245

Mucciarelli A., Lapenna E., Ferraro F. R., Lanzoni B., 2018, ApJ, 859, 75

Murayama T., Taniguchi Y., 1997, PASJ, 49, L13
Murray N., Chang P., 2015, ApJ, 804, 44

Myers P. C., 2009a, ApJ, 700, 1609

Myers P. C., 2009b, ApJ, 706, 1341

Myers P. C., 2012, ApJ, 752, 9

Naiman J. P., Soares-Furtado M., Ramirez-Ruiz E., 2020, MNRAS, 491, 4602

Naranjo-Romero R., Vázquez-Semadeni E., Loughnane R. M., 2015, ApJ, 814, 48

Nardiello D., et al., 2015, MNRAS, 451, 312

Nardiello D., Piotto G., Milone A. P., Rich R. M., Cassisi S., Bedin L. R., Bellini
A., Renzini A., 2019, MNRAS, 485, 3076

Niederhofer F., et al., 2017a, MNRAS, 464, 94

Niederhofer F., et al., 2017b, MNRAS, 465, 4159

Norman C., Silk J., 1980, ApJ, 238, 158

Norris J. E., 2004, ApJ, 612, L25

Norris J., Freeman K. C., 1979, ApJ, 230, L179

Norris J., Pilachowski C. A., 1985, ApJ, 299, 295

Norris J., Cottrell P. L., Freeman K. C., Da Costa G. S., 1981, ApJ, 244, 205

O'Malley E. M., Gilligan C., Chaboyer B., 2017, ApJ, 838, 162

Oey M. S., 2011, ApJ, 739, L46

Offner S. S. R., Arce H. G., 2014, The Astrophysical Journal, 784, 61

Offner S. S. R., Chaban J., 2017, The Astrophysical Journal, 847, 104

Offner S. S. R., McKee C. F., 2011, ApJ, 736, 53

Offner S. S. R., Hansen C. E., Krumholz M. R., 2009, ApJ, 704, L124

Offner S. S. R., Robitaille T. P., Hansen C. E., McKee C. F., Klein R. I., 2012, ApJ, 753, 98

Offner S. S. R., Gaches B. A. L., Holdship J. R., 2019, ApJ, 883, 121

Oh S., Kroupa P., Pflamm-Altenburg J., 2015, ApJ, 805, 92

Osborn W., 1971, The Observatory, 91, 223

Padoan P., Nordlund Å., 1999, ApJ, 526, 279

Padoan P., Haugbølle T., Nordlund Å., 2014, ApJ, 797, 32

Padoan P., Pan L., Juvela M., Haugbølle T., Nordlund Å., 2019, arXiv e-prints,
p. arXiv:1911.04465

Palla F., Stahler S. W., 1991, ApJ, 375, 288

Palla F., Stahler S. W., 2000, ApJ, 540, 255

Palla F., Stahler S. W., 2002, ApJ, 581, 1194

Palouš J., Wünsch R., Martínez-González S., Hueyotl-Zahuantitla F., Silich S., Tenorio-Tagle G., 2013, ApJ, 772, 128

Palouš J., Wünsch R., Tenorio-Tagle G., 2014, ApJ, 792, 105

Pancino E., Carrera R., Rossetti E., Gallart C., 2010, A&A, 511, A56

Pancino E., et al., 2017, A&A, 601, A112

Parker R. J., Alves de Oliveira C., 2017, MNRAS, 468, 4340

Parker R. J., Bouvier J., Goodwin S. P., Moraux E., Allison R. J., Guieu S., Güdel M., 2011, MNRAS, 412, 2489

Parker R. J., Maschberger T., Alves de Oliveira C., 2012, MNRAS, 426, 3079

Parker R. J., Dale J. E., Ercolano B., 2015, MNRAS, 446, 4278

Penston M. V., 1969, MNRAS, 144, 425

Peretto N., et al., 2014, A&A, 561, A83

Perez-Becker D., Chiang E., 2011, ApJ, 735, 8

Peterson R. C., 1980, ApJ, 237, L87

Pineda J. E., Goodman A. A., Arce H. G., Caselli P., Foster J. B., Myers P. C.,
Rosolowsky E. W., 2010, ApJ, 712, L116

Piotto G., et al., 2005, ApJ, 621, 777

Piotto G., et al., 2007, ApJ, 661, L53

Piotto G., et al., 2012, ApJ, 760, 39

Piotto G., et al., 2015, AJ, 149, 91

Pleintinger M. M. M., Siegert T., Diehl R., Fujimoto Y., Greiner J., Krause M. G. H., Krumholz M. R., 2019, A&A, 632, A73

Plummer H. C., 1911, MNRAS, 71, 460

Plunkett A. L., Fernández-López M., Arce H. G., Busquet G., Mardones D., Dunham M. M., 2018, A&A, 615, A9

Portegies Zwart S. F., McMillan S. L. W., 2002, ApJ, 576, 899

Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431

Poveda A., Ruiz J., Allen C., 1967, Boletín de los Observatorios Tonantzintla y Tacubaya, 4, 86

Prantzos N., Charbonnel C., 2006, A&A, 458, 135

Prantzos N., Charbonnel C., Iliadis C., 2007, A&A, 470, 179

Prantzos N., Charbonnel C., Iliadis C., 2017, A&A, 608, A28

Rahner D., Pellegrini E. W., Glover S. C. O., Klessen R. S., 2017, MNRAS, 470, 4453

Rahner D., Pellegrini E. W., Glover S. C. O., Klessen R. S., 2019, MNRAS, 483, 2547

Rappenglück M. A., 2001, Earth Moon and Planets, 85, 391

Rebull L. M., et al., 2014, AJ, 148, 92

Reina-Campos M., Kruijssen J. M. D., Pfeffer J., Bastian N., Crain R. A., 2018, MNRAS, 481, 2851

Reissl S., Klessen R. S., Mac Low M.-M., Pellegrini E. W., 2018, A&A, 611, A70

Renzini A., 2013, MmSAI, 84, 162

Robertson B., Goldreich P., 2012, ApJ, 750, L31

Robitaille T. P., Whitney B. A., Indebetouw R., Wood K., Denzmore P., 2006, ApJS, 167, 256

Rodgers-Lee D., Krause M. G. H., Dale J., Diehl R., 2019, MNRAS, 490, 1894

Rodríguez M. J., Baume G., Feinstein C., 2019, A&A, 626, A35

Rosen A. L., Krumholz M. R., McKee C. F., Klein R. I., 2016, MNRAS, 463, 2553

Roupas Z., Kazanas D., 2019, A&A, 621, L1

Sana H., et al., 2012, Science, 337, 444

Sbordone L., Salaris M., Weiss A., Cassisi S., 2011, A&A, 534, A9

Schaerer D., Charbonnel C., 2011, MNRAS, 413, 2297

Schneider N., Csengeri T., Bontemps S., Motte F., Simon R., Hennebelle P., Federrath C., Klessen R., 2010, A&A, 520, A49

Schober J., Schleicher D., Federrath C., Klessen R., Banerjee R., 2012, Phys. Rev. E, 85, 026303

Schober J., Schleicher D. R. G., Federrath C., Bovino S., Klessen R. S., 2015, Phys. Rev. E, 92, 023010

Semenov D., Wiebe D., Henning T., 2004, A&A, 417, 93

Shepherd D. S., Yu K. C., Bally J., Testi L., 2000, ApJ, 535, 833

Shu F. H., 1977, ApJ, 214, 488

Shu F. H., Adams F. C., Lizano S., 1987, ARA&A, 25, 23

Silich S., Tenorio-Tagle G., Rodríguez-González A., 2004, ApJ, 610, 226

Silich S., Tenorio-Tagle G., Añorve-Zeferino G. A., 2005, ApJ, 635, 1116
Silich S., Bisnovatyi-Kogan G., Tenorio-Tagle G., Martínez-González S., 2011, ApJ, 743, 120

Smith R., Goodwin S., Fellhauer M., Assmann P., 2013, MNRAS, 428, 1303

Sokal K. R., Johnson K. E., Indebetouw R., Massey P., 2016, ApJ, 826, 194

Spitzer L. J., 1958, ApJ, 127, 17

Spitzer L. J., Hart M. H., 1971, ApJ, 164, 399

Stahler S. W., Shu F. H., Taam R. E., 1980, ApJ, 241, 637

Stevens I. R., Hartwell J. M., 2003, MNRAS, 339, 280

Stolte A., Brandner W., Grebel E. K., Lenzen R., Lagrange A.-M., 2005, ApJ, 628, L113

Stolte A., Brandner W., Brandl B., Zinnecker H., 2006, AJ, 132, 253

Stone J. M., Ostriker E. C., Gammie C. F., 1998, ApJ, 508, L99

Sugitani K., et al., 2011, ApJ, 734, 63

Sun N.-C., et al., 2018, ApJ, 858, 31

Sur S., Federrath C., Schleicher D. R. G., Banerjee R., Klessen R. S., 2012, MNRAS, 423, 3148

Sweigart A. V., Mengel J. G., 1979, ApJ, 229, 624

Szécsi D., Wünsch R., 2019, ApJ, 871, 20

Tanaka K. E. I., Tan J. C., Zhang Y., 2017, ApJ, 835, 32

Terebey S., Shu F. H., Cassen P., 1984, ApJ, 286, 529

Terlevich E., 1987, MNRAS, 224, 193

Thévenin F., Charbonnel C., de Freitas Pacheco J. A., Idiart T. P., Jasniewicz G., de Laverny P., Plez B., 2001, A&A, 373, 905

Toomre A., 1964, ApJ, 139, 1217

Tychoniec L., et al., 2018, ApJS, 238, 19

Umebayashi T., Nakano T., 1981, PASJ, 33, 617

Vasiliev E. O., Shchekinov Y. A., Nath B. B., 2017, MNRAS, 468, 2757

Vázquez-Semadeni E., Ballesteros-Paredes J., Klessen R. S., 2003, ApJ, 585, L131

Vázquez-Semadeni E., Ryu D., Passot T., González R. F., Gazol A., 2006, ApJ, 643, 245

Vázquez-Semadeni E., Gómez G. C., Jappsen A. K., Ballesteros-Paredes J., González R. F., Klessen R. S., 2007, ApJ, 657, 870

Vázquez-Semadeni E., Gómez G. C., Jappsen A. K., Ballesteros-Paredes J., Klessen R. S., 2009, ApJ, 707, 1023

Vázquez-Semadeni E., Colín P., Gómez G. C., Ballesteros-Paredes J., Watson A. W., 2010, ApJ, 715, 1302

Vázquez-Semadeni E., Banerjee R., Gómez G. C., Hennebelle P., Duffin D., Klessen R. S., 2011, MNRAS, 414, 2511

Vázquez-Semadeni E., González-Samaniego A., Colín P., 2017, MNRAS, 467, 1313

Vázquez-Semadeni E., Palau A., Ballesteros-Paredes J., Gómez G. C., Zamora-Avilés M., 2019, MNRAS, 490, 3061

Ventura P., D'Antona F., Mazzitelli I., Gratton R., 2001, ApJ, 550, L65

Ventura P., Di Criscienzo M., Carini R., D'Antona F., 2013, MNRAS, 431, 3642

Vink J. S., 2018, A&A, 619, A54

Vishniac E. T., 1994, ApJ, 428, 186

Vorobyov E. I., Basu S., 2005, ApJ, 633, L137

Vorobyov E. I., Basu S., 2006, ApJ, 650, 956

Voss R., Diehl R., Hartmann D. H., Cerviño M., Vink J. S., Meynet G., Limongi M., Chieffi A., 2009, A&A, 504, 531
Walker D. L., Longmore S. N., Bastian N., Kruijssen J. M. D., Rathborne J. M., Jackson J. M., Foster J. B., Contreras Y., 2015, MNRAS, 449, 715

Walker D. L., Longmore S. N., Bastian N., Kruijssen J. M. D., Rathborne J. M., Galván-Madrid R., Liu H. B., 2016, MNRAS, 457, 4536

Wang L., 2020, MNRAS, (arXiv:1911.05077), p. 2759

Wang Q. D., Dong H., Lang C., 2006, MNRAS, 371, 38

Wang Y., Primas F., Charbonnel C., Van der Swaelmen M., Bono G., Chantereau W., Zhao G., 2016, A&A, 592, A66

Ward J. L., Kruijssen J. M. D., 2018, MNRAS, 475, 5659

Ward J. L., Kruijssen J. M. D., Rix H.-W., 2019, MNRAS submitted, p. arXiv:1910.06974

Wareing C. J., Falle S. A. E. G., Pittard J. M., 2019, MNRAS, 485, 4686

Washabaugh P. C., Bregman J. N., 2013, ApJ, 762, 1

Webb J. J., Reina-Campos M., Kruijssen J. M. D., 2019, MNRAS, 486, 5879

Weiss A., Denissenkov P. A., Charbonnel C., 2000, A&A, 356, 181

Whitworth A., Summers D., 1985, MNRAS, 214, 1

Wielen R., 1985, in Goodman J., Hut P., eds, IAU Symp. 113: Dynamics of Star Clusters. pp 449–460

Wilson R. W., Jefferts K. B., Penzias A. A., 1970, ApJ, 161, L43

Wright N. J., Mamajek E. E., 2018, MNRAS, 476, 381

Wright N. J., Parker R. J., Goodwin S. P., Drake J. J., 2014, MNRAS, 438, 639

Wright N. J., et al., 2019, MNRAS, 486, 2477

Wünsch R., Silich S., Palouš J., Tenorio-Tagle G., Muñoz-Tuñón C., 2011, ApJ, 740, 75

Wünsch R., Palouš J., Tenorio-Tagle G., Ehlerová S., 2017, ApJ, 835, 60

Wyrowski F., et al., 2016, A&A, 585, A149

Yong D., Grundahl F., Lambert D. L., Nissen P. E., Shetrone M. D., 2003, A&A, 402, 985

Yong D., Grundahl F., Norris J. E., 2015, MNRAS, 446, 3319

Yuan J., et al., 2018, ApJ, 852, 12

Yusef-Zadeh F., Law C., Wardle M., Wang Q. D., Fruscione A., Lang C. C., Cotera A., 2002, ApJ, 570, 665

Zamora-Avilés M., Vázquez-Semadeni E., 2014, ApJ, 793, 84

Zamora-Avilés M., Vázquez-Semadeni E., Colín P., 2012, ApJ, 751, 77

Zamora-Avilés M., Ballesteros-Paredes J., Hernández J., Román-Zúñiga C., Lora V., Kounkel M., 2019, MNRAS, 488, 3406

Zennaro M., Milone A. P., Marino A. F., Cordoni G., Lagioia E. P., Tailo M., 2019, MNRAS, 487, 3239

Zhang D., Thompson T. A., Murray N., Quataert E., 2014, ApJ, 784, 93

Zhang Y., et al., 2016, The Astrophysical Journal, 832, 158

Zhu Z., Hartmann L., Gammie C., 2009a, ApJ, 694, 1045

Zhu Z., Hartmann L., Gammie C., McKinney J. C., 2009b, ApJ, 701, 620

Zinnecker H., Yorke H. W., 2007, ARA&A, 45, 481

Zinnecker H., Keable C. J., Dunlop J. S., Cannon R. D., Griffiths W. K., 1988, in Grindlay J. E., Philip A. G. D., eds, IAU Symposium Vol. 126, The Harlow-Shapley Symposium on Globular Cluster Systems in Galaxies. p. 603

Zuckerman B., Evans N. J. I., 1974, ApJ, 192, L149

Zuckerman B., Palmer P., 1974, ARA&A, 12, 279

de Grijs R., Gilmore G. F., Johnson R. A., Mackey A. D., 2002, MNRAS, 331, 245
de Mink S. E., Pols O. R., Langer N., Izzard R. G., 2009, A&A, 507, L1