Role of supergiants in the formation of globular clusters

Multiple stellar populations are observed in almost all globular-clusters, but the origin of this phenomenon is still debated. We investigate the role cool supergiants may have played. To do this, we combine two investigative methods: state-of-the-art massive stellar evolution and calculations of the hydrodynamic structure of the cluster-gas. This approach allows us to study how star-formation in young massive clusters depends on the energy- and mass-input of the first-generation of stars, while predicting the chemical composition of the second-generation. We find that the presence of massive (9-500 M$_{\odot}$) metal-poor supergiants in the young cluster leads to a star-formation episode within the first 4 Myr of the cluster's lifetime, that is, before the first core-collapse supernovae explode or the gas is expelled. The stellar winds accumulate in the cluster center, forming the second-generation there. Its composition is predicted to show variations in O&Na abundances, consistently with observations. The abundance of helium is, similarly to other scenarios involving massive stars, higher than what is referred from observations. Supposing dynamical removal of stars from the outskirts of the cluster, or applying a top-heavy initial-mass-function, we can predict a number ratio of the second-generation as high as 20-80%. The effect of metallicity is shown to be important, as the most luminous supergiants are only predicted at low-metallicity, thus limiting $-$ but not excluding $-$ the extent of a polluted second-generation at high-metallicity. These massive stars becoming black-holes suggests globular-clusters hosting gravitational-wave progenitors. Our scenario predicts a correlation between the mass of the cluster and the extent of the multiple population phenomenon.


