Star clusters stand at the intersection of much of modern astrophysics the ISM gravitational dynamics stellar evolution and cosmology. Here we review observations and theoretical models for the formation evolution and eventual disruption of star clusters. Current literature suggests a picture of this life cycle including the following several phases Clusters form in hierarchically structured accreting molecular clouds that convert gas into stars at a low rate per dynamical time until feedback disperses the gas. The densest parts of the hierarchy resist gas removal long enough to reach high starformation efficiency becoming dynamically relaxed and well mixed. These remain bound after gas removal. In the first 100 Myr after gas removal clusters disperse moderately fast through a combination of mass loss and tidal shocks by dense molecular structures in the starforming environment. After 100 Myr clusters lose mass via twobody relaxation and shocks by giant molecular clouds processes that preferentially affect lowmass clusters and cause a turnover in the cluster mass function to appear on 110Gyr timescales. Even after dispersal some clusters remain coherent and thus detectable in chemical or action space for multiple galactic orbits.In the next decade a new generation of space and adaptive opticsassisted groundbased telescopes will enable us to test and refine this picture.
['Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Solar and Stellar Astrophysics']
Star Clusters Across Cosmic Time
{'Mark R. Krumholz 1 , 2 , Christopher F. McKee 3 , and Joss Bland-Hawthorn 2 , 4 , 5': "1 Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia; email: mark.krumholz@anu.edu.au 2 Centre of Excellence for All Sky Astrophysics in Three Dimensions (ASTRO-3D), Australia 3 Departments of Astronomy and of Physics, University of California, Berkeley, CA 94720, USA 4 Sydney Institute for Astronomy, School of Physics A28, University of Sydney, NSW 2006, Australia 5 Miller Professor, Miller Institute, University of California, Berkeley, CA 94720, USA  \nXxxx. Xxx. Xxx. Xxx. 2019. 59:1-76  \nThis article's doi: 10.1146/((please add article doi))  \nCopyright c © 2019 by Annual Reviews. All rights reserved", 'Keywords': 'star clusters, star formation, globular clusters, open clusters and associations, stellar abundances', 'Abstract': 'Star clusters stand at the intersection of much of modern astrophysics: the interstellar medium, gravitational dynamics, stellar evolution, and cosmology. Here we review observations and theoretical models for the formation, evolution, and eventual disruption of star clusters. Current literature suggests a picture of this life cycle with several phases:  \n- · Clusters form in hierarchically-structured, accreting molecular clouds that convert gas into stars at a low rate per dynamical time until feedback disperses the gas.\n- · The densest parts of the hierarchy resist gas removal long enough to reach high star formation efficiency, becoming dynamicallyrelaxed and well-mixed. These remain bound after gas removal.\n- · In the first ∼ 100 Myr after gas removal, clusters disperse moderately fast, through a combination of mass loss and tidal shocks by dense molecular structures in the star-forming environment.\n- · After ∼ 100 Myr, clusters lose mass via two-body relaxation and shocks by giant molecular clouds, processes that preferentially affect low-mass clusters and cause a turnover in the cluster mass function to appear on ∼ 1 -10 Gyr timescales.\n- · Even after dispersal, some clusters remain coherent and thus detectable in chemical or action space for multiple galactic orbits.  \nIn the next decade a new generation of space- and AO-assisted ground-based telescopes will enable us to test and refine this picture.', '1. INTRODUCTION': 'The galaxy is in fact nothing but a congeries of innumerable stars grouped together in clusters. . .  \nSiderius Nuncius (Sidereal Messenger), Galileo Galilei, 1610', '1.1. Historical background and motivation': "After four centuries since Galileo's first observations, star clusters remain mysterious. What are the unifying principles that govern the formation of globular and open clusters, groups, associations and superassociations? As illustrated in Figure 1 1 , clusters cover a huge range of mass, size, and density scales. The clusters shown span ages from ∼ 1 Myr to > 10 Gyr, and masses from ∼ 10 2 -10 6 M glyph[circledot] . Some are so  \nONC, T < 3 Myr  \nArches,  \nT  \n2  \n3 Myr  \nNGC 265,  \nT  \n250 Myr  \nColl 261, T ≈ 7 Gyr  \n<!-- image -->  \n≈  \n-  \nNGC 6535, T > 10 Gyr  \n47 Tuc, T > 10 Gyr  \nNGC 1252: Fake News  \n<!-- image -->", 'Figure 1': 'Images of a range of star clusters, along with NGC 1252, an object previously classified as a cluster but now known to be an asterism. The field of view in all frames is 3 pc × 3 pc, and North is up; angular sizes are indicated by scale bars. Image sources are: Orion Nebular Cluster (ONC) - Robberto et al. (2013); Arches cluster - NASA & ESA; NGC 265 - NASA & ESA; Davide De Martin (ESA/Hubble) and Edward W. Olszewski (University of Arizona, USA); Hyades - NASA, ESA, & STScI; Collinder 261 - ESO/Digitized Sky Survey; NGC 6535 ESA/Hubble & NASA, Gilles Chapdelaine; 47 Tuc - NASA, ESA, and the Hubble Heritage  \n(STScI/AURA)-ESA/Hubble Collaboration; J. Mack (STScI) and G. Piotto (University of Padova, Italy); NGC 1252 - WEBDA database, https://www.univie.ac.at/webda/ .  \ncompact and rich that stars become lost in confusion in the 3 pc-frames shown, others are so sparse and extended that most cluster members are outside the frame. Some are classified as open clusters (Arches, NGC 265, Hyades, Coll 261), some as globular clusters (NGC 6535, 47 Tuc). The Orion Nebula Cluster (ONC) is still forming and, depending on the author, might not even be classified as a cluster at all. NGC 1252 had been classified as an open cluster since 1888, but in 2018 was shown to be merely an asterism (Kos et al. 2018b). Cluster formation is central to the star formation process. Conceivably, all stars formed in groups, clusters, or hierarchies, although, for this to be true, most clusters must have dissolved into the Galactic background soon after formation. However, our understanding of when, how, and why stars cluster remains primitive (see reviews by Krumholz et al. 2014, Renaud 2018, and Adamo & Bastian 2018).  \nA review of clusters is timely for several reasons. One is the emerging overlap between the traditionally separate communities focused on globular clusters and modern-day star clusters and star formation. These have been separate in the past because star clusters in the disk of the Milky Way, generally classified as open clusters (OCs), and those in the halo, generally classified as globular clusters (GCs), appear to occupy distinct loci in the parameter space of mass, age, and metallicity. New data have begun to blur this distinction: even in the Milky Way there is substantial overlap in both metallicity and density between OCs  \n≈  \nHyades, T ≈ 625 Myr  \nand GCs, and in rapidly star-forming local galaxies such as the Antennae (Zhang & Fall 1999) and M82 (McCrady & Graham 2007) there is overlap in mass as well. While there is still a bimodality in the star cluster age and spatial distribution, this may simply reflect the combined effects of cluster dissolution and the history of galaxy assembly. In this view, toward which we shall mostly tend here (see also Kruijssen 2014 and Forbes et al. 2018), GCs are no different than any other type of star cluster in their formation and internal dynamics; only their cosmological history is different. Their origins are therefore best addressed in the context of general models for the formation and evolution of star clusters, hence the need for this review. 2  \nAn additional motivation for this review is to bring together the results of modern studies of stellar kinematics and abundances with the formation and evolution of star clusters. As with GCs and OCs, these have traditionally been somewhat separate fields, with the stellar kinematics and abundances community focused more on the long-term dynamical evolution of stars in and around the Galaxy, and the star formation community paying more attention to the dynamics of interstellar gas and the effects of star formation and stellar feedback on it. This separation is no longer viable. In the era of Gaia and massive spectroscopic surveys, kinematic and abundance data are now becoming accurate enough that it should be possible to trace stars that are now part of the field back to their birth places, or at least to reconstruct some of the now-dissolved structures in which they were born. To guide such reconstructions, however, will require theoretical and observational input from the star formation community.  \nWith these motivations in mind, we first review current observational constraints on star cluster populations, and then discuss our current understanding, or lack therefore, of the processes by which clusters form, evolve, and eventually disperse.', '1.2. Prelude: what is a star cluster?': '1.2.1. The need for a definition. The first question that any discussion of star clusters must face, going back at least to Trumpler (1930), is how to define a cluster and distinguish clusters from multiple systems. In their recent review, Portegies Zwart, McMillan & Gieles (2010) define a star cluster as a set of stars that are gravitationally bound to one another, while the earlier review by Lada & Lada (2003) defines a cluster as a collection of stars with a mass density large enough ( glyph[greaterorsimilar] 1 M glyph[circledot] pc -3 ) to resist tidal disruption in Solar neighborhood conditions, and numerous enough to avoid N -body evaporation for at least 100 Myr. One might also consider defining clusters as concentrations of stars whose mass density significantly exceeds the mean in their galactic neighborhood ( ≈ 0 . 1 M glyph[circledot] pc -3 near the Sun - McKee, Parravano & Hollenbach 2015). For older stellar populations these definitions are in practice nearly identical. A group of stars of mean density ρ ∗ and radius r is unbound only if its velocity dispersion σ glyph[greaterorsimilar] r √ Gρ ∗ , so the time required for it to disperse is of order the crossing time t cr = r/σ glyph[lessorsimilar] 1 / √ Gρ ∗ . At the density threshold proposed by Lada & Lada (2003), this is glyph[lessorsimilar] 10 Myr, so overdensities older than this must be either held together by self-gravity or external forces, or be short-lived transients.  \nIn this review, however, we are interested in the formation of star clusters as well as their evolution, and thus we cannot limit ourselves to stellar populations that are older than their internal dynamical time (c.f. Portegies Zwart, McMillan & Gieles 2010; Gieles & Portegies Zwart 2011). However, younger populations often have complex internal structures such that there is no single way to decompose them into clusters. We illustrate this point with the example of the ONC in Figure 2, where the left two panels show the results of applying two different cluster identification methods that show little agreement: a Bayesian  \nSurveys referenced in this review:', 'PHAT:': "Panchromatic Hubble Andromeda Treasury (Dalcanton et al. 2012)  \nLEGUS: Legacy Extragalactic UV Survey (Calzetti et al. 2015)  \nGOALS: Great Observatories All-sky LIRG Survey (Armus et al. 2009)  \nGRS: Galactic Ring Survey (Jackson et al. 2006)  \nATLASGAL: APEX  \nTelescope Large Area Survey of the Galaxy (Schuller et al. 2009)  \nGALAH: Galactic Archaeology with HERMES (De Silva et al. 2015)  \nRAVE: Radial Velocity Experiment (Steinmetz et al. 2006)  \nFigure 2  \n<!-- image -->  \nSample decompositions of the stars in and around the Orion Nebula Cluster (ONC). In the left two panels, gray dots show stars from the catalog of Da Rio et al. (2009); the axes indicate offset from the position of θ 1 c Ori, the most luminous star in the cluster. In the left panel, ellipses mark the clusters identified by Kuhn et al. (2014) based on a Bayesian decomposition into isothermal ellipsoids. The middle panel shows a minimum spanning tree (MST) decomposition, performed on the Da Rio et al. (2009) catalog using a pruning length of 0 . ' 3; different colours indicate the different sub-clusters identified by MST. The right panel is a synthetic UBV image produced by placing the ONC to a distance of 1 Mpc (assuming that the true distance is 389 pc - Kounkel et al. 2018), and assuming diffraction-limited imaging with a 2.4 m-diameter telescope such as the Hubble Space Telescope . The region shown in the right panel is the same as shown in the two other panels. Colors are scaled logarithmically to give a dynamic range of 100 in flux for all channels.  \ndecomposition into isothermal ellipsoids (Kuhn et al. 2014, left), and a minimum spanning tree (middle). Numerous other decompositions are possible as well (Schmeja 2011). It is not clear which of these, if any, correspond to physically meaningful concepts such as a distinction between structures with positive or negative total energy. Moreover, none of the clustering methods used on resolved stars in the Milky Way are likely to correspond to the clusters picked out by most extragalactic surveys. As the right panel of Figure 2 shows, the ONC would certainly be seen as a single cluster (were it detected at all) in a galaxy at a distance of ≈ 1 Mpc, as targeted by the PHAT survey (Dalcanton et al. 2012, see Margin Note), or ∼ 3 -10 Mpc, typical of the LEGUS survey (Calzetti et al. 2015). We return to the question of how to handle the ambiguity in defining clusters at young ages in § 2.7.  \n1.2.2. Working definition. Given the discussion above, we shall follow Trumpler (1930) and define clusters very generally, so that they can be studied throughout their lives. For our purposes a star cluster is a group of at least 12 stars (to clearly distinguish it from a multiple star system - Trumpler 1930) with a mean density that (i) is at least a factor of a few times the background density (similar to the constraint adopted by Lada & Lada 2003), with the factor being larger in larger galaxies so that the overdensity is statistically significant, and (ii) is much greater than the local density of dark matter. Such groups of stars are physically associated while they form, but there is no constraint on the formation timescale or lifetime; in contrast to Lada & Lada (2003), we do not require the cluster to be denser than the gas cloud out of which it forms, nor do we require that it have enough members to avoid evaporation for 100 Myr. Clusters generally dissipate in time ( § 4), and once they no longer satisfy our criteria we refer to them as 'dissolved clusters' ( § 5).  \nStar clusters can be gravitationally bound or unbound, and young clusters can be a mixture. Unbound  \nclusters are termed associations; bound clusters can be categorized as OCs or GCs. The most massive and young OCs are sometimes referred to as young massive clusters (YMCs) or super-star clusters (SSCs); they are sometimes also called young GCs, although it is unknown if they will in fact evolve into clusters like the GCs we observe today. YMCs have traditionally been defined by mass > 10 4 M glyph[circledot] and age < 100 Myr, but the lines are fuzzy and there is no particular physical significance to them. Similarly, OCs and GCs do not have precise definitions. In the Milky Way, most OCs have masses glyph[lessorsimilar] 5000 M glyph[circledot] and ages glyph[lessorsimilar] 6 Gyr, whereas most GCs have masses glyph[greaterorsimilar] 10 4 M glyph[circledot] and ages glyph[greaterorsimilar] 6 Gyr (Kharchenko et al. 2013), but in other galaxies these catagories have more overlap. The metallicity of OCs overlaps that of the thick-disk population of Galactic GCs, which has [Fe/H] > -0 . 8 (Zinn 1985). Clusters with ages glyph[lessorsimilar] few hundred Myr are often organized hierarchically, with OCs and YMCs being embedded in associations with sizes up to a few hundred pc, associations grouped together in complexes at ∼ kpc scales, etc.; there is no physical division between these scales for the unbound structures (Gouliermis 2018). Despite the wide range of physical properties of these various types of star clusters, there is no evidence that they do not form via the same basic mechanism (Elmegreen & Efremov 1997; § 3). The only clear intrinsic difference among these different types of clusters is that old GCs exhibit anti-correlations in some light element abundances and have multiple stellar populations, hallmarks that are never observed in OCs ( § 2.6). Since these phenomena, which are not understood, have been reviewed recently by Renzini (2013) and Bastian & Lardo (2018), we shall not discuss them further here.", '2. DEMOGRAPHICS: OBSERVATIONAL CONSTRAINTS': "The physical processes that govern the formation and evolution of star clusters encode themselves in the distributions of star clusters' properties - mass, age, size, etc. Any review of star clusters must therefore begin with the observational constraints and how they are obtained. In this section and the subsequent ones we will define a number of symbols, which we summarize for convenience in Table 1.", '2.1. Methods': "Methods for determining the demographics of star clusters can be divided into those that operate on resolved stellar populations, and those that operate on unresolved populations.  \n2.1.1. Resolved Stellar Populations. For star clusters in the Milky Way, the Magellanic Clouds, and M31 (using Hubble Space Telescope resolution), it is possible to resolve individual main sequence stars down to masses of a few M glyph[circledot] or less. This statement oversimplifies the situation somewhat; the luminosity or mass limit down to which it is possible to resolve stars is a complex function of the stellar density and luminosity distribution in the region being studied. It is possible to resolve massive main sequence stars in relatively sparse regions even beyond the Local Group (e.g., Larsen et al. 2011), while confusion is a problem for the most crowded regions even in the Milky Way (e.g., Ascenso, Alves & Lago 2009). There exists a broad range of parameter space where stars are partially resolved, i.e., their separations are comparable to the observational point spread function (psf); the statistical and photometric techniques used to analyze such fields are beyond the scope of this review.  \nWhen individual stars can be resolved, placing them on a color-magnitude diagram (CMD) provides the most direct method of determining the properties of the parent star cluster. In principle, the CMD allows one to read off the masses of the individual stars almost directly and, if the main sequence turn-off is resolved, to infer the age of the population as well. By adding spectroscopy to photometry, one can also determine the stars' compositions. While this method is direct, it does encounter complications. In young clusters where the resolved stars are pre-main sequence, the choice of evolutionary tracks can induce  \nPDF: probability distribution function  \nInitial mass function (IMF): PDF of initial mass distribution of stars in a stellar ensemble  \nPresent day mass function (PDMF): PDF of observed mass distribution of stars in a stellar ensemble  \nInitial cluster mass function (ICMF): PDF of initial mass distribution of star clusters", 'Cluster mass': "function (CMF):  \nPDF of observed mass distribution of star clusters  \nCluster age function (CAF): PDF of observed age distribution of star clusters  \nTable 1 Symbols used in this review  \n§  \nsignificant uncertainties in the final properties (e.g., Da Rio et al. 2012), and a non-negligible range of ages may be present, so that the system cannot be described by a single age (e.g., Reggiani et al. 2011; Getman et al. 2018). If the observations do not resolve down to the peak of the initial mass function (IMF), one may need to extrapolate using an assumed IMF to account for the mass of unresolved stars (e.g., Lada & Lada 2003; Beerman et al. 2012).  \nDespite these limitations, however, two other uncertainties are more important. One is the error introduced by the need to correct demographics for clusters that are too small or too old to make it into the catalog. For resolved observations, catalog completeness is invariably a function not just of the mass and age of the stellar population, but also its concentration, the exact masses and ages of its most luminous few stars, and the density of the background (Silva-Villa & Larsen 2011; Johnson et al. 2015b). The other source of uncertainty is in assignment of stars to clusters. As already discussed in § 1.2, for populations that are not old enough to have relaxed, there may be no unique way to decompose a collection of stars into clusters, and variations in how one defines a cluster can lead to order of magnitude differences in the inferred number of clusters and their properties (Bressert et al. 2010). However, even for somewhat older populations, catalogs constructed using different methods can easily differ by ∼ 10 -20% even when the underlying data are identical (Johnson et al. 2015b). These errors are generally larger than the uncertainties on the properties of individual clusters.  \nFigure 3  \n<!-- image -->  \n-  \n-  \n-  \nColor evolution of a simple stellar population with a fully sampled IMF, including nebular emission. In each panel the line shows the intrinsic color (in AB mag) as a function of population age. Arrows, placed at intervals of 0 . 25 dex in age starting at log( T/ yr) = 6 . 25, show observed color for a visual extinction A V = 0 . 5 mag for the stars, A V = 1 . 05 mag for nebular light. This calculation uses SLUG stellar population synthesis code (da Silva, Fumagalli & Krumholz 2012; Krumholz et al. 2015b), run with a Chabrier (2005) IMF, Solar metallicity, MIST rotating stellar evolution models (Choi et al. 2016), starburst99 atmosphere models (Leitherer et al. 1999), a Calzetti et al. (2000) starburst dust attenuation curve, and HST/WFC3 filters. However, other choices for these parameters yield qualitatively similar results.  \n2.1.2. Unresolved Populations. While CMDs offer the most reliable method for determining star cluster properties, at distances glyph[greaterorsimilar] 1 Mpc it is not possible to resolve any but the most luminous individual stars. Consequently, studies that probe a wide range of environments must use unresolved populations. In an unresolved observation, a star cluster is distinguished from a single star or other point source either by the fact that its light is more extended than the observational psf (for targets at distances below tens of Mpc) or that it is too luminous to be a single massive star (for more distant targets). In the earliest work with unresolved populations authors simply analyzed luminosity functions (e.g., Whitmore et al. 1993, 1999), but one can also infer physical properties using simple stellar population (SSP) models. In practice this procedure can be carried out in several different ways. The simplest is to generate theoretical luminosities as a function of mass and age by assuming a stellar population that fully samples the IMF (including models for nebular emission and dust extinction) and find the mass and age that best match the observations using χ 2 or a similar goodness-of-fit statistic (e.g., Zhang & Fall 1999; Adamo et al. 2010). More complex methods take into account the stochastic fluctuations in light output that occur as a result of partial sampling of the IMF in clusters smaller than ∼ 10 4 M glyph[circledot] (e.g., Fouesneau & Lan¸con 2010; Popescu, Hanson & Elmegreen 2012; de Meulenaer et al. 2013; Krumholz et al. 2015a) or add narrowband photometry targeting nebular lines (usually H α ) to provide better age discrimination (e.g., Fouesneau et al. 2012; Bastian et al. 2012; Chandar et al. 2016).  \nFitting star cluster properties from photometry leads to substantially larger errors than doing so from resolved measurements of individual stars, with the magnitude of the error depending on the locus in color space. We illustrate this point in Figure 3, which shows the trajectory that an SSP traces out in color over its lifetime. The process of inferring cluster age and extinction essentially amounts to placing an observed point on a multidimensional color-color plot such as Figure 3, and finding the location on an extinction-adjusted SSP track such as the one shown that gets as close as possible to the observed point in every dimension. As is apparent from the figure, for ages below ≈ 5 Myr or above ≈ 50 Myr, in at least some color combinations, SSP tracks are relatively straight and not parallel to the extinction vector. Consequently, the matching process is straightforward and there is little ambiguity in the determination of ages; comparing photometric ages to those inferred from CMDs, Elson & Fall (1985, 1988) estimate that  \nExclusive catalog: a star cluster catalog that excludes objects that do not pass morphological tests for compactness and roundness  \nInclusive catalog: a star cluster catalog without additional morphological requirements  \nFigure 4  \n<!-- image -->  \nExamples of unresolved clusters of different morphological classes. Each panel shows a three-color UBV image of a star cluster in NGC 628, from Adamo et al. (2017). The ring shows a radius of 0 . '' 28, approximately 13 . 4 pc at the distance of NGC 628. The morphologies are classified as, from left to right, compact and symmetric, compact and asymmetric, and multiply-peaked. An 'exclusive' catalog, in the sense used in this review, would include the left two objects but exclude the one on the right, while an 'inclusive' catalog would include all three. A comparison to Figure 2 suggests that the ONC might well be excluded from an 'exclusive' catalog.  \nphotometric ages for massive clusters are accurate to ≈ 0 . 3 dex. At ages of ≈ 5 -50 Myr, however, SSP colors oscillate rapidly and tracks in color space frequently cross themselves, making age determinations quite uncertain. For example, notice that SSPs with ages of ≈ 8 Myr (dark green in the figure) have very similar colors to SSPs with ages of ≈ 30 Myr (light blue in the figure), and that the separation between them is mostly along the extinction vector, making these two possibilities difficult to disentangle, particularly once stochastic color variation due to finite IMF sampling is taken into account (also see Ma'ız Apell'aniz 2009 and Krumholz et al. 2015a).  \nOn top of these errors that apply to individual clusters, catalog construction and completeness correction induce additional errors on measurements of population demographics. Completeness correction for unresolved observations faces essentially the same challenges as for resolved data. Recently Krumholz et al. (2018) proposed a method to carry out full forward-modelling of an entire photometric catalog, simultaneously accounting for completeness and stochasticity, but this method is thus far experimental. Catalog construction is the more serious issue, because it raises the question of how one should handle sources where the light distribution is multiply-peaked or contains significant color gradients, as illustrated in Figure 4. Some authors, who limit the definition of cluster to dynamically-relaxed structures (c.f. § 1.2), discard such sources on the grounds that they are not relaxed (e.g., Bastian et al. 2011, 2012), while others adopt a broader definition and thus retain them (e.g., Chandar et al. 2010, 2011, 2014); some split the difference by including all objects but then reporting separate analyses for objects of differing morphology (e.g., Adamo et al. 2017). These different methods of catalog construction have historically led to extreme confusion in interpretation of observations, since inclusive catalogs tend to have ≈ 30 -50% more members in the youngest age bin than exclusive ones, with much smaller differences at older ages (Chandar et al. 2014; Messa et al. 2018b). In the discussion that follows, we will minimize this confusion by identifying whether a particular catalog was constructed using criteria that are 'exclusive' (i.e., the catalog excludes nonsymmetric objects) or 'inclusive', and pointing out when the choice of one method or the other leads to systematic differences. However, we caution that all existing cluster catalogs for galaxies at distances glyph[greaterorsimilar] 20 Mpc are necessarily inclusive, because limited resolution makes morphological measurement impossible in  \nsuch distant samples.  \nAuthors who build exclusive catalogs generally assume that morphology can be used as a proxy for boundedness. This likely holds in one direction: since star-forming regions are morphologically complex, if the observed stars have a round, compact morphology, it is likely that they have relaxed into it and thus are bound. However, the converse need not be true, i.e., there is no reason to assume that a population that is too young to have dispersed is unbound simply because it is not round. The implication is that, for young stellar populations, an exclusive approach probably omits everything that is not bound, but also discards some unknown number of bound systems. By contrast, an inclusive catalog captures all structures above some luminosity threshold without regard to their dynamical state.", '2.2. Mass distribution': "The most basic property of a star cluster is its mass, and thus the most basic distribution for star clusters within a galaxy is the observed cluster mass function (CMF). 3 Figure 5 summarises recent measurements for the CMFs of star clusters in the disks of Local Group galaxies. All observed disk CMFs appear to be reasonably well-described by a power-law dN/dM ∝ M α M , with values of α M = -2 ± 0 . 2. 4 A slope of α M = -2 corresponds to equal mass per logarithmic bin, and thus is the expected slope for a completely scale-free distribution. In normal star-forming galaxies, the mass function slope is generally measured over the range from ≈ 10 3 -10 5 M glyph[circledot] . The lower limit on this range is entirely a function of observational limitations - only in a few cases do we have cluster samples that are complete enough below ≈ 10 3 M glyph[circledot] to enable measurement of a CMF, and in those cases the data are consistent with α M ≈ -2 down to the completeness limit. Indeed, some authors have claimed that the power law continues to masses as small as ≈ 20 M glyph[circledot] , corresponding to individual massive stars (e.g., Lamb et al. 2010).  \nThe nature of the upper limit is less clear. In most galaxies, the total number of clusters is such that very few with masses glyph[greaterorsimilar] 10 5 M glyph[circledot] would be expected even for a pure power-law mass function, and in galaxies with cluster populations large enough that such clusters would be expected - the Antennae, M82, and various luminous infrared galaxies (LIRGs) - they are found. Consistent with this, the luminosity of the brightest cluster in a galaxy is well correlated with the luminosity of the galaxy and its star formation rate (Larsen 2002). However, some authors report evidence that the observed CMF or luminosity function is better fit by power-law that is truncated at a mass M c (either a hard truncation or a Schechter function, dN/dM ∝ M α M e -M/M c ) than a pure power-law (e.g., Bastian 2008; Larsen 2009; Bastian et al. 2012; Adamo et al. 2015, 2017); others question the statistical robustness of this conclusion (Mok, Chandar & Fall 2018). Any truncation would need to vary systematically with some other galaxy property (e.g., the Jeans mass in the disk, § 3.1.1) in order to explain the presence of more massive clusters in more luminous galaxies.  \nThe question of whether there is a truncation has been difficult to settle due to the challenge of deriving realistic confidence intervals from observations where the dominant errors are catalog construction and translation from photometry to cluster masses; some authors using nearly identical data sets nonetheless reach differing conclusions about whether they provide evidence for a truncation (e.g., Adamo et al. 2015 versus Sun et al. 2016 on M83), and the results can be sensitive to the statistical power of the fitting method  \nFigure 5  \n<!-- image -->  \n/circledot  \nCompilation of measured CMF slopes α M . In the main panel, colored points represent clusters in the disks of nearby galaxies, and solid and dashed black lines are evolved Schechter function fits to the Milky Way's GCs and to an average of GCs in Virgo, respectively. The horizontal error bars show the range in cluster mass over which the measurement was made, with open points indicating measurements from exclusive catalogs and filled points indicating inclusive catalogs. Points are displaced slightly from the centers of their corresponding mass ranges to minimize confusion. A straight dotted line on the high mass end of the horizontal error bar indicates that the observation did not detect a truncation to the CMF, while a downwards curve indicates an observation that reported a truncation. For the GC mass functions, α M is not constant in mass, and the line shows α M versus mass for the best-fitting evolved Schechter function. In the main panel only whole-galaxy average measurements are shown. The histogram on the right shows the distribution of measured CMF slopes for disk clusters; blue is the distribution of whole-galaxy average measurements, while hatched shows measurements in which galaxies are broken into independent sub-samples. In the legend, N is short for NGC, and the galaxies listed as ∗ , † , and GOALS are UGC09618NED02, IRAS20351+2521, and an average of 22 LIRGs from the GOALS sample, respectively. References listed in the legend are as follows: MG07 = McCrady & Graham (2007), J07 = Jord'an et al. (2007), FC12 = Fall & Chandar (2012), B13 = Baumgardt et al. (2013), R14 = Ryon et al. (2014), C15 = Chandar, Fall & Whitmore (2015), A15 = Adamo et al. (2015), C16 = Chandar et al. (2016), H16 = Hollyhead et al. (2016), M16 = Mulia, Chandar & Whitmore (2016), J17 = Johnson et al. (2017), A17 = Adamo et al. (2017), L17 = Linden et al. (2017), M18 = Messa et al. (2018b). The hatched histogram on the right also contains sub-samples from Messa et al. (2018a).  \n(Messa et al. 2018a). The most convincing case for a deviation from a pure power-law mass function is in M31, where Johnson et al. (2017) robustly identify a truncation mass of M c ≈ 8 . 5 × 10 3 M glyph[circledot] based on a sample of ≈ 1000 clusters with masses determined from CMDs rather than unresolved photometry. However, this survey does not cover some of the more actively star-forming parts of M31, and it is possible these might host clusters larger than M c . There is also indirect evidence for a truncated CMF in the Milky Way, in the form of an observed truncation in the luminosity function for radio recombination line emission from H ii regions (McKee & Williams 1997), although translating this into a truncation mass is not straightforward.  \nValues of α M tend to be slightly larger for starbursting galaxies than for quiescent spirals. However, it is unclear if this represents a real physical difference or an observational artifact: the starburst galaxies are located at systematically larger distances, and thus observations of them have systematically lower spatial resolution, which will tend to flatten the mass function by blurring multiple small objects into a single large one. This matter will be 'resolved' in the next decade by AO-assisted imagers and spectrographs on  \nextremely large telescopes.  \nThe mass function for the GCs found in the halos of galaxies is substantially different. It has a clear turnover mass below which the number of clusters is either flat or increasing with mass. The distribution can be reasonably well described as a lognormal with a characteristic mass of ≈ 10 5 . 3 M glyph[circledot] (Brodie & Strader 2006). However, a more physically-motivated description is the evolved Schechter function (Fall & Zhang 2001; Jord'an et al. 2007), which is simply a standard Schechter function modified by the assumption that all clusters have lost a fixed amount of mass ∆; we discuss the physical motivation for this model in § 4.2.2. This gives a distribution dN/dM ∝ ( M +∆) α M,i exp[ -( M +∆) /M c ], where α M,i is the slope at formation. The corresponding present-day slope is  \nα M = d d log M ( log dN dM ) = -M ( M -α M,i M c +∆) M c ( M +∆) , (1)  \nIn Figure 5, we show α M for the Milky Way GC system and for the average of GC systems in the Virgo cluster from Jord'an et al. (2007); 5 the best-fit M c and ∆ values are both in the range 10 5 . 3 -10 6 . 6 M glyph[circledot] , spanning the range of truncation masses inferred for disk clusters. However, the functional forms of the GC and disk cluster mass functions are clearly very different, and the relationship between the two is a topic to which we shall return below.", '2.3. Age distribution': 'A second distribution of equal importance to the CMF is the age distribution of star clusters, the cluster age function (CAF). For Milky Way GCs this distribution is centred at ages ∼ 10 Gyr with a width of a few Gyr (Bastian & Lardo 2018). For clusters in disks, observers generally fit the CAF as a power-law dN/dT ∝ T α T , where T is the cluster age 6 ; some observers use a single α T for all clusters, while others consider more complex forms where α T can be different for different age or mass ranges. However, in practice most extragalactic samples cover a small enough dynamic range in mass and age that one can define at most two or three bins with potentially differing α T values.  \nIn a galaxy with a constant cluster formation rate, α T is a measure of how long clusters survive. Consider a simple model, in the spirit of Boutloukos & Lamers (2003), where the probability per unit time that a cluster of age T is destroyed is 1 /aT for some constant factor a . If we consider a cohort of clusters at a single age, our assumed destruction rate implies that their number will decrease with time as dN/dT = -N/aT . Integrating, we find that the number that survive to age T is N = N 0 ( T/T 0 ) -1 /a , where N 0 is the number of clusters at age T 0 . This corresponds to a CAF with α T = -1 /a . Thus α T → 0 corresponds to the case a glyph[greatermuch] 1, i.e., clusters survive for many times their current age, while α T = -1 corresponds to a = 1, i.e., the expected time required to destroy a cluster is equal to its current age.  \nWhile the interpretation of α T is straightforward, divining a consistent value of it from the published literature is much more difficult. In part this is because different groups tend to subdivide the sample to be fit in different ways, and to reach divergent conclusions both about the value of α T and whether it is different at different ranges of cluster age or mass. Figure 6 shows a summary of recent observational determinations for extragalactic systems, along with our fit to the completeness-corrected catalog of Piskunov et al. (2018)  \nα  \nFigure 6 shows that there is a systematic difference in slopes between inclusive and exclusive catalogs at young ages. The fits from exclusive catalogs at T glyph[lessorsimilar] 10 8 yr have α T ≈ -0 . 3 to -0 . 2, while the inclusive catalogs cluster around α T ≈ -1 to -0 . 7. Thus exclusive catalogs suggest that clusters have long survival times compared to their current ages, while inclusive catalogs imply survival times comparable to cluster ages. In addition, although this is not shown in the figure, analysis of exclusive catalogs suggests that more massive clusters have larger values of α T than less massive ones (e.g., Silva-Villa et al. 2014; Adamo et al. 2017; Messa et al. 2018b), implying that cluster survival time is mass-dependent. A second, more subtle dependence worth noting is that starburst galaxies seem to have a slight bias toward slopes near α T = -1 compared to normal spirals and dwarfs - for example, the GOALS galaxies and the Antennae all have α T glyph[lessorsimilar] -0 . 9, while, even for inclusive catalogs, the mean slope for spirals is closer to α T ≈ -0 . 7. A similar increase in α T is seen in the centers of massive spirals compared to their outskirts (e.g., SilvaVilla et al. 2014; Messa et al. 2018a). While such an environmental dependence is in accord with some theoretical expectations (e.g., Kruijssen et al. 2011), one should be cautious in putting too much weight on the observations. The regions with steeper α T values are on average more crowded and distant than the  \n<!-- image -->', 'Figure 6': 'Observational determinations of the CAF in galaxy disks. The horizontal range indicates the age interval over which the measurement was made, while the value plotted in the vertical direction is the index α T ( dN/dT ∝ T α T ) with error bars. For extragalactic data, open points correspond to exclusive catalogs (left panel), filled points to inclusive ones (middle panel). The galaxy and source are as indicated in the legend (using the same abbreviations as in Figure 5), and multiple points for a single galaxy correspond to subsamples of clusters in that galaxy, separated either by mass or by location; we have added a small random displacement to the horizontal positions of points to prevent them from overlapping. The black line shown in both the left and middle panels is a smoothed basis spline fit to the catalog of Piskunov et al. (2018) for clusters within 2 kpc in the Sun. The right panel shows a histogram of α T values for the extragalactic catalogs, with inclusive catalogs and exclusive catalogs for ages below 10 8 yr shown separately. References are as follows: F14 = Fouesneau et al. (2014), SV14 = Silva-Villa et al. (2014), R14 = Ryon et al. (2014), C17 = Chandar et al. (2017), A17 = Adamo et al. (2017), L17 = Linden et al. (2017), M18 = Messa et al. (2018a,b), P18 = Piskunov et al. (2018).', 'for Milky Way clusters within ≈ 2 kpc of the Sun. 7': 'Figure 7  \n<!-- image -->  \nCAF for star clusters within ≈ 2 kpc of the Sun. Points with error bars show measurements from Piskunov et al. (2018, P18). Lines show χ 2 fits to several functional forms, with parameters indicated in the legend: broken power-law (orange; dN/dT ∝ T α T, 1 for log T < log T b , and dN/dT ∝ T α T, 2 for log T > log T b ), Schechter function (green; dN/dT ∝ T α T e -T/T ∗ ), single power-law (faded red; dN/dT ∝ T α T ), and exponential (faded purple; dN/dT ∝ e -T/T ∗ ). The broken power-law and Schechter forms both have reasonable χ 2 values, while the other fits are poor.  \nsystems with shallower α T , and both crowding and distance could create a systematic bias in measured color that would manifest as a difference in the inferred age distribution.  \nMilky Way clusters show a complex dependence of α T on age, which we highlight in Figure 7. In the Milky Way, α T ≈ -0 . 5 at ages glyph[lessorsimilar] 10 9 yr, indicating moderately strong cluster destruction. The age distribution steepens sharply above ≈ 10 9 yr, indicating much more rapid disruption. This change in slope is not seen in the extragalactic samples, which are generally limited to clusters younger than ≈ 10 8 yr for reasons of sensitivity. However, there are also two additional cautions to be mentioned in comparing the Milky Way and extragalactic samples. First, because the Milky Way sample is limited to glyph[lessorsimilar] 2 kpc from the Sun, it consists entirely of low-mass clusters ( M glyph[lessorsimilar] 10 3 M glyph[circledot] ), while the extragalactic sample is for much larger masses, M glyph[greaterorsimilar] 10 3 . 5 M glyph[circledot] . Second, the extragalactic sample ages are derived from photometry, while the Milky Way ages are based on CMDs, which may lead to systematic differences. CMDs are more reliable in general, but this might not hold for the Milky Way sample because many of the clusters in it contain only a small number of stars bright enough to allow placement on the CMD, leading to large age uncertainties.', '2.4. Bound mass fraction': "A third basic statistic for clusters is the fraction of the total stellar mass in gravitationally-bound clusters, denoted Γ. Since this fraction changes with stellar age unless α T = 0, Γ in general is a function of T . To measure Γ one must determine the total mass in bound clusters and the total mass in all stars within the same age interval. Measuring the former invariably involves some degree of extrapolation to account for the mass of clusters too small to detect, but for α M ≈ -2 the extrapolation is fairly modest. For the latter quantity, at ages above ∼ 10 Myr the most reliable means of determining total stellar mass is from CMDs of field stars. However, these are only available for relatively nearby sources. A second-best option is to estimate the total stellar mass by multiplying the star formation rate by the length of the age interval, assuming the star formation rate has been constant. Before proceeding, we note that several authors  \nFigure 8  \n<!-- image -->  \n∗  \n/circledot  \nObservational estimates of the bound mass fraction Γ. The top row shows values of Γ taken from exclusive catalogs, while the bottom row shows inclusive catalogs. In the left column we show Γ as a function of cluster age T , while in the right column we show it as a function of galaxy star formation rate ˙ Σ ∗ . Measurements of Γ in multiple age ranges for the same galaxy are connected by dashed lines. Faded points indicate measurements that are potentially dubious: either measurements of Γ at ages < 10 Myr where it is not feasible to separate bound from unbound systems, or measurements of Γ for ages > 10 Myr in starburst systems where the star formation rate may have fluctuated on this time scale. In the legend, ♣ , ♠ , and ♦ stand for SBS0335-052E, ESO338-IG04, and ESO185-IG13, respectively. References are as follows: G10 = Goddard, Bastian & Kennicutt (2010), A11 = Adamo, Ostlin & Zackrisson (2011), R14 = Ryon et al. (2014), A15 = Adamo et al. (2015), LL15 = Lim & Lee (2015), J16 = Johnson et al. (2016), C17 = Chandar et al. (2017), M18 = Messa et al. (2018b), G18 = Ginsburg & Kruijssen (2018).  \ndescribe Γ as the 'cluster formation efficiency' (e.g., Larsen & Richtler 2000; Bastian 2008; Goddard, Bastian & Kennicutt 2010; Kruijssen 2012), meaning the fraction of stars formed in bound clusters. This definition implicitly assumes that one one can cleanly distinguish bound from unbound at all cluster ages, and that α T ≈ 0 so that Γ is the same when measured for any age range. We have argued neither of these assumptions is strictly correct, but if one makes them, our more general definition of Γ reduces to theirs.  \nAs with α T , measurements of Γ in the literature span a wide range of values. We summarize the currently available set of measurements in Figure 8. These data must be examined with extreme caution, even in comparison to compilations of α M and α T . The most secure measurements, for steadily star-forming galaxies at ages above ≈ 10 Myr, generally show Γ ≈ 0 . 01 -0 . 1, consistent with the observation that the majority of early B stars are found outside clusters, as is most of the UV flux (Pellerin et al. 2007). The low value of Γ at ages glyph[greaterorsimilar] 10 Myr is sometimes referred to as 'infant mortality': essentially all stars are born in clusters (using our expansive definition of cluster that includes both bound and unbound systems), but a majority of the stars are not gravitationally bound, so by ages of a few tens of Myr, most stars are not members of clusters any more.  \nAs shown in Figure 8, a number of authors have reported measurements of Γ at younger ages and in galaxies undergoing starbursts, and have attempted to deduce trends in Γ with either age (e.g., Chandar  \net al. 2017) or star formation activity (e.g., Goddard, Bastian & Kennicutt 2010; Adamo, Ostlin & Zackrisson 2011). We regard such claims as dubious given the methodological uncertainties, and in Figure 8 we fade the points for which significant methodological concerns exist. With regard to age trends, as pointed out by Kruijssen & Bastian (2016), estimates of Γ from inclusive catalogs at ages below ≈ 10 Myr are likely to be heavily contaminated by the presence of non-bound structures that have simply not yet had time to disperse. Conversely, however, exclusive catalogs may be missing a substantial population of bound but not yet relaxed clusters at similar ages. Thus Γ estimates at ages below 10 Myr derived from inclusive and exclusive catalogs should be viewed as upper and lower limits, respectively. An exception is the point from Ginsburg & Kruijssen (2018) for the Sgr B2 region in the Milky Way's Central Molecular Zone, for which the authors check boundedness directly using stellar velocities estimates from the radio recombination lines of the ultracompact H ii regions around each massive star.  \nPossible trends in Γ with star formation activity are also potentially contaminated by bias. While Γ measurements in low-star formation rate (SFR) galaxies tend to be made at ages glyph[greaterorsimilar] 10 Myr to avoid the issues with inclusive versus exclusive catalogs, measurements for high-SFR galaxies are almost exclusively at younger ages. This is both because high SFR galaxies also tend to be merging systems, for which measuring Γ at older ages is problematic because the total SFR may be variable, and because high SFR galaxies are rare and thus tend to be distant, resulting in a magnitude limit that precludes measurements of Γ at ages glyph[greaterorsimilar] 10 Myr (Chandar et al. 2017). This means that there is a degeneracy between two possibilities: high values of Γ in starbursts could be because Γ is in fact larger at higher SFR, but it could equally well be that Γ is higher at younger ages. The Sgr B2 point from Ginsburg & Kruijssen (2018) does not suffer from this methodological problem, but also unfortunately does not help break the degeneracy either, because it has both high SFR and young age. It is unclear to which of these factors its high value of Γ should be attributed. Consequently, we conclude that there is not yet convincing evidence that Γ varies with either age or SFR.", 'IMPROVING PHOTOMETRIC DETERMINATIONS OF CLUSTER PROPERTIES': "The preceding sections should make clear that one of the significant uncertainties for CMFs, CAFs, and Γ are the difficulties of assigning masses and ages to clusters based on photometry. Progress will require an extensive effort 'ground-truthing' photometric property determinations against measurements either CMDs or spectroscopy. Such a systematic comparison has not been performed since the work of Elson & Fall (1985, 1988) for the Magellanic Clouds using ground-based data. Much more accurate work with space-based CMDs should now be possible in both the Clouds and M31 (Johnson et al. 2016). Similarly, spectroscopy can be used to break the degeneracy between ages of ≈ 5 -10 and ≈ 50 Myr from photometry; spectroscopic measurements of a small sample of ambiguous clusters would at a minimum provide a useful prior for Bayesian analysis methods that generally return bimodal posterior PDFs when applied to clusters in the ambiguous parts of color space (Fouesneau et al. 2014; Krumholz et al. 2015a). Unfortunately the improved infrared photometry provided by the James Webb Space Telescope is likely to be of limited use, since clusters' colors in the IR are nearly constant at ages glyph[greaterorsimilar] 6 Myr, when they become dominated by red giants and supergiants (Gazak et al. 2014).", '2.5. Size and density': "In addition to masses and ages, we can also study the physical sizes of star clusters, or equivalently their densities. The data on the distribution of cluster sizes are more limited than those available for mass or age, because measuring a size is obviously only possible if the object in question is at least marginally resolved. While the most physically-meaningful characterization of cluster size and density would be a half-mass radius or similar description of how mass is distributed, the limited resolution available for galaxies beyond the Milky Way and its closest neighbors generally precludes measuring anything other than the projected half-light radius, defined as the projected radius containing half the cluster light. Even this quantity is difficult to estimate for targets more distant than M31, because many clusters have relatively shallow light profiles where much of the light is in an extended halo (see § 4.1) that tends to be lost against the background. For resolved stellar populations, on the other hand, it is usually possible to extract multiple structural parameters from fits of stellar number counts to theoretically-motivated profiles; see § 4.1 and Section 1.3.2 of Portegies Zwart, McMillan & Gieles (2010) for a summary of the various size parameters that can be defined for star clusters. Here we focus on projected half-mass radii (or half-number radii if masses are not available) from the resolved samples, as these are the closest to what is available from the photometric size measurements.  \nTable 2 Summary of masses and radii in current cluster samples  \na The first number gives the sample median, while plus and minus values indicate the 16th to 84th percentile range; b Σ h is the surface density at the half-mass radius, Σ h = M/ 2 πr 2 h ; c References are as follows: K13 = Kharchenko et al. (2013), J12 = Johnson et al. (2012), F14 = Fouesneau et al. (2014), R15 = Ryon et al. (2015), R17 = Ryon et al. (2017), C16 = Chandar et al. (2016), B07 = Barmby et al. (2007), BH18 = Baumgardt & Hilker (2018), MG07 = McCrady & Graham (2007), L18 = Leroy et al. (2018).  \nIn Figure 9 we show the mass-radius relation for cluster samples in nearby galaxies, and we summarize the properties of the cluster samples shown in the Figure in Table 2; for the Milky Way, we have updated Portegies Zwart, McMillan & Gieles (2010)'s compilation of YMCs in Table 3 8 . The plot reveals a few notable features. First, there is a mass-radius relation, but it is weak in terms of both slope and tightness. Cluster radius increases with mass slightly more slowly than r h ∝ M 1 / 3 , so that more massive clusters have mildly higher average density. We refrain from giving a formal fit to the log r h -log M relation,  \n5  \nFigure 9  \n<!-- image -->  \n/circledot  \nMass-radius relation for star clusters in nearby galaxies. For Milky Way and globular clusters, we plot half-mass radii where available, half-number radii otherwise; for all other data sets, we plot half-light radii. Dashed black lines indicate loci of constant density within r h , with the mass density ρ = 3 M/ 8 πr 3 h (in units of M glyph[circledot] pc -3 ) and the corresponding number density n (units of H nuclei cm -3 ) as indicated; dotted lines show loci of constant surface density Σ = M/ 2 πr 2 h (units of M glyph[circledot] pc -2 ). Hexagonal density plots show the log of the density of clusters in the ( M,r h ) plane for Milky Way within 2 kpc of the Sun (data from Kharchenko et al. 2013, K13, with masses estimated from tidal radii using the method of Piskunov et al. 2007), clusters in M31 from PHAT (radii from Johnson et al. 2012 and masses from Fouesneau et al. 2014), and clusters in NGC 628, NGC 1313, and NGC 5236 (galaxies grouped together; Ryon et al. 2015, 2017, R15, R17). The low-mass edges in the extragalactic data are a result of observational limitations, not a physical truncation. Points show individual clusters from smaller samples: clusters in the disk of M51 (Chandar et al. 2016, C16), GCs in M31 (Barmby et al. 2007, B07) and the Milky Way (Baumgardt & Hilker 2018, BH18), SSCs in M82 (McCrady & Graham 2007, MG07) and NGC 253 (Leroy et al. 2018, L18), and YMCs from our own compilation (Table 3).  \nbecause the results would be highly sensitive to how we weighted the various data sets. Regardless of the index of the mass-radius relationship, it is clearly very broad, in the sense that, at any given mass, one can find clusters whose radii vary by a factor of 5; at least some of the spread is undoubtedly the result of measurement errors. There is no obvious difference in either the breadth or slope of the mass-radius relation between the inclusive and exclusive catalogs, although a comparison is not straightforward since, unlike for measurements of α T and Γ, there is little overlap in Figure 9 between measurements carried out using exclusive versus inclusive catalogs. Second, the distribution is bounded at both large and small radii. Cluster densities are bounded below by ρ ∼ 1 M glyph[circledot] pc -3 (or equivalently n ∼ 30 H cm -3 ); GCs in galaxy halos go to only slightly lower densities. On the small radius side, the data appear to be limited to the region Σ glyph[lessorsimilar] 10 5 M glyph[circledot] pc -2 , as noted previously by Hopkins et al. (2010). Third, there is no separation between the loci occupied by clusters in galaxy disks and GCs; the two distributions blend continuously into one another.  \n5  \n8  \nTable 3 Young massive ( M > 10 4 M glyph[circledot] ) clusters (YMCs) in the Milky Way  \n±  \na Error estimates are given when provided by the referenced source; b References are as follows: FMM99 = Figer, McLean & Morris (1999), F06 = Figer et al. (2006), D07 = Davies et al. (2007), HEM08 = Harayama, Eisenhauer & Martins (2008), B08 = Borissova et al. (2008), C09 = Clark et al. (2009), R10 = Rochau et al. (2010), N10 = Negueruela et al. (2010), N11 = Negueruela et al. (2011), G11 = Gennaro et al. (2011), GFN12 = Gonz'alez-Fern'andez & Negueruela (2012), B12 = Borissova et al. (2012), C12 = Clarkson et al. (2012), Z17 = Zeidler et al. (2017); c Cited reference does not list a radius; the value we quote is a by-eye fit based on their published images, and should only be considered accurate to the factor of ∼ 2 level.", '2.6. Elemental abundances': "In systems close enough to permit spectroscopy of individual low-mass stars, we can study the composition of clusters. The stars in all OCs in the Galaxy are observed to be highly homogeneous (∆[Fe/H] glyph[lessorsimilar] 0.05 dex) in their abundances of all elements, while no old GCs are fully homogeneous, at least not in their light elements (Li, C, N, O, F, Na, Mg, and/or Al). Essentially all massive GCs also show an anti-correlation in certain light elements (e.g., Na-O, Mg-Al; Carretta et al. 2010). Like multiple stellar populations observed in CMDs, this has become a distinguishing feature of globulars, including for globulars in the Magellanic Clouds, that is never observed in OCs (e.g., Bragaglia et al. 2017). This difference has led to the view that GC history is more complex than that of OCs (Gratton, Carretta & Bragaglia 2012). However, there is a remarkable exception to the rule: the 12-Gyr old, 10 4 . 8 M glyph[circledot] globular Ruprecht 106 appears to be a true 'single stellar population' with homogeneous abundances and no light element anti-correlations (Villanova et al. 2013); unusually, it is enriched in r-process and s-process elements but does not show the enhanced [ α /Fe] signature of the oldest stars, a characteristic shared by all other Galactic GCs. The globulars Terzan 7 and Palomar 12 also do not have multiple populations and, intriguingly, may be extragalactic through their association with the disrupting Sgr dwarf. Bragaglia et al. (2017, Table 9) catalog all modern attempts to detect multiple populations in open and globular clusters.  \nTypically, the heavy elements (Fe and beyond) in globulars exhibit much less scatter than for the light elements (e.g. Gratton, Carretta & Bragaglia 2012), except in a few systems (e.g., ω Cen, Terzan 5, M22, NGC 5824) that have unusual properties compared to the general family of clusters. Unusually, NGC 5824 has ∆[Fe/H] variations of up to 0.3 dex but it also has a remarkable stellar halo with no tidal truncation detected out to 400 pc in radius, far beyond other clusters (Da Costa, Held & Saviane 2014). There are numerous other oddities (e.g., Marino et al. 2018): NGC 2419 has a large scatter in K-Mg; some globulars show s-process variations (e.g., M22, ω Cen), others (e.g. M15) exhibit r-process variations. With globulars, there are exceptions to every rule. In one particularly deep study of the globular NGC 6752,  \nYong et al. (2013) demonstrate weak inhomogeneities in all elements, even the Fe peak. The existence of inhomogeneities in globulars may be related to the presence of multiple populations in these systems but there is no clear association at the present time (Bragaglia et al. 2014). Multiple stellar populations in globulars are simply not understood (Gratton, Carretta & Bragaglia 2012; Renzini 2013; Bastian & Lardo 2018).  \nAbundance studies in OCs are more limited. The young OC Hyades is nominally homogeneous (De Silva et al. 2006, 2011) although a deep study reveals weak variations at the level of 0.02 -0.03 dex (Liu et al. 2016b). The older OC M67 exhibits slightly larger variations (still glyph[lessorsimilar] 0 . 05 dex) from a differential analysis of two fortuitous solar twins (Liu et al. 2016a). To date, all OCs have been shown to be single stellar populations with a few rare counter claims (e.g., Geisler et al. 2012) refuted in follow-up analysis (e.g., Bragaglia et al. 2014). The Hyades and M67 clusters have been subject to the most sensitive searches for abundance variations, so it is plausible that the majority of OCs will be shown to have low-level inhomogeneities in future studies. Weak metallicity variations may reflect atomic diffusion in stellar atmospheres (Dotter et al. 2017; Souto et al. 2018), planetary infall (Mel'endez et al. 2017), or simply intrinsic variations that existed at the onset of star formation (Feng & Krumholz 2014). We return to this point in § 3.6.", '2.7. The limits of demography at young ages': 'Thus far we have discussed the demographics of star clusters without too much worry about the issues raised in § 1.2 about how one defines a cluster, though our discussion about the differences between inclusive and exclusive catalogs should already hint at potential problems. We now turn directly to the question: at what stellar population age does it make sense to treat as meaningful concepts such as the CMF, CAF, and all the other demographic distributions discussed above?  \nTo answer this question, it is helpful to conduct a thought experiment: consider an annulus within the Milky Way centered on the Solar Circle, and for every star of age T ± ∆ T , find its N th nearest neighbor (either in 3D or in projection) in the same age range, and plot the distribution of nearest neighbor distances. Qualitatively, what does this distribution look like? For stars at an age T ∼ 1 Gyr, the answer is fairly straightforward: almost all stars at this age are field stars, and after 1 Gyr the orbits of these stars through the Galaxy are thoroughly phase-mixed. Separations therefore follow a Poisson distribution, perhaps with some overall structure at ∼ kpc scales as a result of spiral arms. A small fraction of 1 Gyr-old stars are members of open clusters, and for these the neighbor distances are much smaller. Consequently, the overall separation distribution is a Poissonian with an excess at very small separation. Indeed, one sees exactly this sort of distribution if one plots the distribution of distances between molecular clouds and young star clusters in external galaxies (e.g., Kawamura et al. 2009). This makes it straightforward to define old clusters: they are the structures made up of the stars in the small-separation excess.  \nThe results are very different for pre-main sequence stars (selected based on optical colours, infrared spectral energy distributions, or other methods) with ages T glyph[lessorsimilar] 10 Myr. The angular correlation function for these stars is well described as a power-law ξ ( θ ) ∝ θ -p with p ≈ 0 . 5 -1 . 5 on scales from ∼ 0 . 01 pc to ∼ 30 pc (e.g., Hennekemper et al. 2008; Kraus & Hillenbrand 2008; Schmeja, Gouliermis & Klessen 2009; Bressert et al. 2010). The cluster-cluster angular correlation function is a power-law up to ∼ kpc scales (e.g., Elmegreen & Elmegreen 2001; Bastian et al. 2007; Gouliermis et al. 2017; Grasha et al. 2017a,b) see Gouliermis (2018) for a comprehensive review. For such distributions, there is no obvious break or scale that could be used to delineate between cluster and field populations. This presents a challenge for describing the demographics of star clusters, since it is not obvious how to go about defining the clusters to begin with.  \nThis situation leads to two main options. One is a Press-Schechter-like approach, whereby one considers the stellar field smoothed with some characteristic window that defines a size scale. However, there is an  \nimportant difference between our situation and the cosmological one: in cosmology one can compute from first principles the critical overdensity of ≈ 200 that separates bound and unbound structures, while no such calculation is available for star clusters (though there are attempts - see § 3.3). Alternately one can simply appeal to the empirical result that older clusters have a characteristic size ∼ 1 -10 pc, and adopt a similar averaging scale at younger ages. This is what many authors do implicitly by using clustering algorithms such as MST that use a free parameter to set the scale. A second approach is to abandon discussion of cluster demographics entirely for stellar populations glyph[lessorsimilar] 10 Myr old. One can instead discuss clustering -i.e., the non-Poissonian distribution of stellar separations, and the correlations between separation in space and separation in other quantities (age, velocity, elemental abundance, etc.) - rather than discuss the properties of discrete clusters, an approach previously suggested by Zinnecker (2010). This removes the need to impose a size scale on a distribution where none is immediately apparent.  \nAs stars age, their separation distribution must evolve from the power-law form seen at young ages to the Poisson plus small-separation excess form that prevails at older ages. Observations suggest that much of this evolution happens relatively quickly, over about a crossing time of the relevant size scale (e.g., Elmegreen 2000; Grasha et al. 2017b). Thus the type of demographic analysis we have discussed to this point is sensible as applied to structures more than about a crossing time old. However, from a theoretical standpoint the question that any successful theory of cluster formation must answer is how the power-law, scale-free distribution that prevails among young clusters turns into the bimodal distribution that prevails at older ages, and what determines the demographics of the clustered portion of that bimodal distribution. In other words, how does stellar clustering turn into star clusters ?', '3. BIRTH': "In this section we examine clusters' early life, when gas dominates the mass; § 4 covers the gas-free phase. Here we first review the gaseous initial conditions for star cluster formation ( § 3.1). We then discuss observational constraints and theoretical models for the rates at which these structures form stars ( § 3.2) the stellar feedback processes that inhibit this transformation ( § 3.3), and the star formation histories that result ( § 3.4). We conclude with an attempt to synthesise this material into a coherent model for the physical properties ( § 3.5) and elemental abundances ( § 3.6) of the final gas-free cluster population.", '3.1. Initial conditions for star cluster formation': "3.1.1. Giant molecular clouds. Clusters are born in molecular clouds, primarily giant molecular clouds (GMCs). There is considerable evidence that GMCs obey several relations discovered by Larson (1981): (1) The 1D velocity dispersion, σ , is supersonic and varies with size L as σ ∝ L p , where p ≈ 0 . 5 in the Galaxy (Solomon et al. 1987; Falgarone, Pety & Hily-Blant 2009); (2) GMCs have roughly equal kinetic and gravitational potential energies, indicating that they are bound or nearly so; (3) the column density N ≈ nL (or, equivalently, the mass surface density Σ = M/πR 2 ) is approximately constant within a given galaxy. As Larson pointed out, only two of these are independent. To see this, note that the degree of gravitational binding can be measured by the virial parameter, α vir ≡ 5 σ 2 R/GM , which is unity for a sphere of constant density in virial equilibrium with no surface pressure; such a sphere is bound if it has α vir < 2 (Bertoldi & McKee 1992, who also considered non-spherical clouds). This definition can be rewritten as σ = √ ( π/ 5) Gα vir Σ R . Thus bound clouds ( α vir ∼ 1) with the same surface density, Σ, automatically satisfy σ ∝ R 1 / 2 . In terms of the virial parameter, the ratio of the crossing time t cr = R/σ to the free-fall time is t cr /t ff ≈ 2 /α 1 / 2 vir .  \nObservations of molecular gas since Larson (1981)'s work have generally confirmed his findings. For the linewidth-size relation, Solomon et al. (1987) found σ = 0 . 72 R 0 . 5 0 km s -1 , where R 0 = R/ (1 pc), in their  \nsurvey of GMCs in the first Galactic quadrant. Falgarone, Pety & Hily-Blant (2009) compiled observations of molecular gas in the Milky Way and found σ ∝ R 0 . 5 within a factor 3 over the size range 0.1 to > 100 pc. Subsequent observations of Galactic (Rice et al. 2016; Miville-Deschˆenes, Murray & Lee 2017) and extragalactic (Bolatto et al. 2008; Faesi, Lada & Forbrich 2018) molecular gas are generally consistent with the Solomon et al. (1987) result. Faesi, Lada & Forbrich (2018) attribute the absence of a clear linewidthsize relation in some extragalactic observations to low spatial resolution. Finally, we note that the value p ≈ 0 . 5 is expected for supersonic turbulence on theoretical grounds and has been seen in many simulations (e.g., Padoan et al. 2014).  \nOne can determine if a cloud is gravitationally bound by comparing the 'turbulence parameter' C ≡ σ/R 1 / 2 to the surface density Σ (Keto & Myers 1986):  \nα vir = 3 . 70 C 2 Σ 2 , (2)  \nwhere C is in units of km s -1 pc -1 / 2 and Σ 2 = Σ / (100 M glyph[circledot] pc -2 ) (c.f. Bolatto et al. 2008). The results of Solomon et al. (1987) for GMCs in the inner Galaxy ( C = 0 . 72 km s -1 and Σ 2 = 1 . 7) imply α vir ≈ 1 . 1, and they concluded that these GMCs are bound. Heyer et al. (2009) and Roman-Duval et al. (2010) confirmed this, the former by plotting C vs. Σ from 13 CO observations of the inner Galaxy. Sun et al. (2018) were able to determine the virial parameter for the clouds in the different galaxies by measuring surface density and line width at a fixed linear scale. They found α vir on the scale R = 60 pc between 1.3 and 3.2 in all galaxies but M31 and M33 (where they were larger), showing that the clouds are bound or nearly bound. Visual inspection of their results shows that α vir is almost always close to unity for the clouds with the highest surface densities. For a cloud to be bound, it must also be able to avoid tidal disruption. For a flat rotation curve, this requires that the mean cloud density exceed twice the mean density of the galaxy inside the orbit of the cloud (Chernoff & Weinberg 1990), or ¯ ρ > 5 . 4 ( v c, 220 /R kpc ) 2 M glyph[circledot] pc -3 , where v c, 220 is the circular velocity in units of 220 km s -1 , and R kpc is the distance to the center of the galaxy in kpc.  \nLarson's third relation, the constancy of the surface density of GMCs in the inner Galaxy (excluding the Galactic Center), was confirmed by Solomon et al. (1987), who found 〈 Σ 〉 = 170 M glyph[circledot] pc -2 and by Roman-Duval et al. (2010), who concluded that 〈 Σ 〉 = 140 M glyph[circledot] pc -3 . (However, clouds with M < ∼ 10 4 M glyph[circledot] can have lower Σ - Heyer, Carpenter & Snell 2001 - and there can be systematic large-scale variations of the surface density within GMCs - Schneider et al. 2015.) Different galaxies have different surface densities: in a study of 15 disk galaxies, Sun et al. (2018) found that median surface density at 120 pc resolution ranged from about 10 to 200 M glyph[circledot] pc -2 excluding the Antennae, for which it is 2300 M glyph[circledot] pc -2 . The variation in Σ for the normal galaxies could reflect variations in the beam filling factor of GMCs as well as variations in the intrinsic properties of GMCs (Lada et al. 2013). The dispersion in these values within individual galaxies is about ± 0 . 4 dex.  \nThe surface density of a GMC is directly related to both its internal pressure and the pressure exerted on its surface, and insofar as Σ is constant within a galaxy, these pressures tend to be also. The turbulent pressure within a self-gravitating cloud is  \n¯ P turb ≡ ¯ ρσ 2 ≡ 3 π 20 α vir G Σ 2 GMC , (3)  \nfor spherical clouds (e.g., McKee & Tan 2003), so that ¯ P turb /k B = 1 . 0 × 10 5 α vir Σ 2 2 K cm -3 . This should exceed the external pressure, P ext , and indeed Hughes et al. (2013) found that ¯ P turb ∼ (1 -8) P ext in the eight galaxies they studied. For a typical GMC in the Milky Way, with Σ CO ∼ 140 M glyph[circledot] pc -2 , one finds that the surface pressure on the CO cloud is comparable to the mean turbulent pressure inside the cloud. On the other hand, for GMCs in high pressure environments, such as galactic nuclei or starbursts, the  \nsurface density must be large in order for the cloud to be self-gravitating: the virial theorem implies that α vir ∼ (1 -P ext / ¯ P ) -1 in the absence of strong magnetic fields (Bertoldi & McKee 1992), so the requirement that α vir ∼ 1 implies Σ > ∼ 2( P ext /G ) 1 / 2 . High pressures also imply high volume densities: the mean density in a cloud is  \n¯ ρ GMC = 3 4 ( π Σ 3 2 M 6 ) 1 / 2 M glyph[circledot] pc -3 = 179 α 3 / 4 vir M 1 / 2 6 ( ¯ P turb /k B 10 8 K cm -3 ) 3 / 4 M glyph[circledot] pc -3 (4)  \nwhere the final step follows because for a non-spherical cloud Σ ≈ M/V 2 / 3 has the value as for a spherical cloud of the same volume, V , and P turb > ∼ 2 P ext for bound clouds.  \n3.1.2. Clumps and clusters. Clusters are born in density concentrations (clumps) within GMCs. Both the GMCs and the clumps within them are supersonically turbulent. In the absence of self-gravity and for an approximately isothermal equation of state, this turbulence leads to a log-normal density distribution (e.g., McKee & Ostriker 2007) and a corresponding log-normal distribution of surface densities (Brunt, Federrath & Price 2010). Self-gravity leads to the formation of a power-law tail in the PDF for dense gas (Collins et al. 2012). This has been observed in the PDF of surface densities in a number of GMCs, dP/d log Σ ∝ Σ -(2 ± 0 . 7) for clouds with active star formation; this is often associated with an overall powerlaw density distribution in the cloud (Lombardi, Alves & Lada 2015; Schneider et al. 2015). The clumps in GMCs are formed by the combined action of supersonic turbulence and self-gravity. Equation 2 implies that the typical surface density of a bound clump within a bound GMC (each of which therefore has α vir ≈ 1) is about the same as that of the GMC, or Σ clump ∼ 100 M glyph[circledot] pc -2 in the Galaxy. Observational studies of clumps have been carried out only for the Galaxy. Analyzing the results of the APEX 870 µ m survey of the Galactic plane (excluding the central ± 5 · in longitude), which was complete for Σ > 700 M glyph[circledot] pc -2 and M > 1000 M glyph[circledot] within 20 kpc, Urquhart et al. (2018) found about 10 7 M glyph[circledot] in 8000 dense clumps. A by-eye fit to their data on the surface densities of the clumps gives dN clump /d log Σ ∝ Σ -1 . 6 over the range Σ = 600 -6000 M glyph[circledot] pc -2 ; the median surface density is about 700 M glyph[circledot] pc -2 and the maximum is about 3 × 10 4 M glyph[circledot] pc -2 . The clumps are gravitationally bound with α vir ∼ 0 . 1 -1; Kauffmann, Pillai & Goldsmith (2013) and Tan et al. (2013) have suggested that magnetic fields could provide support for sub-virial clumps. At high masses, the clump mass distribution has α M, clump ≈ -2 (indeed, Heithausen et al. 1998 found α M, clump = -1 . 84 over 5 decades in mass) and has an upper limit of about 10 5 M glyph[circledot] . In the Galaxy, the GMCs out of which the clumps form have a power-law mass with α M, GMC ≈ -1 . 6, so that most of the mass of GMCs is in massive clouds. This distribution is truncated at about 10 7 M glyph[circledot] (Rice et al. 2016 and references therein, but see Miville-Deschˆenes, Murray & Lee 2017 for a different view), as expected theoretically from the Jeans mass of the galactic disk (i.e., the Jeans mass expressed in terms of the gas surface density, Σ gal ; Kim & Ostriker 2001),  \nM J = σ 4 G 2 Σ gal = 10 7 ( σ 7 km s -1 ) 4 ( 13 M glyph[circledot] pc -2 Σ gal ) M glyph[circledot] . (5)  \nIf the Toomre Q -parameter, Q = κσ/ ( πG Σ) ∼ 1, where κ is the epicyclic frequency, this is about the maximum mass expected for a bound cloud, so in the Galaxy the most massive dense clump is about 1% of the mass of the most massive GMC.  \nObservations of GMCs in external galaxies show that in some, but not all, cases the mass distribution can be described by a truncated power law with most of the mass in the most massive clouds, as observed in the Galaxy: this describes the GMCs in NGC 300 (Faesi, Lada & Forbrich 2018) and those in the inner region and in the spiral density waves of M51, but not the GMCs in the rest of the galaxy, where α M, GMC < -2, so that low-mass clouds contain most of the mass (Colombo et al. 2014). In the regions  \nof M51 in which α M, GMC > -2, the upper limit on the mass distribution is similar to the Galactic value, ∼ 10 7 M glyph[circledot] .  \nThe conditions required to produce massive clusters, such as globular clusters, occur in interacting galaxies such as the Antennae (e.g., Whitmore et al. 2014), dwarf irregular galaxies, some of which are starbursts and many of which show signs of interaction (Billett, Hunter & Elmegreen 2002), and nuclear starbursts in spiral galaxies such as NGC 253 (e.g., Leroy et al. 2018). Such conditions are rare now, but were common at redshifts z > ∼ 2 (Kruijssen 2014). A detailed survey of the GMCs in the Antennae galaxies has yet to be made, but Wilson et al. (2000) carried out a low-resolution ( ∼ 300 × 500 pc) survey of this system and identified five large molecular cloud complexes with masses in excess of 10 8 M glyph[circledot] ; the largest one not associated with one of the galactic nuclei had a mass of 6 × 10 8 M glyph[circledot] . By contrast, Wilson et al. (2003) pointed out that the most massive molecular cloud complex in M51 had a mass about 10 times less. They suggested that the interaction of the galaxies produced regions of less shear, which enables the existence of more massive GMCs; Billett, Hunter & Elmegreen (2002) emphasized the importance of low shear in producing more massive clusters in dwarf irregular galaxies. For disk galaxies, we can express the Jeans mass (Equation 5) in terms of the Toomre Q -parameter, so that  \nM J = π 4 G 2 Σ 3 Q 4 κ 4 = 2500 Σ 3 2 Q 4 1 . 5 Ω 4 0 M glyph[circledot] , (6)  \nwhere in the second expression a flat rotation curve has been assumed, Q has been normalized to a typical critical value, and Ω 0 is measured in units of Myr -1 (Krumholz & McKee 2005). Hence, if Q self-regulates to be of order unity, a reduction in the epicyclic frequency, κ , by a factor 2 due to the effect of the interaction can increase the maximum GMC mass by more than an order of magnitude. More massive GMCs enable the production of more massive clumps (Harris & Pudritz 1994; Reina-Campos & Kruijssen 2017); Johnson et al. (2015a) have observed a clump in the Antennae with a mass of at least 5 × 10 6 M glyph[circledot] and a radius less than 24 pc, corresponding to Σ > ∼ 3000 M glyph[circledot] pc -2 and ¯ P turb /k B > ∼ 10 8 K cm -3 .  \nThe maximum possible mass of a cloud that forms by gravitational instability in a disk is M GMC , max glyph[similarequal] ( λ T / 2) 2 Σ, where λ T = 4 π 2 G Σ /κ 2 is the Toomre length (Escala & Larson 2008). Disks are stable for Q > ∼ 1 . 5 (Kim & Ostriker 2001), and the Jeans mass (eq. 6) at Q = 1 . 5 is close to M GMC , max . If the disk has a total mass (including stars and dark matter) inside a radius R of M tot ( < R ), a gas mass inside R of M g ( < R ) glyph[similarequal] πR 2 Σ and it is rotationally supported, then  \nM GMC , max M g ( < R ) glyph[similarequal] 4 π ( κ/ Ω) 4 [ M g ( < R ) M tot ( < R ) ] 2 : (7)  \nthe maximum fraction of the gas that can go into a single GMC scales as the square of the gas fraction. As a result, more massive GMCs, and therefore more massive cluster-forming clumps, are expected in gas-rich regions such as galaxies in the process of formation and in starbursts (Escala & Larson 2008).  \nGalaxy interactions also produce shocks that can compress clouds. A collision at a relative velocity v rel produces radiative shocks with a velocity v rel / 2 and a pressure P = ρ 0 ( v rel / 2) 2 , corresponding to P/k B = 3 . 8 × 10 6 n H , 0 ( v rel / 300 km s -1 ) 2 . This leads to localized regions of very high pressure when clouds collide, P/k B ∼ 10 8 K cm -3 (Jog & Solomon 1992), but only ∼ 10 6 K cm -3 for shocks in the intercloud medium, which is not that much larger than the typical pressure in Galactic GMCs.  \nBased on this summary of the observed properties of molecular gas in galaxies, we have the following expectations for the gaseous precursors of star clusters currently forming in the Galaxy: clusters form from clumps with masses ranging from very small values (depending on the star formation efficiency) up to about 10 5 M glyph[circledot] and surface densities from ∼ 100 M glyph[circledot] pc -2 up to about 3 × 10 4 M glyph[circledot] pc -2 . A substantial fraction of the clouds are bound ( α vir < ∼ 2). The typical clump has a column density comparable to that of the  \nt ff : free-fall time, the natural evolutionary timescale for a self-gravitating system  \nglyph[epsilon1] ff : the fraction of a cloud's mass that is transformed into stars per cloud free-fall time  \nglyph[epsilon1] ∗ : the fraction of a cloud's initial mass that is transformed into stars by the time all the initial gas has been been consumed or ejected; sometimes called the star formation efficiency, though we will mostly eschew this term to minimize confusion with glyph[epsilon1] ff  \nt sf : timescale over which conversion to stars occurs: glyph[epsilon1] ∗ ≈ glyph[epsilon1] ff ( t sf /t ff )  \nη : mass loading factor: ratio of mass ejected by star formation feedback to mass converted to stars  \ncloud in which it is embedded. For a clump mass of 10 3 M glyph[circledot] and a typical GMC pressure corresponding to Σ GMC = 140 M glyph[circledot] pc -2 , namely P turb /k B ∼ 2 × 10 5 Kcm -3 , the density is about 50 M glyph[circledot] pc -3 . Clumps have a distribution of surface densities extending to higher values, and the mean density of a clump at a given mass scales as Σ 3 / 2 clump (Equation 4). Conditions in nearby spiral galaxies do not differ substantially from those in the Milky Way. These conditions are consistent with the production of the OCs observed in the local universe, but not of GCs, which are far more massive than the clumps observed in the Galaxy. Molecular clouds that can produce globulars are seen in interacting galaxies like the Antennae, dwarf starbursts such as He 2-10 (Johnson et al. 2018), and nuclear starbursts (Leroy et al. 2018), and are expected in any disk with a high gas fraction (as observed in galaxies at high redshift - Genzel et al. 2011) or a small epicyclic frequency (as observed in interacting and dwarf irregular galaxies).", '3.2. Conversion of gas to stars': "The collapse of molecular clouds eventually produces stars. Here we review only those aspects of this process most relevant to star clusters; for a more general treatment see the reviews by McKee & Ostriker (2007) and Krumholz (2014). The CMF and the production of bound clusters depend strongly on the efficiency with which gas is converted into stars. Efficiency can be defined in multiple ways, which we can illustrate in an idealized example: a cloud of initial mass M and free-fall time t ff forms stars at a rate ˙ M ∗ . At any given time, the remaining gas mass is M g and the mass of stars formed is M ∗ , and we can parametrize the relationship between the star formation rate, current gas mass, and free-fall time as  \n˙ M ∗ = glyph[epsilon1] ff M g t ff , (8)  \nwhere glyph[epsilon1] ff is one of the efficiencies we will define. A region in free-fall collapse with nothing inhibiting star formation has glyph[epsilon1] ff ≈ 1. In addition to gas consumption by star formation, feedback from stars launches a wind at a rate ˙ M w , which we normalize to the star formation rate by defining the mass loading factor η = ˙ M w / ˙ M ∗ . With these definitions, we have  \n˙ M g = -(1 + η ) glyph[epsilon1] ff M g t ff . (9)  \nAll the factors that appear in this equation η , glyph[epsilon1] ff , or t ff - can depend on M g , M ∗ , time, or any number of other variables. Moreover, real clouds are not closed boxes with a fixed initial mass; there is certainly mass inflow occurring simultaneously with star formation (see 'Conveyor-Belt Model for Cluster Formation' in the Supplementary Materials). Indeed, nothing in our formulation even requires that clouds be bound M g may can include both bound and unbound material. For the purposes of illustration, however, we can ignore these complications and treat η , glyph[epsilon1] ff , and t ff as constant, in which case it is trivial to write down expressions for the instantaneous gas and stellar masses:  \nM g M = e -t/t sf M ∗ M = 1 -e -t/t sf 1 + η ≡ glyph[epsilon1] ∗ ( 1 -e -t/t sf ) t sf ≡ glyph[epsilon1] ∗ glyph[epsilon1] ff t ff . (10)  \nThese equations provide another definition of the star formation efficiency: glyph[epsilon1] ∗ , the fraction of the initial cloud mass that has been transformed into stars once all the gas is consumed. In this idealized problem, glyph[epsilon1] ∗ = 1 / (1 + η ). The term t sf we have introduced defines the characteristic timescale for star formation, i.e., M ∗ /M reaches its final value of glyph[epsilon1] ∗ on a timescale t sf . In this section we discuss glyph[epsilon1] ff , and in the next two sections we examine glyph[epsilon1] ∗ and t sf .  \n3.2.1. Observational constraints on glyph[epsilon1] ff . The value of glyph[epsilon1] ff is well constrained by observations. Formally, for any specified region of volume V containing a gas mass M g , we have glyph[epsilon1] ff = ˙ M ∗ / [ M g /t ff (¯ ρ )], where ˙ M ∗ is the instantaneous star formation rate within V , t ff (¯ ρ ) = √ 3 π/ 32 G ¯ ρ is the free-fall time as a function of density, and ¯ ρ = M g /V is the mean density of the region in question. There are several methods available to estimate these quantities. The most direct is to define a column density or extinction threshold, estimate the mass within that threshold, calculate the density and thus the free-fall time from the mass and projected area (assuming the unseen third dimension is comparable in size to the two observed dimensions), and estimate the star formation rate by counting young stellar objects (YSOs) and estimating their masses and the duration of the YSO phase (Krumholz, Dekel & McKee 2012; Federrath 2013; Evans, Heiderman & Vutisalchavakul 2014; Salim, Federrath & Kewley 2015; Heyer et al. 2016; Ochsendorf et al. 2017; also see Heiderman et al. 2010; Lada et al. 2013). All published studies using this method rely on data from Spitzer , which was only sensitive enough to detect ∼ 1 M glyph[circledot] YSOs within a few kpc of the Sun, and massive YSOs out to the distance of the Magellanic Clouds. Consequently, this method is only available out to the Magellanic Clouds, and studies that go beyond the Solar neighborhood must make a correction for the unseen part of the IMF.  \nA second approach, usable throughout the Milky Way and its satellites, is to match catalogs of starforming regions identified by tracers such as infrared or free-free emission with catalogs of molecular clouds identified by CO or dust emission, matching them up if they are sufficiently close in position or velocity, and then using the mass, free-fall time, and star formation rates of the matched clouds and star-forming regions to estimate glyph[epsilon1] ff (Vutisalchavakul, Evans & Heyer 2016; Lee, Miville-Deschˆenes & Murray 2016; Ochsendorf et al. 2017). A third approach, available for extragalactic systems with extensive molecular gas and star formation tracer maps, is simply to pixelize the entire galaxy, and estimate masses, densities, and free-fall times in each pixel (Krumholz, Dekel & McKee 2012; Leroy et al. 2017; Utomo et al. 2018). A fourth method is to observe tracers of dense gas, most commonly HCN, and correlate these with tracers of star formation (e.g., Krumholz & Tan 2007; Garc'ıa-Burillo et al. 2012; Hopkins et al. 2013; Usero et al. 2015; Gallagher et al. 2018). For this method there is no need to spatially resolve the emission, because the molecule itself determines the density and free-fall time; e.g., Onus, Krumholz & Federrath (2018) show that the luminosity-weighted mean density of HCN line-emitting gas is ≈ 10 4 cm -3 . 9  \nWe summarize recent observational constraints on glyph[epsilon1] ff in Figure 10. Comparing the various methods, we see that essentially all studies based on YSO counting or HCN give glyph[epsilon1] ff ≈ 0 . 01, with a study-to-study dispersion of ≈ 0 . 3 dex, and a dispersion of about 0 . 3 -0 . 5 dex within any single study. This dispersion is probably an upper limit, since it includes both the true physical dispersion and the effects of measurement errors. Pixel statistics produce a similar dispersion, but with a median 0 . 3 -0 . 5 dex lower. However, these results are still consistent because there are systematic uncertainties in all the methods at the ∼ 0 . 5 dex level. For pixel statistics one must assume a value for the poorly-constrained line of sight depth through the target galaxy. For YSO counting, one must estimate the duration of the phase during which newly-formed stars would be classified as YSOs, the volume density of the gas seen only in projection, and the IMF correction. For HCN, there is uncertainty in the mean density of the emitting gas and the emissivity per unit mass, which is needed to estimate the total mass. 10  \nThe cloud matching studies disagree with the other three methods, yielding much larger dispersions  \nFigure 10  \n<!-- image -->  \nCompilation of glyph[epsilon1] ff measurements, grouped by method. Each histogram shows the distribution of glyph[epsilon1] ff values measured in a single study. Points and error bars show the median and the 16th to 84th percentile; numerical values are printed in each panel. In panels with two histograms, the colored one is the distribution including non-detections plotted at their 2 σ upper limits, while the gray histogram shows the distribution omitting non-detections; median and percentile values are for the colored distribution. The vertical dashed lines at log glyph[epsilon1] ff = -2 are to help guide the eye. References are as follows: W10 = Wu et al. (2010), L13 = Lada et al. (2013), E14 = Evans, Heiderman & Vutisalchavakul (2014), U15 = Usero et al. (2015), H16 = Heyer et al. (2016), L16 = Lee, Miville-Deschˆenes & Murray (2016), V16 = Vutisalchavakul, Evans & Heyer (2016), S16 = Stephens et al. (2016), L17 = Leroy et al. (2017), O17 = Ochsendorf et al. (2017, appears twice, because this study used two independent methods), G18 = Gallagher et al. (2018), U18 = Utomo et al. (2018). All values of glyph[epsilon1] ff are as reported by the authors of the study, except that we have derived glyph[epsilon1] ff for the Lada et al. (2013) sample following Krumholz (2014), and we have homogenized glyph[epsilon1] ff from the HCN data sets to the calibration of Onus, Krumholz & Federrath (2018).  \nand mostly higher medians. The main difference between cloud matching and the other methods is that cloud matching does not require that the star formation and gas tracers be co-spatial, which makes the results sensitive how one constructs and matches catalogs of clouds and star-forming regions; the 0 . 8 dex difference in median glyph[epsilon1] ff values derived by Vutisalchavakul, Evans & Heyer (2016) and Lee, Miville-Deschˆenes & Murray (2016) is almost entirely due to this sensitivity, since the regions they target and the underlying data sets they use are nearly identical. To gain more insight into why cloud matching differs from other methods, it is helpful to examine the study of Ochsendorf et al. (2017), who use both YSO counting and cloud matching to derive glyph[epsilon1] ff values for the same molecular gas map of the LMC. 11 For YSO counting, they adopt the standard approach of counting YSOs within the same contours used to define clouds, while for cloud matching they decompose an H α map into H ii regions and match clouds and H ii regions if there is any spatial overlap. We show the distribution of glyph[epsilon1] ff values produced by the two methods in the left panel  \nFigure 11  \n<!-- image -->  \n<!-- image -->  \nLeft: cumulative distribution of glyph[epsilon1] ff values measured in the study of Ochsendorf et al. (2017) by YSO counting versus cloud matching, both using the Jameson et al. (2016) H 2 map. The star shows the 30 Doradus region. For the YSO counting method, the blue band indicates the uncertainty range associated with non-detections, with the lower value corresponding to assuming that clouds without YSO detections have star formation rates at the sensitivity limit, and the upper value corresponding to assuming that such regions have star formation rates of zero. Right: zoom-in on the 30 Doradus region, reproduced from Ochsendorf et al. (2017) by permission of the AAS. Grayscale shows the H 2 column density, black contours show CO clouds from the catalog of Wong et al. (2011), red triangles show the positions of massive YSOs, with size indicating inferred mass (Ochsendorf et al. 2016), and blue contours show H ii regions.", 'of Figure 11.': 'Clearly cloud matching yields a broader distribution with a higher median. To understand why, we can look at the example of the 30 Doradus region (indicated by the star in the left panel of Figure 11) for which YSO counting gives glyph[epsilon1] ff = 0 . 065, while cloud matching gives glyph[epsilon1] ff = 0 . 38, a factor of 6 larger. Since the same GMC properties are used in both estimates, the difference arises solely from the imputed star formation rate. Examining the right panel of Figure 11, it is clear how a substantial difference could arise: while the great majority of the YSOs are within the footprint of the cloud, there is significant H α emission that is not, but that partly overlaps it. Treating all this emission as if it comes from the cloud causes the cloud matching method to infer a much larger star formation rate. This emission could indeed be associated with the existing molecular cloud, but it could also be associated with molecular gas that has been dispersed or displaced by stellar feedback.  \nSimilar effects are also apparent in other galaxies. For example, in a survey of NGC 628 at ∼ 50 pc resolution, Kreckel et al. (2018) identify ≈ 1500 H ii regions and ≈ 750 GMCs, but find < 100 overlaps between the two. Naive application of a cloud matching method to this data set would lead one to conclude that most GMCs are inert and have glyph[epsilon1] ff = 0, while most H ii regions formed with glyph[epsilon1] ff glyph[greatermuch] 1, an obviously unphysical result. The data are better explained by the hypothesis that bright H ii regions rapidly disperse or displace the molecular clouds in which they form. In principle it should be possible to account for this effect and extract corrected estimates of glyph[epsilon1] ff from the statistics of the GMC-H ii region correlation (Kruijssen et al. 2018), but this method has not yet been widely applied.  \nThis analysis suggests that cloud matching-based estimates of glyph[epsilon1] ff should be treated with caution. Given', 'THE FUTURE OF YSO COUNTING': "YSO counting is the most reliable method of estimating glyph[epsilon1] ff , but at present it is limited by the sensitivity and angular resolution of Spitzer , the telescope used for all existing YSO counting studies. With Spitzer , the only YSOs we can reliably detect beyond the Milky Way have masses glyph[greaterorsimilar] 8 M glyph[circledot] (Ochsendorf et al. 2016). This sensitivity limit is responsible for the large uncertainty in the shape of the low end of glyph[epsilon1] ff distribution shown in Figure 11. However, the James Webb Space Telescope will be far more sensitive and have significantly better angular resolution. With JWST it will possible to detect ≈ 1 M glyph[circledot] YSOs in the Magellanic Clouds, and to obtain YSO catalogs comparable to those currently available for the LMC and SMC out to M31 at least, possibly farther. This will greatly improve cloud-scale measurements of glyph[epsilon1] ff beyond the Milky Way.  \nthe reasonable consistency in results produced by the other methods, we conclude that the preponderance of the current observational evidence favors glyph[epsilon1] ff ≈ 0 . 01 for regions glyph[greaterorsimilar] 1 pc in size, with a dispersion of glyph[lessorsimilar] 0 . 5 dex and a comparable level of systematic uncertainty.  \n3.2.2. Theory of glyph[epsilon1] ff . The observed value of glyph[epsilon1] ff is surprisingly low, and a number of authors have proposed theoretical models aimed at explaining it. The earliest models relied on magnetic regulation of star formation. Interstellar gas is magnetized, and if the field strong enough it can prevent gas from collapsing until non-ideal magnetohydrodynamic effects reduce the magnetic flux threading the gas. However, more recent measurements of magnetic fields have generally shown that they are too weak to support the gas see Crutcher (2012) for a review. Consequently, attention has focused on three possibilities: glyph[epsilon1] ff may be low because clouds are not bound, because they are turbulent, or due to feedback.  \nUnbound cloud models propose that glyph[epsilon1] ff is low because most of the material observationally-defined as molecular clouds is not in fact self-gravitating (e.g., Dobbs, Burkert & Pringle 2011; Meidt et al. 2018). We have argued in § 3.1.1 that the bulk of the evidence does not support this view, but if it were the case that only ∼ 1% of the mass in a molecular cloud were bound, this would naturally explain why glyph[epsilon1] ff ≈ 1%. While this is a potentially-viable explanation for the low glyph[epsilon1] ff values of GMC as traced by CO, most of the observations compiled in § 3.2.1 are for the much denser gas traced by HCN, or by the YSO counting studies, which is almost certainly bound (e.g., Kauffmann, Pillai & Goldsmith 2013). The unbound cloud hypothesis does not explain why this gas also shows low glyph[epsilon1] ff .  \nTurbulent regulation models propose that glyph[epsilon1] ff is low because star-forming gas is turbulent enough to render most sub-regions of a molecular cloud unbound, even if that cloud is bound on its largest scales. Qualitatively the argument stems from the linewidth-size relation discussed in § 3.1: since velocity dispersion varies with size as σ ∝ R 1 / 2 , then the kinetic energy per unit mass contained within a region of size R scales as e K ∝ R , while the binding energy per unit mass of a region of mass M scales as e G ∝ GM/R . If one chooses a region close to the mean density then M ∝ R 3 and e G ∝ R 2 . Consequently, for a randomly-chosen sub-region of a cloud, the virial ratio obeys α vir ∝ 1 /R , and thus most regions smaller than an entire cloud are unbound. This argument can be made quantitative by integrating over the density PDF. The earliest version of such an argument appeared in Krumholz & McKee (2005), but this has been superseded by numerous later works that model the density PDF and the evolution of its self-gravitating parts with increasing accuracy (Padoan & Nordlund 2011; Hennebelle & Chabrier 2011; Padoan, Haugbølle & Nordlund 2012; Federrath & Klessen 2012; Hopkins 2013; Burkhart 2018; Burkhart & Mocz 2018).  \nThe consensus of recent models is that turbulence does substantially reduce glyph[epsilon1] ff , but not all the way to glyph[epsilon1] ff ≈ 0 . 01 as required by the observations. As time goes by the density PDF of a self-gravitating cloud  \ndeviates develops an increasingly-prominent power-law tail on its high-density end that causes glyph[epsilon1] ff to rise with time (Murray & Chang 2015; Burkhart 2018). This density build-up is most likely counteracted by localized feedback processes that break up high-density regions and suppress the growth of the power-law tail; the most likely mechanism for this is protostellar outflows, since stars begin to launch outflows as soon as they form, and even low-mass stars can produce significant outflow feedback. Modern simulations that include both protostellar outflows and magnetic fields 12 that help couple them to gas are able to achieve values of glyph[epsilon1] ff that, while still a factor of 2 -3 too high compared to the consensus observational value, are approaching agreement with the data (Wang et al. 2010; Myers et al. 2014; Federrath 2015; Cunningham et al. 2018). An alternative possibility is that cloud disruption by feedback is so fast and efficient that typical clouds never have time to develop power-law tails extensive enough to drive glyph[epsilon1] ff to large values (Grudi'c et al. 2018).  \nA final caveat on theoretical explanations for the low value of glyph[epsilon1] ff is that turbulent regulation is effective only if the gas is actually turbulent. There is significant debate in the literature about whether the observed large linewidths of molecular clouds should be interpreted as turbulence, or whether they might represent coherent gravitational collapse (e.g., Heitsch, Ballesteros-Paredes & Hartmann 2009; Ballesteros-Paredes et al. 2011; Zamora-Avil'es, V'azquez-Semadeni & Col'ın 2012; Traficante et al. 2018). While there is significant evidence from both analytic theory (Klessen & Hennebelle 2010; Goldbaum et al. 2011) and numerical simulations (Robertson & Goldreich 2012; Lee & Hennebelle 2016a,b; Birnboim, Federrath & Krumholz 2018) that accretion flows invariably drive turbulence, and recent kinematic measurements with Gaia also appear to disfavor the idea that GMCs are in global collapse ( § 3.4.2), the question is not settled.", '3.3. Feedback and the termination of star formation': "The final star formation efficiency, glyph[epsilon1] ∗ , depends both on glyph[epsilon1] ff and on the ability of stellar feedback to eject gas from nascent star clusters. Unfortunately we cannot measure glyph[epsilon1] ∗ for individual clouds directly, because for a single region we cannot measure both the gas mass at the onset of star formation and the stellar mass at its conclusion. One can attempt to measure the mass-loading factor η , but this is technically very challenging (see Yang et al. 2018 for a recent attempt). The best prospect for measuring glyph[epsilon1] ∗ is likely the use of statistical methods to analyze populations of clouds (Kruijssen et al. 2018), but thus far measurements using this technique are not widely available. For these reasons, estimates of glyph[epsilon1] ∗ come primarily from theory. In the remainder of this section we follow the approach taken by several previous authors (e.g., Matzner 2002; Fall, Krumholz & Matzner 2010; Matzner & Jumper 2015; Rahner et al. 2017) by considering a variety of feedback mechanisms and attempting to determine under what circumstances they become relevant. We summarize these findings in Figure 12.  \n3.3.1. Protostellar outflows. Bally (2016) provides a comprehensive recent review of protostellar outflows, so we focus only on the details most relevant for star cluster formation. Outflows are critical to breaking up dense regions and thus keeping glyph[epsilon1] ff small because they eject ∼ 2 / 3 of the mass from individual protostellar cores (e.g., Matzner & McKee 2000; Offner & Chaban 2017). However, this mass is ejected at relatively low velocities, and thus may not escape the larger-scale protocluster. For this reason, theoretical models (e.g., Matzner & McKee 2000; Matzner & Jumper 2015), numerical simulations (e.g., Wang et al. 2010; Krumholz, Klein & McKee 2012; Murray, Goyal & Chang 2018), and observations (e.g., Nakamura & Li 2014; Li et al. 2015; Plunkett et al. 2015) all suggest that outflows have a limited role in setting final  \nFigure 12  \n<!-- image -->  \n/circledot  \nStar cluster mass-radius relation, showing the same data as Figure 9, with shaded regions indicating where various feedback mechanisms are potentially significant. In constructing this diagram, we have for simplicity considered the case of star-forming clouds where the gaseous and stellar components have comparable masses and radii, but we emphasize that in reality this need not be the case. The feedback mechanisms shown are outflows ( § 3.3.1), supernovae (SNe; § 3.3.6), direct radiation pressure (DR; § 3.3.3), ionization ( § 3.3.2), and indirect radiation pressure (IR; § 3.3.4). Shaded regions fade at masses below which stochastic sampling of the IMF makes them unlikely to be active ( § 3.3.1). For all feedback types except IR, the region where the mechanism is effective is to the left of and above the line.  \nstar formation efficiencies in most clusters. There has yet to be a comprehensive numerical survey of the parameter space, but Matzner & Jumper (2015) estimate analytically that outflows produce glyph[epsilon1] ∗ glyph[lessorsimilar] 0 . 5 only in clusters with escape speeds v esc glyph[lessorsimilar] 1 km s -1 . We show this limit in Figure 12.  \nDespite this limitation, outflows play a unique role, because they do not depend on the presence of massive stars. The momentum per unit mass of stars formed delivered by protostellar outflows to their surroundings is of order the escape speed from a protostellar surface (e.g., Matzner & McKee 2000). Lowmass and high-mass protostars have similar surface escape speeds, so low-mass stars are as effective at providing outflow feedback as massive ones. By contrast, for all the other mechanisms we discuss below, feedback is dominated by massive stars unlikely to be found in small clusters. To quantify this point, we use the SLUG stochastic stellar populations code (da Silva, Fumagalli & Krumholz 2012; Krumholz et al. 2015b) to measure the PDFs for the number of supernovae, and bolometric and ionizing luminosity, in clusters at a range of masses for a Chabrier (2005) IMF. We find that the number of expected supernovae exceeds unity only for clusters with masses glyph[greaterorsimilar] 100 M glyph[circledot] , and that the median (bolometric, ionizing) luminosity of clusters is < 50% of the value for a fully-sampled IMF in clusters with mass glyph[lessorsimilar] (400 , 700) M glyph[circledot] . Consequently, protostellar outflows are likely to be the only feedback mechanism that limits glyph[epsilon1] ∗ for stellar populations smaller than a few hundred M glyph[circledot] . This limit is very blurry due to stochasticity; some clusters with masses of only a few hundred M glyph[circledot] nevertheless have substantial ionizing luminosities (Andrews et al. 2014).  \n3.3.2. Photoionization feedback. A population of stars that fully samples the IMF produces ionizing photons at a rate of Ψ i ≈ 6 × 10 46 s -1 M -1 glyph[circledot] (Murray & Rahman 2010). These can heat gas in star-forming clouds to ≈ 10 4 K, giving it a sound speed of ≈ 10 km s -1 . If not trapped either by the surrounding material or by gravity, this gas will flow outward in a wind known as a champagne flow, directly removing mass from the cloud. The back-reaction created when the ionized gas pushes off the surrounding neutral material, known as the rocket effect, may eject even more mass and exert significant forces, ultimately unbinding clouds entirely. The mass removed in this process can be considerable - the analytic calculations of Williams & McKee (1997) and Matzner (2002) suggest that a 10 4 M glyph[circledot] population of stars can eject ≈ 10 times its own mass over a few Myr, giving glyph[epsilon1] ∗ ∼ 0 . 1. Numerical simulations generally confirm this result, though they show that the value of glyph[epsilon1] ∗ produced by ionization also depends strongly on the initial boundedness of clouds, and on whether or not they are magnetized (e.g., Dale, Ercolano & Bonnell 2012, 2013; Howard, Pudritz & Harris 2017; Geen et al. 2015, 2016; Gavagnin et al. 2017; Dale 2017; Kim, Kim & Ostriker 2018).  \nThe main limitation on photoionization as a feedback mechanism is that the ionized gas must be able to escape, but there is some tension in the simulation literature about how much of a limitation this represents. Dale, Ercolano & Bonnell (2012, 2013) find that photoionization is unable to unbind more than a few percent of the mass for clouds with v esc > 10 km s -1 , while Kim, Kim & Ostriker (2018) find glyph[epsilon1] ∗ glyph[lessorsimilar] 0 . 5 at all escape speeds glyph[lessorsimilar] 20 km s -1 , though it is possible that this is because they include direct radiation pressure as well as photoionization (see § 3.3.3). None of the other published studies explore a range of v esc values. It is plausible that even regions with v esc > 10 km s -1 could lose mass due to photoionization, since freshly-ionized gas rocketing off a dense neutral surface will generally accelerate to small multiples of the ionized gas sound speed. In the analogous problem of photoionized gas driving winds off accretion disks, simulations show that even regions with escape speeds ≈ 30 km s -1 launch significant winds (Woods et al. 1996). Thus escape speed required to suppress the effects of photoionization is of order 10 km s -1 , but its exact value is uncertain; in Figure 12 we adopt 15 km s -1 .  \n3.3.3. Direct radiation pressure. A zero-age population of stars that fully samples the IMF has a light to mass ratio Ψ ≈ 1100 L glyph[circledot] M -1 glyph[circledot] (e.g., Fall, Krumholz & Matzner 2010). The light is emitted mostly at UV wavelengths where the ISM is highly opaque, so that even a modest dust column is sufficient to ensure that essentially all the radiation momentum is deposited in the gas. Whether this affects glyph[epsilon1] ∗ depends on how the associated radiation force compares to pressure and gravitational forces. H ii region gas pressures vary with radius as r -3 / 2 due to ionization balance, while direct radiation pressure varies as r -2 , so gas pressure dominates once H ii regions reach a characteristic size r ch ≈ 0 . 03 S 49 pc, where S 49 is the ionizing luminosity in units of 10 49 s -1 (Krumholz & Matzner 2009). For a stellar population of mass M ∗ with ionizing luminosity S = Ψ i M ∗ , one can show that this radius is larger than the size of the cluster, and thus radiation pressure rather than photoionization dominates mass removal, only if the escape speed v esc glyph[greaterorsimilar] 10 km s -1 . Thus direct radiation pressure becomes dominant roughly where photoionized gas pressure ceases to be effective.  \nThe importance of direct radiation pressure relative to gravity depends on the cloud column density (Fall, Krumholz & Matzner 2010; Murray, Quataert & Thompson 2010). Consider a spherical gas cloud of total (gas plus stellar) mass M and radius R , with a stellar population of mass glyph[epsilon1] ∗ M forming at its center. The column-averaged outward force per unit gas mass supplied by radiation is Ψ glyph[epsilon1] ∗ M/ 4 πR 2 (1 -glyph[epsilon1] ∗ )Σ c , where Σ = M/πR 2 is the mass per unit area. The corresponding column-averaged inward force per unit mass from gravity is ≈ GM/R 2 . Thus the Eddington ratio, which defines the ratio of radiative to gravitational  \nforce, is f Edd = Σ DR glyph[epsilon1] ∗ / (1 -glyph[epsilon1] ∗ ), where  \nΣ DR = Ψ 4 πGc ≈ 340 M glyph[circledot] pc -2 (11)  \nis the surface density below which radiation direct pressure becomes important. Naively we would expect direct radiation pressure to expel gas and thus limit glyph[epsilon1] ∗ only for f Edd glyph[greaterorsimilar] 1, requiring clouds with Σ glyph[lessorsimilar] Σ DR . However, Thompson & Krumholz (2016) point out that in a turbulent medium there will be sightlines with Σ < Σ DR even if Σ > Σ DR on average. Gas can be ejected on these low-column sightlines, leading to a wind. Their models suggest that, in a cloud with glyph[epsilon1] ff ≈ 0 . 01, this wind will remove mass fast enough to produce glyph[epsilon1] ∗ ≈ 0 . 5 at Σ ≈ 10Σ DR ; we show this line in Figure 12.  \nSimulations of the effects of direct radiation pressure have produced mixed results. All simulators find that direct radiation pressure is ineffective at reducing glyph[epsilon1] ∗ at surface densities glyph[greaterorsimilar] 10Σ DR (e.g., Grudi'c et al. 2018), but at the surface densities of Σ ≈ 100 M glyph[circledot] pc -2 characteristic of Milky Way GMCs some simulators find that direct radiation pressure yields glyph[epsilon1] ∗ glyph[lessorsimilar] 0 . 1 (Grudi'c et al. 2018), while others find substantially larger values of ≈ 0 . 3 -0 . 5 (Raskutti, Ostriker & Skinner 2016; Kim, Kim & Ostriker 2018; Howard, Pudritz & Harris 2018). These disagreements may be due to differences in the initial conditions, or to resolution problems in the simulations that can lead to a significant overestimate in the effectiveness of direct radiation pressure in some numerical methods (Krumholz 2018).  \n3.3.4. Indirect radiation pressure. When interstellar dust absorbs radiation, it re-emits the energy at IR wavelengths. IR dust opacities are a factor of ∼ 100 smaller than UV ones, but regions of sufficiently high column density may nonetheless be opaque to the re-radiated IR. In this case the IR photons will be absorbed again, and can deposit additional momentum; in clouds opaque enough to require many cycles of absorption and re-emission before the energy escapes, the resulting force will greatly exceed the direct radiation force (Thompson, Quataert & Murray 2005; Murray, Quataert & Thompson 2010; Thompson et al. 2015).  \nIndirect (i.e., dust-reprocessed) radiation pressure feedback is governed by two processes: frequency diffusion and radiation Rayleigh-Taylor (RRT) instability. Frequency diffusion is important because IR wavelengths are much larger than interstellar grain radii, so dust opacity varies with frequency as κ ∝ ν 2 . Thus as stellar radiation diffuses outward through a dusty envelope and mean photon frequency decreases, dust-radiation coupling weakens. Eventually photons shift to frequencies low enough to escape. This effect greatly limits the importance of dust-reprocessed radiation (Wolfire & Cassinelli 1987; Reissl et al. 2018), and means that the mass-weighted mean opacity of a cloud is a function of the total gas column. Many subgrid models of radiation pressure ignore this effect and simply assume a constant IR opacity (e.g. Rosdahl & Teyssier 2015; Tsang & Milosavljevi'c 2018; Hopkins et al. 2018); given the importance of diffusion, their results should be treated with caution.  \nRRT instability (Jacquet & Krumholz 2011) causes gas that is being accelerated out of a gravity well by indirect radiation to become clumpy, reducing the effectiveness of its coupling to radiation. Early numerical work suggested this would prevent gas ejection entirely (Krumholz & Thompson 2012), but this result proved to be sensitive to the details of the radiation transfer method. Simulations based on flux-limited diffusion (Krumholz & Thompson 2012) or M1 closures (Rosdahl & Teyssier 2015; Skinner & Ostriker 2015) find that RRT prevents gas ejection unless the dust opacity is greatly boosted compared to observed Milky Way values, while those based on variable Eddington tensor (Davis et al. 2014) or implicit Monte Carlo (Tsang & Milosavljevi'c 2015) find that it does not, though it does greatly reduce the ejection velocity. The latter two methods are likely more reliable, so we expect that indirect radiation pressure can limit glyph[epsilon1] ∗ .  \nUnder what conditions will this happen? Krumholz & Thompson (2012) show, and simulations by Davis et al. (2014) confirm, that for any column of gas confined by gravity, there is a maximum radiation flux  \nthat can pass through it without triggering gas ejection. Calculations properly accounting for frequency diffusion show that a gas cloud forming stars at its center will exceed the critical flux before the star formation efficiency glyph[epsilon1] ∗ reaches 0.5 only if its surface density exceeds (Crocker et al. 2018a,b)  \nΣ ∗ glyph[greaterorsimilar] Σ ∗ , IR = 16 √ πGσ SB /c Ψ ( κ 0 T 2 0 ) -1 ≈ 7 × 10 4 ( κ 0 0 . 03 cm 2 g -1 ) -1 M glyph[circledot] pc -2 , (12)  \nwhere T 0 = 10 K and κ 0 is the opacity per unit gas mass to radiation with a color temperature T 0 ; the normalization used above is typical of Milky Way dust. We show this condition in Figure 12. Note that indirect radiation pressure has a minimum surface density at which it becomes effective, while all the feedback mechanisms discussed previously have a maximum surface density (for direct radiation) or escape velocity (for ionization) below which they are effective. The reason for this difference is the frequency diffusion effect: for a fixed radiation flux, a larger gas column traps heat more effectively and thus the amount of indirect radiation force delivered rises superlinearly with the gas column. Consequently, indirect radiation force becomes more effective relative to gravity as the total column density increases.  \n3.3.5. Hot star winds. Stars with surface temperature glyph[greaterorsimilar] 25,000 K and metallicity glyph[greaterorsimilar] 0 . 5 Z glyph[circledot] drive winds with speeds of several thousand km s -1 (e.g., Vink, de Koter & Lamers 2001), producing temperatures glyph[greaterorsimilar] 10 7 K when the wind shocks against the ISM. At this temperature the radiative cooling time is long compared to the dynamical time (unlike for protostellar outflows -§ 3.3.1 - or winds from intermediate-mass stars, Offner & Arce 2015), leading to the formation of hot bubbles that can push on surrounding colder gas and potentially eject it from a cloud. There are a number of analytic models describing this process in one dimension without cooling (see Bisnovatyi-Kogan & Silich 1995 for a review) and with cooling (Koo & McKee 1992), and several authors have published 1D numerical results including cooling (e.g., Krause et al. 2016; Fierlinger et al. 2016; Rahner et al. 2017; Silich & Tenorio-Tagle 2018; Naiman, Ramirez-Ruiz & Lin 2018). The general conclusion from these models is that winds will be the dominant feedback mechanism unless their effects are reduced by two phenomena that cannot be properly captured in 1D: hydrodynamic leakage of hot gas through low-density channels, and sapping of hot gas thermal energy via turbulent mixing at the hot-cold interface. The general conclusion from these models is that winds will be the dominant feedback mechanism unless their effects are reduced by two phenomena that cannot be properly captured in 1D: hydrodynamic leakage of hot gas through low-density channels, and sapping of hot gas thermal energy via turbulent mixing at the hot-cold interface.  \nTo constrain the importance of leakage and mixing, and thus the effectiveness of hot star winds at reducing glyph[epsilon1] ∗ , there are two approaches available: multi-dimensional simulations and observations. Published simulations find that, as a result of efficient leakage (Rogers & Pittard 2013; Calura et al. 2015; Wareing, Pittard & Falle 2017) or mixing-induced losses (Mackey et al. 2015), wind feedback lowers glyph[epsilon1] ∗ by at most tens of percent. However, there has yet to be a comprehensive parameter study. In addition, only Wareing, Pittard & Falle (2017) consider magnetized clouds, and the fields they assume are weak. It is possible that more realistic, stronger magnetic fields might reduce mixing losses (Gentry et al. 2018). Moreover, winds may be important not so much because they remove mass, but because they halt ongoing accretion of low-density gas onto molecular clouds (e.g., Gatto et al. 2017; Haid et al. 2018).  \nX-ray observations probe wind feedback because the pressure in hot gas is directly related to its Xray luminosity. The observed luminosities of diffuse gas around YMCs is rule out the possibility that the wind is trapped (Dunne et al. 2003; Townsley et al. 2003, 2006, 2011), and require significant leakage or mixing-induced radiative loss (Harper-Clark & Murray 2009; Rosen et al. 2014; Rogers & Pittard 2014; Toal'a & Arthur 2018). Indirect diagnostics based on IR line ratios are consistent with this conclusion (Yeh & Matzner 2012). Comparing the hot gas pressure in H ii regions to the pressure exerted by warm ionized  \ngas and radiation, Lopez et al. (2011, 2014) find that hot gas is subdominant in the massive 30 Doradus region, and in a wide range of smaller H ii regions in the LMC and SMC. Pellegrini, Baldwin & Ferland (2011) reach the opposite conclusion for 30 Doradus, but only under the assumption that the hot gas has a small volume filling factor, which Lopez et al. (2014) point out would make it locally significant but unimportant as an agent of gas clearing on large scales. Thus there is currently no observational evidence that hot star winds are important in setting glyph[epsilon1] ∗ . However, as with the simulations, observations cover only a narrow range of parameter space, and it is possible that future work might identify regimes where winds are important.  \n3.3.6. Supernovae. Supernovae (SNe) deliver about the same energy as hot star winds, but all at once rather than gradually over several Myr. Whether this makes them more effective is uncertain. Some simulators find that SNe effectively destroy molecular clouds (Rogers & Pittard 2013; Calura et al. 2015; Wareing, Pittard & Falle 2017), while others find that their effects are modest compared to other feedback mechanisms (Kortgen et al. 2016; Geen et al. 2016; Rey-Raposo et al. 2017). Simulation results appear to be sensitive to both resolution and numerical method (Gentry et al. 2017, 2018), and there has yet to be a systematic survey of parameter space. Moreover, analytic models suggest that photoionization feedback prior to the first SN strongly affects the final outcome (McKee, van Buren & Lazareff 1984; Matzner 2002), but of the published studies, only Geen et al. (2016) include this effect.  \nOne important limitation on the importance of SN feedback is time delay: even the most massive stars require t SN ≈ 3 Myr to explode, and sufficiently dense regions may convert a significant fraction of their gas to stars before the first SN occurs (Fall, Krumholz & Matzner 2010; Kruijssen 2012). A region with free-fall time t ff and star formation rate efficiency glyph[epsilon1] ff and no feedback mechanisms except SNe will convert a fraction glyph[epsilon1] ∗ = glyph[epsilon1] ff t SN /t ff of its mass to stars before the first SN. This yields glyph[epsilon1] ∗ < 0 . 5 only if the density obeys  \nρ < ρ SN ≈ 3 π 128 Gglyph[epsilon1] 2 ff t 2 SN ≈ 5 × 10 5 µ H ( glyph[epsilon1] ff 0 . 01 × t SN 4 Myr ) -2 cm -3 , (13)  \nwhere µ H = 2 . 3 × 10 -24 g is the mean mass per H nucleus for the usual cosmic composition of 27% He by mass. We show this criterion in Figure 12. Moreover, SNe may be unimportant even in regions that obey Equation 13 because other mechanisms might remove gas on a timescale below t SN . There is substantial observational evidence that this is the case in at least some regions: clusters with ages ≈ t SN are observed to be gas free both in the Milky Way (Longmore et al. 2014) and in M83 (Hollyhead et al. 2015), and the spectra of clusters in the process of gas clearing frequently show Wolf-Rayet features, indicating an evolutionary phase that precedes the first SN (Sokal et al. 2016). These observations strongly disfavor SNe as an important mechanism for gas removal. However, much like winds, the main effect of SNe may not be to clear dense gas, but rather to remove low-density gas and thereby prevent star formation from restarting after gas clearing.  \n3.3.7. Summary of feedback mechanisms and expectations on glyph[epsilon1] ∗ . Figure 12 provides an overall picture of stellar feedback. For the smallest clusters, M glyph[lessorsimilar] 100 M glyph[circledot] , protostellar outflows are the only mechanism of interest because the massive stars responsible for other feedback mechanisms are likely to be absent. Moving to slightly higher masses, in the great majority of the parameter space occupied by observed clusters in the disks of modern galaxies, and even for a substantial portion of the globular cluster population, SNe, photoionization, and direct radiation pressure are all at least potentially effective. When all three are capable of operating, both simulations and observations suggest that photoionization will dominate because it has no time delay, and photoionized gas pressure exceeds radiation pressure except in very compact star clusters with large escape speeds. Most star formation in this region should be characterized by efficiencies  \nglyph[epsilon1] ∗ glyph[lessorsimilar] 0 . 3, since simulations suggest that photoionization feedback is quite effective when it operates. This is also in accord with Lada & Lada (2003)'s observation that, for embedded clusters with masses glyph[greaterorsimilar] 10 3 M glyph[circledot] , the stellar mass fraction seems to increase with age, but only to a maximum of ≈ 30%.  \nAs we move to yet higher masses and densities, characteristic of globular clusters, YMCs, and SSCs, these three feedback mechanisms begin to become ineffective (see also Grudi'c et al. 2018). Ionization likely becomes ineffective first, followed by direct radiation pressure and SNe; this ordering is dictated simply by the mechanisms' dependence on mass and radius, as shown in Figure 12. Presumably glyph[epsilon1] ∗ increases in this region. Indirect radiation, unlike all the other feedback mechanisms, only becomes effective once the surface density rises above a threshold value. The observed cluster mass-radius distribution shows a maximum surface density (first pointed out by Hopkins et al. 2010) consistent with this limit, suggesting that indirect radiation pressure either lowers glyph[epsilon1] ∗ at surface densities above the threshold, or causes cloud radii to expand during star formation to drive them below it (Crocker et al. 2018a).  \nInterestingly there is a locus, represented by the unshaded triangle in Figure 12, where no known feedback mechanism is expected to be effective. This locus is occupied by a small population of massive, dense star clusters, generally classified as SSCs. Our analysis suggests that these clusters may form with very high glyph[epsilon1] ∗ values, simply because no feedback mechanism is capable of ejecting much gas from them. However, we caution that it is probably not realistic to assume that protoclusters spend their entire formation history on this locus; as we discuss below, real protoclusters almost certainly grow by gas accretion and form stars simultaneously, and thus will move in Figure 12.", '3.4. Star formation histories and timescales': "3.4.1. Observational constraints. The final ingredient in our simple model of star cluster formation is its duration t sf , or, more generally, the star formation history. Age-dating stars is a fraught problem (Soderblom et al. 2014), so the conclusions we draw must be regarded as somewhat tentative, though the situation has improved considerably in the last few years.  \nBy far the most comprehensively-studied region is the Orion complex. The young stars in Orion are spread over several tens of pc, but molecular gas and ongoing star formation are limited to a smaller ∼ 10 pc region around the ONC, which is ≈ 1 pc in size (Hillenbrand & Hartmann 1998 find r h ≈ 0 . 8 pc) and is centered on θ 1 Ori C, the most massive star in the complex. The ONC is the densest region, and, despite its small volume, represents a significant fraction of the total number of stars in the complex. We illustrate one recent measurement of the age distribution in Orion in Figure 13, which suggests three important conclusions that appear to apply to other star-forming regions as well (e.g., Tan, Krumholz & McKee 2006; Azimlu, Mart'ınez-Galarza & Muench 2015; Getman et al. 2018; Caldwell & Chang 2018). First, the densest regions show age spreads that are substantially larger than their dynamical times. Using the ages shown in Figure 13, the time over which 90% of the ONC stars formed is ≈ 6 Myr (other studies give similar results: Reggiani et al. 2011; Da Rio et al. 2016; Beccari et al. 2017), while the free-fall time is ≈ 0 . 6 Myr (Tan, Krumholz & McKee 2006; Da Rio, Tan & Jaehnig 2014). Consequently, star formation has been ongoing for ≈ 10 free-fall times. While the absolute ages and the magnitude of the age spread depend on the choice of age indicator, the conclusion that there is a substantial age spread is robust against this choice, and against possible contamination due to photometric variability (Messina, Parihar & Distefano 2017) or variable accretion histories (see Sidebar). This conclusion is also bolstered by kinematic observations showing that the ONC is very close to virial equilibrium (Kim et al. 2019), implying that it must have had time to relax dynamically.  \nSecond, while star formation has an extended history, the star formation rate has not been constant over this period; star formation has accelerated (Palla & Stahler 2000; Caldwell & Chang 2018). (In Figure 13, recall that the time axis is logarithmic, so constant star formation rate corresponds to a line of slope unity  \nFigure 13  \n<!-- image -->  \nDistribution of stellar ages in the Orion star-forming complex. In the main panel, points are stars color-coded by age as estimated by Kounkel et al. (2018). The orange and blue contours, respectively, show the density of stars on the sky associated with the components Orion D and Orion A identified in Kounkel et al.'s phase-space analysis; contours are placed at densities of 10%, 50%, and 90% of maximum density. The inset shows a zoom-in on the ONC, defined here as a 1 pc-radius region centered on θ 1 C (shown by the star). The flanking histograms show the age distributions for Orion D, Orion A, and the ONC, as indicated. In the ONC panel, we also show t ff , the free-fall time in the central pc of the ONC as estimated by Da Rio, Tan & Jaehnig (2014), and t 90 ≈ 10 t ff , the time over which 90% of the stellar population formed. The ages we show are those estimated based on spectroscopy (Kounkel et al.'s age HR ) where available, and based on color (Kounkel et al.'s age CMD ) elsewhere. However, in the ONC inset and histogram we show only spectroscopic ages, since the color-based ones are unreliable in regions of high extinction.  \nrather than zero.) Using the ages shown in Figure 13, we find that the time over which 90% of the stars formed is ≈ 10 t ff , but roughly 50% of stars formed in the last ≈ 3 t ff .  \nThird, mean age varies with position. The stars in the ONC are significantly younger than the mean", 'FALSE AGE SPREADS FOR ACCRETING YSOS?': "Some authors have suggested that the dispersion in YSO luminosity at fixed effective temperature, usually interpreted as an age spread, might instead be due to variation in YSO accretion histories (e.g., Baraffe, Chabrier & Gallardo 2009; Baraffe, Vorobyov & Chabrier 2012). If some YSOs were to undergo exclusively 'cold' accretion, meaning that accreting material radiates away all its entropy at the accretion shock, they would descend the Hayashi track more rapidly than YSOs with less efficient shocks, leading to a luminosity spread that could be mistaken for an age spread. However, the current consensus is that observed stellar age spreads are not due to this effect. Although both cold accretion and Kelvin-Helmholtz contraction drive stars down the Hayashi track, they do so at very different rates as a function of stellar mass. As a result, rather than generating a mass-independent age spread, cold accretion would make it appear that stars with effective temperatures glyph[greaterorsimilar] 3500 K are systematically older than cooler stars (Hosokawa, Offner & Krumholz 2011). For these higher T eff stars, cold accretion would also produce a small population that would be extremely luminous and thus appear young (Vorobyov et al. 2017; Jensen & Haugbølle 2018), rendering the overall apparent age distribution bimodal. Neither a bimodal inferred age distribution for hotter stars nor a systematic difference in inferred age as a function of effective temperature are observed, which puts strong limits on how much of the apparent age spread is due to cold accretion.  \nof the larger Orion A region, which is in turn younger than the still-larger Orion D region. Consequently, there is a gradient of increasing mean age and age spread as one moves from denser to less dense regions (Getman, Feigelson & Kuhn 2014; Getman et al. 2018). However, age does not increase as quickly as free-fall time, so while the ONC is ≈ 10 free-fall times old, the lower density regions around it in Orion A are only ∼ 1 -3 free-fall times old (Da Rio, Tan & Jaehnig 2014; Jaehnig, Da Rio & Tan 2015). Thus the outskirts of clusters are older than their centers in an absolute sense, but are dynamically younger. This statement likely applies to Orion D as well, but we cannot be sure because the gas has already been cleared, so the stellar density only provides an upper limit on the free-fall time. 13  \n3.4.2. Proposed scenarios. For an isolated star-forming cloud, the natural evolutionary timescale is t ff , and simulations of clouds that do not include an external environment (or effectively mock one up by using periodic boundary conditions) generally find that all gas is either converted to stars or expelled on such timescales (e.g., Krumholz, Klein & McKee 2012; Kim, Kim & Ostriker 2018; Grudi'c et al. 2018, though for an exception see Wang et al. 2010). Thus the observed spread of ≈ 10 t ff in the ONC is inconsistent with expectations from isolated cloud models. 14 Observed star formation histories also present a second challenge: even with t sf ≈ 10 t ff as in the ONC, since glyph[epsilon1] ff ≈ 0 . 01 (c.f. § 3.2) the ONC should reach a maximum conversion efficiency glyph[epsilon1] ∗ ≈ 0 . 1, probably too low to yield a gravitationally-bound cluster (c.f. § 3.5). It is possible that the ONC will in fact not remain bound (Kroupa, Aarseth & Hurley 2001), though recent Gaia  \nobservations suggest that it will (Kuhn et al. 2018; Kim et al. 2019). Moreover, the general point holds regardless of the fate of the ONC: the observed low median value of glyph[epsilon1] ff means that bound clusters are likely to arise only from regions where glyph[epsilon1] ff is in the upper ∼ 20% of the distribution (plausible, since only ∼ 10% of star formation results in bound clusters), where the free-fall time is short enough for the star formation history to extend over > 10 free-fall times (Kruijssen 2012), or some combination of the two.  \nOne possible explanation for the ONC's star formation history is that it did form in a single free-fall time, but in an extended region of lower density than that the stars currently occupy. This region was in a free-fall collapse that halted when massive stars dispersed the gas. In this scenario, t sf glyph[greatermuch] t ff not because star formation was extended in time, but because t ff decreased after the stars formed (e.g., Zamora-Avil'es & V'azquez-Semadeni 2014; Kuznetsova, Hartmann & Ballesteros-Paredes 2015, 2018; V'azquez-Semadeni, Gonz'alez-Samaniego & Col'ın 2017). This would naturally give rise to an age gradient such as that observed, because stars that formed earlier in the process would retain a higher proportion of their kinetic energy, and thus end up further from the cluster center in the final configuration (Getman et al. 2018).  \nHowever, this scenario encounters serious observational difficulties. First, it would produce a distinct kinematic signature that is not observed. During collapse stellar motions should be radially inward, and once collapse is complete and gas has been expelled, the ∼ 90% of the stars that do not remain bound (since Γ ∼ 10%; § 2.4) should be moving radially outward from the dense center; in either case stellar velocities should point radially toward or away from the densest point (in the example of Figure 13, the ONC). However, Gaia shows that this is not the case in Orion or in other complexes (Ward & Kruijssen 2018; Kuhn et al. 2018; Kounkel et al. 2018). Note that this radial expansion is not seen in the simulations of Kuznetsova, Hartmann & Ballesteros-Paredes (2015, 2018), but this is because the simulations do not include any form of feedback, and thus nearly all the stars they produce remain bound. If the stars in their final configuration were to become unbound somehow, the velocities would become largely radial as soon as the stars expanded significantly compared to their original volume. A second problem is that this scenario would require even more extreme efficiencies, glyph[epsilon1] ff glyph[greaterorsimilar] 0 . 3, to produce bound systems; excluding the discrepant cloud matching results, observed glyph[epsilon1] ff values do not go this high (c.f. § 3.2). A third problem is that this scenario requires star clusters to pass through a dense phase wherein the mass is still largely in the form of gas, so that for any given star cluster, there should be a comparably dense gas cloud. However, surveys have consistently failed to find gas clouds as dense as the densest YMCs (Ginsburg et al. 2012; Longmore et al. 2014; Walker et al. 2016; Urquhart et al. 2018). We illustrate this with a current compilation of gas clouds and clusters in Figure 14. Notice that, for the densest clusters in the disk of the Galaxy (though interestingly not for the Galactic Center), there are no clouds dense and massive enough to be their progenitors if the entire cloud had to be assembled at once. Given these problems, it is unlikely that the extended formation times are illusory.  \nHow then are the long formation times to be explained? One possibility is that existing simulations do not yet have enough physics or enough resolution to produce realistic lifetimes. There would be some precedent for this: for decades simulations struggled to reproduce the small observed values of glyph[epsilon1] ff , but, as noted in § 3.2, simulations that include multiple feedback processes, initial turbulence, and magnetic fields are now approaching this goal. It is possible that the same will happen with the duration of star formation. More likely, however, the problem with current simulations is lack of the larger scale-environment that provides an ongoing mass supply. Observed protoclusters are generally located at the confluence of gaseous filaments, suggesting a 'conveyor belt' picture in which gas accretion onto a protocluster and star formation occur simultaneously (Longmore et al. 2014; Motte, Bontemps & Louvet 2018; see Supplementary Materials). Inflow supplies energy to maintain turbulence and prevent contraction (Klessen & Hennebelle 2010; Goldbaum et al. 2011; Matzner & Jumper 2015; Lee & Hennebelle 2016a,b), which together with local feedback processes keeps glyph[epsilon1] ff low. This state persists as long as the accretion lasts, which can be many times  \nFigure 14  \n<!-- image -->  \nM  \n[  \nM  \n/circledot  \n]  \nComparison of masses and radii of young massive clusters (YMCs) in the disk of the Milky Way and near the Galactic Center, as indicated, and for clumps of molecular gas. The YMCs plotted are those listed in Table 3. For gas clouds in the Galactic disk, we show the sample of Urquhart et al. (2018, U18) from ATLASGAL, which is complete except near the Galactic Center. We supplement this with the Galactic Center clouds reported by Longmore et al. (2012, L12) and Immer et al. (2012, I12). The dashed and dot-dashed lines marked for each cluster show the mass as a function of radius for clouds with density as a function of radius ρ ∝ r -2 and r -1 . 5 , as indicated. These lines indicate the loci below which progenitors clouds for these clusters would have to lie if their densities varied as the indicated power-law, i.e., a gas cloud with density versus radius ρ ∝ r -2 would only be dense enough to form the indicated star cluster if it were to lie below the dashed line emerging from that star cluster. For several disk clusters, there are no gas clouds that satisfy this condition.  \nthe local value of t ff because the accretion flow is dictated by processes occurring on much larger scales with much longer dynamical times. Because glyph[epsilon1] ff is small and gas consumption slow, the star-forming gas mass grows with time, leading to an accelerating star formation rate. Once accretion ceases, and the turbulent energy supply and confining ram pressure of the accretion are removed, any remaining gas is converted to stars or expelled by feedback in a time of order t ff , as happens in the isolated cloud simulations. The final star formation efficiency is ultimately be dictated by the efficiency of feedback at ejecting gas, and if the ratio of mass ejected by feedback to mass converted into stars is glyph[lessorsimilar] 1, then the final value of glyph[epsilon1] ∗ is large even if glyph[epsilon1] ff is always fairly small.  \nThis scenario is consistent with observed values of glyph[epsilon1] ff , naturally produces an accelerating rate of star formation, and also yields a cluster mass-size relation in good agreement with observations (Matzner & Jumper 2015; Lee & Hennebelle 2016b). Because the final mass that goes into a cluster was never assembled in a single cloud all at once, there is no problem with finding clouds dense enough to be the progenitors of the densest clusters. It is not clear if this mechanism is capable of producing the core-halo age gradient observed by Getman, Feigelson & Kuhn (2014) and Getman et al. (2018), and visible in Figure 13 - for the ONC this gradient is ∼ 0 . 5 Myr pc -1 . Such a gradient could arise as a result of earlier-forming stars having longer to undergo dynamical heating in the core of the forming cluster, but this possibility has not been examined quantitatively. A second possible challenge for this picture is that preliminary work suggests that it leads to clusters with significant rotational support Lee & Hennebelle (2016a,b). While some clusters do rotate (e.g. Davies et al. 2011; H'enault-Brunet et al. 2012; Kamann et al. 2018), others do not (Kuhn et al. 2018), and it is not clear if the predicted rotation rate distribution is consistent with observations. However, current theoretical predictions for the rotation rate distribution do not yet take into account the  \npossibility of braking by dynamically-significant magnetic fields.", '3.5. Emergence of the gas-free cluster population': "3.5.1. The initial cluster mass function. The initial cluster mass function (ICMF), defined as the CMF immediately after gas removal, is the starting point for understanding the evolution of clusters. The ICMF is less well defined than its stellar counterpart because the stellar content of a cluster is constantly evolving, particularly for unbound clusters. Adams & Myers (2001) showed that a cluster must contain glyph[greaterorsimilar] 100 stars in order for the dynamical time to exceed the formation time, which they took to be ∼ 1 -2 Myr. Indeed, the ICMF for embedded clusters in the Solar neighborhood turns down below N ∼ 100 stars, corresponding to a cluster mass ∼ 50 M glyph[circledot] (Lada & Lada 2003). As discussed in § 2.2, the CMF for observed clusters usually can be fit with a Schechter function with a power-law index α M ≈ -2 and a truncation mass M c that is related to the environment. Elmegreen (2006) has shown that a slope α M = -2 naturally accounts for the fact that the stellar IMF averaged over a galaxy is very close to the average IMF of individual star clusters.  \n3.5.1.1. The clump mass function. Theoretically, the ICMF is the product of the clump mass function ( § 3.1.2) and the mass-dependent star formation efficiency glyph[epsilon1] ∗ ( M ). This is approximate: in a turbulent medium the clump mass evolves with time, and we have seen in § 3.4 that bound clusters form over multiple dynamical times with continuous accretion, so that the star-forming gas is never assembled into a single monolithic clump. Nonetheless, clumps still provide a useful approximation for modeling the complex process of cluster formation, just as protostellar cores provide a useful basis for modeling the stellar IMF (e.g., Hennebelle & Chabrier 2008). An important factor determining the clump mass function, dN clump /dM ∝ M α M, clump , is that molecular clouds are supersonically turbulent and appear to be selfsimilar from the outer scale of the turbulence (generally encompassing the entire cloud) to the sonic scale ( ∼ 0 . 1 pc in Galactic GMCs). Fleck (1996) was the first to argue that α M, clump = -2 for a self-similar medium in which clumps form hierarchically. Elmegreen & Falgarone (1996) came to a similar conclusion by using the observed fractal dimension of molecular clouds, and for the first time suggested that this would lead to a similar slope for the ICMF. Guszejnov, Hopkins & Grudi'c (2018) have elaborated the argument for α M, clump = -2, which we paraphrase: assume that the clouds form hierarchically, either by fragmentation or by mergers, and that each level of the hierarchy, n , contains N n clumps that have a range of masses ∆ M n ∝ M n . Then the total mass at level n is  \nM tot ( M n ) = M n N n = M n ( N n ∆ M n )( ∆ M n M n ) M n ∝ M 2 n dN dM , (14)  \nwhere the last step follows for large n . If the hierarchy is populated throughout, so that the mass at each level of the hierarchy, M tot ( M n ), is about constant, then α M, clump ≈ -2. This argument is clearly approximate; for example, GMCs have most of their mass at high masses in the Milky Way and at low masses in some other galaxies ( § 3.1.2). In a dynamical system such as a GMC, the clumps at different levels of the hierarchy can also have different lifetimes, thereby changing the observed PDF; for example, if the clumps have mass-independent surface densities and survive for a time proportional to the free-fall time, the observed value of α M, clump is larger than the intrinsic value by 1/4 (Fall, Krumholz & Matzner 2010). This theory of the clump mass function is thus not very precise, but it is broadly consistent with observation.  \n3.5.1.2. From clumps to clusters. As noted above, the ICMF is related to the mass distribution of the natal molecular clumps by the star formation efficiency (SFE), glyph[epsilon1] ∗ ( M ), so that dN ( M ∗ ) ∝ glyph[epsilon1] ∗ ( M ) dN clump ( M ) approximately, where it is possible that more than one cluster forms in a cloud. Photoionization and direct  \nradiation pressure are the most broadly effective feedback mechanisms ( § 3.3). Photoionization removes a mass ∆ M ion ∝ S 4 / 7 t 9 / 7 /ρ 1 / 7 over a time t , where S is the rate of production of ionizing photons (Whitworth 1979). If this is the dominant feedback mechanism, star formation will continue until ∆ M ion ∼ M . In terms of glyph[epsilon1] ∗ , this will occur at the star formation time t sf = glyph[epsilon1] ∗ t ff /glyph[epsilon1] ff , so that  \nglyph[epsilon1] ∗ ∝ glyph[epsilon1] 9 / 13 ff Σ 33 / 52 M 1 / 52 . (15)  \nBound clumps ( α vir ∼ 1) have surface densities that depend only on the turbulence parameter, C = σ/R 1 / 2 (Equation 2), which is observed to be approximately constant in a galaxy (e.g., Sun et al. 2018). For constant glyph[epsilon1] ff ( § 3.2), this implies glyph[epsilon1] ∗ is nearly independent of M , so the ICMF has the same slope as the clump mass function. If on the other hand stars form in a fixed number of free fall times so that glyph[epsilon1] ∗ ∝ glyph[epsilon1] ff , then glyph[epsilon1] ∗ ∝ Σ 33 / 16 M 1 / 16 , very close the results of Kim, Kim & Ostriker (2018) for small glyph[epsilon1] ∗ .  \nDirect radiation pressure injects momentum at a rate proportional to the stellar mass. Stars will form until the injected momentum ∝ glyph[epsilon1] ∗ Mt sf reaches Mv esc , so that  \nglyph[epsilon1] ∗ ∝ ( glyph[epsilon1] ff Σ) 1 / 2 , (16)  \nwith no explicit dependence on M . If t sf ∝ t ff , so that glyph[epsilon1] ∗ ∝ glyph[epsilon1] ff , then glyph[epsilon1] ∗ ∝ Σ, as found by Fall, Krumholz & Matzner (2010). In this case, or if glyph[epsilon1] ff is about constant, the SFE is independent of mass for bound clumps, which have similar surface densities in a given galaxy as noted above. As a result, the slope of the ICMF is the same as that for the clump mass function, α M, clump ≈ α M ≈ -2.  \n3.5.2. Bound clusters. With our broad definition of cluster ( § 1.2.2), essentially all stars are born in clusters (c.f. Lada & Lada 2003). However, most clusters are unbound by gas dispersal and thus disperse rapidly, so that Γ glyph[lessorsimilar] 0 . 1 for T glyph[greaterorsimilar] 10 Myr. In this section we investigate the origin of the low value of Γ after gas dispersal. As discussed in § 3.1.1, most GMCs appear to be bound (i.e., α vir glyph[lessorsimilar] 2), as are the dense clumps within them (Urquhart et al. 2018), so most stars are born in gravitationally bound gas. Even in GMCs that are moderately unbound, α vir ≈ 5, efficient star formation results in most stars being confined to gravitationally bound structures (Clark et al. 2005). Thus the low value of Γ must be related to the star formation efficiency, glyph[epsilon1] ∗ , and to how the gas not converted to stars is ejected.  \nThe first author to estimate Γ considered a spherical cloud in virial equilibrium that turns a fraction glyph[epsilon1] ∗ of its mass into stars and ejects the rest (Hills 1980). If gas removal occurs over many dynamical times, the process is adiabatic, the cluster remains bound for any glyph[epsilon1] ∗ , and its final and initial radii are related by R f /R 0 = 1 /glyph[epsilon1] ∗ . If gas removal is rapid, ejection of more than half the mass ( glyph[epsilon1] ∗ < 0 . 5) unbinds the stars; for glyph[epsilon1] ∗ > 0 . 5 the radius expands, R f /R 0 = glyph[epsilon1] ∗ / (2 glyph[epsilon1] ∗ -1). Even for glyph[epsilon1] ∗ < 0 . 5, however, the low-velocity part of the cluster can remain bound (Lada, Margulis & Dearborn 1984). In simulations, Kroupa, Aarseth & Hurley (2001) find bound remnants remain for glyph[epsilon1] ∗ as small as 0.3.  \nMore realistic starting conditions can modify these conclusions. If the stars are initially subvirial, with Q 0 = 2 T ∗ / | W ∗ | < 1, where T ∗ and W ∗ are the initial kinetic and potential energies of the stars, respectively, then one can show that the minimum glyph[epsilon1] ∗ for the entire star cluster to be bound is glyph[epsilon1] ∗ , min = Q 0 / 2. If a cluster has substructure, then the additional binding energy can help it survive (Allison et al. 2009); whether the cluster remains bound depends unpredictably on the exact spatial distributions of gas and stars and on the timing of gas removal (Smith et al. 2011, 2013). The minimum value of glyph[epsilon1] ∗ for a bound cluster is also reduced if the stars are more concentrated than the gas, so that the local value of glyph[epsilon1] ∗ increases inward. Adams (2000) showed that part of a cluster embedded in a gas distribution approximating an isothermal sphere could survive sudden gas expulsion down to very low values of glyph[epsilon1] ∗ . A cluster that forms in a centrally concentrated gas cloud with a constant value of glyph[epsilon1] ff has a star formation rate ˙ ρ ∗ ∝ ρ 3 / 2 (Parmentier & Pfalzner 2013); in a particular implementation of this model, Shukirgaliyev et al. (2017) found a minimum  \nKrumholz, McKee, & Bland-Hawthorn  \nglobal star formation efficiency glyph[epsilon1] ∗ , min = 0 . 15 for some stars to remain bound. Analysis of a simulation of the formation of a cluster led Kruijssen et al. (2012a) to conclude that most cluster stars form in gas-poor subclusters that have locally high glyph[epsilon1] ∗ and are therefore relatively immune to gas expulsion; however, this work did not include feedback or magnetic fields, and thus its glyph[epsilon1] ∗ values are unrealistically high.  \nReal star-forming regions are hierarchically-structured, containing both dense parts for which mass removal is slow compared to the local dynamical time, and diffuse parts for which it is fast; we see this directly in the Orion star forming complex ( § 3.4). Elmegreen (2006, 2008) was the first to point out that this configuration naturally gives rise to a gravitationally-bound central region and an unbound periphery. Kruijssen (2012) provided a quantitative model for Γ in such an environment, based on three main premises: first, the density PDF is lognormal, as expected in supersonically-turbulent flows; second, glyph[epsilon1] ff ≈ 0 . 01 independent of density, as observations suggestion ( § 3.2); third, gas at all densities is removed rapidly after some fixed timescale t sf set by stellar feedback. Given these ingredients, he calculated a density-dependent star formation efficiency glyph[epsilon1] ∗ ( ρ ) = glyph[epsilon1] ff [ t sf /t ff ( ρ )], and integrated over the density PDF to derive the fraction of mass for which glyph[epsilon1] ∗ glyph[greaterorsimilar] 0 . 5, which he associated with the bound fraction Γ. A generic prediction of this model is that Γ increases with mean ISM density or pressure, although such a prediction is not unique to this model, since many of the feedback mechanisms discussed in § 3.3 become less effective at higher gas surface density.  \nIn this view, the low value of Γ is not so much 'infant mortality' as 'infant weight loss': stars form in hierarchical structures where only a small fraction remains bound after gas dispersal, so rather than Γ = 0 . 1 meaning that 90% of clusters disperse and 10% remain bound, it means that 90% of the stellar mass is in the unbound part of the structure that escapes immediately after gas dispersal. Of course in reality both mortality and weight loss may operate: hierarchical structures may lose much of their stellar mass in low-density regions, and then even the denser regions that survive may disperse soon thereafter if they are insufficiently concentrated to survive stellar evolution and tidal shocking ( § 4).  \nAlthough the case in which stars form in a pre-existing density structure (either uniform or hierarchically-structured) has garnered the most attention in the literature, we have argued that it is more likely star formation and cloud assembly are contemporaneous. One can generalize the star formation model in § 3.2 to allow for gas accretion while the stars form (see 'Conveyor-Belt Model for Cluster Formation' in the Supplementary Materials). However, provided the time over which the cloud accretes is < ∼ t sf (Equation 10) the results are similar to those above. If the accretion stops due to feedback, then it is likely that the subsequent mass loss would be fast, since the same processes stopping the accretion would eject the gas in the clump. On the other hand, if the accretion stops because the supply of new gas runs out, then stars would continue to form and the mass loss would be slow.  \nOne final consideration is that clusters form in tidal gravitational fields ( § 4.1), which sets a minimum density for the bound part of the cluster, ρ ti . After mass loss, the final mean density of a cluster with initial radius R 0 and density ρ 0 is (Mathieu 1983)  \nρ ∗ f ρ 0 = f b glyph[epsilon1] ∗ ( R 0 R f ) 3 →    (2 glyph[epsilon1] ∗ -1) 3 glyph[epsilon1] 2 ∗ fast mass loss, glyph[epsilon1] 4 ∗ slow mass loss, (17)  \nwhere f b is the fraction of the stellar mass that is bound, and we have assumed f b = Q 0 = 1 in the expressions for the limiting cases. 15 Clusters are generally centrally concentrated, and stars in the outer parts of the cluster in which the final stellar density, ρ ∗ f , is less than the tidal limit will be lost. For a simple spherically symmetric model in which (1) the gas density varies as r -k with k glyph[similarequal] 2 (Schneider et al.  \n2015), (2) glyph[epsilon1] ∗ ∝ 1 /t ff (as in the Kruijssen 2012 model), so glyph[epsilon1] ∗ ∝ ρ 1 / 2 ∝ 1 /r , and (3) glyph[epsilon1] ∗ near the center approaches unity (see 'Cluster-Forming Clouds and Clumps' in the Supplementary Materials), the bound fraction is large, f b glyph[greaterorsimilar] 0 . 5, for values of the cloud-averaged glyph[epsilon1] ∗ glyph[greaterorsimilar] 0 . 01; in this model, the mean value of glyph[epsilon1] ∗ for the stars that remain bound is glyph[epsilon1] ∗ b glyph[greaterorsimilar] 0 . 3. 16  \nAltogether, theoretical estimates of the bound fraction are the product of several uncertain factors. First, after the accretion of gas onto the cluster-forming clump has ceased, what fraction of the clumps are bound? Second, what fraction of the stars remain bound after the residual gas is ejected? Observations suggest that star formation is coincident with clump formation, at least in the early stages, but it is not clear whether this process is halted by feedback, leading to rapid mass loss and requiring a relatively high star formation efficiency, or by a lack of gas, which would be consistent with a lower efficiency. Tides due to the local environment can also strip stars from the outer regions of the cluster when mass is lost, a point to which we return in § 4.3.1.  \n3.5.3. Numerical results. A few authors have attempted to simulate the full range of processes - formation of a mass spectrum of clumps, conversion of gas into stars, and gas removal by feedback - that determine the ICMF and Γ. This is very challenging numerically, since capturing the formation of a statisticallymeaningful sample of clusters requires simulating either an entire galaxy or a substantial portion thereof. Modern isolated galaxy or cosmological zoom-in simulations generally reach spatial resolutions no better than ∼ 1 pc, and mass resolutions of ∼ 10 2 -10 3 M glyph[circledot] . Comparing these figures to the star cluster properties shown in Figure 9, even for GCs this corresponds to glyph[lessorsimilar] 10 3 mass resolution elements per cluster, and glyph[lessorsimilar] 3 -10 spatial resolution elements per half-mass radius. Thus simulations must rely on parameterized subgrid models to handle both star formation and feedback.  \nPerhaps not surprisingly, this situation yields little consensus. Renaud, Bournaud & Duc (2015) simulate an Antennae-like galaxy merger, and find that Γ is a few percent pre-merger, but rises to ≈ 10% during the main burst of star formation. In contrast, Li, Gnedin & Gnedin (2018) find Γ values from ≈ 1% to ≈ 50% in their cosmological simulations that run to z ≈ 1 . 5. Renaud, Bournaud & Duc (2015) and Maji et al. (2017) find lognormal ICMFs in simulations of both mergers and quiescent spiral galaxies (contrary to observations), while Dobbs et al. (2017), Li et al. (2017), and Li, Gnedin & Gnedin (2018) find ICMFs that are well-described by power-laws or Schechter functions with α M ≈ -2. Li, Gnedin & Gnedin (2018) test a variety of star formation and feedback prescriptions, all tuned to reproduce the galaxy-scale Kennicutt relation and various other constraints, and find that both Γ and the ICMF are extremely sensitive to the choice of subgrid model.  \nWe conclude that present galaxy simulations have little predictive power when it comes to star clusters; instead, a more productive approach would be to use the observed properties of star clusters as an additional constraint to calibrate the subgrid models. There do, however, appear to be two qualitative results that persist across recipes. First, all simulations find that Γ increases at higher surface densities of star formation, which occur in regions of higher pressure. Second, all authors find that the ICMF extends to higher masses during mergers. This is also consistent with the results of a number of other cosmological simulations showing that early, gas-rich galaxies are capable of producing clusters with initial masses glyph[greaterorsimilar] 10 6 M glyph[circledot] (Kimm et al. 2016; Ricotti, Parry & Gnedin 2016; Kim et al. 2018).", '3.6. Elemental abundance distributions': "There has been surprisingly little theoretical work on the physical origins of abundance homogeneity in star clusters. Murray & Lin (1990) observed that narrow giant branches in GC CMDs implied they were highly uniform in [Fe/H], and suggested that turbulent diffusion could explain why clouds are highly uniform prior to the onset of star formation. This process would mix out abundance inhomogeneities at the scale of a cloud over a cloud crossing time, and smaller-scale inhomogeneities even more quickly. Feng & Krumholz (2014) provide numerical support for this idea, finding in a simulation of star formation in a colliding flow of warm neutral gas that the mixing process washes out homogeneities in the initial colliding streams by a factor of 5 or more. This may be telling us that the small inhomogeneities detected in a few open clusters (e.g. M67; Liu et al. 2016b) reflect metal anomalies (predicted to be ∼ 0 . 1 dex - Krumholz & Ting 2018) in the birth cloud ( § 2.6). Feng & Krumholz also find that even unbound clusters can also be homogeneous, a prediction that has yet to be tested by observations. Armillotta, Krumholz & Fujimoto (2018) revisit the Feng & Krumholz (2014) simulations using more realistic initial conditions, a cloud extracted from a galaxy simulation, and zoom in to follow it collapse into individual stars at a resolution of ≈ 10 -3 pc. They find that the abundance scatter between formed stars decreases on shorter length scales, indicating that turbulent mixing becomes more efficient on smaller and smaller scales ( glyph[lessorsimilar] 5 pc), in agreement with the discussion above. Interestingly, they also find that some stellar clusters are separable from others within the collapsing clouds based on abundances, but not all.  \nA few caveats are in order. First, the simulations of Feng & Krumholz (2014) and Armillotta, Krumholz & Fujimoto (2018) assume that metals are well-coupled to the gas, but this need not be the case for metals in the form of dust grains. Padoan et al. (2006) first raised the prospect that small-scale turbulence could generate low-level variations in the dust to gas ratio. More recently, Hopkins & Lee (2016) and Lee, Hopkins & Squire (2017) find far more dramatic variations in their simulations and even suggested that some stars could form in 'totally metal' regions (Hopkins 2014). But this work has been refuted by Tricco, Price & Laibe (2017) who find no evidence for significant variations; the use of tracer particles in supersonic flows lead to numerical artifacts that exaggerate density contrasts. Second, the cluster formation simulations to date have not considered the effects of magnetic fields, which Sur, Pan & Scannapieco (2014) suggest can suppress the mixing over a wide range in Mach number. Here, the same random stretching and folding that is responsible for turbulent diffusion also amplifies the magnetic energy density by the dynamo mechanism, which in turn can suppress the mixing on small scales.  \nA third caveat is that, while galaxy-scale simulations are starting to consider injection of new metals by stellar evolution (e.g., Emerick, Bryan & Mac Low 2019), simulations that resolve cluster formation have yet to do so. Bland-Hawthorn, Krumholz & Freeman (2010) point out that, even if turbulence homogenizes a cloud, it will only remain homogeneous if star formation ceases before the first SN explodes, since even a single SN will produce measurable abundance variations. For a gas cloud with mass M , surface density Σ, and virial ratio α vir , the crossing time is t cr ≈ ( α vir G ) -1 / 2 M 1 / 4 Σ -3 / 4 (c.f. § 3.1.1). If we adopt α vir = 1 . 5, and the cluster forms over τ sf = t sf /t cr crossing times, then  \nt sf ≈ 3 . 0 ( τ sf 4 )( glyph[epsilon1] ∗ 0 . 2 ) -1 / 4 ( M ∗ 10 4 M glyph[circledot] ) 1 / 4 ( Σ 0 . 3 g cm -2 ) -3 / 4 Myr , (18)  \nwhere glyph[epsilon1] ∗ = M ∗ /M is the star formation efficiency. Thus t sf < t SN , which is conservative given that most SN progenitors have longer lifetimes, holds for OCs (Σ ≈ 0 . 3 M glyph[circledot] pc -2 ) up to 10 5 M glyph[circledot] , and for GCs (Σ ≈ 3 M glyph[circledot] pc -2 ) up to 10 7 M glyph[circledot] . This last estimate is interesting in light of the fact that globulars can show high levels of Fe homogeneity; as we have seen, for other elements, the story is more complicated. How there can be differential scatter across elements is not understood, unless it reflects the state in the collapsing cloud (see below). In Figure 9, there is evidently a lot of scatter but our point is to emphasize the difference between  \nOCs and denser GCs.  \nBoth this calculation and the simulations of Armillotta, Krumholz & Fujimoto (2018) raise the question of how to apply estimates of chemical homogeneity or inhomogeneity in the context of hierarchical starforming systems that do not have clear edges. This topic is currently under investigation using a combination of Gaia data and ground-based, million-star surveys at high spectroscopic resolution (e.g. GALAH; De Silva et al. 2015). The astrometric information can be used to identify every aggregate within dozens of nearby stellar hierarchies (Gouliermis 2018), while the spectra provide abundances for all stars cooler than ≈ 7000 K. This work will provide the first clues for how to map stars between abundance and physical space, i.e., the possibility of determining how closely together in space two stars formed based on the level of similarity or difference in their elemental composition, with major implications for the future of chemical tagging in Galactic archaeology (Freeman & Bland-Hawthorn 2002; Bland-Hawthorn & Freeman 2004). It is also not clear how to apply these models in the context of multiple generations and detectable abundance anomalies in globular clusters ( § 2.6).", '4. LIFE AND DEATH': 'Since the slope of the CAF α T < 0 even for ages older than ≈ 10 Myr ( § 2.3), cluster disruption cannot depend solely on gas expulsion immediately after cluster formation. Instead, there must be processes that destroy star clusters during the gas-free phase of their evolution. It is to these processes that we now turn. We begin with a brief review of star cluster structure and tidal fields in § 4.1, then examine disruption processes driven primarily by mechanisms internal to the cluster ( § 4.2) and those driven primarily by the external environment ( § 4.3), before concluding with attempts to synthesize both into a coherent explanation for the demographics of star clusters § 4.4. Since our focus is on the demographics and evolution of the star cluster population as a whole, we will not discuss processes that alter the internal structure of star clusters while leaving their bulk properties unchanged, e.g., dynamical alteration of the binary fraction or stellar collisions. We refer readers to Portegies Zwart, McMillan & Gieles (2010) and Renaud (2018) for a discussion of these processes.', '4.1. Star cluster structure and tides': "A crossing time after gas is removed, star clusters that remain bound relax to a nearly-spherical, virialized state. From this point up to ages of several hundred Myr, their density profiles are well described by the Elson, Fall & Freeman (1987, EFF hereafter) distribution, which has surface density as a function of radius Σ( r ) = Σ 0 ( 1 + r 2 /a 2 ) -γ/ 2 . Here Σ 0 is the central surface density, a defines the radial scale, and γ describes the fall off in surface density in the outer, power-law region; observed values of γ are generally in the range ≈ 2 . 2 -3 . 2 (Elson, Fall & Freeman 1987; Mackey & Gilmore 2003a,b). The core radius for an EFF profile, defined to match that of a King profile (see below) is r c ≈ a √ 2 2 /γ -1.  \nThe EFF profile is not in energetic equilibrium, and thus clusters observed to have this distribution cannot have reached either internal relaxation or equilibrium with the tidal field of their host galaxy. At ages beyond a few hundred Myr, clusters relax and become well-fit by King (1962, 1966) models where stellar energies are distributed as a 'lowered Maxwellian' of the form f ( E ) ∝ (2 πσ 2 ) -3 / 2 [ exp ( -E/σ 2 ) -1 ] for some constant velocity dispersion parameter σ . Physically, the -1 in square brackets reflects the fact that the velocity distribution cannot be a pure Maxwellian, because the cluster has a finite escape speed, and thus any star that is too far out on the tail of the distribution will leave. Combining this assumed energy distribution with the Poisson equation for self-gravity yields an equation for the radial density distribution and potential. Solutions to this equation form a one-parameter family, with the parameter conventionally specified by setting the dimensionless potential well depth at the cluster center W 0 = ( v esc /σ ) 2 , where  \nFigure 15  \n<!-- image -->  \n<!-- image -->  \nIsosurfaces of the potential φ ( x, y, z ) (Equation 19) for a point-like star cluster embedded in a galaxy, seen in 3d ( left ) and in a cut through the xz plane ( right ). The cluster center is located at the origin, while the galaxy center is located at ( -R g , 0 , 0) for some distance R g much larger than the physical size of the cluster. In the example shown, λ y = 0 and λ z /λ x = -1 / 3. Blue shows the Jacobi surface, defined by φ = E J , while green shows a potential φ = 0 . 95 E J , i.e., with energy 5% greater than the Jacobi energy. Orange shows a sphere of radius r = r ti (Equation 20). The Lagrange points L 1 and L 2 are as indicated. Stars inside the Jacobi surface, shaded blue in the right panel, are bound cluster members if their total energy E < E J . Stars outside the Jacobi surface but at r < r ti (shaded orange), or those inside the Jacobi surface with E ≥ E J , are potentially unbound from the cluster, and can escape through the 'windows' around L 1 and L 2 .  \nv esc is the escape speed from the cluster center. Most observed OCs and GCs are reasonably well-fit by W 0 ≈ 3 -7, though a small number of core-collapsed globulars and some of the most concentrated YMCs have W 0 glyph[greaterorsimilar] 10.  \nFor a King model, one can define the core radius r c = √ 9 σ 2 / 4 πGρ 0 , where ρ 0 is the central density. The core radius is distinct from the half-mass radius r h ( § 2.5), and r h /r c is a function of W 0 , but for W 0 ≈ 3 -7, r h /r c is within a factor of a few of unity. A second useful quantity is the truncation radius r tr , the radius at which the density falls to zero. This can be combined with the core radius to yield the concentration c = log( r tr /r c ). This is a monotonically increasing function of W 0 , so clusters with deeper potential wells are more centrally concentrated. For this reason, clusters are often parametrized in terms of values of c instead of values of W 0 . The observed range W 0 ≈ 3 -7 corresponds to c ≈ 0 . 6 -1 . 5. Note that the EFF model does not have a truncation radius.  \nKing models in principle apply to an isolated star cluster, but it is common to assume that the truncation at r tr is due to the tidal field of the cluster's host galaxy. 17 One can derive an effective potential in the frame  \nco-moving with a cluster orbiting with a galaxy by Taylor-expanding the galactic potential and adding this to the centrifugal acceleration in the rotating reference frame. If one approximates the cluster as point-like, the resulting total potential is (e.g., Renaud, Gieles & Boily 2011)  \nφ ( x, y, z ) = -GM √ x 2 + y 2 + z 2 -λ x 2 ( x 2 + λ y λ x y 2 + λ z λ x z 2 ) , (19)  \nwhere M is the cluster mass, λ x , λ y , and λ z are the three eigenvalues of the tidal plus centrifugal tensor, 18 and the coordinate system is oriented so that the center of the galaxy lies at on the x axis at x < 0, and the tidal plus centrifugal tensor has zero off-diagonal components. For a star cluster on a circular orbit at the midplane of an azimuthally-symmetric galactic potential φ G , adopting a coordinate system where the y direction is spinward, the z direction is out of the galactic plane, and the coordinate system is right-handed, the eigenvalues of the tensor are λ x = -∂ 2 φ G /∂x 2 + ∂ 2 φ G /∂z 2 , λ y = 0, and λ z = -∂ 2 φ G /∂z 2 .  \nThe potential φ has five extrema or saddle points, known as the Lagrange points; of these, by convention the points L 1 and L 2 lie at positions L 1 , 2 = ( ∓ r ti , 0 , 0), where  \nr ti = ( GM λ x ) 1 / 3 (20)  \nis known as the tidal radius. At these points, the acceleration due to the gravity of the cluster is equal to that due to the combination of galactic gravity and the centrifugal force. The value of the potential at L 1 and L 2 is the Jacobi energy E J = -(3 / 2) GM c /r ti , and the locus where φ ( x, y, z ) = E J is called the Jacobi surface (Figure 15).  \nStars with total energy E < E J are unable to cross the Jacobi surface, so any such stars located within it are bound to the cluster. Since the volume within this surface defines the volume of the star cluster, tidal truncation therefore defines a characteristic density ρ ti ≈ 3 M/ 4 πr 3 ti = 3 λ x / 4 πG , which depends only on the tidal tensor and thus on the cluster's orbit and environment, and not on the cluster mass. Stars with energies E ≥ E J are not bound to the cluster, but if they are located inside the Jacobi surface then they can only escape through a 'window' around L 1 and L 2 , with the size of the window depending on E/E J . We show an example escape window for E/E J = 0 . 95 (i.e., for a star with energy 5% greater than the Jacobi energy) in Figure 15. With the exception of those ejected during violent tidal shocks, most stars escaping a cluster have energies only slightly larger than E J , so escape tends to be primarily near the Lagrange points.", '4.2. Internal evolutionary processes': "With these preliminaries understood, we can now ask what mechanisms can alter clusters' properties during the gas-free phase of their evolution.  \n4.2.1. Stellar evolution. At ages from ≈ 3 -40 Myr, stellar populations eject ≈ 20% of their mass in SNe, which produce high-speed gas that almost certainly escapes the parent cluster; moreover, the remnant neutron stars or black holes left behind by SNe receive kicks of several hundred km s -1 (e.g., FaucherGigu'ere & Kaspi 2006), so most escape the star cluster as well. After ≈ 40 Myr, SNe cease and mass loss  \nbecomes dominated by asymptotic giant branch (AGB) stars shedding their envelopes. This process is more gradual, so that the total mass loss reaches ≈ 40% only after ≈ 1 Gyr, and does not reach 50% until ages of order a Hubble time. Moreover, while AGB wind velocities are high enough (tens of km s -1 ) that little of the wind mass is likely to be retained, the white dwarfs left behind receive much smaller kicks, and thus are likely to remain in the cluster (e.g., Kruijssen 2009).  \nBy itself a factor of two decrease in mass would be a small effect. However, the effects of mass loss can be greatly amplified by a tidal field. Adiabatic mass loss makes the cluster radius expand while the tidal radius shrinks; stars that find themselves outside r ti as a result are lost. For clusters that begin their evolution filling their tidal radii, this effect roughly doubles the mass loss rate due to stellar evolution alone (Takahashi & Portegies Zwart 2000; Baumgardt & Makino 2003; Lamers, Baumgardt & Gieles 2010). Moreover, in the presence of a tidal field, for any specified cluster density profile there is a maximum ratio of half-mass radius to tidal radius, r h /r ti , for which virial equilibrium is possible (Chernoff & Weinberg 1990; Fukushige & Heggie 1995). If mass loss drives r h /r ti above this maximum, the stars cannot find a new virial equilibrium, and the entire cluster disrupts. Both the additional mass loss and the propensity to catastrophic disruption are enhanced in mass-segregated clusters, where the stars losing mass carry a disproportionate share of the potential energy (Vesperini, McMillan & Portegies Zwart 2009; Gieles et al. 2010). The importance of stellar evolution depends on a star cluster's initial concentration and ratio of truncation radius to tidal radius r tr /r ti . Clusters with r tr /r ti glyph[lessmuch] 1 or W 0 glyph[greaterorsimilar] 5 ( c glyph[greaterorsimilar] 1) are dense enough to behave qualitatively like clusters that are not in a tidal field. Clusters with low concentration that already fill their tidal radii at the start of gas-free evolution are much more vulnerable.  \nUnfortunately there have been only limited observational studies of c and r tr /r ti for young clusters. Ryon et al. (2015, 2017) survey massive ( > 5000 M glyph[circledot] ) clusters in M83, NGC 628, and NGC 1313, and find that about half have r h /r ti glyph[lessorsimilar] 0 . 15 independent of cluster mass, small enough to render mass loss unimportant. The fraction with larger values of r h /r ti increases with age, as previously observed by Elson, Freeman & Lauer (1989) in the LMC, but too late for stellar mass loss to be effective. However, we caution that accurate measurements of r h require accurate photometry in clusters' extended, low surface brightness halos, which is very difficult to obtain for galaxies at the distances of Ryon et al.'s sample. The result that the concentration does not correlate with mass is more robust, and appears to hold for OCs in the Milky Way (based on our examination of Kharchenko et al. 2013's sample) and the Magellanic Clouds as well (Fall & Chandar 2012, using the measurements of McLaughlin & van der Marel 2005), where measurements are easier. For Milky Way GCs, McLaughlin (2000) reports a weak correlation between concentration and mass, but no such correlation is apparent in the more recent fits of Baumgardt & Hilker (2018). Thus while there are hints that clusters might be born too concentrated for mass loss to be a significant effect, the observational case is far from settled. Moreover, we caution that all observations are potentially subject to both selection bias in favor of more compact clusters and survival bias - a substantial fraction of gas-free clusters at an age of ≈ 10 Myr might be only weakly bound and have W 0 ≈ 2 -3, for example, but these would be dramatically under-represented in OC catalogs because almost none would survive to reach ages more than a few tens of Myr.  \n4.2.2. Relaxation-driven mass loss. Stellar evolution has its greatest effects in the first ∼ 10 -100 Myr of a cluster's lifetime. On much longer time scales, the dominant internal process is mass loss driven by relaxation, sometimes also referred to as dissolution or evaporation. The underlying mechanism is simple: because stars with too much energy escape from clusters, the high-energy tail of a cluster's velocity distribution is underpopulated compared to a Maxwellian. As stars randomly exchange energy via gravitational interactions, the velocity distribution attempts to relax to a Maxwellian, repopulating the tail. Once stars scatter to high enough energy by this process, however, they escape through the windows around L 1 and  \nL 2 , leaving the tail again underpopulated and continuing the cycle. Throughout the Galactic halo, we see direct evidence of this process through the disruption of GCs, e.g., Pal 5 (Odenkirchen et al. 1997, 2002), NGC 5466 (Grillmair & Johnson 2006), and Pal 14 (Sollima et al. 2011). Stars escaping through L 1 are at smaller galactocentric radii and thus move faster in the galaxy than the cluster which gives rise to a leading stream; stars escaping through L 2 are slower than the cluster and therefore trail. The escaping stars give rise to spectacular, cold stellar streams extending up to 60 · across the sky (Odenkirchen et al. 2002; Grillmair & Johnson 2006). Limited evidence for the tidal disruption of OCs has been claimed but is rather less convincing (e.g., Berkeley 17; Bhattacharya et al. 2017).  \nSince this mass loss process is driven by repopulation of the Maxwellian tail as the velocity distribution relaxes, the loss rate scales approximately as ˙ M rlx ∝ M/t rlx , where t rlx ≈ 0 . 1( N/ ln N )( r h /σ ) is the relaxation time at the half-mass radius, M is the cluster mass, N is the number of stars in the cluster, and σ is the velocity dispersion. H'enon (1961, 1965a,b) carried out some of the earliest modeling of this process, finding that ˙ M rlx = f rlx M/t rlx with f rlx ≈ 0 . 05, i.e., clusters lose about 5% of their mass per relaxation time. If we set N = M/m , where m is the mean stellar mass, and σ ≈ √ GM/r h , the mass loss rate scales as ˙ M rlx ∝ √ ρ , where ρ ∝ M/r 3 h and we have dropped the weak mass dependence in ln N . If clusters have density independent of mass (requiring r h vs. M only a slightly steeper than the relation in Figure 9), or if they are tidally limited and thus have ρ ≈ ρ ti independent of M ( § 4.1), then ˙ M rlx is independent of M as well. This process is the justification for the 'evolved Schechter' function form introduced in § 2.2 to fit the globular CMF, and explains how a Schechter-like ICMF could evolve into such a distribution (Fall & Zhang 2001; McLaughlin & Fall 2008).  \nWhile ˙ M rlx ∝ M/t rlx is a reasonable first approximation, N -body simulations suggest a slightly more complex dependence. Fukushige & Heggie (2000), Baumgardt (2001) and Baumgardt & Makino (2003) show that, due to the finite time required for potential escapers to find the L 1 or L 2 windows and pass through them, the mass loss timescale depends on the crossing time t cr = r h /σ as well as the relaxation time, so ˙ M rlx ≈ f rlx M/ ( t p rlx t 1 -p cr ), where p ≈ 0 . 7 -0 . 8 depending on W 0 . The corresponding cluster dissolution timescale is  \nt dis = M ˙ M rlx ≈ { 1 . 6 ρ -1 / 2 2 M 3 Gyr , p = 1 0 . 6 ρ -1 / 2 2 M 0 . 7 3 Gyr , p = 0 . 7 , (21)  \nwhere ρ 2 = ρ/ (100 M glyph[circledot] pc -3 ), M 3 = M/ (10 3 M glyph[circledot] ), and for the numerical evaluation we have used f rlx = 0 . 05, m = 0 . 5 M glyph[circledot] , and ln N = 7 . 6. Subsequent studies have shown that, in addition to W 0 , p depends on the extent to which clusters initially fill their tidal radii (Gieles & Baumgardt 2008; Madrid et al. 2017), the ratios of the eigenvalues of the tidal tensor (Tanikawa & Fukushige 2010; Madrid et al. 2017), and the mass spectrum and degree of mass segregation in the cluster (Hong et al. 2013). Goudfrooij & Fall (2016) argue that p should also depend on cluster mass, and that the observed GC luminosity function implies p ≈ 1 for globulars; however, this analysis does not include possible mass-to-light ratio effects due to correlations between GC metallicity and formation history.  \nA number of authors have proposed analytic models for cluster evolution that attempt to capture the findings of the N -body simulations, though none thus far have been able to capture the full range of numerical results over all parameter space. Lamers et al. (2005) and Lamers, Baumgardt & Gieles (2010) propose a model for cluster evolution incorporating the dependence of mass loss on the crossing time and the initial tidal radius, which gives ˙ M rlx ≈ M 1 -γ /t 0 , with γ ≈ 0 . 65 and t 0 ∝ 1 / Ω, where Ω is the angular velocity of the cluster's orbit through its host galaxy. Gieles, Heggie & Zhao (2011) propose a pair of 'unified evolution equations' that capture the process of expansion and subsequent contraction that occurs when a cluster born with r ti < r tr expands to its tidal radius and then shrinks as mass loss reduces that  \nradius:  \nξ ≡ -˙ N N t rlx = 3 5 ζ t cr t cr , 1 χ ≡ ˙ t cr t cr t rlx = 3 2 ζ ( 1 -t cr t cr , 1 ) , (22)  \nwhere ξ and χ are dimensionless rates of mass loss and cluster expansion and t cr , 1 is the crossing time of the cluster at the point where it has expanded to its tidal radius. The quantity ζ is the fraction of the total energy conducted per t rlx , and has a complex dependence on the stellar mass spectrum and on the cluster core properties that must be calibrated by N -body simulations. The idea that clusters are born with r tr < r ti and expand during their first ≈ 200 Myr of life is consistent with the tentative findings of Ryon et al. (2015, 2017), and with the observation that a higher proportion of the older GC population, though not all, fill their tidal radii (Innanen, Harris & Webbink 1983; Baumgardt et al. 2010).", '4.3. External and environmental processes': "Our discussion of internal processes implicitly assumed a constant tidal field. However, as clusters orbit the tidal field they experience may change due to cluster encounters with massive objects, or because the cluster's orbit is not a simple circle in the galactic plane. Such changes in the tidal potential cause the Jacobi surface to move, and sudden perturbations can also alter stellar energies. Both effects can move stars across the Jacobi surface, thereby unbinding part of a cluster. We now explore processes of this type.  \n4.3.1. Tidal perturbations. For star clusters with near-circular orbits in a galactic plane, the primary source of tidal perturbations is encounters with GMCs (Spitzer 1958). If an encounter takes place over a time shorter than the crossing time of a cluster, one can approximate it as impulsive (e.g., Binney & Tremaine 2008). Formally, this requires that the encounter velocity obey v glyph[greaterorsimilar] b √ GM/r 3 h = 2 . 1 b 0 r -3 / 2 h, 0 M 1 / 2 3 km s -1 , where b is the impact parameter, b 0 = b/ 1 pc, M is the cluster mass, M 3 = M/ 10 3 M glyph[circledot] , and r h is the cluster half-mass radius, r h, 0 = r h / 1 pc. Since the typical encounter velocity in the ISM will be of order the ISM velocity dispersion ( ≈ 10 km s -1 in local galaxies), encounters between GMCs and all but the most massive clusters tend to be impulsive. Impulsive encounters can be further classified into catastrophic ones that completely disrupt the cluster, and diffusive ones that do not. Catastrophic disruption occurs only for perturbers that are denser than the cluster, while most star clusters are at least somewhat denser than Milky Way GMCs, n ≈ 100 cm -3 (c.f. Figure 9), so the diffusive regime is more common.  \nThe specific energy added per tidal shock in the diffusive regime, and the resulting mass loss rate, has been discussed by many authors (e.g., Gnedin & Ostriker 1999; Gieles, Athanassoula & Portegies Zwart 2007; Prieto & Gnedin 2008; Binney & Tremaine 2008). For the case of shocking by GMCs of finite size, Gieles et al. (2006) obtain the approximate relations  \n∆ E E ≈ f ( ˜ b ) ( M GMC M ) 2 ( σ v max ) 2 ∆ M M ≈ f M ∆ E E , (23)  \nwhere E and M are the cluster total energy and mass, ∆ E and ∆ M are the changes in these quantities due to the encounter (with a sign convention whereby ∆ M > 0 corresponds to mass loss), f M ≈ 0 . 25 is the ratio of fractional mass loss to fractional energy gain (less than unity since stars escape with finite energy), σ = √ GM/r h is the stellar velocity dispersion in a cluster, v max is the maximum velocity during the encounter, and f ( ˜ b ) is a numerical factor that depends on the dimensionless impact parameter ˜ b ≡ b/r h as f ( ˜ b ) → const for ˜ b glyph[lessmuch] 1 and f ( ˜ b ) → ˜ b -4 for ˜ b glyph[greatermuch] 1. The time to disrupt a cluster, t dis = M/ ˙ M , can be found by integrating over encounters with a population of GMCs. Second-order energy injection (Kundic & Ostriker 1995) reduces the Gieles et al. (2006) disruption time by a factor C sh glyph[similarequal] 0 . 3 (Kruijssen et al.  \n2011) so that the disruption time for a cluster of mass M becomes  \nt dis ≈ C sh √ π f dis σ ISM G Σ GMC ( ρ ρ GMC ) ≈ 1 . 6 σ 1 Σ -1 2 ( ρ/ρ GMC 1000 ) Gyr (24)  \nwhere σ ISM is the velocity dispersion of the ISM (assumed equal to the velocity dispersion of the clusters), Σ GMC is the surface density of an individual GMC, ρ GMC is the mass density of GMCs in the ISM (not the GMC internal density), ρ = 3 M/ 8 πr 3 h is the cluster density, and f dis is a factor of order unity that depends on gravitational focusing and cluster structure; for Gieles et al. (2006)'s favored parameters, f dis ≈ 0 . 4. We have scaled to Milky Way values (c.f. § 3.1.1) for the other parameters: σ 1 = σ ISM / 10 km s -1 ≈ 1, Σ 2 = Σ GMC / 100 M glyph[circledot] pc -2 ≈ 1, and ρ/ρ GMC ≈ 1000. The 1.6 Gyr disruption timescale is substantially longer than earlier estimates (e.g., Binney & Tremaine 2008), largely because Gieles et al. (2006) find f M ≈ 0 . 25, while some earlier analytic models adopted f M ≈ 1.  \nOne factor that could increase the disruption time is that tidal shocks tend to be self-limiting: a first tidal shock strips off the most weakly bound stars, but the timescale for the cluster to repopulate the high energy orbits these stars occupied is dictated by two-body relaxation ( § 4.2.2), so that subsequent tidal shocks separated in time by less than t rlx remove much less mass (Kruijssen 2015; Gieles & Renaud 2016). On the other hand, the numerical evaluation in Equation 24 could underestimate the importance of tidal heating because it is based on mean conditions in the ISM. However, star clusters are formed in regions of higher than average gas density and structure, i.e., larger ρ GMC and Σ GMC . Star cluster locations remain correlated out to ages ≈ 100 Myr (Grasha et al. 2017a), and during this time they likely experience enhanced tidal perturbations, a phenomenon that Kruijssen et al. (2011) dubbed the 'cruel cradle effect'. Several authors have suggested this effect may explain why cluster number counts decline from 10 -100 Myr, when no other mechanism except stellar mass loss ( § 4.2.1) is expected to operate (Elmegreen & Hunter 2010; Kruijssen 2012; Miholics, Kruijssen & Sills 2017). It is likely to be particularly important in gas-rich environments such as high-redshift galaxies and local mergers, and may explain why globular clusters have a lognormal rather than a power-law mass function: since t dis ∝ ρ , if density increases with mass then low-mass clusters will be preferentially destroyed by tides, while massive clusters will survive (Elmegreen 2010; Kruijssen 2015). However, the importance of this effect depends on the assumed mass-radius relation. Most published models assume mass-independent cluster radii, so ρ ∝ M , which tends to strengthen the preferential destruction of low-mass clusters compared to the mass-radius relation shown in § 2.5.  \nThis situation is somewhat different for clusters in disk-crossing or bulge-crossing orbits, which experience changing tides as they pass through a disk or bulge. During each passage the galaxy exerts a pinch compression on the cluster that acts normal to the disk. The formalism for modelling the effects of a disk passage is largely the same as that for GMC encounters, simply with a different frequency and shock strength, and with some subtle modifications for the fact that GMC tidal shocks are 2D, while disk tidal shocks are 1D (Gnedin & Ostriker 1997). Starting with Ostriker, Spitzer & Chevalier (1972), there is an extensive literature on transiting clusters, with an in-depth review provided by Binney & Tremaine (2008). The timescale for complete destruction of globulars through most of the Milky Way disk is much longer than a Hubble time, but erosion of the outer stars requires only ≈ 10 Gyr. Thus, the stellar escape feeding the cold streams observed in Pal 5 (Odenkirchen et al. 1997, 2002), for example, may have been exacerbated by repeated crossings. Globulars that transit the bulge of a galaxy are subjected to much stronger tides, and may lose substantial fractions of their mass in well under a Hubble time (Madrid et al. 2017).  \nFor open clusters, which are less massive and dense, disk crossing is more perilous. Martinez-Medina et al. (2017, see also Webb et al. 2014) find that the destruction rate is nearly independent of vertical velocity ( ∝ v -0 . 2 z ) in conflict with the v -2 z dependence found by Ostriker, Spitzer & Chevalier (1972); Martinez-Medina et al. find the number of disk crossings is more important than the velocity. Clusters far from the plane transit the disk only rarely, whereas clusters close to the plane experience a weaker gradient  \nwhen they do transit. Consequently, clusters a few hundred parsecs from the plane are the most vulnerable to destruction by tidal shocking. This may explain the clear dichotomy of OCs close to the plane and GCs far from the plane, with very little crossover in phase space. Interestingly, the vast majority of reported OCs more than a few hundred parsecs off the plane are not confirmed in the latest Gaia DR2 analysis (Cantat-Gaudin et al. 2018).  \n4.3.2. Migration and environment. Several mechanisms alter the orbits of star clusters as they evolve, potentially shaping their demographics. Here we discuss processes that move clusters within a single galaxy, deferring discussion of galaxy-galaxy interaction to § 4.4.2. One potentially important mechanism is dynamical friction, which tends to move clusters toward galactic centers by causing them to lose angular momentum. For a cluster of mass M moving with velocity v c through an inertial stellar field described by a distribution function f ( v ) for stars with mass m and velocity v , Chandrasekhar (1943) showed  \nd v c dt = -16 π 2 G ln Λ m ( m + M ) v 3 c ∫ v c 0 v 2 f ( v ) dv v c (25)  \nfor which ln Λ is the Coulomb logarithm; this expression can be simplified considerably if f ( v ) is a Maxwellian. Massive clusters decelerate more rapidly than lower mass ones ( d v c /dt ∝ M for M glyph[greatermuch] m ), because they create larger wakes. This is one of the few mechanisms that encourages the destruction of the most massive clusters: Bekki (2010) finds that dynamical friction causes clusters with masses glyph[greaterorsimilar] 2 × 10 5 M glyph[circledot] in the disks of subL ∗ galaxies to sink to the galactic center in less than a Hubble time, where stronger tides are likely to destroy them (Gerhard 2001). The effect is smaller in larger galaxies like the Milky Way.  \nRadial migration driven by transient spiral density waves, first discussed by Sellwood & Binney (2002) and demonstrated numerically by Roˇskar et al. (2008), can move stars inward or outward. The same mechanism works for clusters, and is non-destructive to them since clusters have internal crossing times much shorter than the time required for them to transit a spiral density wave (Gieles, Athanassoula & Portegies Zwart 2007). Migration occurs because a spiral perturbation with pattern speed Ω P can modify an object's energy E and angular momentum J while conserving its Jacobi integral I J ; in the ( E,J ) plane, objects move along lines of constant I J = E -Ω P J . A single spiral wave near co-rotation can perturb the angular momentum by 20% without significant heating, moving objects from one circular orbit to another, inwards or outwards. Multiple transient spirals lead to a random walk in the ( E,J ) plane at constant I J (Dehnen 2000), so that substantial variations in the angular momentum of an object can occur over Gyr timescales. By contrast, long-term spiral arms produce no net effect, because perturbations as an object enters an arm are cancelled by those as it exits (Lynden-Bell & Kalnajs 1972). Objects that oscillate far from the plane or are on highly radial orbits are much less influenced by the migration mechanism. Minchev & Famaey (2010) and Minchev et al. (2011) propose a related mechanism that does not rely on transience: quasi-chaotic interference between resonances from multiple rotating patterns, e.g., the bar and the spiral arms, with different pattern speeds (c.f. Brunetti, Chiappini & Pfenniger 2011). J'ılkov'a et al. (2012) and Quillen et al. (2018b) investigate this resonance overlap and find that outward migration is possible but relatively inefficient. External influences can also drive radial migration, in particular, radial in-plane orbiting galaxies that come close enough to strongly perturb the disk (Quillen et al. 2009).  \nMinchev & Famaey (2010) provide a comprehensive discussion of the evidence for strong migration of individual stars over billions of years, in particular, the ∼ 1 dex scatter in [Fe/H] at all radii. The high local fraction of metal-rich stars may also be evidence of migration (Kordopatis et al. 2015), since the local ISM has somewhat sub-solar metallicity (Nieva & Przybilla 2012). The case for migrating OCs is less clear. There are at least half a dozen supersolar metallicity OCs within ∼ 1 kpc of the Solar Circle (Quillen et al. 2018b). Moreover, Friel (1995) finds that there is a large scatter in the metallicity and age of OCs at all radii.  \nHowever, given the lower statistics, it is unclear whether this provides compelling evidence for cluster migration within the disk, particularly since there is an alternative possibility: OCs that seem out of place with respect to the average properties of their environment may have been accreted from infalling dwarfs. Carraro & Bensby (2009) argue that the OCs Berkeley 29 and Saurer 1 appear to be associated with the disrupting Sgr dwarf galaxy, and the case for some Milky Way GCs having arrived with Sgr is compelling: Gaia has shown that four GCs (Arp 2, Terzan 7, Terzan 8, M54) confined to the Galaxy's central core are co-moving (q.v. Vasiliev 2018). An OC that is arguably an even stronger candidate for having arrived by infall rather than migration is NGC 6791, an old ( ∼ 8 . 3 Gyr), massive ( ∼ 4000 M glyph[circledot] ), extremely metal-rich ([Fe / H] ≈ +0 . 3) OC along the Solar Circle (Boesgaard, Lum & Deliyannis 2015), 4 kpc from the Sun and 0.8 kpc above the plane (Carraro et al. 2006; Origlia et al. 2006). Its large vertical distance argues against radial migration from the bulge region. While a likely host for infall has not been identified to date, Gaia DR2 data has begun to provide evidence that massive accreted systems like Sgr are dispersed throughout the Galaxy (e.g., Gaia-Enceladus; Helmi et al. 2018b). We anticipate that the more comprehensive Gaia data release (DR3) in 2020 will significantly clarify whether infall can account for all the anomalous OCs, or whether some must be the result of migration.", '4.4. Models of the full cluster population': "Having reviewed the individual processes affecting gas-free clusters, we now discuss attempts to combine them with cluster formation into models for the complete cluster population. We first consider clusters confined to galactic disks, and second examine clusters outside of disks.  \n4.4.1. Disk clusters. Comparison of the observed distribution of cluster masses and ages with theoretical models sheds light on the rate at which clusters are born and on the processes that destroy them. If the mass of a cluster is determined by its initial mass, M i , and its age, T , then the number of surviving clusters in a given age range is the same as that when they were born, so that (Fall & Zhang 2001; Lamers et al. 2005; Fall, Chandar & Whitmore 2009)  \n∂ 2 N ∂M ∂T = C ( T ) M α Mi i ( ∂M i ∂M ) T , (26)  \nwhere C ( T ) is proportional to the birthrate of clusters with a current age T and where we have assumed that the ICMF is a power law in the range of masses considered; as discussed in § 2.2, α Mi ≈ -2. Insofar as the effects of stellar evolution occur prior to relaxation (a good approximation) or tidal effects (a good approximation after the cluster has left its natal complex), one can define M 0 = µ ev ( T ) M i as the mass at the onset of relaxation and tidal effects (Lamers et al. 2005), where µ ev ( T ) drops from unity to ∼ 0 . 25 -0 . 5 during the first few 10 7 years of the cluster's life.  \nTwo classes of models for cluster evolution have been considered, mass-independent destruction (MID; Fall, Chandar & Whitmore 2005; Whitmore, Chandar & Fall 2007; Fall, Chandar & Whitmore 2009) and mass-dependent destruction (MDD; Boutloukos & Lamers 2003; Lamers et al. 2005; Lamers, Baumgardt & Gieles 2010). For the MID case, the cluster distribution function breaks up into the product of a massdependent term and an age-dependent term, which are taken as power laws, so that  \nd 2 N dM dT ∝ M α Mi T α T . (27)  \nSince the effects of stellar evolution are approximately independent of mass as long as initial concentration c does not depend on mass, its effects are naturally included in this formulation. This model is phenomenological, since the rate of cluster destruction is not calculated.  \nFor the MDD case, the destruction time is generally taken to be a power-law function of mass (Boutloukos & Lamers 2003; Lamers et al. 2005), dM/dt = -M/ ( t 0 M γ ), where M is measured in solar masses. If t 0 is constant, then (Lamers et al. 2005)  \n∂ 2 N ∂M ∂T = C ( T ) µ ev ( T ) M α Mi [1 + γT/ ( t 0 M γ )] ( γ +1) /γ , (28)  \nValues of γ estimated from observation (e.g., Boutloukos & Lamers 2003; Gieles 2009) and N -body simulations (Lamers et al. 2005; Lamers, Baumgardt & Gieles 2010) are usually in the range 0 . 6 -0 . 7 (see § 4.2.2). At late times or low masses, the second term in the denominator dominates, so that the distribution function becomes proportional to M α Mi + γ T -( γ +1) /γ . However, because clusters fade with time, it is sometimes the case that in order to be complete, cluster samples must be constructed so that the minimum mass is large enough that the second term is not dominant. It is not clear that any extragalactic sample of clusters has sufficient dynamic range to clearly show the break in the slope of the mass function.  \nThere are two complications with the MDD model. First, t 0 may not be constant: clusters move from a dense environment where they are born to the normal ISM over a period that can extend to 100 Myr (and much longer in galaxy mergers), and over this time t 0 is likely to increase (Elmegreen & Hunter 2010; Kruijssen et al. 2011). Second, the simulations used to calibrate this model (e.g., Baumgardt & Makino 2003; Lamers, Baumgardt & Gieles 2010) have included only relaxation, not tidal losses due to encounters with GMCs, yet in the Galaxy the latter process is likely dominant. The characteristic destruction time for clusters due to GMCs is a few hundred Myr ( § 4.3.1), while for a 10 4 M glyph[circledot] cluster the disruption time due to relaxation (confirmed by simulations - Lamers et al. 2005) is multiple Gyr (Equation 21). The time scale for destruction by encounters with GMCs depends only on the density of the cluster, so agreement between the observed and simulated values of γ occurs only if ρ ∝ M 0 . 6 -0 . 7 - possible given the uncertainties in Figure 9, but by no means certain.  \nObservational evidence has been marshalled on both sides of this debate, which is closely tied to the issue of inclusive versus exclusive catalog construction discussed in § 2.1.2 and § 2.3. Equation 27 predicts a cluster age distribution with constant α T , and thus is 'universal', while Equation 28 predicts a value of α T that smoothly varies from 0, for γT/ ( t 0 M γ ) glyph[lessmuch] 1, to -( γ +1) /γ , for the opposite limit. Thus authors who construct exclusive catalogs, which tend to produce α T close to zero at early ages and then lower values at older ages, generally favor MDD models. Authors who construct inclusive catalogs, which generally have α T close to -1 at all ages, tend to favor MID models.  \nIn the latter category, Whitmore, Chandar & Fall (2007) and Fall, Chandar & Whitmore (2009) show that the data on clusters in the Antennae are well described by a power law of the form in Equation 27 and suggest that this is a 'quasi-universal law' of cluster destruction. Fall & Chandar (2012) extended this analysis to clusters in several other galaxies, including the SMC, the LMC, and M83; in all cases, dN/dM and dN/dT exhibit power law behavior with indexes α M ≈ -2 and α T ≈ -0 . 8. However, for Milky Way OCs, α T glyph[similarequal] -0 . 6 to -0 . 5 over the age range 10 -300 Myr (Figure 6), inconsistent with universality. Fall & Chandar 2012 also search for deviations from power-law behavior in the LMC as predicted by the MDD model, by finding clusters with both γT > t 0 M γ and γT < t 0 M γ . They find that simple MDD models fail for the LMC. However, they did not include stellar mass loss, which is often included in MDD models (e.g., Lamers, Baumgardt & Gieles 2010), allow for a time-dependent destruction time (e.g., Elmegreen & Hunter 2010), or consider the possibility that some of the youngest clusters may be unbound and therefore dissipate rapidly.  \nSupporting the MDD picture, Lamers et al. (2005) show that the cluster population within 600 pc of the Sun (Kharchenko et al. 2005) is consistent with dN/dT calculated from Equation 28 with γ = 0 . 62 and a 1.0 Gyr destruction time for 10 4 M glyph[circledot] clusters; the data cannot be fit by a power law. Bastian et al. (2012)  \ngeneralize the cluster evolution model to include both a mass-dependent and mass-independent destruction term, and find that MID and MDD models give equally good fits to two fields in M83. Adamo et al. (2017) find evidence for mass-dependent destruction of clusters with M < 10 4 M glyph[circledot] in NGC 628, and Messa et al. (2018b) report a similar result for M51. In both cases, however, the evidence comes from the lowest mass bin, so deeper observations or statistical methods for making better use of partially-complete data (e.g., Krumholz et al. 2018) would be beneficial.  \nTheoretically, the key difference between the two models lies in the mass dependence of the cluster disruption time, t dis . As discussed in § 4.3.1, the dominant destruction process for disk clusters is likely tidal shocking by GMCs, yielding disruption time proportional to the cluster density, t dis ∝ ρ . If ρ is approximately independent of mass, then so is t dis and one gets the MID model. Mass-independent densities are in fact expected for tidally limited clusters (Fall & Chandar 2012), but the fact that observed young clusters have EFF rather than King profiles ( § 4.1) and (at least potentially) underfill their tidal radii (Ryon et al. 2015, 2017) suggests they have not yet reached equilibrium. On the other hand, advocates of the MDDmodel (e.g., Lamers et al. 2005) assume ρ ∝ M 0 . 6 -0 . 7 , so that the mass dependence of tidal disruption is similar to that for 2-body relaxation. Their models can therefore implicitly include tidal disruption by GMCs, although not with the correct rate. Unfortunately, cluster densities are difficult to determine ( § 2.5), so to date it has not been possible to use measurements of ρ vs. M to distinguish the two models.  \nMiholics, Kruijssen & Sills (2017) carry out simulations of isolated galaxies using a subgrid model for the formation and evolution of star clusters (Kruijssen et al. 2011, 2012b); the model includes both relaxation and encounters with GMCs. They demonstrate that, in the MDD model, the median age of clusters in a sample is proportional to t 0 , and the median age depends on the surface density and velocity dispersion of the gas. Comparing with observations of M31 (outer galaxy), M51 (three mass ranges), and M83 (inner and outer regions), they find the best fit for the mass dependence of the median age is t med ∝ M 0 . 54 . In all their models, dN/dT decreases with age even though the star formation rates are roughly constant, and they are able to fit the observed median ages for the six cluster samples to within 0.2 dex, ranging from 10 7 . 3 yr for clusters in M51 with masses greater than 10 3 M glyph[circledot] to 10 8 . 3 yr for clusters in the outer parts of M31 and M83. This more complete model for cluster evolution is more complex than the simple MID and MDD models discussed above, and it demonstrates that cluster destruction is environmentally dependent.  \n4.4.2. Non-disk clusters. All star clusters, including the progenitors of GCs, likely began their lives in the disk of some galaxy (though see Mandelker et al. 2018 for a model where this is not true). Thus clusters found outside galaxy disks necessarily have complex histories. Studies of these histories therefore focus on determining possible origin sites for the clusters we see today, on the processes by which those clusters reached their present locations. A great deal of this work has been devoted to modeling the cosmological assembly history of the Milky Way's GC system, and to explaining observed features such as its color bimodality and color-magnitude relationships (the 'blue tilt'; Brodie & Strader 2006). These topics are beyond the scope of our review, which we limit to processes affecting star clusters in general, rather than the specific histories of GCs.  \nSince the transfer of clusters out of disks involves strong tides, simulations of the origin of non-disk clusters must model tidal effects. Since galactic and cosmological simulations cannot resolve cluster interiors ( § 3.5.3), they cannot do so directly, and instead the usual approach is to treat each star cluster as a point mass that is either formed self-consistently in the simulation or inserted following some prescribed initial distribution. As the particle moves along its orbit, one records the tidal potential to which it is subject, and uses this to calculate cluster evolution from an analytic model (e.g., Prieto & Gnedin 2008; Kruijssen et al. 2011; Brockamp et al. 2014; Rossi & Hurley 2015b; Carlberg 2018; Pfeffer et al. 2018) or via a direct N -body simulation (e.g., Renaud & Gieles 2013; Rieder et al. 2013; Rossi & Hurley 2015a; Mamikonyan  \net al. 2017). The latter approach is more accurate, but due to its expense one can only simulate a few examples, and for less than a full Hubble time.  \nOne obvious potential site for the formation of non-disk clusters is in major galaxy mergers such as the Antennae. Renaud & Gieles (2013) and Kruijssen et al. (2012b) both simulate Antennae-like mergers, but reach opposite conclusions: Renaud & Gieles find that the merger has little effect on the cluster population, while Kruijssen et al. find that it transforms the CMF from Schechter-like to lognormal like. The Kruijssen et al. results are likely more reliable. They include gas and form clusters self-consistently during the merger, while Renaud & Gieles carry out pure N -body simulations and study only pre-existing clusters; both differences produce weaker tides in Renaud & Gieles's simulations. Renaud & Gieles also initialize their clusters with much smaller radii, making them more resistant to tidal disruption.  \nDwarf galaxies represent another potential formation site for clusters that end up in the halos of larger galaxies such as the Milky Way: clusters form in the dwarf (and are therefore metal-poor, as observed), and then become part of the halo when the larger galaxy accretes and tidally strips the dwarf. As discussed in § 4.3.2, there are multiple Milky Way clusters that may have originated in this way. Simulations suggest that the accretion process is gentle, so accreted clusters do not experience major tidal losses, and quickly adjust to their new environments and become indistinguishable, except in their abundances, from the pre-existing cluster population (Miholics, Webb & Sills 2014, 2016; Bianchini et al. 2015).  \nExploring the full GC population requires cosmological simulations. Unfortunately, it is not presently possible to use the same simulations to study both the formation of GCs and their subsequent evolution (see Forbes et al. 2018 for more discussion). Simulations that can resolve cluster formation (e.g., Kimm et al. 2016; Ricotti, Parry & Gnedin 2016; Li et al. 2017; Li, Gnedin & Gnedin 2018; Kim et al. 2018) are too expensive to run to redshift zero, so studies that reach the present use potential GC progenitors inserted by hand. Different authors make different assumptions about how this should be done. Prieto & Gnedin (2008) and Carlberg (2018) use dark matter-only simulations; the former initialize cluster particles based on an analytic model of the baryonic mass distribution, while the latter adopt a randomly-oriented disk. They also assume different redshift distributions. Renaud, Agertz & Gieles (2017) and Pfeffer et al. (2018) use dark matter plus baryon cosmological simulations, with Renaud, Agertz & Gieles assuming that clusters follow the distribution of all newly-formed stars and Pfeffer et al. using the analytic models of Kruijssen (2012) and Reina-Campos & Kruijssen (2017) to vary Γ and the ICMF depending on local ISM conditions.  \nThe conclusions are as diverse as the sets of initial conditions: Prieto & Gnedin (2008) and Pfeffer et al. (2018) find that tides and two-body relaxation are able to transform an initially-power-law distribution of cluster masses into a lognormal, consistent with the observed globular CMF. Renaud, Agertz & Gieles (2017) reach the opposite conclusion, finding that their test clusters experience tides too weak to produce such a radical change. The difference is likely because their recipe for cluster formation - clusters simply follow star formation, in contrast to Prieto & Gnedin's prescribed distribution or Pfeffer et al.'s variable Γ and ICMF - places clusters on orbits where they experience weaker tides. Similarly, Carlberg (2018) finds that the effectiveness of tides is sensitive to the assumed distribution of redshifts for cluster formation. As for the ICMF ( § 3.5.3), at present simulation outcomes appear to depend primarily on subgrid recipes, and thus have little predictive power. Fortunately, we are likely to gain far better observational constraints once JWST launches, since it should narrow down the range of plausible models considerably.", '5. AFTERLIFE': 'The clusters we observe today are the surviving remnant of a much larger population that have now dissolved. However, within the Milky Way it is possible to reconstruct some of these clusters via a variety  \nof techniques. Doing so provides important clues about the cluster formation and dispersal process, the Galactic potential and, for clusters in the halo, the nature of dark matter. We therefore focus this final section on cluster reconstruction using kinematic ( § 5.1), action-angle ( § 5.2), and chemical tagging methods ( § 5.3). This field is in its infancy, so the discussion is necessarily more speculative than that up to this point.', '5.1. Kinematics': "The traditional method for cluster reconstruction is to find 'moving groups': clusters of stars in velocity space. 19 Clustering in velocity space is useful because unbound stars disperse in position on a timescale of a crossing time (a few Myr), but only spread out in velocity over timescales of a substantial fraction of their orbital period in the Galaxy (tens of Myr). This technique has a long and checkered history (Griffin 1998) with some early identified groups being questioned in later analysis (q.v. Taylor 2000), reminiscent of the 'faux clusters' identified in Gaia DR2 (Cantat-Gaudin et al. 2018). Nonetheless, it remains a powerful method for finding disrupted clusters.  \nWithin the disk, Dehnen (1998) first used the Hipparcos survey to identify kinematic substructure in the ( U, V ) plane for a very local sample ( D glyph[lessorsimilar] 100 pc); he identified known star clusters, but also discovered the previously-unknown Hercules stream. Bensby et al. (2007) showed that this structure is inhomogeneous in its abundances, and therefore more likely to be stars trapped in a resonance (e.g., Quillen et al. 2018a) than a disrupting cluster. Subsequent wavelet analysis of the Hipparcos survey combined with the RAVE survey identified many more local clumps in the ( U, V ) plane, although they have proven difficult to interpret (Antoja et al. 2012). While resonant streams are interesting, they are 'false positives' from the standpoint of cluster reconstruction. One major challenge for kinematic searches is that the expected rate of such false positives is unknown. The 2D simulations of De Simone, Wu & Tremaine (2004) suggest it is high, but we suspect that it will be lower in a 3D simulation due to the extra degree of freedom. For objects that are real disrupted clusters, one can use dynamical traceback methods to determine the age and other properties of the original cluster (de Zeeuw et al. 1999; Riedel et al. 2017).  \nThe halo has a great deal of substructure, but most known structures are associated with the disrupting Sgr dwarf (Helmi 2008). Exceptions include cold halo streams like those extending from the GC Pal 5 ( § 4.2.2). There are similar linear structures such as GD-1 (Grillmair & Johnson 2006) and the 'Jet' stream (Jethwa et al. 2018) not associated with a surviving progenitor; these are likely clusters that have fully disrupted (e.g., Koposov, Rix & Hogg 2010). Some halo substructures of comparable mass have only been identified kinematically (Kepley et al. 2007) or, for an assumed Galactic potential, by converting from kinematic coordinates to energy-angular momentum (Helmi et al. 2017) or angular momentum-eccentricity space (Helmi et al. 2006). Since Gaia DR2, new studies (q.v. Antoja et al. 2018) reveal many more substructures in the ( U, V ) plane, but follow-up work to identify them has only just begun.", '5.2. Action-angle space': "Kinematic reconstruction fails at ages glyph[greaterorsimilar] 100 Myr because differential acceleration along stars' orbits makes them disperse in velocity space. However, while velocities change as stars orbit, there are conserved quantities that do not, at least to the extent that we can approximate the Galactic potential as smooth and", 'DISRUPTED STAR CLUSTERS AS PROBES': "The coldest streams from disrupting star clusters are proving to be useful probes of the shape and enclosed mass of the Galaxy (Bowden, Belokurov & Evans 2015; Bovy et al. 2016). Recent observations confirm that GD-1 has high-contrast gaps along its length (Price-Whelan & Bonaca 2018). Several authors have argued that these are useful probes of the halo's dark matter substructure, since substructures would cause diffusion of stellar orbits that would blur out sharp features over time (Carlberg 2016; Bovy, Erkal & Sanders 2017).  \ntime-independent. In principle these quantities act as labels that are fixed when stars form, allowing reconstruction of clusters long after they have dispersed in velocity. This insight motivates the idea of cluster reconstruction in action-angle space ( J , θ ), a set of canonical conjugate coordinates defined by (Binney & Tremaine 2008)  \n2 πJ i = ∮ ˙ x i dx i = 1 Ω i ∫ 2 π 0 ( ˙ x ) 2 dθ i (29)  \nwhere x i are Cartesian coordinates and Ω i are the orbital frequencies. For each star, the equations of motion with respect to the Hamiltonian H are  \n˙ θ = ∂ H ∂ J = Ω( J ) ˙ J = -∂ H ∂ θ = 0 (30)  \nAn orbit can be traced with the solution θ i ( t ) = Ω i ( t ) + Ω i, 0 . For the inertial Galactic system, ( R,φ,z ) coordinates 20 are most often used (even for halo studies) such that θ = ( θ R , θ φ , θ z ) and J = ( J R , J φ , J z ) where J φ ≡ L z , i.e., J φ is the angular momentum with respect to the spin axis. A stellar orbit is described by the three actions J . Orbits with J z = 0 lie in the Galactic plane, orbits with J R = 0 are circular, and we choose units so that a circular orbit at the Solar Circle has J φ = 1; circular orbits with J φ < 1 lie inside the Solar Circle. The orbits of all stars in the Galaxy are described by the distribution function f ( J ).  \nTo find clusters in action-angle space one must convert the observed positions and velocities of stars into actions. This requires knowledge of the Galactic potential Φ, which enters the Hamiltonian H . We will not fully review efforts to measure the Galactic potential and the stellar density distribution, which must be determined simultaneously. Arguably, the best efforts to date are from the RAVE survey, which provides an explicit fit for the vertical stellar density distribution (Piffl et al. 2014) as opposed to the parameterized potentials used in earlier work (e.g. McMillan 2011). Moreover, Φ is at least weakly time-dependent. In addition to the outer warp, there is clear evidence near and far of large-scale corrugation waves propagating throughout the disk (e.g. Xu et al. 2015; Antoja et al. 2018). The spiral arms and bar provide additional challenges.  \nEven with the uncertainties, the action-angle approach reveals rich structure throughout the halo and Galactic disk. Sellwood (2010) mapped Hipparcos moving groups to action space and found that this improved the delineation of the Hyades stream, but that other structures were barely evident. With Gaia DR2, the situation has improved and over a much larger volume than Hipparcos . Trick, Coronado & Rix (2018) show that the highest density peaks in the ( U, V ) plane map to the highest density peaks in the action plane, specifically ( √ J R , J φ = L z ), but the mapping is not simple. They trace structure out to 1.5 kpc from the Sun and find much of it has low vertical action, J z , suggestive of disk resonances that operate  \nmost efficiently in the plane (Katz et al. 2018). The Helmi streams were discovered with a dozen or fewer stars using action analysis (Helmi 2008), structures that were invisible against the background in kinematic space. Arguably the most impressive action-based discovery to date is Gaia-Enceladus, the remnants of a massive dwarf spread throughout the inner Galaxy that was accreted long before Sgr (Helmi et al. 2018a). To date, no OCs or GCs (or their streams) have been discovered using actions, but these are early days. There are many clumps in this space that have yet to be followed up. Stars bunched in phase space are typically spread across the sky requiring extended observing programs over 1 -2 years.", '5.3. Chemical tagging': "Dynamical 'invariants' like actions are useful labels for reconstruction, but they are not strictly conserved over a star's lifetime, e.g., when transiting the spiral arm or bar. These time-dependent features of the Galactic potential induce perturbations in the actions of the stars they affect; indeed, non-conservation of J φ due to arms and bars is responsible for stellar migration ( § 4.3.2), and the same structures also likely cause non-conservation of J R (Solway, Sellwood & Schonrich 2012). This ultimately limits how far back in time action-angle reconstructions for clusters in the disk can go. Longer-term reconstructions therefore require another invariant, for which elemental abundances are a natural choice. Stars conserve their abundances over almost their entire lives, and stars in the same cluster have nearly-identical abundances ( § 2.6). This led Freeman & Bland-Hawthorn (2002) to suggest that, with abundance data of sufficient quality and enough elements, clusters that have dispersed even in action space could be detected via clustering in C -space, defined so each star is a point in a high dimensional space, i.e., C ([Fe/H], [ α /Fe], [X 1 /Fe], [X 2 /Fe], . . . ).  \nWeak chemical tagging, defined as identifying distinct components separated by a surface in C -space, has been used to good effect. For example, stars in the Solar neighborhood show a clear bimodality in C ([Fe/H], [ α /Fe]) (e.g., Bensby et al. 2005) and C ([ α /Fe], [(C+N)/Fe], [Al/Fe], [Mg/Mn]) (Hawkins et al. 2015). The Sgr dwarf and stream separate cleanly from all components of the Milky Way in the high-dimensional space C ([(C+N)/Fe], [O/Fe], [Mg/Fe], [Al/Fe], [Mn/Fe], [Ni/Fe]) (Hasselquist et al. 2017). Reconstructing star clusters, however, requires strong chemical tagging, defined as identifying a cluster in C -space against the background (q.v. Ting et al. 2012; Ting, Conroy & Goodman 2015). In searching for dissolved clusters, the choice of elements depends on the target. OCs are highly homogeneous so many elements are prospective candidates, while GCs are inhomogeneous in some or all elements, but have distinctive anti-correlations, e.g., O/Mg vs. Na/Al. Martell et al. (2016) use the latter signature to identify stars in the Galactic halo from long-dissolved GCs. Ting, Conroy & Goodman (2015) search for clustering of abundances in the APOGEE survey of the α -rich disk, but do not identify any previously-unknown clusters. Related techniques applied to the α -poor disk using the GALAH survey have met with more success (Quillen et al. 2015; Kos et al. 2018a). Even when clusters are found, however, there may not be a simple one-to-one mapping between C -space clusters and OCs. For example, De Silva et al. (2015) identify a pair of OCs that have identical radial velocities and overlapping stars in C -space. These were likely formed in close proximity, so abundance analysis does correctly identify real structures, but they are not necessarily the same clusters one would identify in physical space.  \nAs for kinematics, strong chemical tagging suffers from a poorly-constrained false positive rate. Ness et al. (2018) identify pairs of field stars with identical abundance patterns that are unlikely to have been born in the same cloud, but the overall frequency of such chance pairs remains uncertain. The problem should be less severe at low [Fe/H], where there are fewer stars and thus a reduced background of interlopers. This led Bland-Hawthorn et al. (2010) to propose applying chemical tagging to the Milky Way's dwarf satellites rather than in its disk. Prima facie, this appears to work (e.g., Karlsson et al. 2012; Webster, Frebel & Bland-Hawthorn 2016), and if confirmed the patterns identified thus far are evidence of star clusters significantly more metal-poor than the lowest metallicity GCs, which have [Fe/H] ≈ -2 . 4.  \nThe challenges of strong chemical tagging mean that, thus far, it has been most effective when used to characterize groups found by other methods, and to test for potential membership. For example, most stellar streams identified to date appear to be inhomogeneous, implying that they are not disrupting star clusters, and are instead either stars trapped by disc resonances (e.g., Bensby et al. 2007) or the remnants of much more massive ( > 10 8 M glyph[circledot] ) systems (e.g., Hasselquist et al. 2017). However, a few streams have proven to be homogeneous, e.g., the 2-Gyr old moving group HR 1614 (Feltzing & Holmberg 2000; De Silva et al. 2007), first identified by Eggen (1978) and confirmed with Hipparcos (Eggen 1998). For these systems, elemental homogeneity provides strong evidence that the progenitor is a disrupted star cluster. In the future, it seems likely that the most powerful applications of chemical tagging will continue to be in conjunction with other methods (e.g., Quillen et al. 2015) rather than by itself. Kos et al. (2018a) provide an impressive example of how such combinations can work, by combining chemical tagging and Gaia proper motions to identify Pleiades members several degrees from the parent cluster. To date, other identified streams for which chemical tagging may prove useful are faint, so follow-up work awaits multi-object high resolution spectrographs on the ELTs (e.g., MANIFEST on GMT).", '6. CONCLUSIONS AND FUTURE PROSPECTS': 'Star clusters stand at a crossroads of scales. Traditionally star cluster research has focused on clusters as discrete entities whose formation, evolution, and eventual dissolution can be viewed in isolation, at most treating the background galaxy as a source of a passive tidal field. In this paradigm, groups of stars can be neatly classified into open clusters, globular clusters, associations, and field stars, each with distinct properties and formation histories. Such isolation is no longer viable. Observations now reveal that the objects that go on to form bound star clusters are merely the innermost parts of a hierarchy that extends to the scales of galaxies. Once clusters form, the tidal perturbations to which they are subject are affected first by the immediate star-forming environment, and over longer times by the full cosmological history of galaxy assembly. Going to high redshift brings us to a point where the distinction between open clusters and globular clusters dissolves. Consequently any modern understanding of star clusters must bridge from sub-pc to cosmological scales. While many questions remain, the outlines of a model are starting to become clear:', 'SUMMARY POINTS': "- 1. Stars form in giant molecular clouds (GMCs) over a few free-fall times, at a rate per free-fall time glyph[epsilon1] ff ∼ 1%, with a spread of ∼ 0 . 5 dex. At formation, the stars are hierarchically-structured and cannot be separated neatly into clusters. Distinct star clusters emerge only after gas is cleared, which occurs in less than a few million years. After about 30 Myr in Milky Way-like galaxies, ∼ 10% of the stars are in clusters. There are good theoretical arguments that this fraction should vary with galactic environment, but the evidence either for or against variation is thus far unconvincing.\n- 2. The portions of clouds that form bound clusters have at most mildly elevated glyph[epsilon1] ff values, but form stars over multiple free-fall times, allowing them to reach elevated total star formation efficiencies. Extended star formation renders these regions dynamically relaxed and well-mixed in abundance. In contrast, the stars that are not part of the bound regions never dynamically relax, and may or may not be well-mixed. They only gradually drift apart; they are not part of a coherent flow either collapsing towards or expanding away from the dense regions that become bound.\n- 3. Star formation is ultimately terminated by feedback, but which type of feedback is dominant depends on the local environment. Outflows dominate in clouds too small to contain massive stars,  \nand radiation pressure dominates in the most massive and dense systems, but for the bulk of star clusters photoionization is likely to be most important. Feedback sets the timescale for gas clearing, and thus determines the density at which star formation transitions from taking place over one to many dynamical times, marking the boundary between the bound and unbound parts of a hierarchical star-forming cloud. The initial cluster mass function (ICMF) of both bound and unbound clusters is well described by a power-law dN/dM ∝ M α M with α M ≈ -2, likely with an environmentally-dependent high-mass cutoff. Evidence for any stronger variation in cluster formation with environment is weak.  \n- 4. In the first ≈ 10 -100 Myr after a gas-free cluster population emerges, it is subject to mass loss via stellar evolution-driven expansion and tidal shocking by gas in its immediate environment. The rate at which the cluster age function (CAF) declines from ≈ 10 -100 Myr, as parameterized by its slope α T ( dN/dT ∝ T α T ), provides a powerful constraint on the importance of these processes. In the Solar neighborhood α T ≈ -0 . 5, implying moderately strong cluster destruction, but the CAF slope outside the Milky Way is uncertain because present samples are compromised by uncertainties in cluster selection and age assignment. The strength of early mass loss mechanisms depends strongly on the distributions of cluster concentration and ratio of size to tidal radius at the end of gas clearing, which are poorly constrained both theoretically and observationally, though there are observational hints that clusters in Milky Way-like galaxies are born fairly concentrated and under-filling their tidal radii. There is tension between these hints and the Milky Way cluster age distribution, which requires moderately strong early cluster disruption.\n- 5. Over timescales of Gyr, the dominant cluster destruction mechanisms are relaxation and tidal shocking (by GMCs for disk clusters, by disk traversal for halo clusters). Relaxation preferentially destroy low-mass clusters; tidal shocking likely does as well, though this depends on how strongly density increases with mass. These mechanisms thus provide a plausible mechanism for transforming the power-law mass distribution ubiquitously observed for young clusters into the peaked distribution observed for globular clusters. However, enhanced tidal shocking during clusters' first ≈ 100 Myr of life also preferentially destroys low-mass clusters, and potentially provides an alternate explanation. The relative importance of the various processes, and the timescale over which the CMF shape changes, likely depends on the galactic environment.\n- 6. Even after clusters come apart in physical space, the stars remain coherent in kinematic, actionangle, and chemical space for timescales of tens, hundreds, and thousands of Myr, respectively. In principle it should be possible to reconstruct clusters in these spaces. However, efforts to do so are still in their infancy, and most present techniques still suffer from high rates of false positives.", 'FUTURE ISSUES': '- 1. A recurring theme of this review is the uncertainty in inferences of various quantities - cutoffs in the CMF, cluster ages, etc. - based on integrated light. When one cross-checks these results against more accurate methods based on resolved stellar populations and young stellar objects (YSOs), the level of agreement is often worse than a naive interpretation of the stated error bars would suggest. However, integrated light observations will remain a vital tool for the foreseeable future, since even in the era of JWST we will not be able to resolve stellar populations in clusters at even ∼ 100 Mpc distances, let alone at high redshift. Consequently, there is an urgent need for a program of checking integrated light techniques against resolved stellar populations, and developing new methods that  \nreturn realistic error bars.  \n- 2. There is a pressing need to extend extragalactic cluster catalogs to lower mass and larger age. At present the Milky Way and extragalactic samples are almost completely non-overlapping in mass, and the extragalactic samples have very limited dynamic range in mass, making it hard to draw conclusions. The best prospects for improving this situation are in pushing to lower mass clusters in galaxies at distances of a few Mpc. This will become possible with extreme AO-fed instrumentation on the next generation of extremely large telescopes. The same instruments will also make it possible to probe star formation in cosmologically denser regimes compared to the low overdensities in the Local Volume.\n- 3. Both analytic models and simulations of star cluster formation that begin from isolated clouds, ignoring the galactic context, are likely reaching the end of their utility. They consistently fail to reproduce observed star formation histories and stellar kinematic, and are likely incapable of correctly predicting the cluster formation efficiency overall. The future is in simulations that start self-consistently from the galactic or cosmological scale, including all relevant physics, and then zoom in to individual clusters, or that do form and evolve clusters using semi-analytic models calibrated to zoom simulations.\n- 4. The largest theoretical uncertainties in cluster demographics now lie at the transition between gasdominated and gas-free evolution. Pure N -body simulations of the gas-free phase have advanced to the point where it is now possible to simulate a cluster of 10 6 stars directly, but a continuing effort to push these simulations to ever-larger numbers of stars and longer integrations may have hit the point of diminishing returns because the largest uncertainty is the choice of initial conditions, which cannot be addressed without modeling the gas.\n- 5. Observations of 6D phase space from ESA Gaia and ESO Gravity may address the issue of the poorly known initial conditions, since they can simultaneously supply initial conditions to the N -body modelers and provide a benchmark against which gas simulators can test their work. This project will, however, require new statistical tools. Statistics like the CMF are problematic for unrelaxed populations that cannot be uniquely decomposed into clusters, while structure-agnostic statistics such as two point correlation functions are less informative than statistics like the CMF when applied to populations of relaxed, discrete clusters. Observed star-forming complexes such as Orion, however, appear to transition smoothly between smooth and relaxed at high density and fractal and unrelaxed at low density. We need statistical tools that can characterize this transition and work across it.\n- 6. The most powerful future methods for reconstructing disrupted clusters are likely to combine elemental abundance and kinematic data. Each method suffers from significant false positives by itself, but the false positives to which each are subject are quite different, and thus the two methods can be used to cross-check one another. However, this will require the availability of abundance measurements going to fainter magnitudes than is now possible. Ground-layer, AO-corrected wide-field spectrographs planned for the next generation of telescopes will allow detailed abundance analysis on stars to fainter limits ( V ≈ 20) than possible today, greatly aiding this effort.', 'DISCLOSURE STATEMENT': 'The authors are not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.', 'ACKNOWLEDGMENTS': "We thank A. Adamo, N. Bastian, R. Chandar, E. van Dishoeck, B. Elmegreen, E. Falgarone, M. Fall, D. Kruijssen, E. Ostriker, and A. Sternberg for helpful comments on manuscript. MRK acknowledges support from the Australian Research Council's (ARC) Discovery Project and Future Fellowship funding schemes, awards DP160100695 and FT180100375. CFM acknowledges support by NASA through NASA ATP grant NNX13AB84G. JBH acknowledges support through the ARC Discovery Project and Laureate Fellowship funding schemes, awards DP150104667 and FL140100278. JBH acknowledges the hospitality of UC Berkeley and thanks the Miller Foundation for their generous support of the 2018 Visiting Professorship. MRK and JBH acknowledge support from an ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.", 'LITERATURE CITED': "```\nAdamo A, Bastian N. 2018. The Lifecycle of Clusters in Galaxies . In The Birth of Star Clusters , ed. S Stahler, vol. 424 of Astrophysics and Space Science Library 1.1 Adamo A, Kruijssen JMD, Bastian N, Silva-Villa E, Ryon J. 2015. MNRAS 452:246-260 2.2, 5, 8 Adamo A, Ostlin G, Zackrisson E. 2011. MNRAS 417:1904-1912 8, 2.4 Adamo A, Ostlin G, Zackrisson E, Hayes M, Cumming RJ, Micheva G. 2010. MNRAS 407:870-890 2.1.2 Adamo A, Ryon JE, Messa M, Kim H, Grasha K, et al. 2017. ApJ 841:131 4, 2.1.2, 2.2, 5, 6, 2.3, 4.4.1 Adams FC. 2000. ApJ 542:964-973 3.5.2 Adams FC, Myers PC. 2001. ApJ 553:744-753 3.5.1 Allison RJ, Goodwin SP, Parker RJ, de Grijs R, Portegies Zwart SF, Kouwenhoven MBN. 2009. ApJ 700:L99-L103 3.5.2 Andrews JE, Calzetti D, Chandar R, Elmegreen BG, Kennicutt RC, et al. 2014. ApJ 793:4 3.3.1 Antoja T, Helmi A, Bienayme O, Bland-Hawthorn J, Famaey B, et al. 2012. MNRAS 426:L1-L5 5.1 Antoja T, Helmi A, Romero-Gomez M, Katz D, Babusiaux C, et al. 2018. ArXiv e-prints 5.1, 5.2 Armillotta L, Krumholz MR, Fujimoto Y. 2018. MNRAS 481:5000-5013 3.6, 3.6 Armus L, Mazzarella JM, Evans AS, Surace JA, Sanders DB, et al. 2009. PASP 121:559 1.1 Ascenso J, Alves J, Lago MTVT. 2009. A&A 495:147-155 2.1.1 Astropy Collaboration, Robitaille TP, Tollerud EJ, Greenfield P, Droettboom M, et al. 2013. A&A 558:A33 1 Azimlu M, Mart'ınez-Galarza JR, Muench AA. 2015. AJ 150:95 3.4.1 Ballesteros-Paredes J, Hartmann LW, V'azquez-Semadeni E, Heitsch F, Zamora-Avil'es MA. 2011. MNRAS 411:65-70 3.2.2 Bally J. 2016. Annual Review of Astronomy and Astrophysics 54:491-528 3.3.1 Baraffe I, Chabrier G, Gallardo J. 2009. ApJ 702:L27-L31 3.4.1 Baraffe I, Vorobyov E, Chabrier G. 2012. ApJ 756:118 3.4.1 Barmby P, McLaughlin DE, Harris WE, Harris GLH, Forbes DA. 2007. AJ 133:2764-2786 2, 9 Barnes AT, Longmore SN, Battersby C, Bally J, Kruijssen JMD, et al. 2017. MNRAS 469:2263-2285 9 Bastian N. 2008. MNRAS 390:759-768 2.2, 2.4 Bastian N, Adamo A, Gieles M, Lamers HJGLM, Larsen SS, et al. 2011. MNRAS 417:L6-L10 2.1.2 Bastian N, Adamo A, Gieles M, Silva-Villa E, Lamers HJGLM, et al. 2012. MNRAS 419:2606-2622 2.1.2, 2.1.2, 2.2, 4.4.1 Bastian N, Ercolano B, Gieles M, Rosolowsky E, Scheepmaker RA, et al. 2007. MNRAS 379:1302-1312 2.7 Bastian N, Lardo C. 2018. Annual Review of Astronomy and Astrophysics 56:83-136 1.2.2, 2.3, 2.6 Baumgardt H. 2001. MNRAS 325:1323-1331 4.2.2 Baumgardt H, Hilker M. 2018. MNRAS 478:1520-1557 2, 9, 4.2.1 Baumgardt H, Makino J. 2003. MNRAS 340:227-246 4.2.1, 4.2.2, 4.4.1 Baumgardt H, Parmentier G, Anders P, Grebel EK. 2013. MNRAS 430:676-685 5 Baumgardt H, Parmentier G, Gieles M, Vesperini E. 2010. MNRAS 401:1832-1838 4.2.2\n```  \n```\nBeccari G, Petr-Gotzens MG, Boffin HMJ, Romaniello M, Fedele D, et al. 2017. A&A 604:A22 3.4.1 Beerman LC, Johnson LC, Fouesneau M, Dalcanton JJ, Weisz DR, et al. 2012. ApJ 760:104 2.1.1 Bekki K. 2010. MNRAS 401:2753-2762 4.3.2 Bensby T, Feltzing S, Lundstrom I, Ilyin I. 2005. A&A 433:185-203 5.3 Bensby T, Oey MS, Feltzing S, Gustafsson B. 2007. ApJ 655:L89-L92 5.1, 5.3 Bertoldi F, McKee CF. 1992. ApJ 395:140-157 3.1.1, 3.1.1 Bhattacharya S, Mishra I, Vaidya K, Chen WP. 2017. ApJ 847:138 4.2.2 Bianchini P, Renaud F, Gieles M, Varri AL. 2015. MNRAS 447:L40-L44 4.4.2 Billett OH, Hunter DA, Elmegreen BG. 2002. AJ 123:1454-1475 3.1.2 Binney J, Tremaine S. 2008. Galactic Dynamics: Second Edition . Princeton University Press 4.3.1, 4.3.1, 5.2 Birnboim Y, Federrath C, Krumholz M. 2018. MNRAS 473:2144-2159 3.2.2 Bisnovatyi-Kogan GS, Silich SA. 1995. Reviews of Modern Physics 67:661-712 3.3.5 Bland-Hawthorn J, Freeman KC. 2004. PASA 21:110-120 3.6 Bland-Hawthorn J, Gerhard O. 2016. ARA&A 54:529-596 19 Bland-Hawthorn J, Karlsson T, Sharma S, Krumholz M, Silk J. 2010. ApJ 721:582-596 5.3 Bland-Hawthorn J, Krumholz MR, Freeman K. 2010. ApJ 713:166-179 3.6 Boesgaard AM, Lum MG, Deliyannis CP. 2015. ApJ 799:202 4.3.2 Bolatto AD, Leroy AK, Rosolowsky E, Walter F, Blitz L. 2008. ApJ 686:948-965 3.1.1, 3.1.1 Borissova J, Georgiev L, Hanson MM, Clarke JRA, Kurtev R, et al. 2012. A&A 546:A110 3 Borissova J, Ivanov VD, Hanson MM, Georgiev L, Minniti D, et al. 2008. A&A 488:151-159 3 Boutloukos SG, Lamers HJGLM. 2003. MNRAS 338:717-732 2.3, 4.4.1, 4.4.1, 4.4.1 Bovy J, Bahmanyar A, Fritz TK, Kallivayalil N. 2016. ApJ 833:31 5.1 Bovy J, Erkal D, Sanders JL. 2017. MNRAS 466:628-668 5.1 Bowden A, Belokurov V, Evans NW. 2015. MNRAS 449:1391-1400 5.1 Bragaglia A, Carretta E, D'Orazi V, Sollima A, Donati P, et al. 2017. A&A 607:A44 2.6 Bragaglia A, Sneden C, Carretta E, Gratton RG, Lucatello S, et al. 2014. ApJ 796:68 2.6 Bressert E, Bastian N, Gutermuth R, Megeath ST, Allen L, et al. 2010. MNRAS 409:L54-L58 2.1.1, 2.7 Brockamp M, Kupper AHW, Thies I, Baumgardt H, Kroupa P. 2014. MNRAS 441:150-171 4.4.2 Brodie JP, Strader J. 2006. ARA&A 44:193-267 2.2, 4.4.2 Brunetti M, Chiappini C, Pfenniger D. 2011. A&A 534:A75 4.3.2 Brunt CM, Federrath C, Price DJ. 2010. MNRAS 403:1507-1515 3.1.2 Burkhart B. 2018. ApJ 863:118 3.2.2 Burkhart B, Mocz P. 2018. ApJ Submitted, arXiv:1805.11104 3.2.2 Caldwell S, Chang P. 2018. MNRAS 474:4818-4823 3.4.1 Calura F, Few CG, Romano D, D'Ercole A. 2015. ApJ 814:L14 3.3.5, 3.3.6 Calzetti D, Armus L, Bohlin RC, Kinney AL, Koornneef J, Storchi-Bergmann T. 2000. ApJ 533:682-695 3 Calzetti D, Lee JC, Sabbi E, Adamo A, Smith LJ, et al. 2015. AJ 149:51 1.1, 1.2.1 Cantat-Gaudin T, Jordi C, Vallenari A, Bragaglia A, Balaguer-N'u˜nez L, et al. 2018. ArXiv e-prints 4.3.1, 5.1 Carlberg RG. 2016. ApJ 820:45 5.1 Carlberg RG. 2018. ApJ 861:69 4.4.2 Carraro G, Bensby T. 2009. MNRAS 397:L106-L109 4.3.2 Carraro G, Villanova S, Demarque P, McSwain MV, Piotto G, Bedin LR. 2006. ApJ 643:1151-1159 4.3.2 Carretta E, Bragaglia A, Gratton RG, Recio-Blanco A, Lucatello S, et al. 2010. A&A 516:A55 2.6 Chabrier G. 2005. The Initial Mass Function: from Salpeter 1955 to 2005 . In The Initial Mass Function 50 Years Later , eds. E Corbelli, F Palla, H Zinnecker, vol. 327 of Astrophysics and Space Science Library . Dordrecht: Springer 3, 3.3.1 Chandar R, Fall SM, Whitmore BC. 2015. ApJ 810:1 5 Chandar R, Fall SM, Whitmore BC, Mulia AJ. 2017. ApJ 849:128 6, 8, 2.4 Chandar R, Whitmore BC, Calzetti D, Di Nino D, Kennicutt RC, et al. 2011. ApJ 727:88 2.1.2 Chandar R, Whitmore BC, Calzetti D, O'Connell R. 2014. ApJ 787 2.1.2\n```  \n```\nElmegreen BG, Efremov YN. 1997. ApJ 480:235-245 1.2.2 Elmegreen BG, Elmegreen DM. 2001. AJ 121:1507-1511 2.7 Elmegreen BG, Falgarone E. 1996. ApJ 471:816 3.5.1.1 Elmegreen BG, Hunter DA. 2010. ApJ 712:604-623 4.3.1, 4.4.1 Elson RA, Fall SM. 1988. AJ 96:1383 2.1.2, 2.4 Elson RAW, Fall SM. 1985. ApJ 299:211-218 2.1.2, 2.4 Elson RAW, Fall SM, Freeman KC. 1987. ApJ 323:54 1, 4.1 Elson RAW, Freeman KC, Lauer TR. 1989. ApJ 347:L69 4.2.1 Emerick A, Bryan GL, Mac Low MM. 2019. MNRAS 482:1304-1329 3.6 Escala A, Larson RB. 2008. ApJ 685:L31 3.1.2, 3.1.2 Evans II NJ, Heiderman A, Vutisalchavakul N. 2014. ApJ 782:114 3.2.1, 10 Faesi CM, Lada CJ, Forbrich J. 2018. ApJ 857:19 3.1.1, 3.1.2 Falgarone E, Pety J, Hily-Blant P. 2009. A&A 507:355-368 3.1.1 Fall SM, Chandar R. 2012. ApJ 752:96 5, 4.2.1, 4.4.1 Fall SM, Chandar R, Whitmore BC. 2005. ApJ 631:L133-L136 4.4.1 Fall SM, Chandar R, Whitmore BC. 2009. ApJ 704:453-468 4.4.1, 4.4.1, 4.4.1 Fall SM, Krumholz MR, Matzner CD. 2010. ApJ 710:L142-L146 3.3, 3.3.3, 3.3.6, 3.5.1.1, 3.5.1.2 Fall SM, Zhang Q. 2001. ApJ 561:751-765 2.2, 4.2.2, 4.4.1 Faucher-Gigu'ere CA, Kaspi VM. 2006. ApJ 643:332-355 4.2.1 Federrath C. 2013. MNRAS 436:3167-3172 3.2.1 Federrath C. 2015. MNRAS 450:4035-4042 3.2.2 Federrath C, Klessen RS. 2012. ApJ 761:156 3.2.2 Feltzing S, Holmberg J. 2000. A&A 357:153-163 5.3 Feng Y, Krumholz MR. 2014. Nature 513:523-525 2.6, 3.6 Fierlinger KM, Burkert A, Ntormousi E, Fierlinger P, Schartmann M, et al. 2016. MNRAS 456:710-730 3.3.5 Figer DF, MacKenty JW, Robberto M, Smith K, Najarro F, et al. 2006. ApJ 643:1166-1179 3 Figer DF, McLean IS, Morris M. 1999. ApJ 514:202-220 3 Fleck Jr. RC. 1996. ApJ 458:739 3.5.1.1 Forbes DA, Bastian N, Gieles M, Crain RA, Kruijssen JMD, et al. 2018. Proceedings of the Royal Society of London Series A 474:20170616 1.1, 4.4.2 Fouesneau M, Johnson LC, Weisz DR, Dalcanton JJ, Bell EF, et al. 2014. ApJ 786:117 6, 2.4, 2, 9 Fouesneau M, Lan¸con A. 2010. A&A 521:A22 2.1.2 Fouesneau M, Lan¸con A, Chandar R, Whitmore BC. 2012. ApJ 750:60 2.1.2 Freeman K, Bland-Hawthorn J. 2002. ARA&A 40:487-537 3.6, 5.3 Friel ED. 1995. Annual Review of Astronomy and Astrophysics 33:381-414 4.3.2 Fukushige T, Heggie DC. 1995. MNRAS 276:206-218 4.2.1 Fukushige T, Heggie DC. 2000. MNRAS 318:753-761 4.2.2 Gallagher MJ, Leroy AK, Bigiel F, Cormier D, Jim'enez-Donaire MJ, et al. 2018. ApJ 858:90 3.2.1, 10 Gao Y, Solomon PM. 2004. ApJ 606:271-290 10 Garc'ıa-Burillo S, Usero A, Alonso-Herrero A, Graci'a-Carpio J, Pereira-Santaella M, et al. 2012. A&A 539:A8 3.2.1 Gatto A, Walch S, Naab T, Girichidis P, Wunsch R, et al. 2017. MNRAS 466:1903-1924 3.3.5 Gavagnin E, Bleuler A, Rosdahl J, Teyssier R. 2017. MNRAS 472:4155-4172 3.3.2 Gazak JZ, Davies B, Bastian N, Kudritzki R, Bergemann M, et al. 2014. ApJ 787:142 2.4 Geen S, Hennebelle P, Tremblin P, Rosdahl J. 2015. MNRAS 454:4484-4502 3.3.2 Geen S, Hennebelle P, Tremblin P, Rosdahl J. 2016. MNRAS 463:3129-3142 3.3.2, 3.3.6 Geisler D, Villanova S, Carraro G, Pilachowski C, Cummings J, et al. 2012. ApJ 756:L40 2.6 Gennaro M, Brandner W, Stolte A, Henning T. 2011. MNRAS 412:2469-2488 3 Gentry ES, Krumholz MR, Dekel A, Madau P. 2017. MNRAS 465:2471-2488 3.3.6 Gentry ES, Krumholz MR, Madau P, Lupi A. 2018. MNRAS Submitted, arXiv:1802.06860 3.3.5, 3.3.6 Genzel R, Newman S, Jones T, Forster Schreiber NM, Shapiro K, et al. 2011. ApJ 733:101 3.1.2\n```  \n```\nGerhard O. 2001. ApJ 546:L39-L42 4.3.2 Getman KV, Feigelson ED, Kuhn MA. 2014. ApJ 787:109 3.4.1, 3.4.2 Getman KV, Feigelson ED, Kuhn MA, Bate MR, Broos PS, Garmire GP. 2018. MNRAS 476:1213-1223 2.1.1, 3.4.1, 3.4.2, 3.4.2 Gieles M. 2009. MNRAS 394:2113-2126 4.4.1 Gieles M, Athanassoula E, Portegies Zwart SF. 2007. MNRAS 376:809-819 4.3.1, 4.3.2 Gieles M, Baumgardt H. 2008. MNRAS 389:L28-L32 4.2.2 Gieles M, Baumgardt H, Heggie DC, Lamers HJGLM. 2010. MNRAS 408:L16-L20 4.2.1 Gieles M, Heggie DC, Zhao H. 2011. MNRAS 413:2509-2524 4.2.2 Gieles M, Portegies Zwart SF. 2011. MNRAS 410:L6-L7 1.2.1 Gieles M, Portegies Zwart SF, Baumgardt H, Athanassoula E, Lamers HJGLM, et al. 2006. MNRAS 371:793-804 4.3.1, 4.3.1, 4.3.1 Gieles M, Renaud F. 2016. MNRAS 463:L103-L107 4.3.1 Ginsburg A, Bressert E, Bally J, Battersby C. 2012. ApJ 758:L29 3.4.2 Ginsburg A, Kruijssen JMD. 2018. ApJ 864:L17 8, 2.4 Gnedin OY, Ostriker JP. 1997. ApJ 474:223-255 4.3.1 Gnedin OY, Ostriker JP. 1999. ApJ 513:626-637 4.3.1 Goddard QE, Bastian N, Kennicutt RC. 2010. MNRAS 405:857-869 8, 2.4, 2.4 Goldbaum NJ, Krumholz MR, Matzner CD, McKee CF. 2011. ApJ 738:101-+ 3.2.2, 3.4.2 Gonz'alez-Fern'andez C, Negueruela I. 2012. A&A 539:A100 3 Goudfrooij P, Fall SM. 2016. ApJ 833:8 4.2.2 Gouliermis DA. 2018. PASP 130:072001 1.2.2, 2.7, 3.6 Gouliermis DA, Elmegreen BG, Elmegreen DM, Calzetti D, Cignoni M, et al. 2017. MNRAS 468:509-530 2.7 Grasha K, Calzetti D, Adamo A, Kim H, Elmegreen BG, et al. 2017a. ApJ 840:113 2.7, 4.3.1 Grasha K, Elmegreen BG, Calzetti D, Adamo A, Aloisi A, et al. 2017b. ApJ 842:25 2.7 Gratton RG, Carretta E, Bragaglia A. 2012. Asron. Astrop. Rev. 20:50 2.6 Griffin RF. 1998. The Observatory 118:223-225 5.1 Grillmair CJ, Johnson R. 2006. ApJ 639:L17-L20 4.2.2, 5.1 Grudi'c MY, Hopkins PF, Faucher-Gigu'ere CA, Quataert E, Murray N, Kereˇs D. 2018. MNRAS 475:3511-3528 3.2.2, 3.3.3, 3.3.7, 3.4.2, 14 Guszejnov D, Hopkins PF, Grudi'c MY. 2018. MNRAS 477:5139-5149 3.5.1.1 Haid S, Walch S, Seifried D, Wunsch R, Dinnbier F, Naab T. 2018. MNRAS 478:4799-4815 3.3.5 Harayama Y, Eisenhauer F, Martins F. 2008. ApJ 675:1319-1342 3 Harper-Clark E, Murray N. 2009. ApJ 693:1696-1712 3.3.5 Harris WE, Pudritz RE. 1994. ApJ 429:177-191 3.1.2 Hasselquist S, Shetrone M, Smith V, Holtzman J, McWilliam A, et al. 2017. ApJ 845:162 5.3 Hawkins K, Jofr'e P, Masseron T, Gilmore G. 2015. MNRAS 453:758-774 5.3 Heiderman A, Evans II NJ, Allen LE, Huard T, Heyer M. 2010. ApJ 723:1019-1037 3.2.1 Heithausen A, Bensch F, Stutzki J, Falgarone E, Panis JF. 1998. A&A 331:L65-L68 3.1.2 Heitsch F, Ballesteros-Paredes J, Hartmann L. 2009. ApJ 704:1735-1742 3.2.2 Helmi A. 2008. Asron. Astrop. Rev. 15:145-188 5.1, 5.2 Helmi A, Babusiaux C, Koppelman HH, Massari D, Veljanoski J, Brown AGA. 2018a. Nature 563:85-88 5.2 Helmi A, Navarro JF, Nordstrom B, Holmberg J, Abadi MG, Steinmetz M. 2006. MNRAS 365:1309-1323 5.1 Helmi A, van Leeuwen F, McMillan PJ, Massari D, Antoja T, et al. 2018b. A&A 616:A12 4.3.2 Helmi A, Veljanoski J, Breddels MA, Tian H, Sales LV. 2017. A&A 598:A58 5.1 H'enault-Brunet V, Gieles M, Evans CJ, Sana H, Bastian N, et al. 2012. A&A 545:L1 3.4.2 Hennebelle P, Chabrier G. 2008. ApJ 684:395-410 3.5.1.1 Hennebelle P, Chabrier G. 2011. ApJ 743:L29 3.2.2 Hennekemper E, Gouliermis DA, Henning T, Brandner W, Dolphin AE. 2008. ApJ 672:914-929 2.7 H'enon M. 1961. Annales d'Astrophysique 24:369. English translation by F. Renaud available as arXiv:1103.3499\n```", '4.2.2': "```\nKim JG, Kim WT, Ostriker EC. 2018. ApJ 859:68 3.3.2, 3.3.3, 3.4.2, 3.5.1.2 Kim Jh, Ma X, Grudi'c MY, Hopkins PF, Hayward CC, et al. 2018. MNRAS 474:4232-4244 3.5.3, 4.4.2 Kim WT, Ostriker EC. 2001. ApJ 559:70-95 3.1.2, 3.1.2 Kimm T, Cen R, Rosdahl J, Yi SK. 2016. ApJ 823:52 3.5.3, 4.4.2 King IR. 1962. AJ 67:471 4.1 King IR. 1966. AJ 71:64 1, 4.1 Klessen RS, Hennebelle P. 2010. A&A 520:A17+ 3.2.2, 3.4.2 Koo BC, McKee CF. 1992. ApJ 388:103 3.3.5 Koposov SE, Rix HW, Hogg DW. 2010. ApJ 712:260-273 5.1 Kordopatis G, Binney J, Gilmore G, Wyse RFG, Belokurov V, et al. 2015. MNRAS 447:3526-3535 4.3.2 Kortgen B, Seifried D, Banerjee R, V'azquez-Semadeni E, Zamora-Avil'es M. 2016. MNRAS 459:3460-3474 3.3.6 Kos J, Bland-Hawthorn J, Freeman K, Buder S, Traven G, et al. 2018a. MNRAS 473:4612-4633 5.3 Kos J, de Silva G, Buder S, Bland-Hawthorn J, Sharma S, et al. 2018b. MNRAS 480:5242-5259 1.1 Kounkel M, Covey K, Su'arez G, Rom'an-Z'u˜niga C, Hernandez J, et al. 2018. AJ 156:84 2, 13, 3.4.2 Kraus AL, Hillenbrand LA. 2008. ApJ 686:L111-L114 2.7 Krause MGH, Charbonnel C, Bastian N, Diehl R. 2016. A&A 587:A53 3.3.5 Kreckel K, Faesi C, Kruijssen JMD, Schruba A, Groves B, et al. 2018. ApJ 863:L21 3.2.1 Kroupa P, Aarseth S, Hurley J. 2001. MNRAS 321:699-712 3.4.2, 3.5.2 Kruijssen JMD. 2009. A&A 507:1409-1423 4.2.1 Kruijssen JMD. 2012. MNRAS 426:3008-3040 2.4, 3.3.6, 3.4.2, 3.5.2, 3.5.2, 4.3.1, 4.4.2 Kruijssen JMD. 2014. Classical and Quantum Gravity 31:244006 1.1, 3.1.2 Kruijssen JMD. 2015. MNRAS 454:1658-1686 4.3.1 Kruijssen JMD, Bastian N. 2016. MNRAS 457:L24-L28 2.4 Kruijssen JMD, Maschberger T, Moeckel N, Clarke CJ, Bastian N, Bonnell IA. 2012a. MNRAS 419:841-853 3.5.2 Kruijssen JMD, Pelupessy FI, Lamers HJGLM, Portegies Zwart SF, Bastian N, Icke V. 2012b. MNRAS 421:19271941 4.4.1, 4.4.2 Kruijssen JMD, Pelupessy FI, Lamers HJGLM, Portegies Zwart SF, Icke V. 2011. MNRAS 414:1339-1364 2.3, 4.3.1, 4.3.1, 4.4.1, 4.4.2 Kruijssen JMD, Schruba A, Hygate APS, Hu CY, Haydon DT, Longmore SN. 2018. MNRAS 479:1866-1952 3.2.1, 3.3 Krumholz MR. 2014. Phys. Rep. 539:49-134 3.2, 10 Krumholz MR. 2018. MNRAS 480:3468-3482 3.3.3 Krumholz MR, Adamo A, Fumagalli M, Calzetti D. 2018. MNRAS In press, arXiv:1810.10173 2.1.2, 4.4.1 Krumholz MR, Adamo A, Fumagalli M, Wofford A, Calzetti D, et al. 2015a. ApJ 812:147 2.1.2, 2.1.2, 2.4 Krumholz MR, Bate MR, Arce HG, Dale JE, Gutermuth R, et al. 2014. Protostars and Planets VI :243-266 1.1 Krumholz MR, Dekel A, McKee CF. 2012. ApJ 745:69 3.2.1 Krumholz MR, Fumagalli M, da Silva RL, Rendahl T, Parra J. 2015b. MNRAS 452:1447-1467 1, 3, 3.3.1 Krumholz MR, Klein RI, McKee CF. 2012. ApJ 754:71 3.3.1, 12, 3.4.2 Krumholz MR, Matzner CD. 2009. ApJ 703:1352-1362 3.3.3 Krumholz MR, McKee CF. 2005. ApJ 630:250-268 3.1.2, 3.2.2 Krumholz MR, Tan JC. 2007. ApJ 654:304-315 3.2.1 Krumholz MR, Thompson TA. 2012. ApJ 760:155 3.3.4 Krumholz MR, Ting YS. 2018. MNRAS 475:2236-2252 3.6 Kuhn MA, Feigelson ED, Getman KV, Baddeley AJ, Broos PS, et al. 2014. ApJ 787:107 2, 1.2.1 Kuhn MA, Hillenbrand LA, Sills A, Feigelson ED, Getman KV. 2018. ApJ In press, arXiv:1807.02115 3.4.2, 3.4.2 Kundic T, Ostriker JP. 1995. ApJ 438:702-707 4.3.1 Kuznetsova A, Hartmann L, Ballesteros-Paredes J. 2015. ApJ 815:27 3.4.2 Kuznetsova A, Hartmann L, Ballesteros-Paredes J. 2018. MNRAS 473:2372-2377 3.4.2 Lada CJ, Lada EA. 2003. ARA&A 41:57-115 1.2.1, 1.2.2, 2.1.1, 3.3.7, 3.5.1, 3.5.2 Lada CJ, Lombardi M, Roman-Zuniga C, Forbrich J, Alves JF. 2013. ApJ 778:133 3.1.1, 3.2.1, 10\n```  \n```\nLada CJ, Margulis M, Dearborn D. 1984. ApJ 285:141-152 3.5.2, 15 Lamb JB, Oey MS, Werk JK, Ingleby LD. 2010. ApJ 725:1886-1902 2.2 Lamers HJGLM, Baumgardt H, Gieles M. 2010. MNRAS 409:305-328 4.2.1, 4.2.2, 4.4.1, 4.4.1 Lamers HJGLM, Gieles M, Bastian N, Baumgardt H, Kharchenko NV, Portegies Zwart S. 2005. A&A 441:117-129 4.2.2, 4.4.1, 4.4.1, 4.4.1, 4.4.1 Larsen SS. 2002. AJ 124:1393-1409 2.2 Larsen SS. 2009. A&A 494:539-551 2.2 Larsen SS, de Mink SE, Eldridge JJ, Langer N, Bastian N, et al. 2011. A&A 532:A147 2.1.1 Larsen SS, Richtler T. 2000. A&A 354:836-846 2.4 Larson RB. 1981. MNRAS 194:809-826 3.1.1 Lee EJ, Miville-Deschˆenes MA, Murray NW. 2016. ApJ 833:229 3.2.1, 10, 3.2.1 Lee H, Hopkins PF, Squire J. 2017. MNRAS 469:3532-3540 3.6 Lee YN, Hennebelle P. 2016a. A&A 591:A30 3.2.2, 3.4.2 Lee YN, Hennebelle P. 2016b. A&A 591:A31 3.2.2, 3.4.2 Leitherer C, Schaerer D, Goldader JD, Delgado RMG, Robert C, et al. 1999. ApJS 123:3-40 3 Leroy AK, Bolatto AD, Ostriker EC, Walter F, Gorski M, et al. 2018. ApJ In press, arXiv:1804.02083 2, 9, 3.1.2, 3.1.2 Leroy AK, Schinnerer E, Hughes A, Kruijssen JMD, Meidt S, et al. 2017. ApJ 846:71 3.2.1, 10 Li H, Gnedin OY, Gnedin NY. 2018. ApJ 861:107 3.5.3, 4.4.2 Li H, Gnedin OY, Gnedin NY, Meng X, Semenov VA, Kravtsov AV. 2017. ApJ 834:69 3.5.3, 4.4.2 Li H, Li D, Qian L, Xu D, Goldsmith PF, et al. 2015. ApJS 219:20 3.3.1 Lim S, Lee MG. 2015. ApJ 804:123 8 Linden ST, Evans AS, Rich J, Larson KL, Armus L, et al. 2017. ApJ 843:91 5, 6 Liu F, Asplund M, Yong D, Mel'endez J, Ram'ırez I, et al. 2016a. MNRAS 463:696-704 2.6 Liu F, Yong D, Asplund M, Ram'ırez I, Mel'endez J. 2016b. MNRAS 457:3934-3948 2.6, 3.6 Lombardi M, Alves J, Lada CJ. 2015. A&A 576:L1 3.1.2 Longmore SN, Kruijssen JMD, Bastian N, Bally J, Rathborne J, et al. 2014. Protostars and Planets VI :291-314 3.3.6, 3.4.2, 3.4.2 Longmore SN, Rathborne J, Bastian N, Alves J, Ascenso J, et al. 2012. ApJ 746:117 14 Lopez LA, Krumholz MR, Bolatto AD, Prochaska JX, Ramirez-Ruiz E. 2011. ApJ 731:91 3.3.5 Lopez LA, Krumholz MR, Bolatto AD, Prochaska JX, Ramirez-Ruiz E, Castro D. 2014. ApJ 795:121 3.3.5 Lynden-Bell D, Kalnajs AJ. 1972. MNRAS 157:1 4.3.2 Mackey AD, Gilmore GF. 2003a. MNRAS 338:120-130 4.1 Mackey AD, Gilmore GF. 2003b. MNRAS 338:85-119 4.1 Mackey J, Gvaramadze VV, Mohamed S, Langer N. 2015. A&A 573:A10 3.3.5 Madrid JP, Leigh NWC, Hurley JR, Giersz M. 2017. MNRAS 470:1729-1737 4.2.2, 4.3.1 Ma'ız Apell'aniz J. 2009. ApJ 699:1938-1955 2.1.2 Maji M, Zhu Q, Li Y, Charlton J, Hernquist L, Knebe A. 2017. ApJ 844:108 3.5.3 Mamikonyan EN, McMillan SLW, Vesperini E, Mac Low MM. 2017. ApJ 837:70 4.4.2 Mandelker N, van Dokkum PG, Brodie JP, van den Bosch FC, Ceverino D. 2018. ApJ 861:148 4.4.2 Marino AF, Yong D, Milone AP, Piotto G, Lundquist M, et al. 2018. ApJ 859:81 2.6 Martell SL, Shetrone MD, Lucatello S, Schiavon RP, M'esz'aros S, et al. 2016. ApJ 825:146 5.3 Martinez-Medina LA, Pichardo B, Peimbert A, Moreno E. 2017. ApJ 834:58 4.3.1 Mathieu RD. 1983. ApJ 267:L97-L101 3.5.2 Matzner CD. 2002. ApJ 566:302-314 3.3, 3.3.2, 3.3.6 Matzner CD, Jumper PH. 2015. ApJ 815:68 3.3, 3.3.1, 3.4.2 Matzner CD, McKee CF. 2000. ApJ 545:364-378 3.3.1 McCrady N, Graham JR. 2007. ApJ 663:844-856 1.1, 5, 2, 9 McKee CF, Ostriker EC. 2007. ARA&A 45:565-687 3.1.2, 3.2 McKee CF, Parravano A, Hollenbach DJ. 2015. ApJ 814:13 1.2.1\n```"}