Star Clusters Near and Far

Star clusters are fundamental units of stellar feedback and unique tracers of their host galactic properties. In this review, we will first focus on their constituents, i.e. detailed insight into their stellar populations and their surrounding ionised, warm, neutral, and molecular gas. We, then, move beyond the Local Group to review star cluster populations at various evolutionary stages, and in diverse galactic environmental conditions accessible in the local Universe. At high redshift, where conditions for cluster formation and evolution are more extreme, we are only able to observe the integrated light of a handful of objects that we believe will become globular clusters. We therefore discuss how numerical and analytical methods, informed by the observed properties of cluster populations in the local Universe, are used to develop sophisticated simulations potentially capable of disentangling the genetic map of galaxy formation and assembly that is carried by globular cluster populations.


Introduction
Deciding if a grouping of stars is a gravitationally bound star cluster is very challenging when the kinematics of the single stars cannot be traced. To use star clusters as tracers of galaxy evolution, it requires the understanding of the formation and evolution of gravitationally bound stellar systems; i.e. star clusters that are likely to survive for a certain time span within their host systems. A definition of "star cluster" valid at any position in space and time in our Universe is very much a challenge. We will discuss the dynamical conditions of clusters when we can access their stellar kinematic information, i.e. properties of star cluster populations resolved in their stellar components (young star clusters, YSCs, and globular clusters, GCs) in our own Galaxy and local neighbours (Sect. 2). We point out the interested reader that a detailed description of the physics of cluster formation and evolution as single entity can be found in Krause et al. (2020), another review in this series.
The focus of this work is to review cluster formation and evolution in the framework of galaxy assembly. To achieve this goal it requires to move beyond the Local Group. At increasing distances we lose information about the detailed components of star clusters (we will provide a detailed description of the assumptions made) but we acquire a statistical view of the entire cluster population and wider ranges of galactic environments and diverse physical conditions for star and cluster formation. Throughout Sect. 3 to Sect. 9, we will describe the most recent results we have acquired on cluster populations from studies of local galaxies and the lesson we have learned so far about cluster formation and evolution. However we cannot restrict a review on cluster populations to our local Universe. In recent years, we have made huge progresses in probing under which conditions star formation operates in increasingly younger galaxies where the bulk of their stellar mass is assembled. GC populations are witnesses of these key phases of galaxy evolution. We need a way to read the information encoded in their properties as single entities and overall populations. In Sect. 10 and Sect. 11 we will discuss how we can link GC formation and evolution to the galaxy assembly history of our Universe. Analytical implementations of cluster populations (based on the lesson we learn from studying YSC populations in the local Universe) or self-consistent formation of clusters into cosmological simulations of our Universe will be summarised in light of their capabilities to successfully reproduce well established GC properties (mass function, blue tilt, age-metallicity relation, multi-peaked metallicity distributions, etc.). Simulations are not the only accessible way to probe cluster formation at redshift beyond zero. In recent years, thanks to the help of magnification lenses, we have been able to directly detect the formation sites of GCs. These detections are paving the way to new fundamental questions such as what conditions are required for GCs to form, at which redshift GCs start to form, did GCs contribute to reionise our Universe at redshift beyond 6. In the last Section (Sect. 12) we summarise the current observational results, predictions for the redshift formation of GCs and how next generation facilities can open our horizon into the formation of gravitationally bound stellar structures at any redshift.

Definition of a Bound Star Cluster
We all agree what a star cluster is when we see one, but coming up with a quantitative definition is hard. Most people agree that a star cluster is a group of stars, but so are galaxies, so not all stellar groupings are star clusters. Any additional specification brings a risk of over-simplifying things and chopping up continuous distributions in discrete samples. For example, we could try and exclude galaxies from the definition and define a star cluster as a stellar grouping that formed in a single burst, from a well-mixed molecular cloud, implying that all stars have the same age and chemical composition. However, all stellar systems will have an age spread to some degree (∼ 10 5-6 yr), and light-element abundance spreads (HeCNONaMgAl) are the norm in old, massive GCs, leading Carretta et al. (2010) to define a bona fide GC by the presence of a Na-O anti-correlation (more on this in Sect. 1.2). Also, nuclear star clusters display clear spreads in age (Gyrs) and [Fe/H] (∼ 1 dex, see Seth et al. 2019, for a recent review). Alternatively, we could argue that star clusters are dark matter (DM) free, which removes galaxies. However, the first star clusters to form in the Universe may have formed in small DM halos (∼ 10 9 M ) which have since been stripped (e.g., Trenti et al. 2015).
Despite the risks that are associated with defining what a star cluster is, there is a need for a working definition that separates star clusters from associations, where the former are assumed to be gravitationally bound and the latter are usually assumed to be unbound. This is particularly important in extra-galactic samples, where large stellar groupings ( 10 pc) are only marginally spatially resolved and therefore resemble compact star clusters. To test whether a system is bound, one needs kinematics to compare to an independent photometric mass estimate. This is not available for most samples. Gieles and Portegies Zwart (2011) propose a simple proxy, namely an estimate of the dynamical age of the stellar group, defined as the ratio of the age over the crossing timescale: Π ≡ age/τ cr , where τ cr ∝ R 3 eff /(GM), with R eff the half-light radius, G 0.0045 pc 3 M −1 Myr −2 the gravitational constant and M the cluster mass. Stellar groups with Π 1 have evolved for more than a crossing time, suggestive of them being bound. For groups with Π 1 it is not possible to say whether the group is bound, hence by applying this cut at young ages, some bound systems might be excluded. However, because unbound groups expand, they quickly evolve to Π 1, hence at older ages a separation at Π = 1 is safe. We note that a more accurate definition of the crossing time involves the velocity dispersion: τ cr ∝ R eff /σ . For a system in virial equilibrium the two definitions can be used interchangeably, but for unbound systems that are flying apart with high velocities R eff /σ R 3 eff /(GM). Hence using the crossing time based on kinematics -albeit more accurate -could lead to an erroneous classification as bound (Π 1) if virial equilibrium is assumed. Defining τ cr in terms of M and R eff avoids this, and has as the additional advantage that it is straightforward to determine this ratio for large numbers of stellar groups in extra-galactic samples using photometry only. A final word of caution is in place: the distribution of Π for young stellar groupings in the Milky Way and nearby galaxies is continuous (Gieles and Portegies Zwart 2011), making a classification of objects with Π 1 arbitrary and care should be exercised when discussing individual objects near the boundary. The simplicity of the definition and the ease of determining it for a large (extragalactic) sample of stellar groups has made it a popular tool to attempt to separate bound star clusters from unbound associations (e.g., Bastian et al. 2012;Johnson et al. 2016;Ryon et al. 2017a).