INTRODUCTION
Young massive clusters (YMCs) are compact star forming regions with a radius of only a few parsecs (Portegies Zwart et al. 2010;Longmore et al. 2014). Since their projected lifetimes are consistent with those of old globular clusters (GCs, Maíz-Apellániz 2002), they have been suggested to become GC-like objects eventually. In turn, old GCs, observed to populate the bulges and halos of many galaxies including our own, are hypothised to start out as massive clusters (Brodie & Strader 2006;Andersen et al. 2016).
Both YMCs and GCs, as well as their suggested connection, are surrounded by observational puzzles. For example, why do we see multiple stellar populations in practically all GCs (e.g. Yong et al. 2003;Gratton et al. 2004;Harris 2010;Da Costa et al. 2013;Bastian & Lardo 2018), and possibly in other clusters with ages up to 2 Gyr (e.g. Martocchia et al. 2018a)? Since one of the main indications that a cluster harbors multiple populations, is the anomalous ratios of light elements-e.g. the observed ratio of sodium and oxy-gen, which can only be synthesized at temperatures as high as 60−100 MK-it has long been suggested that a first generation of massive or intermediate mass stars is responsible for the formation of an anomalous second generation. But responsible in which sense? What are the conditions under which a second star formation episode can happen that feeds on the material ejected from the first generation? Or, to turn the question around, is the amount of material ejected from the first generation stars enough to produce the observed number of second generation stars? This latter puzzle is usually referred to as the 'mass budget problem', since most scenarios suggested so far do struggle to answer yes.
It all may have something to do with metallicity, as we know that massive stellar evolution strongly depends on this (e.g. Meynet & Maeder 2002;Yoon et al. 2006;Brott et al. 2011a;Georgy et al. 2013;Sanyal et al. 2017;Vink 2018). But how does the metallicity of the cluster influence the second generation star formation? In particular, how do the hydrodynamic structure of the cluster gas depend on the metal-licity of the first stellar generation? How does the composition of the second generation depend on that?
Recently, Szécsi et al. (2018) suggested that cool supergiants may play a role in the formation of the multiple populations in GCs. They investigated a scenario in which the second generation forms in a photoionization-confined shell around such cool supergiants, and speculated that the mass budget problem may be solved by only forming low-mass stars-which, contrarily to some of the other scenarios that assumed the same (e.g. de Mink et al. 2009b;D'Ercole et al. 2010), is better justified in these exotically shaped star forming regions (cf. Sect. 4.7 of Szécsi et al. 2018). Nevertheless, they also suggested that even without such a shell, the wind of supergiant stars may play an important role in GCs that shall be more closely inspected-and this is the objective of present work.
Here we investigate these questions by combining up-todate theories of massive stellar evolution (those that predict cool supergiants) with calculations of the cluster's hydrodynamic structure. Earlier studies involving massive or intermedate mass stellar evolution were able to predict e.g. chemical pollution and element ratios (e.g. Karakas et al. 2006;Decressin et al. 2007b;de Mink et al. 2009b;Denissenkov & Hartwick 2014;Szécsi et al. 2018), but not able to say too much about the hydrodynamic behaviour of the gas in the cluster or, for that matter, under which conditions the formation of the second generation of stars happens. On the other hand, studies of the gas reinserted by massive stars within young clusters and its eventual accumulation leading to secondary star formation (e.g. Silich et al. 2004;Tenorio-Tagle et al. 2005;Wünsch et al. 2011Wünsch et al. , 2017Palouš et al. 2014;Martínez-González et al. 2016;Silich & Tenorio-Tagle 2017) have hardly ever taken into account newly found peculiar aspects of stellar evolution. Combining these two research areas is therefore a viable and auspicious approach.
To that end, we use a hydrodynamical, semi-analytic calculations of the cluster structure (taken from Wünsch et al. 2017) which accounts for the winds of the first stellar generation as an input. We apply two different sets of single stellar evolutionary models for this first generation (taken from Brott et al. 2011a;Köhler et al. 2015;Szécsi et al. 2015). They correspond to metallicities of the Large Magellanic Cloud (LMC) and of the low-metallicity dwarf galaxy I Zwicky 18 (I Zw 18). Thus we are able to investigate both the questions when? and with which composition? the second generation of stars may form in a young massive cluster. Additionally, we are able to study the process' dependence on metallicity, as well as the role that cool supergiants play in it.
This paper is organized as follows. The semi-analytic hydrodynamic code which determines the cluster structure, as well as the stellar evolutionary models, are described in Sect. 2. The synthetic population of stars that we created from the massive stellar models is discussed in Sect. 3. In Sect. 4 we perform calculations of the cluster's hydrodynamic structure applying the synthetic populations. We also discuss the conditions under which a second generation of stars may form. Sect. 5 then investigates the chemical composition of this second generation. Sect. 6 deals with the mass budget. Sect. 7 discusses caveats and future directions, while Sect. 8 summarizes and concludes the work.

Rapidly cooling shocked winds
Young massive clusters include large populations of massive stars concentrated in a rather small volume. Thus their stellar winds are expected to collide with each other and heat up to high (∼ 10 6 − 10 7 K) temperatures. The overpressure of this hot gas then drives a star cluster wind (Chevalier & Clegg 1985). If the cluster is massive and compact, and hence the density of the hot gas within it is high, the gas becomes thermally unstable, cools down to ∼ 10 4 K and forms dense clumps. The clumps are initially warm and ionized due to the radiation of nearby massive stars, but a fraction of them falls into the cluster centre due to the cluster gravity, where the gas accumulates until its column density is high enough to self-shield against the ionising radiation. Then, the gas cools further to lower temperatures and forms new stars. This scenario was explored extensively in a series of papers by Silich et al. (2003Silich et al. ( , 2004; Tenorio-Tagle et al. (2007); Wünsch et al. (2008Wünsch et al. ( , 2011; Palouš et al. (2013Palouš et al. ( , 2014 and others. Wünsch et al. (2017) studied this model of rapidly cooling shocked stellar winds by means of 3D hydrodynamic simulations including gravity (of both stars and gas), radiative cooling of the hot gas and EUV radiation of massive stars. They estimated a fraction of stellar winds that accumulates inside the cluster depending on various cluster parameters.
They compared the results of the 3D simulations to the outcome of a much simpler and much faster 1D semi-analytic code (see below) which is also able to estimate the mass of the second stellar generation. They found a good agreement. As for the first stellar generation, they relied on the predictions of the stellar synthesis code Starburst99 by Leitherer et al. (1999). They found that, using solar metallicity Star-burst99 models of single stars up to M top = 120 M and a standard initial mass function, a substantially massive second generation of stars form only if the heating efficiency 1 , an observationally poorly constrained parameter, is very low. Here we apply the updated version of the semi-analytic code using stellar populations with different underlying physics (includ-1 Heating efficiency means the fraction of the mechanical energy of stellar winds that is transformed into thermal energy of the hot shocked gas inside the cluster. ing different metallicities and M top ). Below we shortly describe how the code works and what initial parameters we assume when running it.
In a gravitationally bound star cluster (SC), it can safely be supposed that properties of the stellar winds vary on a much longer timescale than the cluster wind crossing time. Therefore, one can search for a stationary solution of a set of spherically symmetric (1D) hydrodynamic equations to describe the star cluster wind (Chevalier & Clegg 1985). If the cluster is massive and compact enough, the hot gas is subject to radiative cooling due to its high density, and thus the set of stationary hydrodynamic equations to be solved should include the appropriate cooling term. A code to solve such a set with the assumption of spherical symmetry was developed first by Silich et al. (2004). Here we use a similar code described in (Wünsch et al. 2011(Wünsch et al. , 2017 updated with terms describing the effect of the star cluster gravity on the gas.
The set of stationary spherically symmetric hydrodynamic equations has a form where ρ, u and P are the wind density, velocity and pressure, respectively. The mass and energy input rate densities, q m and q e respectively, represent stellar winds approximated by a spatially smooth source described with a generalised Schuster distribution (∼ [1 + (r/R c ) 2 ] −β for r < R SC ) with parameters R c , R SC and β being the core radius, cut-off radius and the slope of the distribution, respectively. The gravitational potential Ψ includes only a contribution from stars (i.e. the gas gravity is ignored) and it is assumed that the mass is distributed with the same generalised Schuster distribution as the wind sources and that the total mass is M SM . The cooling term has a form Q = n i n e Λ(T, a j ) where n i = n e = ρ/µ H are the ion and electron number densities and Λ(T, a j ) is a cooling function calculated by Schure et al. (2009)  (1)-(3) exists only in a subset of the parameter space. It is possible to define a critical luminosity L crit so that the solution exists only if the total mechanical luminosity of all stellar winds is smaller than the critical value, L SC < L crit . The semi-analytic code determines the critical luminosity iteratively by trying to solve Eqs. (1)-(3) for a given set of parameters: R c , R SC , β and properties of the wind of a stellar population with a given mass, age and chemical composition. It has been shown by Wünsch et al. (2017) that L crit can also be used to estimate the rate of clump formation: the mass of clumps formed over time is the difference betweenṀ SC andṀ crit , the former being the mass deposition rate of the cluster, and the latter the corresponding quantity but taking L SC = L crit . The mass accumulation rate,Ṁ acc is then determined by taking into account only clumps which are formed with initial velocity (the same as the cluster wind velocity u) smaller than the escape velocity u esc ≡ √ 2Ψ As the semi-analytic code models the star cluster using a smooth distribution of mass and energy sources, the calculations cannot represent discrete events such as supernova explosions nor, therefore, the effect of the dust they produce. We discuss why and when it is justified to omit supernovae in our calculations in Sect. 3.5.
In this work, we discuss two cluster models differing in the initial chemical composition of massive stars (cf. Sect. 2.2). We assume that all first generation stars are formed abruptly at t = 0 and we follow the cluster evolution for 10 Myr. The total mass of the first generation stars is M SC = 10 7 M for both models, the stars are assumed to form with the standard Initial Mass Function (IMF, Kroupa 2001), and they are represented by the stellar models described in Sect. 2.2. The second generation stars are represented only as the accumulated mass (M acc ), they do not contribute to our calculations by e.g. their stellar winds. The stellar density profile of the cluster is given by the generalised Schuster distribution with R c = 1 pc, R SC = 3 pc, and β = 1.5. As opposed to Wünsch et al. (2017), here we assume that the all mechanical energy of stellar winds is converted to the thermal energy of the hot gas (i.e. the heating efficiency is unity), and that the mass loading is zero. On the other hand, we use different upper mass limit, M top , of our first generation stellar population, as explained in Sect. 3.3. Additionally, we carry out a parameter space study (details given in Sect. 6.2), where we vary the initial cluster mass, M SC , and the index of the IMF for stars more massive than 1 M .

Stellar evolutionary models
To account for the first generation of massive stars, we apply two sets of models, both computed with the BEC code (see e.g. Heger et al. 2000;Yoon et al. 2012;Szécsi et al. 2015, and references therein). The models with initial composition of the LMC were created by Brott et al. (2011b) and Köhler et al. (2015), representing a subsolar-metallicity environment with ∼0.4 Z (i.e. [Fe/H] ∼ −0.4). Those with initial composition of the dwarf galaxy I Zw 18 were created by Szécsi et al. (2015) and Szécsi et al. (2018) The low-metallicity models between M ini = 10-300 M have all been followed until the end of core-helium-burning (completing on the work of Szécsi et al. 2015, who only followed them until the end of core-hydrogen-burning). From this point on, we mainly refer to the composition of the LMC as high-metallicity and to that of I Zw 18 as low-metallicity.
The most massive models (>70 M in the high-metallicity set and >300 M in the low-metallicity set) have been, however, only computed until core-hydrogen-exhaustion (i.e. terminal-age main-sequence). Therefore, we extrapolate for how much mass they would lose during their remaining evolution if the mass-loss rate was the same as that at the end of the computation. This is clearly a simplistic approach that brings some additional uncertainty into our predictions.
As for rotation, all the models have zero or slow rotation (i.e. 0 or 100 km s −1 initially). They evolve with a distinct core-envelope structure towards lower surface temperatures. It shall be a future task to add models that have more extreme rotation rates (and, for example, evolve chemically-homogeneously). Also, effects of binarity are omitted at this point as we only apply single stellar models.
The wind velocity, v wind , of any given stellar model is calculated from the escape velocity from the stellar surface as v wind = 1.3·v esc and v wind = 2.6·v esc for models below and above a surface temperature of 21 kK, respectively, following the theory of line driven winds (Lamers & Cassinelli 1999). Additionally, following e.g. Leitherer et al. (1992), the wind velocity is corrected for the metallicity of the wind material, Z, by multiplying it by a factor (Z/Z ) 0.13 . Since supergiants' winds are not expected to be line driven, we checked that the outcome of our calculations is not, in fact, sensitive to the actual values of supergiant wind velocity we use, as long as they are below 80 km s −1 , which they indeed are.

STELLAR POPULATIONS
3.1. Comparing the two sets of models Fig. 1 shows the Hertzsprung-Russell diagrams of the two sets of models. The most important difference between the high-metallicity models (LMC) and the low-metallicity mod-els (I Zw 18) is the presence of very massive ( 100 M ) cool supergiants in the latter case, populating the top-right corner of the diagram. The reason the very massive, high-metallicity models do not evolve to the supergiant branch is that their mass-loss rate during the main-sequence is high enough to remove almost the whole envelope, turning them into hot stars such as Luminous Blue Variables or Wolf-Rayet stars. As we show in Sect. 4.1, the presence of very massive supergiants at low metallicity is the key to form a second generation of stars early enough so that they show chemical abundances attributed to a subsequent population.
Other differences between the high-and low-metallicity sets of models are: • both the zero-age main-sequence and the terminal-age main-sequence lie at lower surface temperatures when the metallicity gets higher; • blue supergiants are found in the low-metallicity models with 9-30 M during core-helium-burning (i.e. blue loop); • red supergiants with 9-45 M are found amongst the high-metallicity models, these are also core-heliumburning objects; • the presence of Luminous Blue Variables at M ini ∼ 70-100 M in the high-metallicity models; Luminous Blue Variables are found in the low-metallicity models only at masses above 300 M ; • the presence of Wolf-Rayet stars at M ini 100 M in the high-metallicity models; no Wolf-Rayet stars are found in the low-metallicity models.
Since the very massive supergiants at low-metallicity are the key objects responsible for the multiple population phenomenon in our model, it is important to discuss them here a bit further.

Cool supergiants
Low-metallicity models with M ini 80 M evolve to the supergiant branch even during their core-hydrogen-burning phase due to envelope inflation. This phenomenon has been investigated by Sanyal et al. (2015Sanyal et al. ( , 2017 who found that the reason these stars inflate their envelopes and thus expand is their proximity to the Eddington-limit. What is extremely intriguing in the context of the chemical composition and the multiple populations in star clusters, is that these massive stars becoming cool supergiants during their core-hydrogen burning phase means they have a convective envelope. Convection mixes the material between the core-where the CNO-cycle operates, together with side reactions that can synthesise Na and Al at the expense of O and Mg-and the surface. The surface layers are then removed by the stellar wind, thus polluting the interstellar gas with nuclear ashes that have undergone hot-hydrogen-burning. Interestingly, the convective core does not reach down again into the burning regions during core-helium burning (since these layers are much deeper inside than those of core-hydrogen burning), thus avoiding the ejection of helium-burning products. This makes these supergiants quite ideal to be suggested as potential pollution sources in GCs.
If these stars exist in nature, is a question for future investigations. Envelope inflation has been recently studied by several authors (Gräfener et al. 2012;Grassitelli et al. 2015a,b;Sanyal et al. 2015Sanyal et al. , 2017. In particular, Moriya & Langer (2015) suggested that such supergiants may be responsible for some supernovae of the superluminous type. An additional caveat of using simulations of supergiants is that their mass loss rates are, despite great efforts to constrain them, still quite uncertain. For a comprehensive discussion on the subject, we refer to the book of Levesque (2017) as well as to the relevant literature on our stellar models (Sect. 5 in Szécsi et al. 2015 andSect. 2.1 in Szécsi et al. 2018).

Population synthesis
To create a synthetic population, we suppose that all the first generation stars of the cluster were formed during a single starburst episode, almost instantaneously. We assume a standard piece-wise power-law IMF with three intervals (0.01 − 0.08, 0.08 − 0.5 and above 0.5 M ) with indeces −0.3, −1.3 and −2.3, respectively, suggested by Kroupa (2001). In Sect. 6.2 we explore how the results change if a top-heavy IMF is used, by introducing an additional fourth interval for masses above 1 M with index, α 4 , varying between −2.3 and −1.1.
Since both sets of stellar evolutionary models only contain ten models between 9-500 M , we need to interpolate between these tracks to get two, smoothly changing grids (cf. dashed lines in Fig. 1). We split the whole range of stellar masses (0.1-500 M ) 2 equidistantly in the logarithmic scale into 5000 intervals. For intervals above 9 M (there is ∼1800 of them) we apply the interpolated tracks, while stars less massive than that are taken into account only as mass holders. When interpolating between the stellar models, we use logarithmic scale to interpolate linearly between age, mass, mass-loss rate, surface temperature, radius and surface abundances; while bolometric luminosity, escape velocity and wind energy rate are calculated as follows: L bol = σ SB · T 4 e f f · 4π · R 2 ; v 2 esc = 2G · M/R and L = 0.5 ·Ṁ · v 2 esc . The tracks of this smooth grid are weighted with the IMF and interpolated for a selected age. Then sums of various quantities, such as mass and energy deposition rates, and mass weighted mean abundances, are calculated. This calculation is carried out for 500 points distributed uniformly throughout the followed period of the cluster evolution (10 Myr). The resulting time evolutions of these quantities is then used as the input for the semi-analytic code (described in Sect. 2.1) which computes the structure of the star cluster wind and eventually the rate of the mass accumulation and its chemical composition for each point in time.
We sampled the stellar evolution models so that major steps occur more or less at the same evolutionary stage; in particular, since in our cluster calculations the most important property is the mass-loss rate, we made sure that the interpolated mass-loss rate behaves well around the bi-stability jump (i.e. at T e f f ∼ 21 kK where the mass loss rates change abruptly due to an increase in the line acceleration of Fe III below the sonic point of the stellar wind, cf. Vink et al. 1999;Vink et al. 2000).
As opposed to earlier works on the same subject (such as Wünsch et al. 2017) who used M top = 120 M , here we include very massive stars up to M top = 500 M . This choice is motivated by the findings that low-metallicity stellar models between 150-500 M display significant variations in their surface abundances of light elements ). There are observational implications for the existence of stars up to 315 M in the LMC (Crowther et al. 2010(Crowther et al. , 2016. At low-metallicity, radiation driven winds are less effective, so ideally even more massive stars than that may form. Figure 2 presents the properties of the stellar wind such as its mass loss rate, mechanical luminosity, as well as its velocity for all model stars in the population. The mass loss rate typically becomes higher with stellar mass. While slowly increasing in the first half of the stars' life, it then experiences a local minimum and then a sudden jump. This jump is attributed to certain changes in the wind structure at T eff ∼ 21 kK leading to an increased mass loss. (It is called the bi-stability jump, and concerns the fact that additional iron line transitions become effective in driving the wind under this temperature; see Vink et al. 1999). Soon after this, the stars reach their post-main-sequence phase, during which the mass loss rates are 1-2 orders of magnitude higher than during the main-sequence phase. The models' evolution is computed until core-hydrogen exhaustion, but subsequent evolutionary phases are very short (< 1%) compared to the total lifetimes, so they can be safely omitted from considerations of stellar wind mass loss.

Wind properties of the populations
Models with low-metallicity have typically lower mass loss rates than their high-metallicity counterparts of the same mass. The reason is that the mass loss rate is not only the function of stellar mass, but also of metallicity (and, at some extent, of other stellar paratmeters such as radius and surface composition). The mass loss rate's dependence of metallicity is prescribed asṀ ∼ Z 0.86 in both set of models at every evolutionary phase, while its dependence on the actual stellar mass followsṀ ∼ M 1.13 .
There is one evolutionary phase when the mass loss experienced by low-metallicity models is higher than that experienced by the high-metallicity ones: in the case of the most massive models' late evolution (i.e. above 100 M ). Here the mass loss of the high-metallicity models follow a prescription typical for Wolf-Rayet stars, while that of the low-metallicity ones follow a prescription typical for red supergiants. While the latter does not predict significanty higher mass losses than the former, we have to take into account another effect also playing a role. Namely, that since the mass loss rates during the first half of the main sequence are lower in the case of low-metallicity models, these stars are typically more massive during their later phases than those predicted by highmetallicity models of the same initial mass, and therefore their mass loss rates are now higher. The jumpy mass loss rates visible for example during the late phases of the highmetallicity model with 70 M or the low-metallicity model with 575 M , are due to these models being associated with an LBV phase. Models below 30 M at low-metallicity experience a blue loop during their core-helium-burning phase. The mass loss attributed to this blue supergiant phase is typically lower than what is expected for a red supergiant, leading to a plateau in these models' post-main-sequence mass loss rates.
The behaviour of wind velocity of the models in Fig. 2 can be understood as follows. Initially, the low-metallicity models are hotter than the high-metallicity ones due to their surface opacity being lower. Therefore, these models are typically smaller in radial size. As wind velocity is computed from the escape velocity, it is expected that their wind velocity is higher than their high-metallicity counterparts'. As the evolution progresses however, the low-metallicity models evolve towards lower surface temperatures and larger radii, thus their wind velocity drops. The same is happening with the high-metallicity models. One crucial difference is that while the low-metallicity models become supergiants, the high-metallicity ones, at least those above 70 M , become LBVs or WR stars. These models' wind velocity is high. As for masses below 70 M , both grids predict supergiantsthe low-metallicity grid blue supergiants below 40 M with a loop in the wind velocity.
Energy flux of the wind is closely related to both the mass loss and the wind velocity, as it is computed by these two as L wind = 1/2Ṁ v 2 wind . From this we can calculate the mechanical luminosity inserted into the cluster by the first generation of stars, by properly considering every individual stellar model's contribution to the population (i.e. weighting with the IMF). Thus we arrive to the value L SC presented in Fig. 3.

Supernova explosions
Massive stars end their lives in various ways, depending on mass and metallicity. If supernova explosions happen, they may contribute to the cluster's subsequent evolution significantly. Core-collapse supernova explosions in particular may increase the cluster's iron content. For this to happen however, a mechanism is needed that traps the (possibly very energetic) supernova ejecta inside the cluster's potential well. Such mechanism was suggested e.g. by Tenorio-Tagle et al. (2013) and further studied by Martínez-González et al. (2018), however, its discussion is out of the scope of this work.
Our semi-analytic model does not include this effect. Therefore, we need to discuss when and what kind of supernova explosions we expect from our massive stellar models (if any) so that we can carefully evaluate the validity of our semi-analytical approach. The majority of globular clusters show no variation in their iron content, so to account for them with our model, we need to play extra attention to the time periods when supernovae's contribution enter the picture. Some globular clusters such as e.g. ω Cen and M 54 are peculiar in this regard, displaying variations in iron and even in their total sum of C, N and O content. Still, their formation may have happened fundamentally differently from the average cluster (ω Cen in particular was suggested to start out as a dwarf galaxy, cf. Schiavon et al. 2017). Therefore we only account for the majority of clusters here, those that do not exhibit star-to-star variations in iron abundance.
Connecting stellar models to supernova types is somewhat uncertain. It depends on both the model in question, and the assumptions about the nature of the explosion. Here we simply rely on the work of Heger et al. (2003) and establish that, in the case of our low-metallicity population, the first supernova that may pollute the cluster with iron (which is a 40 M star), explodes at a cluster age of 4.5 Myr. The reasons are the following.
In the low-metallicity population, massive stars up to an initial mass of 25 M are expected to explode as corecollapse supernovae of type II P. The total lifetime of a 25 M star is 7 Myr; those at lower masses live even longer. We expect them to form neutron stars as remnants. As for stars between an initial mass of 25−40 M , they are expected to explode as weak II P supernovae and form black holes as remnants. The total lifetime of a 40 M star is about 4.5 Myr.
Above this mass but below 140 M , it is expected that the metal-poor models do not explode but fall into a black hole directly due to their immense self-gravity at the moment of iron-core collapse. An initially 140 M model has a total lifetime of 2.4 Myr. This means that stars with lifetimes between 2.4 and 4.5 Myr do not explode as supernovae; athough they contribute to the (stellar-mass) black hole content of their clusters. The same fate awaits those stars that have an initial mass above 260 M -that means many of our metal-poor models form black holes without an explosion.
Between 140−260 M (i.e. total lifetimes of 1.9−2.4 Myr), the models in the low-metallicity set are again predicted to explode. This time though, it is not due to iron-core collapse but another effect: pair-creation 3 . During their oxygenburning phase, the creation of electron−positron pairs disturbs their hydrostatic stability, and makes them explode without leaving a remnant. Such a pair-instability supernova does not pollute the cluster with iron, as nuclear fusion has not yet produced an iron-core. The core that explodes con-tains mainly carbon and oxygen, which undergo explosive burning during the supernova event, producing a unique nucleosynthetic signature (Burbidge et al. 1957;Langer 1991;Heger et al. 2003;Langer et al. 2007;Kozyreva et al. 2014).
A pair-instability supernova is expected only for a small number of all massive stars in our low-metallicity population (for a 10 7 M cluster with the standard IMF, we expect about 1 such supernova in every 270 years). For our present purposes, we do not investigate how these supernovae may contribute to our cluster's structure or chemical composition, but consider all stars above an initial mass of 40 M (that is, an age of 4.5 Myr) not to pollute the cluster with iron. We also do not investigate the effect of pulsational pair-instability (Woosley et al. 2007;Moriya & Langer 2015), pointing out that this process may also play some, yet to be investigated, role in polluting the cluster.
As for the high-metallicity set of models, the situation is quite different. Here we also expect type II P supernovae from models between 10−25 M initial mass, with stars between 25−40 M ending up as a type II L/b supernovae. But above that limit, instead of falling into black hole due to self-gravity or undergoing pair-creation induced instability, stars are expected to explode as supernovae of type I b/cthat is, due to iron-core collapse. The reason for the highmetallicity models exploding as core-collapse supernovae instead of some more exotic scenario, is that they lose so much mass during their lifetimes that their final mass is in the range where neither self-gravity is strong for direct black hole formation, nor is pair-creation playing a role. So, all our highmetallicity stellar models undergo a core-collapse induced supernova explosion that may, if the ejecta is trapped in the cluster, pollute the gas with iron. The first supernova, that of our most massive model (initial mass 500 M but final mass only 34 M ) explodes at the age of 2.1 Myr.
When discussing our results in the next sessions, we always point out where and when supernovae are expected. The extent of which the supernova ejecta stays trapped to mix with the gas, remains a question. It may as well depend on the nature of the explosions themselves, if they are weak ('failed') or strong ('successful') supernovae, a question currently undergoing some investigation (MacFadyen & Woosley 1999;O'Connor & Ott 2011;Smartt 2015). Also, it remains a question how the supernova feedback influences star formation of the second generation stars. It may enhance it or stop it; with our current method, we have no ways to know. So, all these questions around supernovae are left to be investigated in future work.

CLUSTER WIND AND SECONDARY STAR FORMATION
As explained in Sect. 2.1, mass accumulation happens when the cluster luminosity L SC is lower than a certain criti-   Figure 3. Time evolution of the cluster wind. Whenever the cluster wind luminosity LSC exceeds the critical luminosity Lcrit, wind mass is supposed to accumulate into the cluster center (marked by the shaded regions). Mass loss rate of the stellar populations is also shown (blue lines, values on the rigth axis). In the I Zw 18 population, mass accumulation happens early on: before the first supernova explosions happen at 4 Myr. Thus we conclude that this early mass accumulation episode provides a 'window' for the undisturbed formation of a second generation of stars.  Figure 4. Time evolution of the mass lost from massive stars (i.e. inserted into the cluster wind by stellar winds, Mins and the mass which is accumulated into the cluster center (Macc), both on a cumulative scale (left axis). Also shown is the helium mass fraction (right axis). The periods when mass is actually accumulated (cf. the shaded regions in Fig. 3) are indicated with thick lines, the four colors (black, blue, brown, white) corresponding to every 25% of the total accumulated mass.
cal luminosity L crit . The former is simply the combined mechanical energy of all stellar winds, while the latter is determined by the hydrostatic structure of the cluster (computed in our semi-analytical model) and takes into account various effect such as cooling and the stability of the cluster wind. Mass accumulation means that mass is removed from the wind and supposed to accrete onto forming proto-stars of a secondary generation-given of course that its velocity, inherited by the gas' velocity, is lower than the escape velocity from the cluster's gravitational potential. Figure 3 shows how the cluster luminosity L SC relates to the critical luminosity L crit over the calculations. The phases when L SC exceeds L crit are are shaded in the figure; however, the size of the shaded area should not be taken as indicative of anything, especially not the efficiency of star formation; it only serves to help us see when and how long the star formation is going on.

Mass accumulation at low-metallicity
In the case of our low-metallicity stellar population, mass accumulation starts at ∼1.6 Myr and lasts until ∼3.9 Myr, resulting in a star formation episode that is ∼2.3 Myr long. Figure 3 also shows the total mass loss rate of the clusterthat is, the mass loss rate of all stars in the population combined. The mass accumulation episode coincides with a pronounced peak in the mass loss rate; but this is not really a coincidence, as both effects are caused by the presence of cool supergiants in the population. Indeed, as discussed in Sect. 3.1, massive and very massive supergiants are present in the low-metallicity population (but are absent in the highmetallicity population). These stars have low surface temperature (thus a slow stellar wind) and a high mass loss rate. Therefore, they facilitate star formation by not heating the gas too much since they eject material with a small velocity and thus keep the cluster wind velocity also rather low. Additionally they deposit a huge amount of mass into the cluster in just the right time for it to accumulate in the center.
As discussed in Sect. 3.5, the first core-collapse supernova explodes when the low-metallicity cluster is about 4.5 Myr old. Mass accumulation in this cluster happens before that age. Thus we conclude that in our low-metallicity calculations, the proto-stars of the second generation will have been already formed out of the stellar wind material before iron is deposited into the gas via supernova explosions.
Another caveat that this cluster avoids due to accumulating mass early, is that observationally, YMCs typically remove their gas by ∼4 Myr (Hollyhead et al. 2015). Figure 4 shows the mass that is lost in stellar winds (i.e. inserted in the cluster by stars, M ins ) as well as that accumulated in the center (M acc ). The mass starts to accumulate when L crit first exceeds L SC (i.e. at ∼1.6 Myr). Almost all the mass that is inserted into the cluster wind accumulates in the center. When the process ends (at ∼3.9 Myr), the total accumulated mass is almost 10 5 M . This is the mass budget from which a second generation of stars form.

Mass accumulation at high-metallicity
In the case of our high-metallicity stellar population in Fig. 3, mass accumulation starts at a later time than at lowmetallicity, at 4.2 Myr. After this time, L SC exceeds L crit and keeps exceeding it until the end of our calculation.
In this case however, we do not imply that a second generation of stars should be expected to form out of the accumulated mass. As mentioned in Sect. 3.5, this population experiences supernova explosions starting at the age of 2.1 Myr. Thus, our calculation should not be taken on face value after 2.1 Myr. We can nonetheless draw some interesting conclusions from it.
The reason this population does not produce a mass accumulation episode as early as the low-metallicity population, is that the very luminous supergiants are practically absent. In this population, stars above 40 M become LBVs or WR stars (as discussed in Sect. 3.1) which are hot stars with fast winds. So, while the mass deposited into the gas from the stellar winds is high due to the high mass loss rates of LBVs and WR stars (higher than in the low-metallicity population at any given point in time), the cluster wind luminosity does not exceed the critical value until the first red supergiants, those with ∼45-40 M and below, appear. Indeed, the to-tal lifetime of a 45 M model is 4.2 Myr, marking the point where L crit > L SC and the mass accumulation starts.
That is, if we disregard supernova explosions. It falls outside the scope of current work to investigate what happens to the supernova ejecta under the conditions in this cluster: whether it gets shocked and cools staying and mixing with the gas, or leaves the cluster; and whether it enhances or stops star formation. What we can conclude from our calculation nonetheless, is that mass accumulation starts ∼2.5 Myr later at high-metallicity (i.e. in the absence of cool supergiants) than at low-metallicity (when their contribution dominates). In the latter case, we expect that the accumulated mass forms a second generation of stars; while in the former case, we cannot be certain if a second generation forms without conducting 3D hydrodynamic simulations of the supernova ejecta and its contribution to star formation in the cluster.

Light element variations
Light element abundance variations are a well-established observational fact for practically all GCs which have been extensively studied spectroscopically (see e.g. Bastian & Lardo 2018, for a recent review). In particular, Na overabundance is always observed together with O depletion; while the sum of C, N and O is constant, suggesting that the CNO-cycle is operating. In some GCs, Al and Mg display variations as well. This implies that at least one population of low-mass stars (usually referred to as the second generation) is made up of material that has previously undergone hot-hydrogen burning. This process is known to be active only in massive and intermediate mass stars because low-mass stars' core temperatures are not high enough for that.
We do not study the contribution of intermediate mass stars (i.e. asymptotic giant branch, AGB, stars) to our cluster, since their contribution happens after ∼30 Myr (while our computation ends at 10 Myr). Nonetheless, AGB stars have been used to account for the observed light-element variations, see e.g. Cottrell & Da Costa (1981); Karakas et al. (2006) Here we focus only on massive stars, but point out that AGB stars may contribute to the cluster gas' composition at later ages. Figure 5 shows the observed anticorrelation between two pairs of light elements, O vs. Na and Mg vs. Al. The data is taken from Carretta et al. (2009) which is a FLAMES-UVES survey of ∼200 red giant stars in 17 GCs; as well as from Pancino et al. (2017, and private comm.) which contains ∼150 red giant stars in 9 GCs with Gaia/ESO-UVES abundance measurements of all four elements in question. Since the abundance scales used by these two surveys differ, we took into account a small (compared to the internal spreads)   . Observational error is typically between 0.05−0.12 dex. The composition of the pristine gas, that is, the composition we attribute to the first stellar generation in our calculation is marked by the filled circle (labelled '1st Gen.'). Our calculation predicts that the mass accumulated in the cluster center (cf. Fig. 3) have a composition shown by the thick line. The four colored stripes overplotted on top of the thick line (black, blue, brown and white) mark the four quadrants of the total mass (i.e. every 25%, starting with the black and ending with the white; cf. the same color coding in Fig. 4). Thin lines show the surface composition of the original stellar models during all their evolution; for the meaning of the colors, see Fig. 1. We expect that the accumulated mass will mix with the original gas in the center. Thus, if a star forms out of it (as we suppose is the case for the low-metallicity cluster, cf. Sect. 4.1) its composition will be a mixture of the pristine and the accumulated composition.
offset of ∼0.10−0.15 dex when plotting the two data sets next to each other, as suggested by Pancino et al. (2017, and private comm.). We account for low-metallicity and high-metallicity clusters in a simple way, by dividing the observational samples into two categories: one with [Fe/H] > −0.9 and another with < −0.9. The choice of this value is motivated by Fig. 3 of Harris (2010) in which the metallicity histogram of a large catalogue of GCs seems to show two peaks, with an arbitrary division at around −0.9. In future work this has to be refined by investigating a range of different metallicities, as discussed in Sect.7.2. According to Fig. 5, both anticorrelations are observed to be much more pronounced at low-metallicity than at high; with no significant Mg depletion found amongst high-metallicity clusters whatsoever in these data sets.
When comparing our theoretical predictions to the observed spreads, we suppose that our low-metallicity cluster (with [Fe/H] = −1.7) is representative for all GCs below −0.9 and the same for our high-metallicity cluster (with −0.4) for above −0.9. As the stellar models do not use an α-enhanced mixture (as suggested for GC stars by e.g. Decressin et al. 2007b, see their Table 3) but a mixture suitable for certain galaxies, therefore when creating Fig. 5, the initial O, Na, Mg, and Al abundances of our models are scaled to match the composition of the unpolluted red giants. Below we discuss what our calculations predict for the composition of the second generation of stars (or more precisely, as explained in Sect. 4.2, for that of the mass accumulated in the cluster center, out of which a second generation forms at lowmetallicity).

Na & O at low-metallicity
Our calculation of a low-metallicity cluster predicts that the mass accumulated in the center has a high sodium value and a large range of oxygen values. In fact, about half of the accumulated mass has extremely low oxygen abundance with high sodium (cf. black and blue stripes), while another half is spread out in oxygen. This is related to when the mass is accumulated: if it is accumulated early, its composition is dominated by the winds of the most massive supergiants (which evolve to the supergiant branch earlier); while if late, it is dominated by the less massive ones.
It is expected that in the center, this material mixes with the pristine gas (out of which the first generation of stars formed). Thus, the mixture of the accumulated mass and the original gas can possibly produce the whole observed range of Na-O abundances in stars of low-metallicity clusters.
Such dilution of polluted gas with pristine gas is typically also invoked by other scenarios (e.g. in Eq. 7 of Decressin et al. 2007a). One caveat here is that observationally, YMCs have removed their gas by ∼4 Myr (Hollyhead et al. 2015), which presents a challenge for all models that form the second generation after this age. But in our low-metallicity cluster governed by the presence of supergiants, this caveat is avoided by the 'window' for starformation opening between 1.6−3.9 Myr (as shown in Fig. 3).

Na & O at high-metallicity
Metal-rich clusters show a smaller extent of both Na and O variations than metal-poor ones, with the lowest observed Na value being about 0.4 dex higher. As explained by Carretta et al. (2009) this is because the plateau of minimum Na established by (a previous generation) of supernova nucleosynthesis is a function of the metallicity.
Although we cannot infer from our calculation of a highmetallicity cluster that a second generation of stars would form from the accumulated mass (due to the uncertainties associated with supernova explosions, cf. Sect. 3.5), it is nonetheless interesting to compare our predictions to observations of high-metallicity cluster stars in Fig. 5. As opposed to the low-metallicity case discussed above, now the most massive stars' mass loss have no contribution to the accumulated mass, as mass accumulation starts when these are already dead. The composition in this case is dominated by supergiants of initial mass below 40 M . They do produce some Na by destroying some O, but much less of an extent than higher mass stars do, especially when it comes to destroying oxygen. This is in accordance with the observational data in Fig. 5.

Al & Mg at low-metallicity
This is the case where our calculation struggles to account for the whole extent of observations. While red giants in lowmetallicity GCs display a broad spread in both Mg and Al, our prediction is that the second generation of stars would have an Al spread only about 0.6 dex with almost no Mg being destroyed. This is so even though our most massive models do lose material with a very low Mg abundance (cf. yellow and black thin lines corresponding to stellar models with 257 and 575 M ).
The reason for this lies in the specifics of the star formation episode. We accumulate mass in our calculation when the hydrodynamical conditions in the cluster are just right (cf. Sect. 4.1 and Fig. 3). This means that we have a star formation episode that lasts from 1.6 Myr to 3.9 Myr, during which the material lost by massive stars is integrated together (with properly weighting by the IMF) to produce the stripes in Fig. 5. Our very massive models gradually lose their outer layers via their stellar winds, starting with those layers that are less effected by hot-hydrogen burning. The first stars of the second generation form out of this material (black stripe). When these very massive stars are already losing their deeper, magnesium-deficient layers, lowermass stars have evolved to the supergiant branch, contributing to the total composition significantly. The material of their surface layers is, therefore, accumulated together with the deeper layers of the very massive stars, resulting in some slight decrease in the combined abundance (blue stripe) but clearly not enough to account for the whole observed spread in Mg. The later phases of star formation, when even lower mass supergiants dominate the composition, produces a second generation with an even higher Mg (yellow and white stripes).
Nonetheless, all the extremely low Mg values were observed in the same cluster, NGC 2808 (in both data sets). As discussed by Carretta et al. (2009), there is a quite significant cluster-to-cluster variation when it comes to Mg and Al (see their Fig. 6), with some clusters displaying large, and some displaying small, spread of these elements. Indeed, the phenomenon of magnesium depletion is not a common feature among GCs, rather an exception, with extremely low magnesium abundances only observed in a handful of clusters (e.g. NGC 2419, NGC 2808) and even there the situation is further complicated by the phenomenon's apparent dependence on mass and metallicity (cf. Pancino et al. 2017).
Our calculation is way to simplistic to account for this cluster-to-cluster variation, as we use only two sets of single stellar models at two given metallicity values, together with some-reasonable, but certainly improvable-assumptions about the secondary star formation. We discuss ways to improve our theory in Sect. 7.2. It is, for example, quite conceivable that some stars do form out of the pure material (that is, mixed neither with the other stars' ejecta nor with the pristine gas) of the most massive supergiants. Such a scenario was suggested by Szécsi et al. (2018) to possibly operate in the case of some very massive supergiants.

Al & Mg at high-metallicity
The high-metallicity models do not show the Mg-Al anticorrelation. Not even in the most massive case (up to 500 M ). The reason for this is related to their core temperatures. To destroy magnesium and produce aluminium, the 24 Mg(p,g) 25 Al chain should be active, which happens at a core temperature of ∼80−100 MK, according to Ventura et al. (2011). Below and above this temperature range, the reaction rate of the 24 Mg(p,g) 25 Al chain is too small for it to produce any effect in stellar models.
But having the correct core temperature is not enough. In order to destroy a lots of Mg, the right thermodynamical conditions should last for a long time (because the chain first creates 24 Mg → 25 Mg, and then slowly destroys 25 Mg too; see Fig. 1 of Ventura et al. 2011). So we need stars that not only pass through the correct temperature range while, for example, collapsing or re-structuring, but keep burning their fuel with exactly the right temperature for a long time. The longer the time the core has the right temperature, the more Mg is destroyed and converted into Al.
This differentiates between the LMC models and the I Zw 18 models. For example, our LMC model with M ini =260 M has T c < 55 MK during almost all its main sequence lifetime. It only reaches the range 80−100 MK when the core contracts to ignite helium, but this is a rather short phase in the star's life ( 10 kyr), after which the temperature increases way above 100 MK. The I Zw 18 model with 257 M , on the other hand, already starts its evolution with 60 MK and then slowly increases to 75 MK. Although the literature cites 80 MK as the nominal lower limit for the reaction to be effective, we find in our models that already at >65 MK there is significant Mg depletion and Al production if this temperature lasts for a long time (in our 257 M model, for ∼0.6 Myr). This is in accordance with observations (Fig. 5), which show that high-metallicity clusters have no significant variation in Mg.

C, N, O and He abundances
Observations of GCs typically show that carbon, nitrogen and oxygen abundances are in accordance with the CNOcycle's equilibrium values. It means that the sum of these three atoms is constant; they simply act as catalysts in the cycle. Nonetheless, C and O drops and N increases in later populations due to the CNO-equlibrium values being different from the abundances of the original gas. This is confirmed by our theoretical calculations. Our massive stellar populations (both at low-and at high-metalicity) do conserve the sum of C, N and O, with their respective abundances being consistent with the CNO-equilibrium values. However, as the available C and N data for the first generation stars are sparse and the interpretation of C and N variations is complicated by some evolutionary effects in the red giant branch phase (e.g. Boothroyd & Sackmann 1999;Gratton et al. 2000), we refrain from fitting the C & N anticorrelation here.
The helium mass fraction of the accumulated mass in our calculation is shown in Fig. 4. It reaches a much higher value than what is inferred from observations of any GC Bastian & Lardo 2018;Milone et al. 2017). We discuss the implications of this finding in Sect. 7.1. Colored stripes in Fig. 4 represent the same as in Fig. 5, that is, the four quadrants of the mass accumulated in the cluster center. Comparing these two figures, we find that the most extreme oxygen depletion is produced together with a helium mass fraction, Y, ranging between 0.5−0.7 (i.e. black and blue stripes). This composition is a mixture of the material ejected from our most massive stars down to 150 M . However, a less extreme oxygen depletion is possible to reach with a helium mass fraction between 0.3−0.5 (white stripe). This happens at end of the star formation episode when the last 25% of the mass is accumulated. This mass is made of the winds of supergiants with an initial mass 40−60 M . For further discussion on this issue, we refer to Sect. 7.1.
No helium-burning products, nor products of later burning phases, are found in our calculations at any metallicity. 6. MASS BUDGET 6.1. On the mass budget problem and dynamical removal of stars The fraction of stars with anomalous chemical composition, i.e. the second generation, varies in the range of ∼30−90% among GCs, with a mean value around 67 % (Milone et al. 2017;Bastian & Lardo 2018). A difficulty of most models to predict such a high fraction of second generation stars is called the mass budget problem.
In our calculation of a low-metallicity cluster with a total initial mass of 10 7 M , we find that about 10 5 M is available to form second generation stars. We emphasise that this is the mass that is ejected from massive stars via their winds and accumulated inside the cluster center-both processes accounted for in our calculations-and not, for example, the total mass in these massive stars. We took into account a standard IMF with an upper mass limit of 500 M . As opposed to earlier works accounting for the mass accumulation process (e.g. Wünsch et al. 2017), we did not include additional ad-hoc parameters such as mass-loading or heating 7 7 6.5 6.5 6.5 6.5 6.5 6.5  Figure 6. Results of our parameter space study. We vary the initial mass of the cluster, MSC, shown on the X axis, and the index of our IMF shown on the Y axis. Contours present the outcome of our calculation in terms of the amount of mass that is accumulated in the cluster center.
Since the amount of accumulated mass (as well as the number of second generation stars, cf. Fig. 7) is an increasing function of the initial mass, we suggest that this may explain why we observe a correlation between today's GC masses and the extent of a polluted second generation.
efficiency, but found that the presence of supergiants already facilitates the process. Globular clusters today have a typical total mass of a few times 10 5 M . If they indeed used to be massive clusters born with 10 7 M , they must have lost 90% of their mass during their lives. It has been suggested (D'Ercole et al. 2008;Decressin et al. 2010;Vesperini et al. 2010;Khalaj & Baumgardt 2015) that clusters may lose a huge fraction of their stars via dynamical 'evaporation'. This process mainly effects stars located in the outskirts of the cluster.
In our calculations, mass accumulation happens in the center of the cluster 4 . This means that even if mass is removed from the cluster very efficiently by dynamical evaporation, the ∼10 5 M proto-stars of the second generation would be very hard to eject via this process, due to them being centrally located. For example, N-body simulations of the cluster's long-term dynamical evolution carried out by Khalaj & Baumgardt (2015) show that it is indeed possible to explain present day observations by a cluster that contained a second generation number fraction of 10% initially, on the conditions that a substantial amount of gas is kept after the formation of the second generation, and that this gas is then expelled on a very short timescale. If such a process takes place in the cluster leading to the loss of almost only first generation stars, our low-metallicity model presented in Fig. 4 is able to fulfill the mass budget, having already accumulated and converted ∼10 5 M into second generation stars.
However, it is not conclusively established that the dynamical evolution leads to the loss of only first generation stars. A recent study (Reina-Campos et al. 2018) suggests that while GCs have indeed lost 90−95% of their initial masses, the present-day ratios of first vs. second generation reflect the initial values. Another argument comes from Kruijssen (2015) who suggests that if mostly first generation stars were lost due to tidal interactions with the host galaxy, we should expect to observe GCs with increasing mass-loss towards smaller galactocentric radii, with higher gas pressures at birth and with higher cluster metallicities (cf. Sect. 2.1.2 and especially argument (v) on page 1661 of the cited paper). While it could be insisted that the measurements of these quantities are significantly impacted by uncertainties, it is nonetheless clear that there are far too many open questions regarding the dynamical evolution during which YMCs become GCs for it to be called an established theory.
We have no means of solving any of these open questions here, due to us only focusing on the relatively short term phase of starformation. Indeed, our calculation only involves the first 10 Myr of the cluster's life, as we are mainly interested in the mass accumulation process. We predict that the mass is accumulated in the cluster center and the second generation of stars form there. Still, we cannot directly quantify which fraction of the first generation would be lost over the subsequent lifetime; nor, therefore, the final ratio of first vs. second generation we may expect in today's GCs after them having been undergone several gigayears of dynamical evolution.
We can nonetheless provide upper limits and predict some trends, by studying how our calculation is affected if we vary initial conditions such as the total mass of the cluster at birth or the IMF. This is done in what follows. We emphasize that from now on, we do not suppose that the dynamical ejection process prefers the first generation. If it does, it helps our case; but there are other options too to alleviate the mass budget problem.

Parameter space study
We repeated our low-metallicity calculation with varying two of the input parameters, the total initial mass of the cluster, M SC , and the high-mass index of the IMF, α 4 (our IMF is explained in Sect. 3.3). We vary M SC between 10 5 ..10 8 M and α 4 between −1.1..−2.3. The results are summarized in Fig. 6 (and Sect. 6.3) for the amount of accumulated mass and in Fig. 7 (and Sect. 6.4) for the number ratio of the second generation stars vs. the total, N 2 /(N 1 + N 2 ). Our motivation for testing the high-mass index of the IMF is that recently, some attention has been paid to measuring this value in various star forming regions with various methods. Some of these works report the finding of a top-heavy IMF (Kalari et al. 2018;Schneider et al. 2018a;Zhang et al. 2018), although the debate seems not to be over (Bastian et al. 2010;Khorrami et al. 2016;De Masi et al. 2018;Farr & Mandel 2018;Hopkins 2018;Schneider et al. 2018b).

Mass accumulation as a function of initial cluster mass and IMF
We find that the amount of accumulated mass is an increasing function of both the total initial mass and the number of massive stars in the population. Clusters in the top left corner of Fig. 6 do not accumulate a noteworthy amount of mass; in the case of a cluster with M SC = 10 5 M and a classical IMF index, no mass accumulation is taking place. This suggests that at least the lowest mass GCs with multiple populations must have started out as more massive clusters and then gotten rid of some of their mass over their lives.
In the most extreme, probably hypothetical case of a 10 8 M cluster with an extemely top-heavy IMF, as much as 10 7 M is available to form the second generation. The consequences of this are discussed further in Sect. 6.6.
As for the chemical composition, we find the following. The total spreads predicted in Na and O, as well as those in Mg and Al, together with the high helium abundance (Sect. 5.2), are practically unchanged when varying the initial parameters. Nonetheless, more massive clusters produce more polluted gas, though with a less extreme average composition. To be able to quantitatively compare these results to observations, we would need to model the mixing with primordial gas, which will be done in future work. Here we only establish a correlation between the initial cluster mass and the amount of (polluted) mass accumulating in the cluster center. What this correlation means in terms of number of stars, is discussed in the next sections.
6.4. The ratio of second vs. first generation stars Figure 7 shows the same parameter study as before, but here we include assumptions about the number of stars formed out of the available mass, as follows. The number of low mass stars in the first generation, N 1 , is taken to be stars with masses between 0.08 and 0.8 M , as these are expected to be still alive after ∼ 12 Gyr of cluster evolution, and thus to be observed today. To calculate N 1 then, we assume a first generation with mass M SC and a given IMF (depending on α 4 , with a mass range of 0.01 − 500 M ). The number of low-mass stars in the second generation, N 2 , is estimated from the amount of accumulated mass, M acc , by assuming that all this mass is used to form the second generation with a standard IMF. Similarly as before, we apply an IMF between 0.01 − 500 M , and take the number of stars in the mass range 0.08 − 0.8 M to be our N 2 . The ratio of generations plotted in Fig. 7 is then N 2 /(N 1 + N 2 ).
We find that, depending on the IMF index, the ratio of the generations can be anything between 1 and 80%. We emphasize again that we do not suppose a second generation forming stars only up to 0.8 M (as done in e.g. de Mink et al. 2009b), but apply a regular IMF which includes all stars, even massive and very massive stars. Additionally, the dynamical evolution of the star cluster should be taken into account: it may increase the N 2 /(N 1 + N 2 ) ratio even further, as the second generation is expected to be more centrally concentrated (see the discussion in Sect. 6.1). For observational evidence of a more centrally concentrated second generation, see e.g. Milone et al. (2012); Dalessandro et al. (2016).
Thus we conclude that either a top-heavy IMF or the dynamical removal of the non-centrally located stars-or a combination of these two effects-is our suggested explanation for the observations of the second generation being as populous as ∼30−90%.

On the correlation of second generation and cluster mass
It has been reported by several authors that the fraction of enriched stars in GCs strongly correlates with the cluster mass. See e.g. Figs. 14 and 16 of Carretta et al. (2010), Fig. 20 of Milone et al. (2017) and Fig. 7 of Carretta & Bragaglia (2018). To directly confirm these observations with our calculations, we would need to test the effects of two processes: the mixing with pristine gas (cf. Sect. 6.3) and the dynamical removal of stars (cf. Sect. 6.1). These tests fall outside the scope of the present study. However, we point out that both the total initial mass and the accumulated (i.e. enriched) mass do correlate with the number ratio of second generation stars in our calculations. Thus, unless the processes mentioned above cancel out this correlation for some reason (e.g. ejecting only centrally located stars; or removing more mass from the more massive clusters, thus making the initial cluster mass strongly anticorrelate with the present day Number ratio of 2G stars Figure 7. The same as Fig. 6, but for the number ratio of low-mass stars in the second generation (i.e. N2/(N1 + N1), as explained in Sect. 6.4). Note that we do not suppose that only low-mass stars form in the second generation, but apply a regular IMF with an index −2.3.
cluster mass), we predict this trend to be indeed imprinted on GCs observed today.
6.6. More than two generations?
In our model, we allow the second star formation episode to also create massive stars (see the discussion in Sect. 6.4). Thus, it is possible in principle that the whole scenario repeats another time. As seen in Fig. 6, our most extreme cluster with M SC = 10 8 M can accumulate as much as 10 6 − 10 7 M in the center-which may again form massive stars and thus repeat the cycle again. Without testing this possibility quantitatively here, we speculate that it may explain the fact that some clusters are observed to have more than two stellar generations (e.g. NGC 2808, cf. D' Antona et al. 2005;Milone et al. 2015). 7. DISCUSSION 7.1. On the helium problem Direct helium abundance measurements are difficult to obtain spectroscopically, since the relevant photospheric transitions need a much higher surface temperature than what red giant stars (those we have extensive spectroscopic studies of when it comes to light elements; cf. the dataset plotted in Fig. 5) have. Helium abundance is usually inferred indirectly from photometric measurements via fitting stellar models (isochrones) in the observed color-magnitude diagrams. These studies suggest that almost all GCs show variation in helium abundance, albeit rarely up to Y obs = 0.35.
In Sect. 5.2 we show that our models overpredict helium at an extent (Y ∼ 0.5−0.7) that is not reconcilable with observations. However, as pointed out by several authors already (Lochhaas & Thompson 2017;Bastian & Lardo 2018;Szécsi et al. 2018), this is a generic problem in the field. The current understanding of the nucleosynthetic origin of the sec-ond generation is that in order to reach the extremely low levels of e.g. oxygen (together with the other extreme values of light elements in Fig. 5), hot-hydrogen burning is needed because it has the relevant side reactions (namely, the Mg-Al chain and the Ne-Na chain) that produce the observed light element ratios. But of course by burning hydrogen, helium is created.
These problems indicate that our general understanding of the multiple population phenomenon is still far from perfect. For example, we may not have the complete picture of how the Mg-Al chain and the Ne-Na chain operate in massive and very massive stars. This would not be surprising, given the scarcity of observations when it comes to massive stars especially at low-metallicity (e.g. Sect. 1 in Kubátová et al. 2018).
Alternatively, there may be a way to separate Na, O, Mg and Al from He either inside the first generation of massive stars or in the second generation of low-mass stars (gravitational settling of elements may play some role?). That said, there may be a missing ingredient in our theory of the cluster gas dynamics (e.g. the mixing of the pristine gas to the ejecta of massive stars may happen in an unexpected way; turbulence in particular, cf. Hopkins 2013, may play a role?) or even that of the star formation process itself. Also, the method of measuring He variations indirectly from photometry and isochrone fitting, although currently our most reliable method to infer He dispersions (Cassisi et al. 2017), may in the future undergo some now not yet seen developments, leading to the revision of what we know about helium in GCs.
And finally, although some observational features are well explained by the hypothesis that a first generation of massive stars synthesised the observed element ratios, the problem around He may still mean that the massive star hypothesis needs to be looked beyond, and completely new hypotheses need to be suggested (see also the discussion in Sect. 2.2.2 of Bastian & Lardo 2018).

Directions for future work
We combined two research fields, massive stellar evolution and cluster gas dynamics; both have studied GCs so far mostly from their own perspectives. Combining them now opens up new pathways for investigation.
For example, we only applied single stellar models. But massive stars have a very high binary fraction (Sana et al. 2012). In a close binary, the interaction may lead to rather interesting outcomes such as the formation of gravitationalwave progenitors (e.g. de Mink et al. 2009a;Belczynski et al. 2016;de Mink & Mandel 2016;Marchant et al. 2016Marchant et al. , 2017Szécsi 2017a,b;Vigna-Gómez et al. 2018). In fact, interacting massive binaries have been suggested as the source of anomalous light element ratios by de Mink et al. (2009b) and investigated further by Bastian et al. (2013) and Elmegreen (2017). We note that our single stellar models have been computed with very similar physical assumptions as the binary system in de Mink et al. (2009b). Therefore we expect our stellar models to, upon putting them next to a companion in such a system, pollute the cluster at some extent via the binary channel as well. This however needs to be quantitatively investigated in the future, preferably by computing a set of detailed binary models. This would be a long-awaited completion of the work by de Mink et al. (2009b) who only computed one such system and extrapolated therefrom. Alternatively, although less elaborate, our set of single star models may be applied in a synthetic binary population, similarly to what we did here for a synthetic single star population, but with some necessary assumptions about the interaction of the companions. This binary population then can be applied as input for the cluster gas dynamics simulations, potentially allowing us to test uncertain parameters of stellar and binary evolution (e.g. mass transfer efficiency, nuclear reaction rates etc.).
Another way to make use of our new method, is to test the predictions of different pollution scenarios against each other. For example, the 'fast spinning star' scenario (Decressin et al. 2007b) can be tested against ours, or the interactive binary scenario against these etc., in terms of facilitating star formation and with which chemical composition. Alternatively, the role of rotation may be explored by using not only slow-rotating models of Szécsi et al. (2015) but moderate, and even fast rotators that predict chemically homogeneously evolution. These so-called 'TWUIN stars' (Kubátová et al. 2018) are needed for the star forming shell scenario of Szécsi et al. (2018) to operate. So if we want to explain out our present work's insufficient accounting for the high extent of Mg depletion by the shell-scenario (as done in Sect. 5.1.3), we need to make sure that the contribution of hot TWUIN stars do not brake down our 'window' of star formation by heating up the gas too much.
The role of metallicity should also be studied in more detail, as this may indeed be important especially when it comes to the most fragile elements (e.g. magnesium). Also, the very massive supergiants, those responsible for accommodating the secondary starformation, are only predicted in our lowmetallicity population and not in our high-metallicity one; but the exact metallicity value below which supergiants start to play their role in the formation of GCs, should be specified further in future work. We suggest to apply several sets of stellar models covering the range in metallicities between, and even above and below, of the two values investigated here.
Our low-metallicity supergiants are theoretical predictions; if they indeed exist, their observational discovery is yet to be carried out. They are of course not expected to be found in globular clusters, because even if they used to be there, they are long dead for now. As our models show, they are only expected to be present in clusters as young as between ∼2−4 Myr. Another caveat for their discovery is that their natal cloud needs to be sufficiently metal-poor and sufficiently massive to be able to form them at all. Nonetheless, from the models we know that they may be extremely bright objects. Szécsi et al. (2015, see their Sect. 5) explains that at the distance of 18 Mpc, they would appear with a visual magnitude of 19 mag. It has been suggested that brightness variations due to pulsations (with periods of the order of months to years, see also Moriya & Langer 2015) may reveal them as stars rather than star clusters in photometric multi-epoch observations.
When it comes to the calculations of the cluster gas dynamics, there are ways to improve here too. We suggest for example to apply our stellar models (both the low-and highmetallicity ones) as input for 3D radiation-hydrodynamic simulations. There are basically two ways to do this. One is to compute the so-called smooth source hydrodynamics of the cluster (following Wünsch et al. 2017), that is, to suppose a population of several hundred massive stars providing energy and mass which is inserted smoothly into the cluster. Another, more elaborate but also computationally more expensive way is to model the cluster evolution applying individual sources in the 3D simulations.
Less concerned with stellar evolution or cluster gas dynamics, ways to improve our general understanding of the multiple population phenomenon has been pointed out throughout the text. To summarize these: the effect of supernovae on our star formation episode should be assessed (Sect. 3.5), the conditions for dynamical mass removal including its effect on changing the ratio of first vs. second generation stars should be quantified (Sect. 6.1), and the extent of mixing the enriched material with the pristine gas in the cluster center should be investigated (Sect. 5.1.1). And last but not least, the issues around helium should be resolved (Sect. 7.1).
An interesting observational conundrum that arose recently, is that there seems to be a cut in the occurrence of the multiple population phenomenon at ∼2 Gyr. Some clusters with this age, like NGC 1978, do show multiple populations, whereas slightly younger clusters with comparable presentday masses do not (Mucciarelli et al. 2008(Mucciarelli et al. , 2014Martocchia et al. 2018a,b). If this proves to be an established fact in the future, any theoretical models should be able to account for it. It will need to be done, however, by including the investigation of the cluster's longterm dynamical evolution. As we have no means to do that yet (see also our discussion in Sect. 6.1), we leave this question open for now.

Gravitational-waves
With direct detections of gravitational waves (Abbott et al. 2016b,a;Bagoly et al. 2016;The LIGO Scientific Collaboration et al. 2017;Szécsi 2017a,b), many authors suggested globular clusters to be the host of these events (Rodriguez et al. 2015;Antonini et al. 2016;Belczynski et al. 2016;Askar et al. 2017). In Sect. 3.5, discussing the final fate and remnants of our supergiants, we pointed out that many of our stellar models are expected to form black holes or neutron stars after they explode; only those that explode as pair-instability supernova are not expected to leave any remnant. Thus, without quantitatively investigating this question, we point out that our scenario qualitatively predicts a significant number of these compact objects to be present in the young massive cluster, and thus to potentially merge over the long lifetimes of these clusters via dynamical interactions. Our work is therefore an important motivation to look for gravitational-wave emission, as well as its compact object progenitors, in globular clusters. The same is true for shortduration gamma-ray bursts, the origin of which have been associated with gravitational-wave emitting compact object mergers (Abbott et al. 2017).

CONCLUSIONS
We realized a novel synergy between two research fields, massive stellar theory and cluster gas dynamics. In particular, we explored whether the model of rapidly cooling shocked stellar winds combined with state-of-the-art stellar evolution models can contribute to the explanation of multiple stellar populations observed in globular clusters.
The model of rapidly cooling shocked stellar winds predicts that the hot gas within star clusters can become thermally unstable and form warm clumps. These clumps fall into the cluster centre where they cool further and form a second generation of stars. The new stars are necessarily enriched by the nuclear ashes synthesised in the first generation massive stars.
We apply stellar evolutionary models as input for the calculations of the cluster gas dynamics. The models are computed for two chemical compositions: for low-metallicity corresponding to [Fe/H] ∼ −1.7, and for a higher, but still subsolar metallicity corresponding to [Fe/H] ∼ −0.4. By applying these two sets of models, we evaluate the impact of metallicity on the secondary star formation.
We find that at low-metallicity, cool supergiant starspredicted to have very high mass loss rates and, in the same time, a low wind velocity-help to make the hot gas thermally unstable very early on. Their winds include products of hot hydrogen-burning, thus making them a suitable candidate for explaining the multiple population phenomenon.
Our calculations are run for the initial 10 Myr life of the clusters, and predict the amount of mass accumulated inside the cluster center, as well as its chemical composition depending on the cluster mass, slope of the initial mass function and metallicity. We draw the following conclusions: 1. A 'window' for undisturbed star formation. At low-metallicity, mass accummulation starts early (at ∼1.6 Myr), and a significant amount of mass is available for star formation before the first supernovae start to explode (at around 4.5 Myr) or before the gas is expelled from young massive clusters (typically observed around 4 Myr). This is thanks to the slow winds of massive supergiants. At high-metallicity however, mass accumulation starts later, after ∼4 Myr, thus limiting-but not necessarily excluding-our scenario to work.
2. Agreement with light-element abundance ratios.
Our calculations reproduce the Na-O spread sufficiently well. Also, only hydrogen-burning products are ejected (i.e. the sum of the C, N and O atoms are preserved, with C & O being depleted and N enhanced), but no products of helium-burning or those of subsequent burning stages. The spread our calculations predict in Mg-Al is lower than observed; although the stellar models do provide the right chemical composition (i.e. heavily depleted in Mg and enriched in Al), it is not captured by our 'window' of star formation. This, together with predicting higher than observed helium abundances, points to future directions of improvement.
3. Cluster center captures all the mass of stellar winds. Our metal-poor clusters with initial mass larger than several 10 6 M capture almost all the mass ejected by stellar winds. This 10 4 -10 5 M material accumulates in the center of the cluster and forms new stars there. Thus, under the assumption that the massive cluster is evolving to become a globular cluster by dynamically removing its not too centrally located (i.e. mainly first generation) stars over several gigayears, we are able to consistently explain the mass budget of present day globular clusters by applying a normal initial mass function for both the first, and the second, stellar generations.

4.
A top-heavy IMF helps though. With a normal initial mass function (α 4 = −2.3) our massive clusters form a second generation as populous as 1% of the total. But applying a top-heavy initial mass function with α 4 = −1.5 raises this fraction to close to 50%, and a very extreme index of α 4 = −1.1 up to 80%. This means that even without supposing dynamical removal of the first generation (see above), a top-heavy IMF can explain the observed high ratio of second vs. first generation of stars. If both effects contribute however, we suggest that a moderate amount of dynamical ejection together with a moderately top-heavy IMF is enough to account for the mass budget of present day globular clusters.
5. More than two generations? In principle our scenario is able to produce more than two stellar generations since we allow the newly formed generation to contain massive stars, thus possibly repeating the cycle.
6. Globular clusters as hosts of gravitational wave emission. In our scenario, the second generation is formed from the winds of massive stars. The massive stars themselves are predicted to end up mostly as compact objects, supporting the hypothesis that gravitational-waves should be expected from globular clusters.
7. Fraction of second generation correlates with GC mass. We predict that both the amount of accumulated mass and the total initial mass of the cluster correlate with the number of low-mass stars in the second generation. This may lead to the observed correlation between the mass of the GC and the extent of the multiple population phenomenon.
Our way of investigating the multiple population phenomenon in GCs by combining stellar evolutionary models with calculations of cluster wind dynamics, should be considered an important method of testing stellar theories in the future-especially those that are very hard to find observational evidence for. Such theories include metal-poor massive stars, both in single and binary systems. The potential of these systems in explaining exotic explosions (gravitationalwaves, gamma-ray bursts, several types of supernovae and superluminous supernovae) is quite high; nonetheless, many of the existing theories awaits future tests. We suggest that studying GCs by combining stellar models with cluster wind dynamics is a viable new approach by which these tests could be done.