The Rich Chemical Footprints of GCs and Viable Channels to Their Origin
As discussed in the previous section, the dynamical state of a stellar grouping can be used as a classification of star clusters and associations. Apart from the clarity it provides when discussing observational samples, this classification might also provide insight into open questions regarding their formation mechanism. Do bound star clusters form in a fundamentally different way than unbound associations and fields stars? Or are star clusters simply the high density tail of the density distribution of the interstellar gas from which they form?
In this section we discuss the additional clues that come from the stellar populations. The (very small, or even absent) star-to-star spread in [Fe/H] and α-and s-process elements in open clusters in the Milky Way is consistent with the measurement uncertainty of highresolution spectra ( 0.05 dex, e.g., De Silva et al. 2007a), suggesting that all stars formed from a chemically homogeneous (i.e. well-mixed) molecular cloud. Somewhat surprisingly, the same level of homogeneity was found for unbound moving groups (De Silva et al. 2007b). These moving groups tends to be older (∼ 50 Myr) than associations (∼ 10 Myr). If they are the same objects in different evolutionary phases, it suggests that both bound clusters and unbound associations form in a similar way.
The picture is completely different for old Galactic clusters: nearly all GCs in the Milky Way ( 10 Gyr) display star-to-star light-element abundance variations (He, C, N, O, Na, Mg, Al). In fact, GC stars that have field-like chemical composition (usually called first population, 1P) are only a minority (10 to 30%, depending on the GC), while the chemically anomalous stars (or second population, 2P) dominate (e.g., Bastian and Lardo 2018). The abundance variations display particular anti-correlations, that are the tell-tale of hot-hydrogen burning (Denisenkov and Denisenkova 1990;Prantzos et al. 2007Prantzos et al. , 2017 20 MK), 45 MK) and 70 MK) anti-correlations. Additional support for the hot-hydrogen burning hypothesis comes from the fact that some GCs display broadened or split main sequences in the colour magnitude diagram (Anderson 2002;Bedin et al. 2004;Piotto et al. 2007Piotto et al. , 2012Han et al. 2009), which has been attributed to He spreads (Norris 2004;D'Antona et al. 2005;Charbonnel 2016;Lagioia et al. 2019). More recently, the so-called chromosome map was introduced, which is based on specific combinations of HST filters that are sensitive to He and N (Milone et al. 2015); it turns out to be a very powerful photometric tool to detect the presence of multiple stellar populations in Galactic and extra-galactic GCs (Milone et al. 2017;Zennaro et al. 2019;Saracino et al. 2019). In the last few years, evidence for N enhancement has been found in younger clusters ( 2 Gyr) in the Magellanic Clouds (e.g., Hollyhead et al. 2017;Martocchia et al. 2018a,b). Although a N spread is not necessarily pointing at the same origin as the variations in ONaMgAl seen in the older halo GCs, it is at the moment something that can not be explained by (single) stellar evolution models. As mentioned earlier, the O-Na anti-correlation is so ubiquitous that it has been suggested to be the unique identifying property of a genuine GC (Carretta et al. 2010). The ubiquity does not mean that all GCs are similar; quite the contrary, the details of the multiple populations are different in every GC (Milone et al. 2017). Some note-worthy trends with GC properties have been identified, that might provide clues to the origin of these multiple populations. For example, the Mg-Al anti-correlation is only found in the most massive and metal-poor GCs (Carretta et al. 2009a;Pancino et al. 2017;. In these clusters, also the minimum(maximum) O(Na) abundance is lower(higher) (Carretta et al. 2009b). Both findings are expected if these anti-correlations are the result of hot-H burning and the temperature of the polluter was higher in more massive and metal-poor GCs. In addition, both the He spread inferred from the main sequence broadening and the fraction of stars with anomalous abundances correlate with GC mass (Milone et al. 2014(Milone et al. , 2017, implying that more polluted material is produced per unit of cluster mass in massive GCs. More details on the observational signatures can be found in Bastian and Lardo (2018) and Gratton et al. (2019).
It is broadly accepted that anomalous (2P) GC stars formed out of original proto-cluster gas mixed with the H-burning yields ejected by massive (M > 5 M ) and short-lived 1P stellar polluters (e.g., Prantzos et al. 2017, and references therein). However, there is no consensus on the nature of the polluter(s) and the pollution mechanism (e.g., Renzini et al. 2015;Charbonnel 2016;Bastian and Lardo 2018). The vast majority of GCs show no spread in iron abundance, which suggests that (self-)enrichment by supernovae plays no role (e.g., Simmerer et al. 2013;Marino et al. 2015, 2018, andreferences therein). GCs with clear [Fe/H] spreads, such as Omega Centauri (ω Cen) and M54, are among the most massive clusters and are (former?) nuclear clusters. Several possible polluters have been proposed, which all reach the required high temperatures at some stage of their evolution: Asymptotic Giant Branch (AGB) stars (∼ 5-6.5 M , Ventura et al. 2001), massive stars ( 20 M , Maeder andMeynet 2006;Prantzos and Charbonnel 2006;Decressin et al. 2007a,b), massive binaries ( 10 M , de Mink et al. 2009;, and supermassive stars (SMSs, 10 3 M , Denissenkov and Hartwick 2014; . The different aspects and the pros and cons of the corresponding scenarios are described in Krause et al. (2020), a review of this series. Here we would like to conclude on the fact that observations of proto-GC at high redshift will help to discriminate between the different options in the near future. As we present in Sect. 10, proto-GC candidates have already been detected with the aid of gravitational lensing. Their number detections will significantly increase in the James Webb Space Telescope (JWST) and the European Extremely Large Telescope (E-ELT) era and spectroscopic sampling of their light will certainly help to confirm the nature of the stellar populations hosted in these systems. On the theoretical side, evolutionary synthesis models are developed to predict the characteristics of proto-GCs in the early Universe (Renzini 2017;Pozzetti et al. 2019). In this context, Martins et al. (2020) developed the first synthetic models of proto-GCs hosting multiple stellar populations and a SMS. They compute theoretical combined spectra and synthetic photometry in UV, optical, and near-IR bands over a wide range of redshift (1 to 10), and predict the detectability of cool SMS in proto-GCs through deep imaging with JWST NIR-CAM camera.

Young Star Cluster Populations Within the Local Group
The properties of YSCs in the local, resolvable Universe span a wide parameter space. Their masses range from low-mass associations, such as the Orion Nebular Cloud (ONC, e.g., Robberto et al. 2013), the Upper Scorpius association (e.g., Preibisch and Mamajek 2008), or the Pleiades, to the massive super star clusters (SSCs, defined in this review as clusters with stellar masses above 10 5 M ) like R136 (Hunter et al. 1995; in the Large Magellanic Cloud (LMC) or in other galaxies of the Local group (e.g., Hunter 2001;Sabbi et al. 2008). While it is still possible to resolve the most luminous stars in nearby galaxies Sacchi et al. 2018), only star clusters located in the Milky Way or the Magellanic Clouds are close enough that a major fraction of their stellar population can be resolved with existing telescopes, such the Hubble Space Telescope (HST) or from the ground using adoptive optic (AO) systems to correct for atmospheric turbulence (e.g., the Very Large Telescope (VLT), Keck, or Gemini).
Observing the low-mass stellar populations is crucial to understand cluster formation and evolution. While most of the energy is emitted by the most massive stars, the majority of the mass budget is bound in the low-mass stars, which influence the gravitational potential of their host systems' long-term evolution.
New instruments like large integral field units (IFUs), such as the Multi Object Spectrographic Explorer (MUSE, Bacon et al. 2010) mounted at the VLT, the Gaia satellite (Prusti et al. 2016;Brown et al. 2018), and long baseline photometric observations allow us, for the first time, to study the detailed 3D dynamics of the majority of stars in these resolved star clusters, including the dynamics of the gas (e.g., Kamann et al. 2013;McLeod et al. 2015;Zeidler et al. 2018;Lennon et al. 2018;Wright and Mamajek 2018;Ward and Kruijssen 2018;Ward et al. 2019;Getman et al. 2019;Zari et al. 2019). This provides insights into the star cluster formation modes: Do star clusters form hierarchically, following the structure of the giant molecular cloud (GMC) (e.g., Kruijssen et al. 2012b;Parker et al. 2014;Longmore et al. 2014;Walker et al. 2015Walker et al. , 2016Barnes et al. 2019;Ward et al. 2019), or do they form in monolithic, central starburst-like events (e.g., Lada et al. 1984;Bastian and Goodwin 2006;Banerjee and Kroupa 2015)? Future missions and telescopes, such as JWST, the E-ELT, or the Thirty Meter Telescope (TMT), will provide the necessary angular resolution and wavelength ranges to further investigate the low-mass end of the initial mass function (IMF) 1 and the embedded objects in the surrounding HII regions. This will lead to a better understanding of the star formation and feedback processes in these HII regions under the influence of a large central population of massive, luminous OB stars, eventually shedding light on the formation process of the stellar populations within GCs.
The scope of this section is not to describe the detailed, individual parameters of each star-forming region but rather give a more general overview over the observed parameter space provided by star clusters that are close enough to be resolved into their components. When looking at YSCs in local galaxies outside the Local group, the rich information contained in each of these star-forming regions and very young clusters will be collapsed into a handful amount of pixels. We give up on their single physical components and look at them in a statistical meaningful way.

Young Massive Star Clusters in the Milky Way and Magellanic Clouds
Compared to other, more distant galaxies (e.g., Gascoigne and Kron 1952;Hodge 1961a,b), the Milky Way host relatively few massive YSCs, none of which are expected to survive a significant fraction of a Hubble time (Krumholz and McKee 2019, and references therein). Yet, together with the Magellanic Clouds, these are the only places where individual stars can be resolved down to the hydrogen burning limit or below, even in the dense star clusters (e.g., Stolte et al. 2006;Zeidler et al. 2015). Milky Way and Magellanic Cloud star clusters are crucial to understand the first few million years, during which star formation is on-going and a significant amount of gas is still present. Stellar and gas dynamics and interactions, feedback processes, and possible secondary triggered star formation in the surrounding HII regions is still poorly understood due to the lack of sufficient, large scale observations. To trace star cluster evolution over a longer time scale sophisticated simulations are necessary, yet these are only as good as their initial conditions. To understand cluster formation in a more distant Universe, unresolved with current telescopes, the local star cluster observations must suffice to deepen the knowledge about these cluster initial conditions.
With ∼ 5 × 10 4 M (Gennaro et al. 2017, and references therein), Westerlund 1 (Wd1) is the most massive YSC in the Milky Way. With an age of ∼ 5 Myr it has already undergone several supernova explosions and is dynamically more evolved than other, younger massive Milky Way star clusters (m > 10 4 M ), such as Westerlund 2 (Wd2, Westerlund 1961) or NGC 3603. The Milky Way also hosts two YSCs in extreme environments, the Arches and Quintuplet Cluster (e.g., Figer et al. 1999;Stolte et al. 2010Stolte et al. , 2015. Being only ∼ 60 pc away from the Galactic center , the star-forming regions in this environment are characterised by high stellar and gas densities (Walker et al. 2015), highly compressive tidal fields (Kruijssen et al. 2019c), turbulent motion (Oka et al. 2001;Henshaw et al. 2016), and they are located in a very steep gravitational potential. Their observations however are challenging because they are not only one of the densest and most efficient star-forming regions but are affected by a visual extinction exceeding 20 mag. The LMC hosts the only nearby SSC. With an estimated age of 1-2 Myr, R136, located in 30 Doradus (30 Dor) or the Tarantula Nebula, hosts the most massive and luminous stars known with masses up to ∼ 300 M and a spectral types of O2-3V .
While the young Milky Way star clusters mainly have Solar metallicity, star clusters in the LMC (distance: 50 kpc, A V ≈ 0.3, Schaefer 2008;Imara and Blitz 2007) and the Small Magellanic Cloud (SMC; distance: 62 kpc Hilditch et al. 2005) are located in a more metal poor environment, corresponding to the properties at higher redshifts (z). The typical metallicity in the LMC is 0.5 Z with a dust-to-gas ratio of ∼ 1/3 of the Milky Way and the SMC has even lower values (0.25 Z , dust-to-gas ratio: ∼ 1/6 of the Milky Way, e.g., Russell and Dopita 1992;Rolleston et al. 1999;Lee et al. 2005;Roman-Duval et al. 2014). This allows star and star cluster formation and evolution studies in lower metallicity environments, where effective stellar temperatures and luminosities are higher, resulting in faster stellar evolution and lower mass-loss rates (Kudritzki and Puls 2000). An increased stellar temperature leads to higher far-ultra-violet (FUV) fluxes emitted by the most massive Oand B-stars. Additionally, the low extinction toward the Magellanic Clouds allows detailed observations in the ultra-violet (UV) and FUV.

Observing Young Star-Forming Regions
Because of the large dynamical range of the physical processes in star clusters, observations across the full electromagnetic spectrum are necessary to fully understand these systems.
Stars and star clusters form due to the gravitational collapse of (parts of) GMCs (Kennicutt and Evans 2012; . These GMCs mostly contain neutral (HI) and molecular (H 2 ) hydrogen at temperatures of 10-50 K. These temperatures make it necessary to observe star formation sites with radio telescopes, such as the Atacama Large Millimeter Array (ALMA) to understand the GMC dynamics by observing various CO transition lines of the cold interstellar medium (ISM, e.g., Yonekura et al. 2005;Furukawa et al. 2009;Heyer et al. 2009a;Sun et al. 2018;Tsuge et al. 2019). The unprecedented high spatial resolution of radio telescopes also allows for the direct observation of protoplanetary and debris disks around young stars in the later stages of the star formation process (e.g., Brogan et al. 2015;. Observations from the infrared (far to near) to the FUV are necessary to observe the stellar IMF and the resulting wide range of stellar masses (from 80-100 M OB stars to the faint 0.1 M dwarf stars, see Fig. 1), high (differential) extinction, still deeply embedded objects, YSOs, and disk-bearing objects. Class 0-III YSOs and active star formation may still be present in remaining gas and dust clouds (Carlson et al. 2007). Mid and far-infrared space telescopes, such as Spitzer, the Stratospheric Observatory for Infrared Astronomy (SOFIA), and Herschel, are able to look more deeply into the dense gas clouds and observe and classify these objects (e.g., Gaczkowski et al. 2013).
Combining various optical and NIR broad-band filters, such as UBV I J H K-filters allow the construction of various two-color diagrams and color-magnitude diagrams (CMDs), which can be used together with model isochrone fitting to determine age, distance, and extinction of the stellar population as well as the individual stellar masses. This method is widely applied to very young and open clusters in the Milky Way and in general to star clusters in both the Magellanic Clouds and M31 (e.g., Zinnecker and Yorke 2007;Glatt et al. 2010;Johnson et al. 2015).
Stars that are still in their pre-main sequence phase show strong excess emission in the NIR due to their circumstellar disks, Balmer line emission in the optical, and X-Ray emission due to large hot stellar coronae, magnetic coupling of the disks to the stellar surface, and flaring of the central stars. While X-ray observations require space missions (e.g., Chandra or XMM-Newton), the optical and NIR observations can also be obtained from the ground, and with extreme AO even at similar spatial resolutions as HST provides. Combining these broad-band observations with narrow-band observations such as the Hα or Paβ filter allows to detect pre-main sequence stars with active mass accretion (e.g., De Marchi et al. 2011;Beccari et al. 2015;Zeidler et al. 2016a;Kalari 2019) revealing protoplanetary disks. NUV and FUV photometry and spectroscopy from space is necessary to study and classify the most massive OB-stars. Their spectral energy distribution (SED) peaks in the (F)UV and most of their parameters are degenerate in optical CMDs. With these data stellar winds and the FUV flux budget can be measured , which is responsible for accelerated disk dispersal of the lower mass stars in the close vicinity (Clarke 2007), creating photo-dissociated regions (PDRs) in the surrounding gas cloud, as well as triggering secondary generation star formation. NUV photometry also allows to directly probe the UV extinction curve via the stellar color excess.
Large star-forming regions such as 30 Dor in the LMC (Sabbi et al. 2012) or the Carina Nebula Complex (Smith and Brooks 2008;Zeidler et al. 2016b) are highly substructured and show a multitude of individual star clusters and associations of various masses and ages (e.g., Trumpler 14, 16, NGC 3324, and the Treasure Chest in the Carina Nebula Complex or R136, NGC 2070, and Hodge 301 in 30 Dor). These regions are dominated by feedback processes from massive stars, supernova explosions, and the formation of new stars in the surrounding GMCs. The individual clusters within these star-forming regions span wide mass and age ranges, i.e., Grebel and Chu (2000) derived an age of 10-25 Myr for Hodge 301, while R136 contains stars as young as 0.5 Myr (Walborn and Blades 1997). Individual, isolated star clusters or star clusters within these larger star-forming regions may themselves show sub clustering and highly complex structures (see Fig. 2 and e.g., Kuhn et al. 2014).
The by far best studied young star-forming region is the Orion Nebula (Messier 42 and NGC 1976) and its associated star cluster, the ONC. Although not very massive (4.6 × 10 3 M , Hillenbrand and Hartmann 1998), its close proximity (∼ 440 pc, O'Dell and Henney 2008) makes it a perfect target to study the stellar and gas content. The ONC revealed the first direct imaging of protoplanetary disks (proplyds, McCaughrean and O'dell 1996;. A recent 3D kinematic study by Zari et al. (2019) using Gaia DR2 data (Prusti et al. 2016;Brown et al. 2018) showed that the ONC is highly sub-structured, with stellar ensembles of different ages in several kinematic groups, mixed in 3D space, which are overlapping in projection. Jerabkova et al. (2019) suggested that these YSCs may harbor multiple populations. Using OmegaCAM photometry they identified three populations with an age difference of 3 Myr between the oldest and the youngest sequence. These sequences cannot be described with binary or triplet systems alone leading to the conclusion that they are real, which is in agreement with the above findings by Zari et al. (2019) and suggesting star formation happens sequentially possibly triggered by the luminous OB stars.
Although multiple populations have not been detected in any other YSCs, mainly due to observational limitations, the majority of clusters and star-forming regions is still highly sub-structured showing multiple smaller clumps and are far from a spherical shape. In Wd2, Zeidler et al. (2018) recently discovered that the cluster stellar population shows multiple velocity components using MUSE observation to measure stellar radial velocities (RVs). These components appear to be spatially correlated with its two coeval subclumps (Hur et al. 2015;Zeidler et al. 2015), suggesting that they are, given the young age (∼ 1 Myr, Zeidler et al. 2015), an imprint of the formation history of the cluster. Other clusters, such as NGC 346 in the SMC are constituted even more complicated and show more than 16 individual sub clusters ). Wd1 is with ∼ 5 Myr older and dynamically more evolved, is elongated, which is probably a product of the past merging of former subclumps (Crowther et al. 2006;Gennaro et al. 2011). Other young massive star clusters in a similar mass and age range, such as NGC 3603 (Stolte et al. 2006;Pang et al. 2013) or R136  (Hunter et al. 1995;Sung and Bessell 2004;Sabbi et al. 2012; do not show present sub clustering. Their spherical shape may be explained through a spherical, burst-like single star formation event or through dynamical evolution suggesting that both hierarchical and in-situ cluster formation may be possible.

The Stellar Mass Function
The IMF is a key parameter that affects almost all the fields of astrophysics from stellar populations up to formation of first galaxies and galaxy evolution in general. Empirical studies in the Milky Way and the Magellanic Clouds have revealed a remarkably constant IMF, regardless of location, age, or metallicity (e.g., Chabrier 2003; Offner et al. 2014). This lead to the idea that the IMF is constant across the Universe, which implies constant and somehow regulated star formation processes. More recent observational studies in more extreme, extra galactic star-forming regions, low-metallicity star clusters, and likelihood studies have started to challenge this view (van Dokkum and Conroy 2010; Dib 2014; Kalari et al. 2018).
While the high-mass end of the IMF is relatively well-studied via simple star counts, this is more difficult for the low-mass end of the IMF due to observational limitations. Therefore, the shape of the IMF below a critical mass remains uncertain. Studying the IMF in the Milky Way and Magellanic Cloud clusters showed that massive stars appear to be over-abundant showing O-star candidates with a high probability to be ejected from the cluster core. Objects coloured in red are those with a high ejection probability while the two stars marked in blue have a trajectory consistent with a reduced probability to originate from the star cluster. The lines mark the direction of movement. This figure was published as Fig. 7 in Drew et al. (2019) compared to the expectations from a standard Salpeter (1955) slope resulting in a slightly top-heavy IMF. This holds for star-forming regions in varying environments (e.g., Zeidler et al. 2017;Schneider et al. 2018;Kalari et al. 2018;Hosek et al. 2019). Most of the studies on the cluster IMF assume that if the cluster is massive and young enough < 2-3 Myr the upper main sequence is fully populated. Yet, even if this assumption holds (even the most massive stars have lifetimes of a few million years) most observations have a limited survey area and, therefore, fast runaway stars may have left the immediate vicinity of the cluster and the survey area. Recent studies of several massive clusters showed that a significant fraction of O and B stars may have been ejected from the cluster center within the last million years (see e.g., Fig. 3 and Lennon et al. 2018;Drew et al. 2018Drew et al. , 2019. Although the reasons for the ejection are not clear these studies show that a significant number of stars can be missed using the traditional method of star counts in the closer cluster region. This argument though would lead to an even more top-heavy IMF. The majority of YSCs are highly mass segregated. Mass segregation describes the overabundance of high-mass stars relative to low-mass stars in the cluster center, compared to the outer regions of a star cluster. Mass segregation can have a significant influence on the cluster evolution and survival. The majority of massive stars will go supernova with the first ∼ 5 Myr. These stars are located deeper in the cluster's gravitational potential well and in the case of remaining gas within the cluster, these supernova explosions and the resulting abrupt mass ejection may disrupt the cluster faster. Additionally, the low mass stars that are on larger orbits around the center can be stripped away more easily while moving through the ISM. Both effects lead to an accelerated cluster dispersal. The origin of mass segregation is assumed to be two-fold: 1) primordial mass segregation, where more massive stars have formed in the cluster center (suggested for e.g., Wd2, Zeidler et al. 2017, or NGC 346, Sabbi et al. 2008 and 2) dynamical mass segregation, where the high-mass stars migrated inwards due to two-body relaxation driving the cluster towards (but never fully reaching, see e.g., Trenti and van der Marel 2013, Bianchini et al. 2016, Parker et al. 2016, Spera et al. 2016) energy equipartition (e.g., the Arches cluster, Habibi et al. 2013, or NGC 3603, Sung andBessell 2004). The origin of mass segregation for a specific cluster is often difficult to establish due to their dynamical evolution. Sub-clustering and the nonespherical distribution of stars highly influences the determination of the cluster mass, as well as the crossing time (e.g., Binney and Tremaine 1987), both on which the mass segregation time scale depends. McMillan et al. (2007) argue that merging sub-clusters keep their mass segregation imprint, the individual sub-clusters are less massive, and therefore, have shorter dynamical mass-segregation time scales. Conclusively, the origin of mass segregation of a star cluster that formed through hierarchical merging may be dynamical although the system as a whole suggests otherwise. Yet, the low-number high-mass stars in each individual subclump makes difficult to observe this effect.

Feedback Processes and Stellar and Gas Dynamics
Stellar feedback dramatically modifies the appearance of the region where star clusters form. During the first ∼ 10s of Myr, feedback from massive stars (mass > 8 M ), in the form of photoionisation and mechanical feedback (radiation pressure, stellar winds, SN explosions) ionises the left-over gas in the region and it imparts energy and momentum on the dust and gas out of which the star clusters form. We refer the interested reader to Dale (2015) and to Krause et al. (2020), another review in this series, for theoretical and numerical reviews of the stellar feedback from young star clusters. We summarise here some key observations of local massive star-forming regions in the Milky Way and Magellanic clouds carried out in recent years.
The combination of ground-based, AO supported telescopes using optical and NIR IFUs and space-based photometry and FUV spectroscopy provide astonishing insights into the feedback processes of these YSCs (e.g., see Fig. 4). FUV, NUV, and optical spectroscopy of massive OB stars allow accurate spectral typing, stellar wind strength measurements, and FUV flux determinations, providing information about the ionizing flux budget of the central star cluster (e.g. Smith et al. 2016). The early feedback originating within the star cluster could be the main driver for possible, triggered secondary star formation events (McLeod et al. 2018;Zeidler et al. 2018) in the shell of gas and dust still surrounding the cluster, which emphasized the necessity of a detailed analysis of the gas.
Mapping optical strong and forbidden line ratios yield insightful information on the properties of the ionised gas. Typical Hydrogen Balmer line decrements (e.g., Hα/Hβ) provide extinction information of the star-forming regions, leading to the reconstruction of the relative 3D location of the stars within the cluster, the individual PDRs, and gas and dust rims (McLeod et al. 2015). Studies show (e.g., McLeod et al. 2016a) that gas pillars need a minimum density to survive a given, local ionizing radiation level, confirming existing, theoretical mass-loss rate models. Combining multiple optical gas emission lines (i.e., the Balmer lines or [SII] 6717,6731) lead to the detection of embedded objects in the gas, like Herbig-Haro (HH) jets and bipolar outflows. Recently this has become feasible also outside the Milky Way, namely in the Magellanic Clouds (McLeod et al. 2016b(McLeod et al. , 2018 ratios), will correspond to an increase in the optical depth of the HII region and, therefore, smaller chances of for the ionising radiation to escape the region. In the LMC, Pellegrini et al. (2012) show the power of this technique, by tracing regions that are ionisation bound (optically thick to their ionising radiation) and/or density bound (optically thin, therefore leaking ionizing radiation). It is quite complex to relate the 2D picture of a star-forming region provided by IFU observations to the real status of the system. Numerical simulations show that high ionisation channels might form in low gas density regions created by turbulence in the gas . Recent, highresolution cosmological simulations have reported to find that most leaked ionizing photons are from star-forming regions that usually contain a feedback-driven kpc-scale superbubble surrounded by a dense shell. Young stars in the bubble and near the edge of the shell can fully ionize some low-column-density paths pre-cleared by feedback, allowing a large fraction of their ionizing photons to escape (Ma et al. , 2020.
Another leap forward will be made with JWST, both for detecting and probing accreating protostars and for investigating the effect of stellar feedback in the photo-dissociated (PDRs) and most dense gas regions that remains inaccessible at optical wavelengths. JWST will have the necessary spatial resolution in the NIR and MIR to detect YSOs in the gas rim, to study their properties in detail, and to determine to which extent the central ionizing cluster drives star formation into the surrounding gas cloud as seen e.g., for NGC 602 (Carlson et al. 2007). It will also give insight into the evolution and distribution of disk-bearing objects throughout the cluster. Observations (e.g., De Marchi et al. 2011;Zeidler et al. 2016a) hint that the close proximity to the OB-star population accelerates mass accretion processes in protostellar disks, leading to a faster disk dispersal, and eventually hinders planet formation, confirming various theoretical studies (e.g., Clarke 2007;Anderson et al. 2013;Winter et al. 2020). JWST will also allow us to map the gradual evolution of the gas and dust within the star-forming region as a function of its ionising stellar population. Polycyclic aromatic hydrocarbon (PAH) emissions, combined with photoionisation line emissions in the MIR, will provide an extinction-free view of the earliest phases of interaction between the source  Kharchenko et al. (2013a) of feedback and the surrounding ISM. These studies are currently limited to star-forming regions in the Milky Way and Magellanic Clouds (e.g., Chevance et al. 2016Chevance et al. , 2020c, but these will be extended beyond the Local Group, enabling a much more complete understanding of the early phases of star formation in a large variety of galactic environments and physical properties.

The Young Star Cluster Population of the Milky Way; Properties of Open Clusters
Although the number of massive YSCs in the Milky Way are limited, our Galaxy hosts numerous open clusters. The detection of those clusters can be challenging due to their potential lower stellar density (surface densities are not much higher than the those of the field stars) and the lack of gas. Extensive all-sky star catalogs, such as the ASCC-2.5 bright star catalog (Kharchenko 2001 Fig. 5). This catalog was further extended by another 202 clusters by Schmeja et al. (2014) and Scholz et al. (2015) leading to a total number of 3061 open and 147 GCs. The analyzed open clusters cover a wide range of ages (6.0 ≤ log(t) ≤ 9.8, Kharchenko et al. 2016) with absolute integrated K S -band magnitudes between −11.7 mag and 0.6 mag. Kharchenko et al. (2016) also analyzed the cluster luminosity function (CLF) with respect to cluster age and distance to the Galactic center using 2MASS photometry. The slope of the CLF appears to decrease with increasing age up to an age of log(t) ≈ 7.2 slightly increase for 8.3 < log(t) < 8.8 and decrease again for older ages. This behaviour may be explained by stellar evolution, changing the relative number of red giant stars in the individual clusters, which dominate the luminosity in the NIR. Additionally, Kharchenko et al. (2016) found that the CLF slope increases from the inner to the outer Galactic disk, which may indicate that massive clusters tend to be located more in the inner disk. One needs to caution a possible bias due to the limited depth of the survey toward the fainter end of the CLF, especially in the direction of the Galactic center, where extinction is higher.
The Gaia DR2 (Prusti et al. 2016;Brown et al. 2018), subsequent data releases, and machine learning techniques will increase the sample of open clusters by the analysis of the 5D phase space (l, b, , μ α , μ δ   Many of these newly detected clusters are located closer than 2 kpc from the Sun and, although, these studies probably consist of overlapping samples it shows that the statement of having compiled an almost complete sample (out to 1.8 kpc from the Sun) by Kharchenko et al. (2016) has to be handled with care. It also shows that with the further data releases of the Gaia mission and an introduction of new techniques in machine learning and neural networks using big data may reveal many more clusters in the near future.

Star Cluster Populations in the Local Universe; A Statistical View of Their Formation and Evolution
As we move away from the Local Group, we lose resolution on the single components of star clusters but gain a viewpoint into whole cluster populations forming into a much larger spectrum of galactic environments than offered by the Local Group. YSCs form in the densest regions of GMCs. Turbulent energy regulates the density structure and distribution of the cold gas. When gravitational fragmentation takes over, the densest regions in a cloud start to collapse. The interplay between turbulence and gravity results in a clustered and hierarchically distributed star formation. However, only gravitationally bound stellar systems, with stellar densities sufficient to overcome the tidal field of the galaxy and the destabilising gravitational pull of the remaining gas (Elmegreen and Hunter 2010;Kruijssen et al. 2011) will move away from their birth place and survive for a certain time span within their host galaxies.
We will introduce the main properties of GMC populations and, in general, the conditions of dense gas in local galaxies (Sect. 4). We will then focus on the statistical properties of cluster populations. We will show how YSCs trace the clustering properties of star formation and the largest coherent regions of star formation in different galaxies (Sect. 5). We will discuss the cluster size-mass relation as determined from measurement of cluster populations in local galaxies and its implications for cluster formation and evolution (Sect. 6).
To link GMC to YSC populations it is necessary to take into account that only a fraction of the dense gas forming stars will result in bound stellar systems. To date, contrasting evidence are reported in the literature both in favour and against a variation of the fraction of stars forming in bound clusters as a function of galactic environment. We will combine results available in the literature, with recent multi-band imaging survey of a large spectrum of local galaxies. Namely, we will refer to the analyses based on the HST treasury program Legacy ExtraGalactic UV Survey (LEGUS, Adamo et al. 2017;Cook et al. 2019). LEGUS sampled typical star-forming galaxies within 4 and 18 Mpc. The other project is Hubble imaging probe of extreme conditions and clusters (HiPEEC, Adamo et al. 2020a to be subm.). HiPEEC consists of 6 starburst/merger systems observed and analysed in a similar fashion as the LEGUS targets. The distances of these systems are between 30 and 80 Mpc and SFRs are above 10 M yr −1 , i.e. this program extend the cluster analysis from the LEGUS galaxy spectrum to highly efficiently star-forming systems. We will summarise the current status of the field from the observational, theoretical and numerical point of view (Sect. 7). After formation, the masses of the newborn YSCs follow mass distributions which have a power-law shape of index close to −2. However, the description of the YSC mass function requires the addition of an upper mass truncation (M c ) that we will discuss more in detail in Sect. 8. Finally, to fully describe YSC populations in local galaxies we also need to account for their dissolution rates. We will provide a description of the possible models put forward and how they are reflected in the literature in Sect. 9.

Properties and Conditions of the Molecular Gas in Local Galaxies
To understand the conditions under which YSCs form, it is important to understand how the properties and spatial distributions of GMCs depend on the environment (i.e. galaxy structure and dynamics, ISM pressure, etc.) and how these are linked to the properties of YSCs. From early Milky Way observations, GMCs seem to have relatively uniform properties and follow the relations described by Larson (1981), showing correlations between their sizes, line-widths and luminosities (e.g., Solomon et al. 1987;Heyer et al. 2009b). These relations describe clouds as having constant surface densities, being in virial equilibrium and following a size-line width relation. However, the universality of GMC properties and of Larson's relations has been questioned since. Early theoretical works already predicted an environmental dependence of cloud properties, such as their surface density, velocity dispersion, mass, and size distributions (e.g., Elmegreen et al. 1989). However, until recently, it has been challenging to extend GMC observations to other galaxies. To probe cloud properties and their dynamical state requires surveys of the molecular gas in nearby galaxies at high sensitivity and resolution, with a large coverage to provide a statistically significant sample, in a large variety of environments. Recent sub-millimetre facilities such as the Plateau de Bure Interferometer (PdBI), ALMA and the NOrthern Extended Millimeter Array (NOEMA) have now made this possible.
Several observational studies have since then revealed significant variations of the molecular cloud properties in nearby star-forming galaxies (e.g Koda et al. 2009;Hughes et al. 2013;Donovan Meyer et al. 2013;Leroy et al. 2013;Colombo et al. 2014;Leroy et al. 2016;Sun et al. 2018;Schruba et al. 2019), as well as in starburst galaxies and merging systems (Downes and Solomon 1998;Wilson et al. 2003;Wei et al. 2012) and it is now clear that there exists an environmental dependence of cloud properties. For instance, in the nearby spiral galaxy M51, Colombo et al. (2014) highlight the change of GMC properties in different regions of the galaxy (e.g. bar, bulge, disk, spiral arms, interarm regions). In particular, the GMC mass spectrum is found to vary (in terms of slope, normalisation and maximum mass; see Fig. 6) between arms and inter-arm regions: the population of clouds in the inter-arm regions is dominated by lower-mass objects (with a power law slope of the mass spectrum steeper than −2), than the population located in the arms Cumulative GMC mass spectra normalised by the observed area for different galaxies (left; from ) and for different regions of the spiral galaxy M51 (right; from Colombo et al. 2014). The differences in slopes and maximum masses in different galaxies and for different kinematic environments of a given galaxy suggest an environmental dependence of the GMC mass distribution (with a slope shallower than −2). This suggests that different mechanisms are regulating the growth and destruction of GMCs in different regions. Differences between the mass spectra of GMCs have also been observed between galaxies. In particular,  show that the mass distribution is flatter in denser and more massive galaxies (e.g. M51) compared to lower mass galaxies (such as M33; see Fig. 6).
The physical state of GMCs in extragalactic environments has also been compared to that seen in the Milky Way. Clouds in extragalactic environments seem in general to follow the Milky Way size-linewidth relation relatively well (Bolatto et al. 2008;Faesi et al. 2018). In addition, molecular clouds in other galaxies are typically observed to be bound or marginally bound structures, with a viral parameter close to unity (e.g., , and references therein), although there are exceptions, especially in low surface density galaxies. This agrees with observations of GMCs in our Galaxy. Sun et al. (2018) observed relatively universal virial parameters throughout a sample of nearby galaxies (α vir = 1-3, excluding M31 and M33; see Fig. 7), suggesting that molecular clouds are close to virial energy balance. However, they find a wide range of different turbulent pressures, with ranges of ∼ 1-2 dex within galaxies and a variation across the sample of four orders of magnitude. In particular, in the gas-rich, turbulent environment of the Antennae galaxies, which is the nearest major merger, the internal pressure of the gas is considerably elevated by the merging process compared to disc galaxies. However, this does not seem to significantly affect the dynamical state of the gas -the measured scaling relation between the CO line width σ , and the gas surface density Σ in the Antennae galaxies follows the average relation observed in the discs of star-forming galaxies. This example is particularly interesting, because it possibly forms a bridge to the extreme environmental conditions of high-redshift galaxies. As shown by Dessauges-Zavadsky et al. (2019), the GMC population detected in a typical starforming galaxy at z ∼ 1 has physical properties similar to those detected in local starbursts. Schruba et al. (2019) extend this result by showing that, statistically speaking, GMCs are in ambient pressure-balance virial equilibrium. Clouds are near energy equipartition in highpressure (molecular-dominated) environments (α vir ∼ 1-2, considering self-gravity only)  Sun et al. (2018). The top panel shows the relation between the CO line width σ and the gas surface density Σ at a common resolution of 120 pc for the discs of a sample of 15 nearby galaxies. The bottom left panel presents the mass-weighted distribution of the virial parameter α vir and the bottom right panel the distribution of turbulent pressure P turb for the disc (circles) and center (star symbols) of all galaxies. The spread in molecular gas dynamical state and internal turbulent pressure is clearly visible within and between galaxies, in particular when comparing normal star-forming disc galaxies with a merger system such as the Antennae, or more quiescent galaxies such as M31 and pressure confined by the diffuse ambient medium in low-pressure (atomic-dominated) environments, leading to higher viral parameters (α vir ∼ 3-20).
The environmental dependence of ISM structure and molecular cloud properties also affects the process of star formation and feedback. Recent work by Chevance et al. (2020a) analysing a sample of galaxies from the PHANGS survey (Physics at High Angular resolution in Nearby Galaxies; Sun et al. 2018, Schinnerer et al. 2019Leroy et al. in prep.) shows that the molecular cloud lifetime is not constant between and within galaxies, suggesting that the cloud lifecycle, star formation, and feedback are regulated by different physical processes in different galaxies. Specifically, the lifetimes of CO clouds sitting in environments of high global (kpc-scale) molecular gas surface density (≥ 8 M pc −2 ) are regulated by galactic processes (in particular the gravitational free-fall of the mid-plane ISM and shear, as predicted by Jeffreson and Kruijssen 2018 and in agreement with theoretical predictions by Pringle 2013 andRey-Raposo et al. 2017). By contrast, CO clouds in environ-ments of low global molecular gas surface density (≤ 8 M pc −2 ) decouple from the galactic environment and live for a free-fall time or a crossing time, i.e. their lifetime is regulated by internal dynamics. More details about the lifecycle of molecular clouds can be found in Chevance et al. (2020b), another review in this series. After the onset of massive star formation, the rate at which YSCs will destroy their parent molecular cloud through feedback is also likely to be environmentally dependent. The duration of this feedback phase has been shown to be relatively short (a few Myr after the formation of the first massive stars; e.g., Kawamura et al. 2009;Hollyhead et al. 2015;Grasha et al. 2018;Kruijssen et al. 2019d;Chevance et al. 2020a), suggesting a rapid cycling of the gas, with a low integrated efficiency of star formation per formation event (Kruijssen et al. 2019d;Chevance et al. 2020a). Future, multi-wavelength, high-resolution observations of the gas during the early phases of star formation and feedback with recent and coming facilities such as MUSE and the JWST in a large variety of environments will help understand how the properties of YSCs are affected by the properties of their natal molecular cloud and by the large-scale galactic environment.

The Clustering Properties of Young Star Clusters
Star formation is clustered, carrying the imprint of the gas from which stars form (Lada and Lada 2003;Ward et al. 2019). The gas in galaxies is hierarchically distributed, with power law mass distributions measured for molecular clouds (e.g., Elmegreen and Falgarone 1996;Roman-Duval et al. 2010), for the gas within both molecular (Lombardi et al. 2015) and pre-molecular clouds (Miville-Deschênes et al. 2010), and for dense cores (Stanke et al. 2006;Alves et al. 2007) and young stellar objects (e.g., Schmeja et al. 2008). In star-forming regions, the ISM fragments into smaller and smaller substructures, driven by supersonic turbulence aided by gravity (Elmegreen and Scalo 2004). At the smallest scales of the hierarchy are the stars, which also form fractal, scale-free structures of increasing density and decreasing scale from large star-forming complexes to bound star clusters (Elmegreen 2011). Observations reveal that young stellar populations, associations, and clusters are in fact clustered (e.g., Bastian et al. 2009; de la Fuente Marcos and de la Fuente Marcos 2009; Gouliermis et al. 2015;Grasha et al. 2015Grasha et al. , 2017a. YSCs trace the densest peaks of the hierarchy, and can be used to trace the clustering and its relation to the hierarchy of the gas. Pre-supernova feedback from massive stars, in the form of stellar winds and photoionisation, exposes stellar clusters (Hollyhead et al. 2015;Smith et al. 2016) and even disperses molecular clouds within the first 1-5 Myr (Kim et al. 2018;Rahner et al. 2019;Kruijssen et al. 2019d;Chevance et al. 2020a), i.e., well before secular and bulk motions act on star clusters to disperse them out of the parent environment. As a result, emerged star clusters with ages below a few Myr are closely associated with their parent cloud: the median age of the clusters whose location is projected within the area occupied by a molecular cloud is about 4 Myr and 2 Myr in the two galaxies M51 and NGC 7793, respectively (Grasha et al. 2018(Grasha et al. , 2019. If the clusters are bound and survive, they tend to migrate away from their parent cloud as they age; in the same two galaxies, clusters that are more than 4 times separated from the closest molecular cloud are about 12 and 3.5 times older, respectively, than those that are coincident with the cloud's footprint. Although the median ages of the star clusters not associated with the molecular clouds are drastically different for the two galaxies, ∼ 50 Myr in M51 and ∼ 7 Myr in NGC 7793, the differences disappear when the ages are normalized by the median age of the entire young (< 200 Myr) cluster population. The result is that the amount of time a cluster takes to migrate away from the parent molecular cloud is a fixed fraction, ∼ 1.1-1.3, of the median age of the cluster population. This result is likely related to the measured sizes of the molecular clouds that host YSCs, as well as to other effects (e.g., the dispersion velocity of the cluster population): the median radius of the clouds is ∼ 10 pc in NGC 7793 and ∼ 40 pc in M51, which is reflected in the footprints that are used for the association between clouds and clusters.
The two-point correlation function (TPCF) is a standard tool for measuring the clustering of a population, by quantifying how much the distribution of pairs deviates from a random distribution as a function of the pairs' separation. According to this metric, Zhang et al. (2001) report that younger star clusters in the Antennae are more clustered and more associated to longer wavelength tracers of star formation. Recently, Grasha et al. (2018Grasha et al. ( , 2019 find that molecular clouds are randomly or almost randomly distributed in M51 and NGC 7793, but massive clouds are clustered, in agreement with findings that massive clouds are preferentially located in the spiral arms and other galactic structures (Koda et al. 2009;Colombo et al. 2014). This is matched by the TPCF of the youngest, < 10 Myr, star clusters which are as strongly clustered as the massive clouds which they likely originated from.
As the star clusters age, they disperse or migrate within the galaxy, and this trend is also reflected in their TPCF (Grasha et al. 2015(Grasha et al. , 2017a. In general, clusters younger than about 40 Myr have a TPCF that is best described as a power law with exponent ∼ −0.6 to − 0.8, while star clusters older than ∼ 40-60 Myr are consistent with a random distribution. The age difference between any two pairs of clusters increases for increasing separation, according to a power law (Age) ∼ (Sep) α with α ∼ 0.3-0.6 (Efremov and Elmegreen 1998; de la Fuente Marcos and de la Fuente Marcos 2009; Grasha et al. 2017b). For reference, in turbulent-driven star formation α = 0.5 (Elmegreen and Efremov 1996). The power law is truncated at separations between 200 pc and 1 kpc, depending on the galaxy; this maximum separation marks the maximum 'cell of coherent star formation' present in galaxies (also see Kruijssen et al. 2019d). The size and age difference at the truncation point define a velocity, which is likely related to the average speed at which turbulence moves through the 'cell of coherence'. This speed is a constant multiple, about a factor 2-3, of the velocity difference imparted by shear in each of the galaxies (Grasha et al. 2017b). Thus, while turbulence is likely responsible for the age-separation relation, the maximum size of the cell of coherent star formation in a galaxy could be determined by its shear. A recent analysis of the nearby flocculent spiral galaxy NGC 300 shows that the cell size closely matches the gas disc scale height, suggesting that in this galaxy the cell size is set by stellar feedback breaking out of the host galaxy disc, rather than shear (Kruijssen et al. 2019d). It remains an important open question which physical mechanisms set the length scale for the independent building blocks of galaxies as a function of the galactic environment.

The Cluster Mass-Radius Relation; Insights Into the Dynamical State of Young Star Clusters
The radius of a star cluster is usually expressed in the effective radius (R eff ), defined as the radius containing either half the cluster light (for unresolved clusters) or half the number of observed stars (for resolved clusters). The mass-radius relation of cluster populations at various evolutionary stages provides insight in cluster formation and evolution. From early HST observations of young massive clusters in NGC 3256, Zepf et al. (1999) reported a surprisingly shallow size-luminosity relation: R eff ∝ L 0.07 , i.e. a nearly constant radius. Larsen (2004) found a similar shallow slope between R eff and cluster mass for young clusters ( 100 Myr) in several spiral galaxies, with a typical R eff ∼ 3 pc. A near constant radius was also found for clusters in several other galaxies (e.g., Scheepmaker et al. 2007;Ryon et al. 2015Ryon et al. , 2017b. The near constant radius implies that massive clusters are denser than low-mass clusters and it is not clear whether this relation is the result of nature or nurture. These findings are surprising, because molecular clouds -from which clusters form -have a constant surface density (i.e. a radius increasing as the square-root of the mass). But a word of caution is in place, because for these extra-galactic samples the resolution limit imposes a constant lower limit to the values of R eff that can be resolved, possibly biasing the mass-radius relation to a constant value. In addition, in most of these samples, clusters with a range of ages are included, making it difficult to separate formation from evolution effects.
Both the resolution and age effect can be addressed by looking at the youngest Galactic clusters. For Galactic embedded clusters with 10-100 stars, Adams et al. (2006) find a steep mass-radius relation of the form R eff ∝ N 1/2 * , where N * is the number of stars, i.e. a constant surface density. Because these clusters still have gas associated with them, this is likely as close we can get to observing the initial mass-radius relation of star clusters. The selection procedure of these clusters likely puts a lower limit on the observable surface density, possibly biasing the results to this steep relation. The slightly older Galactic open clusters in the catalogue of Kharchenko et al. (2013b), with masses derived from the tidal radii by Piskunov et al. (2007) show a near constant volume density (i.e. R eff ∝ M 1/3 ). It is important to realise that the masses and radii are simultaneously determined from King (1962)  However, because of the increasing spatial resolution limit with galaxy distance, it is difficult to make definite statements about this relation.
By splitting in age, Portegies Zwart et al. (2010) showed that a sample with clusters younger than 10 Myr contains clusters with M 10 5 M and R eff 1 pc, which are not found in the older sample 10-100 Myr. This may point at an expansion, something that was also noticed from the expansion with age of the core radii of extra-galactic clusters ) and the radii of Galactic clusters and OB associations (Pfalzner 2009). This expansion could be the result of residual gas expulsion (Banerjee and Kroupa 2017) or internal two-body relaxation ). In Sect. 6 of Krause et al. (2020), another review in this series, we discuss in details the dynamics of stars within a cluster (we refer the interested reader to that review for more information). It is here important to point out that two-body relaxation leads to a faster expansion of low-mass clusters, potentially erasing a mass-radius correlation, or even inverting it to an anti-correlation.

Observational Constraints
After more than a decade, it has yet not been possible to reach a final agreement on whether or not the fraction of stars that form in bound stellar clusters will depend on the intensity of the star formation event and on the general physical properties of the galactic environment where clusters form. Already from the very beginning, thanks to the HST high-spatial resolution optical/UV imaging, it was observed that SSCs preferentially formed in galaxies experiencing starburst events, like merger systems (e.g., Meurer et al. 1995;Whitmore and Schweizer 1995) or in dwarf galaxies (Billett et al. 2002). However, it was soon recognised that galaxies with higher SFR would likely host more massive (luminous) star clusters, simply because a larger number of clusters are formed and, therefore, the likelihood of sampling the cluster mass function at the high-mass end increases (Whitmore 2000;Larsen 2002). These relations simply describe a "size-of-sample effect". On the other hand, an increase in the fraction of stars forming in bound clusters implies a change in the clustering nature of star formation and in the efficiency at which bound stellar structures can be formed. We will quantify this process defining the cluster formation efficiency (hereafter CFE or Γ ) as the fraction of total stellar mass formed in clusters per unit time over a given age interval (cluster formation rate, CFR in units of M yr −1 ) divided by the SFR of the galaxy or region of the galaxy where clusters have been detected (e.g., Bastian 2008).
The pioneering work by Goddard et al. (2010) suggested that the CFE would steadily increase in galaxies with higher SFR per unit area. Since then, several observational works in the literature have extended this positive correlation both at high and low Σ SFR and galactic and sub-galactic scales (e.g., Adamo et al. 2011;Cook et al. 2012;Ryon et al. 2014;Johnson et al. 2016;Ginsburg and Kruijssen 2018, among many others, see references in Fig. 8). Kruijssen (2012a) derived an analytical model that reproduces the positive correlation between the two physical quantities (Γ and Σ SFR ). In this theoretical framework, bound star clusters form in the high-density peaks of the hierarchically organised ISM, where the free-fall time is shorter and the star formation efficiency higher. Additionally, it includes a prescription for how tidal perturbations caused by the encounters with dense GMCs set a minimum limit for the formation of bound star clusters. Overall, the model predicts the Γ vs. Σ SFR relation given three observable galactic properties, the gas surface density Σ gas , the Toomre parameter Q, and the angular velocity Ω, by converting Σ gas into Σ SFR with a star formation relation (e.g. the Schmidt-Kennicutt relation or the Bigiel et al. 2008 formulation, see e.g. Kennicutt and Evans 2012).
However, it is important to note that not all the data reported in the literature support the Γ vs. Σ SFR relation (e.g., Chandar et al. 2017;Fensch et al. 2019). Chandar et al. (2017) raised one important point regarding the observed Γ vs. Σ SFR relation. The data at high Σ SFR have Γ preferentially estimated using short time scales (e.g., 1-10 Myr), while data at low Σ SFR are estimated over a longer time range (e.g., up to 100 Myr). In their work they report to find a constant CFE at formation (over an age range of 1-10 Myr) close to 24%. The CFE constantly and rapidly declines to few percents in the age range 10-100 and 100-400 Myr, because of rapid cluster disruption, equally affecting the overall cluster populations of their sample. Therefore, they conclude that the observed Γ vs. Σ SFR relation is driven simply by mixing data in the literature that have CFR derived over different time ranges.
We take now this discussion a step further. These contrasting observational results may be understood in light of limitations and assumptions that go into the estimate of the CFE and Σ SFR . The estimates of Γ relies on: 1. a significant fraction of cluster candidates used to estimate the CFR being gravitationally bound; 2. reliable cluster age and mass determinations and detection limits; 3. a SFR tracer that is sensitive to the same age interval as the cluster population.  Adamo et al. (2020b, to be subm.). In both plots, empty symbols are used for plotting Γ values detected on sub-regions of galaxies. See main text for discussion Hence, the challenge to estimate Γ resides in the difficulty to create homogeneous reliable cluster catalogues combined with a star formation tracer that is sensitive to variations on time scales of tens of million years. In a recent review on star clusters,  carefully discuss the data available in the literature in light of the different assumptions made to estimate the CFE. As pointed out in their review, different works take different steps in constructing their cluster catalogues. For instance, the Chandar et al. (2017) results are based on catalogues considered "inclusive", following the terminology of . This means that catalogues are constructed automatically. Some steps are taken to remove potential stellar systems and interlopers, but no human visual inspection takes place. The latter task is undertaken by several studies in the literature as a necessary step to clean the cluster catalogues by fake cluster candidates, i.e. systems that may appear to have a light spread function larger than a stellar one, but simply because of interposition chances along the line of sight. This is a clear problem for detection of cluster candidates, as it is very likely to have these spurious detections in regions with higher and clustered stellar densities, i.e. typical star-forming regions in local galaxies that are also the place where bound clusters are formed (Adamo et al. 2017). Visual inspection has then become an important step in mitigating the contamination of unbound or spurious systems (see for example Bastian et al. 2012;Ryon et al. 2014;Johnson et al. 2015;Adamo et al. , 2017Cook et al. 2019). These cluster catalogues are referred to as "exclusive" in the terminology of . Analyses on the boundness criterion (Π , Gieles and Portegies Zwart 2011, and discussed in Sect. 1.1), confirm that the majority of candidate clusters contained in the exclusive catalogues satisfy the Π > 1 boundness condition (e.g., Bastian et al. 2012;Ryon et al. 2014Ryon et al. , 2017aJohnson et al. 2015). Independent approaches, like the one used by , confirm indeed that the fraction of massive stars forming in clustered star-forming regions are constant with SFR and that the time scales for dissolution of these regions is of ∼ 10 Myr. This implies that gravitationally unbound associations will heavily contaminate the fraction of bound star formation measured in 'inclusive' catalogues at early age ranges and more so in galaxies with lower Σ SFR , where the overall CFE is only a few percent and star formation is therefore dominated by unbound systems. At high Σ SFR , the bound clusters instead dominate the SFR, such that the inclusion of unbound systems only has a relatively small effect on Γ .
In recent years, to overcome the "human intervention and subjectivity" on the cluster catalogues and improve upon the reproducibility of these catalogues some attempts have been done to introduce supervised training of neural network algorithms capable of classify sources into cluster candidates or contaminants (e.g., Messa et al. 2018a;Grasha et al. 2019;Wei et al. 2019). These first attempts report ∼ 70% agreement between machine learning and human morphological classifications, similar to the agreement reported among several human classifiers according to Grasha et al. (2019). In the near future, improvements on the cluster catalogues by better training sets and recognition algorithms will certainly open the way to huge advancements in our understanding of cluster formation and evolution.
Finally, before presenting the observed Γ vs. Σ SFR relation, by using all data available in the literature to date, it is also important to note that at distances beyond 20 Mpc, clusters become point-like sources even at HST resolution. It requires different assumptions and approaches to produce cluster candidate catalogues (e.g., Adamo et al. 2010;Goddard et al. 2010;Linden et al. 2017). We assume that the light spread function of the compact object is dominated by the stellar cluster within the region. However, clusters are never formed in isolation but in a star-forming region with elevated stellar clustering. The reader can think of the 30 Doradus region (described in Sect. 2). 30 Doradus hosts a very massive cluster R136 of about 10 5 M only a few Myr old (Zinnecker and Yorke 2007) with its light dominated by very massive stars (Crowther 2019). Other, less massive clusters have been found within distances of tens of pc. For example, Hodge 301 is significantly older, containing red super giant stars. At a distance of 80 Mpc the entire region will fit within a few HST pixels. In UV and optical bands the light of the region will be dominated by the O stars within R136. At the NIR, the red super-giant stars in Hodge 301 will dominate the integrated flux of the region. With increasing galactic distance the approximation of having a single cluster within the compact source becomes weaker and weaker (Randriamanakoto et al. 2013) and eventually even the approximation of single stellar population fails. At distance beyond 80 Mpc starforming regions with the size of 30 Doradus become unresolved, we enter the domain of the so-called stellar clumps (e.g., Messa et al. 2019) studied up to redshift z ∼ 6 (e.g., Shibuya et al. 2016).
In Fig. 8, we collect data on the CFE available in the literature. We divide them into two groups accordingly to the age range used to estimate Γ . Galaxies for which Γ is estimated within 1-10 Myr range are showed in the left side plot. On the right side, we plot all the Γ estimates obtained using a longer age interval. The interval 10-100 Myr is used in values extracted in the LMC and SMC by Goddard et al. (2010) (Adamo et al. 2020b, to be subm.). In both plots we use filled symbols to indicate that values have been estimated over a large fraction (or the whole) galaxy, while empty symbols show Γ values estimated in sub-regions of a given galaxy. We do not separate the sample into inclusive vs. exclusive. However, we note that most of the data showed in the left plot for Γ (1-10 Myr) are obtained with inclusive cluster catalogues except for the M83 data points ) and the Sagittarius B2 complex in our own galaxy (Ginsburg and Kruijssen 2018). Except for Sagittarius B2, data with log(Σ SFR [M yr −1 kpc −2 ]) −1 have all been estimated in galaxies with distances 20 and lower than 80 Mpc. On the right side plot, most of datapoints have been derived with exclusive cluster catalogues except for the estimates by Overall, we observe that in spite of the large scatter (in part introduced by the different approaches to sample definition), both age intervals are statistically consistent with a positive correlation between Γ and Σ SFR . We do not see a drastic decline in Γ in the longer age interval (right plot) confirming the results reported in the literature that Γ values estimated at different age range are similar within uncertainties Johnson et al. 2016;Messa et al. 2018b). This latter result also suggests that in galaxies with log(Σ SFR [M yr −1 kpc −2 ]) −1 cluster disruption is not significantly affecting Γ estimates within the age range 1-100 Myr. We also notice that Γ estimated in sub-regions of galaxies (empty symbols) tend to occupy the region of the Kruijssen model (dashed and dotted line) obtained by using the Bigiel formulation of the Kennicutt-Schidmit relation, derived for kiloparsec-size regions within galaxies.

Comparison to Numerical Simulations
So far, we have discussed the evidence in support of and limitations in the data that can affect the observed Γ vs. Σ SFR relation. Another powerful tool to test whether such a relation arises from the physical properties of the star formation process is to use numerical approaches. In recent years, the increasing computational power combined with improved numerical recipes that account for sub-galactic scale physics such as stellar feedback and self-consistently evolving multi-phase ISM, has made possible to follow cluster formation and evolution in combination with galaxy evolution (e.g., Kruijssen et al. 2011Renaud et al. 2015;Li et al. 2017;Choksi et al. 2018;Lahén et al. 2019;Li and Gnedin 2019). However, computational resources are not unlimited and the initial generation of works modelling YSC populations focus on isolated or merging galaxies in idealized, non-cosmological simulations (e.g., Kruijssen et al. 2011;Renaud et al. 2015;Lahén et al. 2019). These setups were designed for developing methodologies and numerical methods, but they lack the cosmological context that determines and drives galaxy assembly history. The dependence of star cluster formation on galactic-scale properties means that modelling the formation of realistic star cluster populations also requires modelling the formation and evolution of galaxies and their environment. We will provide a more detailed description of different approaches to simulate cluster populations in a cosmological context in Sects. 10 and 11 and we refer the interested reader to Forbes et al. (2018) for a careful review of the field.
In this section, we compare observational data compiled from the literature to two different sets of cosmological simulations, that use radically different approaches and numerical recipes. In Fig. 9 we show the Γ vs. Σ SFR space. We combine all the observational data (black dots) included in the two plots of Fig. 8. Chandar et al. (2017) reports several Γ estimates using different age ranges for each galaxy in their sample. In this figure, we report Γ Fig. 9 Γ vs. Σ SFR plane. We combine the observed data available in the literature (filled grey symbols), i.e, the data are not separated by age ranges used to estimate Γ (see Sect. 7 and Fig. 8 for description of the observed data). We also include median and 25 and 75% quartiles of the observed data (filled black symbols) estimated in Σ SFR bins, of size indicated by the horizontal bars. Left: We include the resulting Γ vs. Σ SFR extracted from the simulation sets by . The curves show the median and quartiles of different sets of simulations using different feedback prescriptions. Right: We include the extracted Γ and Σ SFR for the E-MOSAICS simulations (filled square symbols) by . The values have been extracted from galaxies at redshift 0 and are plotted color-coded by the galaxy stellar mass in the same snapshot estimated at the youngest age interval for each galaxy to avoid contamination by disruption. We report in Table 1 the median and quartiles of the observed Γ values estimated in Σ SFR intervals. These values are also included in Fig. 9 to facilitate the comparison with simulations. It is important to notice that in order to estimate the median trends of the observational results, we have combined heterogeneous datasets that does not uniformly sample the Σ SFR space and, therefore, these trends should not be over-interpreted. Overplotted on the left side of Fig. 9 are the resulting Γ vs Σ SFR obtained by the galaxy simulations of . On the right side, we overplot the Γ and Σ SFR reported by  obtained in the E-MOSAICS simulations (MOdelling Star cluster population Assembly In Cosmological Simulations within EAGLE Pfeffer et al. 2018;. Despite the fact that both sets of simulations use vastly different numerical approaches, the resulting fraction of stars forming in clusters increases with Σ SFR in both works. If we look into the details of each simulation we can try to understand how they can inform us on the origin and physical meaning of the Γ vs. Σ SFR relation.  simulate an isolated overdense region within a box of 4 comoving Mpc across, using a Eulerian gas dynamics and N-body adaptive refinement tree (ART) code. They incorporate many state-of-the-art physical processes, such as non-equilibrium chemical networks and radiative transfer at very high spatial resolution (about 5 pc at the redshift range of the run).
Due to the computational cost, the simulation runs included in the left plot of Fig. 9 stop at redshift z ∼ 1.5. Each run has different assumptions for the star formation efficiency within each grid cell (see , for more details). Clusters are not analytically implemented in the run but "form" in the cells with the highest gas densities. Only a fraction of the formed stars belong to a cluster, and this fraction is determined locally by considering the star formation efficiency and gas condition within each cell. The Γ and Σ SFR are determined at each snapshot within the region of the simulation that is gravitationally part of the central galaxy.  note that different prescriptions for the SFE have no significant impact on the SFH of the central galaxy, which also follow the Schmidt-Kennicutt law, but it significantly affects the normalisation of the positive relation between Γ and Σ SFR . This means that models with higher star formation efficiency have higher fraction of stars belonging to bound clusters. Therefore, as already suggested by the observational data, the change in the gas conditions (reflected by the observational quantity Σ SFR ) will change the clustering efficiency of the star formation. A direct comparison between the literature data and the simulations by  shows that simulations with low SFE (SFE10) cannot explain the currently observed Γ and Σ SFR . Yet, from observations it is not clear how clustered star formation proceeds at the very low Σ SFR ranges (< −3.5 in log scale). On the other hand, the remaining sets of simulations reproduce quite well the space occupied by the observations. They suggest that the observed scatter in the data might be the result of varying SFE (as predicted by Kruijssen 2012a).
The simulations (square symbols color-coded accordingly to their host galaxy stellar mass at z = 0) in the right plot of Fig. 9 are taken from the E-MOSAICS simulations and were presented in . The E-MOSAICS project ) couples the cosmological, hydrodynamical simulations EAGLE  to cluster formation and evolution via analytical implementations (Kruijssen et al. 2011;Kruijssen 2012a;Reina-Campos and Kruijssen 2017). They utilize a sub-grid model where a fraction of the stellar particles formed by the simulation are converted into a cluster population. The analytical model that predicts CFE by Kruijssen (2012a) is used to set the fraction of the stellar mass formed in clusters. The fiducial runs (showed in Fig. 9) use the analytical model of Reina-Campos and  to predict the maximum mass that a cluster population can have and sample the cluster mass distributions according to a Schechter (1976) initial cluster mass function. Once formed, clusters evolve together with their host systems. The simulations also account for cluster mass loss and disruption because of interactions with the tidal field of their host galaxy, and internal processes such as stellar evolution and tidal evaporation.
The subgrid cluster formation models used in the E-MOSAICS simulations have an environmental dependence on the local gas conditions (density, pressure, dynamical state) within each galaxy. In particular, the data included in the plot here are published in  and are obtained by a suite of zoom-in re-simulations of 25 galaxies selected at redshift z = 0 to have Milky Way-mass halos. The CFE is determined locally in each simulation based on the gas conditions and used to determine what fraction of stellar mass "forms" in clusters. The values reported in the right plot of Fig. 9 are estimated in galaxies at z = 0, using clusters with ages younger than 300 Myr. In general, we observe that the majority of the E-MOSAICS data cluster around a CFE between 1 and 10%, quite close to the values reported for spiral galaxies like the Milky Way and the Andromeda galaxy. A fraction of the galaxies at z = 0 have also more elevated CFE, tracing the median trend reproduced from the observed data. Due to the choice of simulated galaxies, the simulations do not significantly cover the higher Σ SFR values. Nonetheless, the resulting Γ vs. Σ SFR relation from the E-MOSAICS simulations has a similar normalisation and slope as recovered in the observed data, thereby confirming the dependence of the CFE on the physical parameters used in the model, i.e. Σ gas , Q, and Ω. The inclusion of a third parameter, i.e., the galaxy stellar mass, M * , can help us to make a few more considerations. Galaxies with low stellar masses (< 10 9 M e.g., dwarf systems with sub-solar metallicities) rarely reach CFE of 10%. This implies that cluster formation in these systems is highly stochastic, as indeed reported by Cook et al. (2012Cook et al. ( , 2019. On the other hand, CFE can change from a few to very high rates in galaxies with larger stellar masses, indicating that the positions of galaxies in the Γ vs. Σ SFR relation will change as a function of the evolutionary phase that the galaxy experiences. This was indeed recently outlined also in the simulations analysed by Lahén et al. (2019).

The Mass Distributions of Young Star Clusters in Local Galaxies
The mass distribution is a fundamental observable of a star cluster population because it links the formation of bound stellar systems to the star formation process. As already discussed in Sect. 4, GMC populations have mass distributions described by a power-law slope close to α = −2 (Kennicutt and Evans 2012). This slope is similar to the slope recovered for the luminosity function of HII regions and star-forming clumps in the local Universe (e.g., Thilker et al. 2002;Bradley et al. 2006;Elmegreen et al. 2006) and mass distributions of stellar clumps at redshift ∼ 2 (Dessauges-Zavadsky and Adamo 2018).
From the very beginning, YSC luminosity and mass functions appeared to be consistent with a power law function of slope −2 (see review by Adamo and Bastian 2018). The power-law shape that characterises the distributions of masses and luminosities, from largest coherent star-forming regions to the densest and compact structures such as star clusters, down to clumps of proto-stars inside embedded clusters, is very likely the result of fragmentation produced by the balance between gravitational collapse and turbulence compression (Elmegreen 2011).
Typically, the shape of the YSC (and OC in our Milky Way) mass function is thought to closely follow the initial intrinsic mass function at formation. Evolved systems, like GCs, have mass distributions that are better described by a log-normal function, peaked at a luminosity (mass) values that remains almost unchanged across the local Universe (Brodie and Strader 2006). Whether the GC mass function is the result of cluster mass loss and dissolution, remains still under debate.

Observations of the Cluster Mass Function from Local Galaxies
It is widely accepted that the shape of the cluster mass function can be describe to first order by a power law. However what has emerged in the last decade is that a pure power-law distribution is not be sufficient to explain the dearth of very massive clusters both in our own and in local galaxies (Larsen 2006(Larsen , 2009. A power law with a truncation at the high-mass end of the cluster mass distribution, in the shape of an exponential cut-off above a characteristic mass, M c (typically referred to as a Schechter function) might be a more realistic representation of the true cluster mass function. First evidence of a possible truncation at the high-mass end was suggested by Gieles et al. (2006) and Gieles (2009). It was also suggested that the M c would change as a function of galactic environment, i.e. higher M c would be found in galaxies experiencing elevated star formation events (like the Antennae, see e.g., Larsen 2009;Portegies Zwart et al. 2010).
In the literature, however, discrepant conclusions are reached by different authors. From the cluster mass distribution analyses in two grand-design spiral galaxies, M51 and M83, Chandar et al. ( , 2016 report that there is no clear evidence of the presence of a truncation in the cluster mass function and that a pure power-law shape with slope consistent with −2 is statistically the preferred solution. Similar conclusions are reached by Mok et al. (2019) on a small sample of local galaxies. On the other hand, Bastian et al. (2012), , Messa et al. (2018a), report to have determined a M c ∼ 10 5 M for both M83 and M51 cluster populations. Variations in the recovered M c have been found in M83 , where M c decreases at increasing galactocentric distance. No significant variations have been found across M51 cluster disk population (Messa et al. 2018b). Johnson et al. (2017a) reported the lowest yet determined M c = 8.5 × 10 3 M in the cluster population of M31.
The exact description of the cluster mass function has paramount implications not only for our understanding of cluster formation and evolution, but also for stellar feedback. Numerical simulations that aim at understanding the link between the multi-phase ISM and the stellar feedback have shown that the clustering of massive stars and SNe, like in massive star clusters, is a key player in the star-formation cycle of galaxies (Krause et al. 2013;Gentry et al. 2017;Kim et al. 2017). Because star clusters host a large fraction of the massive stars forming in the host galaxy, especially massive clusters are fundamental units to maintain a multi-phase ISM and regulate the star formation process. The difference between a pure power-law mass function and a Schechter-type function provide very different predictions on the number of massive clusters that will form in a galaxy (e.g., see discussion in Johnson et al. 2017a;Adamo et al. 2017).
To understand why contrasting conclusions are reached on the intrinsic shape of the cluster mass function, it is important to make some considerations both on the size of the sample typically analysed and the methods used. As mentioned above, cluster formation is a stochastic process, which implies that galaxies with higher SFR will be able to form more numerous clusters and therefore better sample the mass function at the high-mass end. The implications are twofold.
1. Very often the number of detected clusters in the field of view barely reach 100 objects (e.g., see Cook et al. 2019, or Adamo et al. 2020c to be subm.) either because the galaxies intrinsically have small cluster populations or because the HST imaging coverage is limited to a portion of the disk of the system. Small numbers implies a poor sampling of the mass function and thus a degeneracy when fitting for two parameters (M c and α) instead of one (α, see appendix in Messa et al. 2018b). 2. Until very recently, the fit is done on binned distributions. However, as pointed out in Adamo et al. (2017), in case of small number statistics a binned distribution always bias the fit against the upper mass-end of the distribution. Equal size bins will always have higher error bars at the high-mass end because of the small number of objects they contained and therefore weight less significantly on the resulting slope. Equal-number of object bins have also a biased impact because they tend to become very large at the high-mass end, washing out the presence of a truncation.
To overcome these limitations, different methods have recently been implemented. For example, by fitting cumulative mass distributions , using maximumlikelihood methods applied to cumulative distributions (Adamo et al. 2017;Messa et al. 2018a), and a maximum-likelihood analysis combined with a Markov Chain Monte Carlo (MCMC) technique to sample the posterior probability distributions of the Schechter and power-law mass function parameters (e.g., see Johnson et al. 2017a;Messa et al. 2018b, for more details). In the latter case, the analysis does not depend on the binning technique used or the functional distribution applied, therefore, overcomes several limitations discussed  Table 2 and obtained by combining the cluster catalogues of 8 LEGUS spiral galaxies with SFR < 0.5 M yr −1 (lowSFR spirals, black cross), 6 spiral galaxies with SFR > 0.5 M yr −1 (highSFR spirals, black square), 6 merger systems belonging to the HiPEEC sample (SFR > 10 M yr −1 , black diamond). The horizontal bar associated with each black symbol shows the Σ SFR interval that each sample spans (see main text). The range of Σ SFR of the dwarf sample is plotted as a dot-dashed horizontal line. Right: Power-law slope of the cluster mass function vs. star formation rate surface density. Values plotted in this plot are reported by Johnson et al. (2017a), and Adamo et al. (2020a, c, to be subm.) above. However, as pointed out in Johnson et al. (2017a) and Messa et al. (2018b) the convergence of the fit will still be affected by small number statistics, i.e. if the cluster sample is small, the M c parameter remains unconstrained and single power-law fit is preferred. In Fig. 10, we compile a sample of statistically significant recovered M c and α by fitting cluster catalogues available in the literature with the method introduced by Johnson et al. (2017a). For the sake of homogeneity and consistency, we include in these plots only data that have been fitted with the same method. The M c best values are sampled out of the marginalized posterior probability distribution function (PDF) for each of the Schechter function parameters, accompanied by a 1σ confidence interval defined by the 16th to 84th percentile range. For α PL we report the median and 1σ confidence interval for the single parameter, α (see respective source papers Johnson et al. 2017a;Messa et al. 2018b, Adamo et al. 2020a for details on the cluster mass function analysis). In the left plot of Fig. 10, we plot M c vs. Σ SFR . The Antennae (orange diamond) and the M83 data (magenta circles) are taken from the literature and included in the plot for completeness because they were used by Johnson et al. (2017a) to derive an analytical formulation of the positive correlation between M c and Σ SFR . As reported in Adamo et al. (2020a), of the 6 HiPEEC galaxies analysed only 4 have a statistically meaningful constraint on M c . Those values are included in the plot along with the M c value (black diamond) obtained by fitting the combined cluster catalogues of the 6 merger systems (we report in Table 2 the SFR range of the galaxies, the age range of the combined clusters, the resulting number after a conservative mass limit cut at M = 5 × 10 4 M has been applied to mitigate the effect of incompleteness due to detec- Table 2 Maximum-likelihood fit outputs of the cluster mass function in diverse galaxy environments (see Sect. 8). Clusters populations of similar galaxy types (see column 1) have been combined to increase the statistical significance of the fit. The cluster catalogues of dwarf and spirals systems have been obtained from the LEGUS survey. The merger/starburst system catalogues from the HiPEEC survey. The columns list, in order, the galaxy SFR range, the cluster age range used, the total number of clusters included in the fit, the determined M c and α Sch , the slope of the power-law fit, α PL . The last column summarises whether the two parameters (M c and α Sch ) describing the Schechter function are uniquely determined (therefore, statistically significant) or cluster mass distribution is better described by a power-law function (Sch or PL, respectively).  Sch tion limits). Within the LEGUS survey we have analysed in total 31 galaxies. We report here the results obtained by using only systems classified as class 1 and 2 (compact cluster candidates), with masses M > 5000 M and ages younger than 200 Myr. First, the marginalized posterior probability distribution analysis of the Schechter function parameters of the cluster population of each of the 17 dwarf systems published by Cook et al. (2019) does not produce a statistically meaningful constraint on the M c . It is also important to note that the fit to the combined cluster catalogues of the 17 LEGUS dwarf galaxies does not provide a tight constraint on M c , as reported in Table 2 and in Cook et al. (2019). In total 14 LEGUS spiral galaxies have complete cluster catalogues. Of those only 3 galaxies (M51, NGC 628, NGC 1566) have converging marginalised PDF for the M c . These values are reported in the plot as blue diamonds. The spiral sample is then divided in two groups based on their SFR (determined within the HST field of view where the cluster population has been analysed). The combination of cluster samples has the advantage of increasing the number of clusters for a given SFR range so that it mitigates the size-of-sample effects mentioned above. Indeed, the convergence of the marginalized posteriors of the Schechter function parameters in both samples confirms that the number statistics is the main problem in determining these parameters in local galaxies as also pointed out by Larsen (2009), Adamo et al. (2017, and Elmegreen (2018). The M c retrieved for both sub-samples of spirals is included in the left plot of Fig. 10 as black square (high SFR spiral sample) and cross (low SFR spiral sample). The addition of new datapoints confirm the overall positive correlation between M c and the Σ SFR of the host galaxy as proposed by Johnson et al. (2017a). If we consider the Σ SFR as a tracer of gas surface density (and pressure), the relation suggests that galaxies experiencing elevated star formation episodes, high gas density (pressure) have higher probability to form massive star clusters. From Table 2, we can see that dwarf galaxies have modest SFR values, which explains the small number of clusters that are available for the fit. However, dwarf galaxies have also very compact star-forming regions, resulting in Σ SFR overlapping with populous spiral galaxy systems. Extrapolating the M c vs. Σ SFR relations for the starburst dwarfs, like NGC 4449, NGC 5253, NGC 4656, with the highest Σ SFR , we expect them to have M c ∼ a few times 10 5 M . However, their combined population counts ∼ 200 clusters (as opposed to the almost 2500 clusters for the high SFR spirals), i.e. too small to provide a good constraint on the M c . So cluster formation in dwarf galaxies is either fundamentally different than in spiral and merger systems, or simply more affected by stochastic effects and low number statistics which prevent to make definitive conclusions. On the right side of Fig. 10 we collect the recovered α PL obtained by marginalising the posterior distributions assuming a power-law function of 17 LEGUS galaxies and 6 HiPEEC systems. We remind the reader that this is in many galaxies not the best description of their cluster mass distributions. We assume a power-law function because for a Schechter function the M c and α are degenerate parameters, therefore, we want to prevent mixing slopes derived by assuming different functions. A similar plot was presented by  using data published in the literature. There, they noticed that the slope of the cluster mass function appears to fattens in environments with increasing Σ SFR . The advantage, here, is that we adopt the same method across a large galaxy spectrum and control the lower mass limit of each sample against completeness issues. As one can see from the values listed in Table 2, and showed by numerical exercises Johnson et al. 2017a;Messa et al. 2018b), α PL is systematically steeper than α Sch obtained for the same sample. Therefore, the trend observed in the right plot of Fig. 10, would not disappear even if we could ideally determine both M c and α for the entire sample. Our homogeneous sample confirms the trend outlined by . The change in the star formation conditions (higher Σ SFR corresponds to increasingly elevated Σ gas and therefore pressures) are reflected in shallower slopes and a larger number of massive clusters.
It is interesting to compare the properties of the cluster mass function to those of the GMC mass spectrum reported in Sect. 4. There is not a one-to-one correspondence between a GMC and the formation of a star cluster. The fragmentation process combined with the star formation efficiency within each collapsing core could result in one bound cluster, several, or none. Irrespectively, some fraction of these dense gas regions will be able to form star clusters. The mass distributions and other properties of cluster populations (like the presence of a truncation mass) correlate in some way with the variations observed in the GMC mass function, as a function of galaxy scale dynamics, gas content, and Σ SFR . For completeness, we report that different conclusions have been recently reached by Mok et al. (2020), who suggest that GMC and YSC populations are unrelated to the global properties of the galaxies where they form. Clearly, the field has not yet converged on a single interpretation and definitive answers may require larger samples subjected to a homogeneous analysis.
Finally, the presence of a truncation M c has also been observed in GC populations of galaxies belonging to the Virgo Galaxy cluster (Jordán et al. 2007). The GC M c appears to positively scale with galaxy stellar mass. It is not yet clear whether the GC M c is a property of the initial cluster mass function of the progenitor clusters or the results of evolution. Observationally both Johnson et al. (2017a) and Adamo et al. (2020c, to be subm.) report no significant evolution within the uncertainties on the recovered M c as a function of age bins. Numerical works agree that significant fractions of massive clusters that could be considered GC progenitors are formed during the merger events that determine the assembly history of the hosts (e.g., Forbes and Bridges 2010;). However, numerical simulations also agree that in spite of their large masses these systems have very slim chances to survive if they remain trapped in the dense gas environments where they form (e.g., Kruijssen et al. 2011Kruijssen et al. , 2012cRenaud et al. 2015). Based on this evidence Lamers et al. (2017) argue that the recovered M c of GC populations may be the result of nurture (i.e. cluster disruption).

The Origin of the M c ; From Observations to Theory
To understand what regulates the formation of the most massive clusters and the origin of the M c different theoretical models have been put forward.
In the model initially proposed by Kruijssen (2014) the maximum cloud mass and the stellar M c might have a common origin and it may correspond to the Toomre mass. The latter depends on the maximum size of a region that will overcome the shear of the disk and the kinetic pressure of the gas and start collapsing (Toomre 1964). In this model for a given gas surface density, this length-scale directly provides the Toomre mass. This shear-dependent model was subsequently refined by Reina-Campos and Kruijssen (2017). In addition to the shear-limited maximum mass model they take into account that feedback from young stars might disrupt the cloud before the global collapse of the shear-limited area is completed. If the feedback time (i.e. the time it takes for the stellar feedback to destroy the cloud) is smaller than the free-fall time of the shear-limited region the resulting collapsed mass is smaller than the shear-limited Toomre mass. The model is able to reproduce the observed trend of maximum GMC mass and maximum cluster mass as a function of increasing galactocentric distance in M31 (maximum GMC and cluster mass peak at the star-forming ring), M83 (declining maximum GMC and M c as a function of distance from the centre Freeman et al. 2017), and M51 (only small variations as a function of distance from the centre Messa et al. 2018b). The Reina-Campos and  model has analytically been implemented in the E-MOSAICS simulations ). In  the authors recover similar M c and α vs. Σ SFR relations as showed here for observations in Fig. 10. They suggest that not only gas surface density and pressure (via the dependence from Σ SFR ) play a role in determining the maximum GMC mass and maximum cluster mass, but also the angular velocity of the host galaxy from which the Toorme size depends.
Recently, Trujillo-Gomez et al. (2019) have extended the Reina-Campos and Kruijssen (2017) model to investigate the effect of feedback on the formation of low mass clusters, hence the intrinsic shape of the cluster mass function at the low mass end. The model evaluates which parts of the star-forming region remain bound given the time-scales for gravitational collapse, star formation, and stellar feedback that also determine the upper mass distribution. In this model, galaxies like the dwarf Fornax in the Local Group, with high specific frequency of GCs, might have never formed the lower mass counterpart expected for a power-law cluster mass distribution.
Based on observational evidence, Elmegreen (2018) derives an analytical model that is based on the minimum SFR and Σ SFR necessary to form a GC progenitor, i.e. a star cluster of ∼ 10 6 M . A SFR of 1 M yr −1 is necessary in a region to sample the cluster mass function to the high-mass end. A minimum Σ SFR of 1 M yr −1 kpc −2 will ensure that the gas density and therefore pressure is high enough to form bound stellar systems with stellar densities similar to GCs. It is therefore the gas pressure in the ambient medium that determine M c . If the SFR is elevated but spread over a large area (low Σ SFR ) then the pressure-limit mass, that correspond to the M c , is never reached and the most massive cluster is simply determined by size-of-sample effects. If both conditions are satisfied, galaxies will be able to sample their cluster mass function up to the M c , which should then be traced in the cluster mass function analysis. The predicted trend of the Elmegreen (2018) models agree with the observed distribution of observed datapoints in the M c vs. Σ SFR plane of Fig. 10 reinforcing the idea that pressure could be a driving parameter for the change of M c .
Finally it is interesting to note that, at difference of the E-MOSAICS simulations where the cluster mass function is analytically implemented, the cosmological simulations of  find steeper cluster mass functions than the canonical −2 slope, which they interpret as a possible presence of a truncation, though quantitatively it remains unconstrained.

The Survival Rates of Young Star Clusters Forming in Local Galaxies
Observationally, cluster disruption has traditionally been traced using the demographics of the cluster population, such as their combined age and mass distributions or cluster frequencies per age dex (e.g., Lamers et al. 2005b,a;Whitmore et al. 2007;Bastian et al. 2012;Fouesneau et al. 2014;Chandar et al. , 2017Messa et al. 2018b;Cook et al. 2019). As summarised in recent reviews (e.g., Portegies Zwart et al. 2010;Adamo and Bastian 2018;, there is no observational consensus on the way clusters disrupt. Observationally, different teams reach different conclusions by counting numbers of clusters in age and mass bins. The debate has focused on whether clusters simply disrupt because gas expulsion has destabilised their potential (e.g., Fall et al. 2005;Chandar et al. 2010) and therefore disrupt independently of their mass and environment, or mainly because they lose mass via interaction with their galactic environment (e.g., Baumgardt and Makino 2003;Lamers et al. 2005b,a;Gieles et al. 2006, resulting in a gradual mass loss, otherwise known as 'mass dependent cluster disruption'. While observationally there is no consensus yet, Adamo and Bastian (2018) noted, based on a compilation of results available in the literature, that the cluster disruption rate may change as a function of environment, becoming progressively higher in galaxies with higher gas surface densities (and thus higher external pressures), as suggested by Elmegreen and Hunter (2010) and Kruijssen et al. (2011). Indeed, Messa et al. (2018b) analysing the number of clusters per age bin and as a function of position in the disc of M51, find a positive correlation between the increasing number rate and the gas surface density of the region considered (also see Miholics et al. 2017). We are therefore moving away from the initial dichotomy (massdependent as apposed to mass-independent disruption) into an environmental dependence of the disruption rate of clusters, that goes from considerably longer dissolution timescales in low-pressure galaxies to very fast disruption rates in highly star-forming galaxies like the Antennae.
From a theoretical perspective, after their formation clusters lose mass due to stellar evolution, two-body relaxation-driven evaporation in the host galaxy's tidal field, and tidal perturbations ('shocks') resulting from GMC passages, spiral arm passages, pericentre passages, and disc crossings (e.g., Spitzer 1958;Gnedin and Ostriker 1997;Baumgardt and Makino 2003;Gieles et al. 2006Kruijssen et al. 2011;Krause et al. 2020). Recent analyses of Gaia's second data release reveal tidal tails emerging from clusters in the Galactic disc (e.g., Fürnkranz et al. 2019;Röser et al. 2019). After carrying out a completeness correction (which is non-trivial), the detailed phase space information of these tails can potentially be used to directly measure dynamical mass loss rates. A rich variety of theoretical work has investigated dynamical cluster mass loss. Next to predictions for cluster mass loss rates from analytical theory (e.g., Gnedin and Ostriker 1997;Lamers et al. 2005a;Kruijssen 2015;Gieles and Renaud 2016), a wide variety of numerical approaches has been used, ranging from direct N -body simulations (e.g., Baumgardt and Makino 2003;Gieles and Baumgardt 2008;Renaud and Gieles 2013;Rieder et al. 2013;Miholics et al. 2016;Webb et al. 2019) to hydrodynamical simulations that link cluster mass loss to the formation and evolution of the host galaxy (e.g., Kruijssen et al. 2011Miholics et al. 2017;Reina-Campos et al. 2018Li and Gnedin 2019).
Modeling cluster disruption is critical for reconstructing the properties of cluster populations at the time of their formation from the observed cluster demographics. This holds for young cluster populations, which can experience significant attraction due to tidal shockdriven disruption on timescales as short as 50-100 Myr (Gieles et al. 2006;Miholics et al. 2017;), but even more so for old GC populations, which have undergone ∼ 10 Gyr of evolution before attaining their present-day properties. The physical mechanisms driving cluster disruption are discussed in more detail in Krause et al. (2020), another review in this series. As stated above, these are a combination of stellar evolution, tidal evaporation, and tidal shocks. Numerical models need to include all of these processes, including their environmental dependence. This means that a complete model for cluster formation and destruction requires self-consistent simulations of galaxy formation and evolution, including descriptions for star formation, feedback, chemical enrichment, cluster formation, cluster disruption, and stellar evolution. Recent years have seen major steps in this direction (e.g., Li and Gnedin 2019), and we will discuss these efforts in more detail in Sects. 10 and 11.

Current and Future Observations and Simulations of Massive Cluster Formation Across Cosmic Time
Direct detections of proto-GC candidates are reported in the literature by Vanzella et al. (2017), Bouwens et al. (2017), Johnson et al. (2017b), Vanzella et al. (2019), and Vanzella et al. (2020). All these detections rely on the aid of gravitational lensing, which can magnify several times the fluxes of these very young systems in their UV and blue optical rest-frames. Indeed because of their compact sizes these massive and compact star-forming regions have larger chances to be magnified and detected and might become the only signposts of their otherwise faint host galaxies (Zick et al. 2018). The intrinsic sizes of these systems are very often uncertain because they rely on lensing models and are limited to upper limits of a few tens of parsec, i.e. the sizes of entire star-forming regions. Spectroscopic studies of a handful of these lensed proto-GCs at redshift ∼ 2-3, show that their FUV light is dominated by spectral signatures of very massive stars (with lifetimes shorter than 4-5 Myr), WR HeII broad wind emission, and direct evidence of Lyman continuum escape radiation (Rivera-Thorsen et al. 2019;Vanzella et al. 2020). These initial results put proto-GCs as potential contributors to the reionisation era, or towards the end of it ( 7). Indirect evidence of proto-GCs formation at high redshift is produced by the physical properties of stellar clumps. Stellar clumps dominate the UV rest-frame of their host systems (e.g., Elmegreen et al. 2013;Shibuya et al. 2016;Messa et al. 2019) and have very high SRF surface densities, making them a natural sight for very massive cluster formation (Elmegreen 2018). Indeed, assuming that proto-GCs form in stellar clumps following a Schechter mass function distribution and accounting for cluster disruption, it produces a first order calculation of the surving GC populations at redshift z = 0 that are comparable to the number and mass functions of the GCs detected around galaxies like the Milky Way and M31 (e.g., Shapiro et al. 2010;Adamo et al. 2013). At least in part motivated by the impending launch of JWST (and facilitated by technical improvements in hardware and modelling), a wide variety of recent studies has made predictions of the properties of young GC populations at high redshift. These range from back-of-the-envelope estimates, assuming simple scaling relations (e.g., Katz and Ricotti 2013;Renzini 2017;Boylan-Kolchin 2018), to expectations based on observations of nearby GC populations (e.g., Zick et al. 2018) and high-redshift star formation (e.g., Vanzella et al. 2017), and predictions from numerical simulations of galaxy formation, either by 'painting on' GCs using an ad-hoc model (e.g., Renaud et al. 2017;Halbesma et al. 2019;Madau et al. 2019;Phipps et al. 2019) or by including a physical description for GC formation and evolution (e.g., . The resulting predictions for when GCs formed and what their corresponding detectability is with future observations with JWST or 30-m class telescopes vary greatly between these different approaches. The reason is two-fold. 1. Different models assume different formation scenarios for GCs. Recently, numerous works suggest that GCs are the natural byproduct of 'normal' star formation under the high-pressure conditions in gas-rich high-redshift galaxies (e.g. Shapiro et al. 2010;Kruijssen 2015Kruijssen , 2019Li et al. 2017;El-Badry et al. 2019;Keller et al. 2020). However, additional formation scenarios that have been considered are mergers of galaxies or dark matter substructure (e.g. Kim et al. 2018;Li and Gnedin 2019;Madau et al. 2019), or special conditions during reionisation (e.g. Katz and Ricotti 2014;Trenti et al. 2015;Kimm et al. 2016;Ricotti et al. 2016;Creasey et al. 2019). All of these scenarios predict increased GC formation towards high redshift, because of an increase in gas pressure, merger rate, or specific early-Universe conditions, but the specific redshift range in which GCs are predicted to have formed still varies, from slightly preceding (but largely tracing) the cosmic star formation history (e.g.  to extending deeply into the epoch of reionisation (e.g. Katz and Ricotti 2014). In the latter case, it is possible that GCs may have played an important role in reionisation (e.g. Ricotti 2002;Ricotti 2013, 2014;He et al. 2020). 2. Even within families of models adopting a similar formation scenario for GCs, there exist considerable differences in terms of GC formation redshifts and luminosity functions. For instance, numerical simulations of galaxy formation that describe GCs as the natural outcome of regular star formation roughly fall into two categories. 2 The first category of models takes a 'normal' galaxy formation simulation and uses some set of conditions (e.g. cuts in age, metallicity, or halo mass) to associate GCs to star particles (e.g., Renaud et al. 2017;Halbesma et al. 2019;Madau et al. 2019;Phipps et al. 2019). This 'particle tagging' technique has had difficulties to reproduce the observed demographics of GCs, such as their total number per unit galaxy mass and their age or metallicity distribution. If an age cut is made, this obviously directly defines the redshift range where GCs are expected to be observed. The second category of models employs a sub-grid model for stellar cluster formation and disruption within galaxy formation simulations (e.g., Li et al. 2017;. In general, these models predict a lower formation redshift of GCs (e.g., , as well as lower initial numbers and masses of proto-GCs than in other models (e.g., ).
Unfortunately, age measurements of old GCs in the local Universe have uncertainties of ∼ 1 Gyr, which limits their diagnostic power to distinguish between these ideas. However, the observation that GC formation extends well down to z ∼ 1 (Marín-Franch et al. 2009;Forbes and Bridges 2010;Leaman et al. 2013) and still takes place in the present-day Universe (e.g., Ashman and Zepf 1992;Elmegreen and Efremov 1997;Kruijssen 2014) suggests that regular star formation under high-pressure conditions can indeed lead to GC formation. While it is currently under discussion as to whether or not the requirement of a high gas pressure still makes GCs analogues to young massive clusters in the local Universe (Renaud 2019), a consensus is emerging that there is no evidence for two 'modes' of cluster formation, but GCs rather seem to form at the extreme end of a continuum. Occam's Razor thus suggests that the models describing GCs as the outcome of regular star formation provide the most accurate description of GC formation, without a major contribution by mergerinduced cluster formation (Keller et al. 2020). While some of the more exotic formation mechanisms may have contributed a handful of GCs, it is unlikely that they dominated the formation of the GC population in galaxies (Keller et al. 2020). Fortunately, the great variety of predictions for the initial demographics of the GC population, and in particular their UV luminosities (compare e.g. Renzini 2017 and) and the redshift range of their formation (see the above discussion), means that future observations of young GCs in the early Universe with JWST, the ELT, the TMT, and the GMT will be able to distinguish between the currently considered formation scenarios.

Cluster Populations as Tracers of Galaxy Assembly
The discussion in this review so far has emphasised the prominent environmental dependence of both cluster formation and their dynamical disruption. This environmental dependence is critical in the context of galaxy formation. Firstly, it implies that cluster formation and evolution cannot be described in isolation, but require the multi-scale coupling to the galaxy formation context. Secondly, it implies that (especially old, i.e. globular) cluster populations can be used as tracers of galaxy formation and evolution, provided that a comprehensive model for GC formation and evolution can be constructed.
The ultimate goal of comprehensive models of GC formation and evolution in the context of galaxy formation is to simultaneously reproduce the demographics of the evolved GC populations observed in the local Universe and those of the initial GC populations that will be observed with JWST and 30-m class telescopes (see Sect. 10). Thanks to considerable efforts aimed at integrating GC formation and evolution in cosmological simulations, this goal is now starting to come within reach, and at a very opportune moment, given the impending arrival of the next generation of major facilities. However, despite significant progress, current numerical simulations still face a number of shortcomings.
As stated, GC formation and evolution is one of the greatest multi-scale problems in astrophysics, spanning the scales of black hole binaries (< 0.1 pc) to galaxy formation (∼ 1 Mpc). This large dynamic range implies a great computational cost when trying to resolve the dynamics of individual stars within GCs, and doing so self-consistently in relation to galaxy formation and evolution for the entire GC population of a galaxy throughout its history. This cost is so prohibitive, that numerical studies need to choose to either model the formation and evolution of individual clusters at high resolution (e.g., Kim et al. 2018;, or to employ a sub-grid prescription for the formation and evolution of the entire cluster population (e.g., Li and Gnedin 2019). Both approaches are entirely complementary, but they suffer from different problems.
1. Simulations that focus on maximally resolving the relevant physics obviously provide the most fundamental descriptions of GC formation, but they cover a limited range in terms of the number of clusters and the redshift interval, prohibiting predictions for the demographics of the GC population at any redshift, such that the z = 0 GC population is completely out of reach. In addition, these studies often need to focus on a single formation environment, which means that the dependence of GC formation and evolution on the formation and assembly history of the host galaxy is not accounted for. Therefore, the results may not be statistically representative. 2. Simulations that instead employ a sub-grid model for cluster formation and evolution can predict the demographics of the cluster population at all redshifts (Li et al. 2017Li and Gnedin 2019;. If a sufficiently large number of galaxies is modelled, they can even do so as a function of the galaxy formation and assembly history ). However, these simulations are fundamentally limited in terms of the (unresolved) physical process that are described by the sub-grid physics. Most critically, GC formation and evolution relies on the interaction with the cold ISM, either because it sets the gas pressure from which clusters form (e.g., Elmegreen and Efremov 1997;Kruijssen 2012a;Johnson et al. 2016Johnson et al. , 2017a, or because it sets the disruption rate of the resulting clusters (e.g., Gieles et al. 2006;Kruijssen 2015). To date, simultaneously modelling the physics of the cold ISM and GC formation down to z = 0 has been too expensive to enable multiple galaxies to be simulated (e.g., Li and Gnedin 2019). Even the most sophisticated suites of cosmological simulations modelling large numbers of galaxies and their GC populations currently lack a cold ISM , which most prominently limits their ability to accurately describe cluster disruption, especially at high metallicities (see  for a discussion). The next generation of galaxy formation simulations will overcome this problem (Reina-Campos et al. in prep.).
Despite the need for further development, the current generation of galaxy formation simulations describing the formation and evolution of the GC population that include a large number of galaxies (developed as part of the E-MOSAICS project, ) has already revealed a number of quantitative connections between the present-day properties of the GC population and the assembly history of the host galaxy. In many ways, this is the most explicit realisation to date of the promise that GCs can be used as tracers of galaxy formation and assembly. Many of these predicted correlations rely on the diagnostic power of the age-metallicity distribution of GCs (e.g., , also see Choksi et al. 2018 for an earlier discussion on the importance of the GC age-metallicity distribution in this context). Specific examples of strong correlations are: 1. The total number of mergers experienced by the host galaxy is traced by the slope of the GC age-metallicity distribution and the total number of GCs. 2. The fraction of these mergers taking place before z = 2 can be estimated by measuring the median age of the GC population in a galaxy. 3. The median age of the GC population also probes the redshift by which certain fractions of the host galaxy mass have been assembled, thus tracing its assembly rate.
These examples are non-exhaustive, and upcoming works will extend the range of correlations beyond the GC age-metallicity distribution, also including GC kinematic and spatial distributions (Trujillo-Gomez et al. 2020, Reina-Campos et al. in prep.). Most importantly, these results are not affected by the current omission of a cold ISM in these simulations, because the global GC demographics trace the assembly history of the host galaxy even if GC disruption is not sufficiently efficient.  applied the insights from E-MOSAICS to the GC population of the Milky Way, with the goal of constraining our galaxy's assembly history and reconstructing its merger tree. The analysis showed that the Milky Way experienced a total of ∼ 15 mergers throughout its history, with the last major merger having taken place at z > 4, thus sharpening previous limits of z > 2 (e.g., Wyse 2001;Hammer et al. 2007;Stewart et al. 2008). Improved determinations of the ages of extragalactic GCs may enable the application of these correlations to galaxies beyond the Milky Way in the near future (Usher et al. 2019). In addition, the analysis presented by  predicted the existence of the satellite Gaia-Enceladus/Sausage (GES), which was discovered only weeks later in the data from Gaia's second data release (e.g., Myeong et al. 2018;Helmi et al. 2018), as well as the enigmatic galaxy Kraken, which together with the GES accretion event forms the two most massive galaxies ever accreted by the Milky Way. Recent analyses of the phase space distribution of GCs in the Gaia data by Massari et al. (2019) suggest that the relics of Kraken have been found ).
An extremely promising avenue of research is the association of GCs in the Galactic halo with fossil stellar streams from disrupted satellite dwarf galaxies that were accreted by the Milky Way. Gaia now enables the identification of these systems, and galaxy formation simulations including a model for the GC population provide quantitative predictions for the properties and demographics of such streams and their GCs (e.g., Hughes et al. 2019). In the long run, extremely sensitive, wide-field photometry with LSST will aid the efforts to link the GC population to accretion events, not just in the Milky Way, but also in other galaxies.
In summary, the falsification and validation of GC formation models in the context of galaxy formation and evolution are now becoming possible. The combination of revolutionary observational data from Gaia, JWST, LSST, and 30-m class telescopes will provide a comprehensive picture of GC populations in the local Universe and their progenitors at high redshift, for the first time enabling models to be tested simultaneously at the time of GC formation and after nearly a Hubble time of evolution. In addition, observations spanning a wide redshift range will provide direct tests of the connection between GC demographics and the evolution of the host galaxy. It is to be expected that such observations will cement the use of GCs as quantitative tracers of galaxy formation and assembly.

Open Questions, Outlook, and Future Steps
The development of new observational techniques such as large IFUs and the increased low-cost computational capacities have made a large impact on the studies of resolved starforming regions. Over the past decade is has become feasible to study the multi dimensional phase space using 3D positions and velocities of stars to measure velocity dispersions to an accuracy of a few km s −1 and refining stellar distances and cluster membership probabilities. Also it allows us to study the star-gas interactions and feedback processes in unprecedented details. Yet, many aspects of the earliest stages of star cluster formation and evolution still remain hidden due to observational limitations leaving unanswered questions: How is the low-mass and brown dwarf mass function constituted? How do massive stars, binary, and higher-order systems influence planet formation? To what extend does star formation progress into the HII regions around a massive star clusters with a high ionizing potential?; to name a few.
To address these unsolved questions, larger and more sensitive telescopes are necessary. Over the next few years some revolutionizing missions will go online: JWST (spring 2021) a new generation 6.5 m class space telescope will observe from the Lagrange point 2. This joint mission between NASA, ESA, and CSA will allow us to observe (photometry and spectroscopy) in the near and mid infrared at a spatial resolution rivaling that of HST. The E-ELT is scheduled to have first light in 2025 observing in the optical and NIR. This 39 m ground-based telescope will be equipped with state-of-the-art adaptive optics systems for high-resolution photometry and spectroscopy. The Wide-Field Infrared Space Telescope (WFIRST) will be a HST-class ultra wide field surveyor designed to settle essential questions in the areas of dark energy, exoplanets, and infrared astrophysics. These new telescopes in combination with the existing ones will allow us to study the lowest-mass star and planet formation processes and track the earliest stages of star formation (Class 0 -III) inside their parental molecular clouds. This will be even pushed further with telescopes planned for the late 2020s and 2030s, such as the TMT and future flagship space missions like the Large UV/Optical/IR Surveyor (LUVOIR), a 9-16 m class space telescope combining the capabilities of HST and JWST at an unprecedented precision, the Origins Space Telescope (OST), the Habitable Exoplanet Observatory (HabEx), or the Lynx X-ray observatory.
Our current understanding of cluster formation and evolution remains still fragmented even if we limit ourself to the local Universe. Thanks to the multi-wavelength and highspatial resolution coverage offered by the HST we begin to witness what unique carrier of information cluster populations are. YSCs are potentially the link between small scale events like the formation of stars and the multiphase-ISM necessary to regulate the star formation process throughout the galaxy assembly history. During their formation and early evolution they can be considered fundamental units of stellar feedback, they are the nurse of very massive stars and, thus, carriers of radiative and mechanical energy. While aging they will keep records of the gas conditions where they formed provided we understand how they evolve. GC populations have encoded in them the history of their host galaxies, yet we are not able to fully map their genetic sequence.
The question is how to move forward. While HST reaches its endpoint, we are gaining a large FUV and optical archive that will allow us to reinforce our studies of cluster populations across all the star-forming galactic spectra available in the local Universe. We will loose for more than a decade our FUV window to the local Universe. We will, however, gain, with the JWST and E-ELT, access to the same window at redshift 2 to 10, where we expect to be able to detect proto-GCs.
For the local Universe, we gain a new wavelength regime at spatial resolutions never achieved before on such large field of views. The NIR and MIR spectral wavelength carry nearly dust-free information into the onset of cluster formation and provide direct statistical constraints on the time scales a cluster requires to emerge from its natal cloud. These time scales can be obtained across a vast sample of galaxies (we are no more limited by the Spitzer resolution to our Local Group), with improved cluster ages (by the use of NIR and MIR indicators) and improved cluster masses (by breaking the extinction-age-metallicity degeneracy plaguing cluster analysis). Multi-object spectroscopy and IFU capabilities of JWST and E-ELT will allow to get spectra of a significant fraction of the cluster population of a galaxy on a time scale comparable to what imaging requires nowadays. Such information is still highly distance dependent and remains limited to a handful of bright YSCs in other galaxies. As we prepare for the next decades of new observatories and instruments, we also need to take major steps in 1. improving stellar population models to account for binary fractions and stochastic sampling of the IMF; 2. introduce cognitive algorithms at all phases of the cluster analysis, from detection to constraining of their physical parameters.
Observations of cluster populations in local and high redshift galaxies are vital to inform theoretical and numerical modeling of galaxy formation and evolution that are challenged to link parsec scale physics to inter-galactic scale processes.
From the theoretical perspective, the past years have seen a major effort in developing numerical simulations of galaxy formation and evolution that simultaneously describe the formation and evolution of GCs at high redshift and of YSCs in the local Universe. The impending observational revolution with the arrival of JWST, LSST, and 30-m class telescopes, combined with the recent observations of molecular gas with ALMA and galactic archaeological relics with Gaia, will for the first time allow the demographics of the cluster population predicted by these state-of-the-art numerical simulations to be tested as a function of cosmic time. This comprehensive range of tests is unprecedented and will lead to the most transformational change in our understanding of cluster formation across cosmic history that the field has ever experienced. In order to realise this major step, a number of key developments are necessary, both on the theoretical and observational fronts.
1. Numerical simulations of galaxy formation must include an accurate description of the cold ISM. The descriptions of cluster formation and evolution must account for the expected environmental dependences of these processes, and include a model for cluster disruption by tidal shocks. Suites of such simulations must cover a sufficiently large number of galaxies to probe how these processes and the resulting cluster demographics depend on the formation and assembly history of the host galaxy. Finally, in order to fully develop our understanding of tidal-shock driven cluster disruption, systematic suites of direct N -body simulations resolving the collisional stellar dynamics are needed. These should fully map a comprehensive parameter space of shock durations, strengths, and successions. 2. Observational programmes targeted at mapping the properties of proto-GC populations at high redshift with the next generation of telescopes should be allocated sufficient observing time to obtain a comprehensive census of the demographics of the cluster population, measuring the luminosity function, peak masses, and redshift distribution of the young GC population. Simultaneously, deep observations with 30-m class telescopes will be able to detect the brightest GCs at intermediate redshifts, revealing how the correlations between GC population properties (e.g. their numbers and the luminosities of the brightest GCs) and the host galaxy properties (e.g. mass, radius, metallicity) evolve with cosmic time. Finally, galactic archaeological surveys in the Milky Way and other galaxies must develop ways of statistically associating GCs with the enormous richness of fossil stellar streams that will be discovered, and thereby reconstructing the properties (e.g. masses and chemical compositions) of the disrupted galaxies that brought GCs into the halo of the central galaxy.
By connecting the above developments, the field will achieve a step-change in our understanding of cluster formation and evolution throughout cosmic history. In particular: 1. We will have a complete picture of the physics driving the cluster demographics as a function of cosmic time and environment. 2. We will be considerably closer to uncovering the intricate connection between stellar cluster populations and their host galaxies. 3. We will be able to use cluster populations as the fossils revealing the formation and assembly histories of galaxies.
HST brought about a true revolution in studies of stellar cluster populations, both young and old. While a much wider range of telescopes is necessary to make the next major step, and this will require a much closer connection to theoretical models than has been common practice in the past, this holistic approach will hopefully make the impact on our understanding of cluster formation and evolution even greater. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.