# Rebounding Cores to Build Star Cluster Multiple Populations

G. Parmentier<EMAIL_ADDRESS>Astronomisches Rechen-Institut,
Zentrum für Astronomie der Universität Heidelberg, Mönchhofstr. 12-14, D-69120
Heidelberg, Germany  A. Pasquali  Astronomisches Rechen-Institut, Zentrum für
Astronomie der Universität Heidelberg, Mönchhofstr. 12-14, D-69120 Heidelberg,
Germany

(Accepted 2021 October 23)

###### Abstract

We present a novel approach to the riddle of star cluster multiple
populations. Stars form from molecular cores. But not all cores form stars.
Following their initial compression, such ’failed’ cores re-expand, rather
than collapsing. We propose that their formation and subsequent dispersal
regulate the gas density of cluster-forming clumps and, therefore, their core
and star formation rates. Clumps for which failed cores are the dominant core
type experience star formation histories with peaks and troughs. In contrast,
too few failed cores results in smoothly decreasing star formation rates. We
identify three main parameters shaping the star formation history of a clump:
the star and core formation efficiencies per free-fall time, and the time-
scale on which failed cores return to the clump gas. The clump mass acts as a
scaling factor. We use our model to constrain the density and mass of the
Orion Nebula Cluster progenitor clump, and to caution that the star formation
histories of starburst clusters may contain close-by peaks concealed by
stellar age uncertainties. Our model generates a great variety of star
formation histories. Intriguingly, the chromosome maps and O-Na anti-
correlations of old globular clusters also present diverse morphologies. This
prompts us to discuss our model in the context of globular cluster multiple
stellar populations. More massive globular clusters exhibit stronger multiple
stellar population patterns, which our model can explain if the formation of
the polluting stars requires a given stellar mass threshold.

galaxies: star clusters: general — stars: formation — ISM: clouds

##  1 Introduction

The presence of multiple stellar populations in star clusters remains a great
riddle. Increasingly complex patterns have been revealed in the color-
magnitude diagrams of Galactic globular clusters, which were once thought to
epitomise simple stellar populations (see e.g. Milone et al., 2008; Anderson
et al., 2009; Milone et al., 2015). These patterns stem from star-to-star
abundance variations in He (Piotto et al., 2007), and in C, N and O (Milone et
al., 2015). CNO abundance variations in individual globular clusters had
already been identified by earlier spectroscopic studies (e.g. Kraft, 1994),
and Na-Mg-Al abundance variations are also common (Kraft, 1994; Carretta et
al., 2006; Charbonnel, 2016, and references there in). These abundances are
either correlated (Na-N-He) or anti-correlated (Na-O, N-C), which is the
signature of hot-hydrogen burning (i.e. hydrogen-burning through the CNO-cycle
and the NeNa- and MgAl-chains, Langer et al., 1993). The current consensus is
that the low-mass stars enriched in Na-N-He and depleted in C and O formed out
of gas polluted by the ejecta of stars massive enough to host hot-hydrogen
burning. The nature of the polluters has remained debated, however. Candidate
polluters include asymptotic giant branch stars (mass $\simeq 4-8\,M_{\odot}$,
D’Antona et al., 2002; D’Ercole et al., 2008; D’Antona et al., 2016), fast-
rotating massive stars (mass $>20\,M_{\odot}$, Prantzos & Charbonnel, 2006;
Decressin et al., 2007), and super-massive stars (mass $\simeq
10^{4}\,M_{\odot}$, Denissenkov & Hartwick, 2014). Given their very different
masses, measurements of stellar age spreads in globular clusters could help
pinpoint the most likely scenario. Yet, this remains a challenging task,
hindered by the old age of globular clusters (Oliveira et al., 2020).

To complexify the picture further, the color-magnitude diagram of some younger
clusters presents extended main sequence turn-offs, as well as split main
sequences (for ages younger than 2 Gyr and 700 Myr, respectively). Observed in
the Large Magellanic Cloud clusters and in open clusters of the Galactic disk,
such features stem from different stellar rotation velocities (Dupree et al.,
2017; Marino et al., 2018; Cordoni et al., 2018).

Last but not least, some forming clusters present stellar age distributions
with peaks and troughs, suggesting that their star formation activity
periodically recedes, to be re-invigorated later on. Examples are the Orion
Nebular Cluster (ONC) in the Galactic disk (Beccari et al., 2017), and the
young cluster R136 at the heart of the 30 Dor HII region in the Large
Magellanic Cloud (Schneider et al., 2018). Interestingly, the age-span between
successive peaks is of order 1 Myr for both the ONC and R136. It remains
unclear whether and how these complex star formation histories relate to the
multiple stellar populations of old globular clusters.

To explain the bumpy star formation history of the ONC, Kroupa et al. (2018)
devise a scenario in which clusters with masses in the range
$600-20,000\,M_{\odot}$ form at the converging point of molecular gas inflows.
Newly-born O-stars shut off the inflow and star formation until they get
ejected out of the cluster via stellar interactions. The gas inflow then
resumes and builds the next stellar population.

We will explore another avenue, one in which the gas density of a cluster-
forming clump is regulated by the formation and dispersal of ’failed’ cores.
Molecular cores are often considered the precursors of individual stars and
small stellar systems (e.g. André et al., 2014; Ballesteros-Paredes et al.,
2020, and references therein). Not all cores are prestellar, however. Some,
possibly the majority, fail to form stars, ’bounce back’ and disperse in the
clump gas out of which they initially form (Taylor et al., 1996; Vázquez-
Semadeni et al., 2005; Clark & Bonnell, 2005; Gómez et al., 2007; Lada et al.,
2008; Belloche et al., 2011a). The cycle of their formation and dispersal must
therefore impact the mass and density of the clump gas, hence the core and
star formation rates. We consider these ’failed’/’rebounding’ /’dispersing’
cores as a source of gas replenishment possibly leading to another star
formation episode.

We will investigate whether different regions of the model parameter space can
render different star formation histories, from sequences of fairly discrete
star formation episodes to smoothly evolving star formation rates. Such a
finding could be a major gateway to understand the large variance observed in
the multiple stellar populations of old globular clusters. Consider the
distributions of stars along their O-Na anti-correlation (Carretta, 2015;
Charbonnel, 2016) and in their chromosome map (Milone et al., 2017). They
differ from cluster to cluster, being sometimes more clumpy, sometimes more
continuous. Are the successive star formation episodes of some forming
clusters (e.g. the ONC) the low-mass counterpart of the clumpy multiple
stellar populations in some old globular clusters?

For a new star formation episode to rise, the previous one must first recede.
We therefore opt for a model able to predict a decreasing star formation
activity in cluster-forming clumps. We adopt the volume-density driven
cluster-formation model of Parmentier & Pfalzner (2013), which predicts a star
formation rate declining with time under the assumptions of constant clump
mass and size. This decline is driven by two factors: (i) the decrease with
time of the available gas mass due to ongoing star formation in an isolated
constant-size clump, and (ii) the corresponding increase of the gas free-fall
time (see Parmentier et al., 2014, for a detailed discussion). Evidence for
decreasing star formation histories have been found for the Chamaeleon I
molecular cloud (Luhman, 2007; Belloche et al., 2011b) and, possibly, for the
Perseus molecular cloud (Hatchell et al., 2005) (see also Preibisch, 2012).
Additionally, the amplitude of the three star formation episodes in the ONC is
also a declining function of time (Beccari et al., 2017).

Parmentier & Pfalzner (2013) assume that the clump gas is directly converted
into stars. Yet, stars form out of molecular cores, which are themselves
condensations forming out of the clump gas. With this in mind, we add two
major elements to their model: (1) the formation of cores constitute an
intermediate stage in the conversion of the clump gas into stars, and (2) some
cores - rather than collapsing and forming stars - re-expand into the clump
gas reservoir, thereby re-supplying it in gas. In this preliminary study, we
focus on cluster-forming clumps with uniform density.

The outline of the paper is as follow. In Sec. 2, we introduce a preliminary
model, out of which we develop a criterion for the rise of a second star
formation episode. Section 3 reviews the evidence for failed cores in
simulations and observations of star-forming regions. In Sec. 4, we expand the
model of Sec. 2 by differentiating between star-forming and failed cores, and
we explore the parameter space in Sec. 5. We further upgrade the model in Sec.
6, where we relate the dissolution time-scale of the failed cores to the free-
fall time of the host clump, using the Orion Nebula Cluster as a test-bed. In
Sec. 7, we discuss how our model star formation histories could be related to
multiple stellar populations in old globular clusters, the limitations of our
current model, and how an initial gas density gradient will impact our current
results. We present our conclusions in Sec. 8.

##  2 Model Fundamentals

The model of Parmentier & Pfalzner (2013) predicts that the star formation
rate of cluster-forming clumps is a decreasing function of time. This behavior
stems from their model clumps having a constant radius and being isolated
(i.e. no inflow, no outflow). Declining star formation rates such as those
shown in Fig. 3 in Parmentier et al. (2014) and Fig. 4 in Parmentier (2019)
define the first star formation episode of our model clumps. We aim at
pinpointing the conditions leading to a second star formation episode, while
retaining the assumptions of an isolated clump with a fixed size. The pivotal
point here is: to form stars, a molecular clump must first form starless
molecular cores. Some of these cores are bound, gravitationally unstable,
collapse and form a star or a small stellar group (i.e. binary or triple). Not
all cores are prestellar, however. Others remain unbound and eventually bounce
back to the clump gas reservoir out of which they formed (e.g. Taylor et al.,
1996; Vázquez-Semadeni et al., 2005; Clark & Bonnell, 2005; Lada et al., 2008;
Belloche et al., 2011a, see also Sec. 3 for a discussion). The formation and
dispersal of these ’failed’ cores equates with removing some gas from the
clump gas reservoir (core formation), and restoring it at a later time (core
dispersal). When failed cores are abundant enough, their dispersal increases
the leftover gas mass of a clump, renewing its core and star formation. To
take into account the lifecycle of molecular cores - especially the non-
collapsing ones - will therefore yield star formation histories different from
those shown in Parmentier et al. (2014) and Parmentier (2019). In what
follows, we coin return time the time-span needed for non-collapsing cores to
disperse. The imprint of this parameter is the time-span between two
successive peaks of a star formation history.

###  2.1 Building Bumpy Core Formation Histories

To model the core formation history of a cluster-forming clump, we use Eq. 18
of Parmentier & Pfalzner (2013), substituting the core formation rate per unit
volume, $d\rho_{\rm core}(t)/dt$, for the star formation rate per unit volume
$d\rho_{\rm star}(t)/dt$, and the core formation efficiency per free-fall
time, $\epsilon_{\rm ff,co}$ for the star formation efficiency per free-fall
time. We obtain

$\frac{d\rho_{\rm core}(t)}{dt}=\frac{\epsilon_{\rm ff,co}}{\tau_{\rm
ff}(t)}\rho_{\rm gas}(t)=\sqrt{\frac{32G}{3\pi}}\epsilon_{\rm ff,co}\rho_{\rm
gas}(t)^{3/2}\;.$ (1)

In this equation, $\rho_{\rm gas}(t)$ is the volume density of the clump gas
at time $t$ after core formation onset, $\tau_{ff}(t)$ is the clump gas free-
fall time and $G$ is the gravitational constant. For a clear understanding of
how the gas density impacts our model, we consider clumps with a uniform
density. The evolution of clumps with an initial gas density gradient will be
investigated in a future paper.

Equation 1 shows that the core formation rate per unit volume at a given time
$t$ depends on the corresponding gas density, hence on the gas mass inside a
constant-radius clump. The gas mass equates with the clump total mass minus
the mass in cores and stars. At first, ongoing core formation depletes the gas
reservoir, thereby decreasing the core formation rate. The higher the initial
density of the clump gas (i.e. the shorter its initial free-fall time
$\tau_{\rm ff,0}$) and/or the greater the core formation efficiency per free-
fall time $\epsilon_{\rm ff,co}$, the higher the initial core formation rate,
but also the faster the gas density and core formation rate dwindle. We note
that a gas density gradient would hasten the system evolution as well (see
Fig. 4 in Parmentier, 2019).

For a second core formation episode to occur, the core formation rate must
stop declining and start rising. Eq. 1 shows that this happens if the clump
gas reservoir is re-supplied in gas. In our model, this gas comes from cores
formed at an earlier time, which now get dispersed following their failure to
form stars. Consider cores formed at time $t\gtrsim 0$ (shortly after core
formation onset). After a return time $t_{return}$, a fraction $f_{return}$ of
the mass in cores returns to the clump gas. If the gas supply from dispersing
cores overcomes the gas depletion due to ongoing core formation, the clump gas
mass and density are on the rise. The return fraction $f_{return}$ is
obviously an important parameter. Here we can consider two extreme scenarios.
If cores consist exclusively of non-collapsing ones, the return fraction is
unity, $f_{return}=1$. If only pre-stellar cores form, then the return
fraction is their envelope mass fraction and $f_{return}\simeq 2/3$ (i.e. the
star formation efficiency of a single core is $\simeq 1/3$; Matzner & McKee,
2000; Alves et al., 2007). We therefore adopt an intermediate value
$f_{return}=0.85$, and we keep it constant. For now, we assume that star-
forming cores release their envelope and failed cores disperse on the same
time-scale $t_{return}$. We will lift this hypothesis in Sec. 4. Armed with
such a model, we can already anticipate the possibility of generating not just
two, but multiple star formation episodes. Indeed, if the second core
formation episode is significant enough, it will feed the emergence of a third
one.

Having detailed how our preliminary model unfolds, we now consider a spherical
clump of mass $m_{clump}$ and radius $r_{clump}$, and we infer the evolution
with time of its gas density and core formation rate. The initial free-fall
time of the clump gas is

$\tau_{\rm ff,0}=\sqrt{\frac{3\pi}{32G\rho_{0}}}$ (2)

with $\rho_{0}$ the clump volume density, equivalently the initial gas volume
density:

$\rho_{0}=\rho_{gas}(t=0)=\frac{m_{clump}}{\frac{4}{3}\pi r_{clump}^{3}}\;.$
(3)

The top panel of Fig. 1 shows as the solid red line the evolution of the gas
density of a clump of mass $m_{clump}=10^{4}\,M_{\odot}$, radius
$r_{clump}=1\,$pc and core formation efficiency per free-fall time
$\epsilon_{\rm ff,co}=0.30$. The return time and the return fraction of the
core gas are set to $t_{return}=1$ Myr and $f_{return}=0.85$. The return time
is shown as the vertical dotted black line. These parameter values define the
fiducial model of Fig. 1. As long as $t<t_{return}$, the evolution of the gas
density is as predicted by Parmentier & Pfalzner (2013, their Eqs 17-18,
written for cores rather than for stars). It steadily decreases due to the
formation of new cores

$\frac{d\rho_{\rm gas}(t)}{dt}=-\frac{d\rho_{\rm core}(t)}{dt}\;.$ (4)

At $t\gtrsim t_{return}$, the gas density increases. This shows that the
dispersal of cores returns enough mass to the clump gas to compensate for the
gas used by ongoing core formation. The evolution of $\rho_{\rm gas}(t)$
already tells us that the core formation history will present two well
distinct peaks followed by a third weaker one and a long trail of decreasing
star formation activity.

How does modifying the fiducial model affect the time evolution of the gas
density? We consider the following modifications (one at a time): a 10-times
more massive clump ($m_{clump}=10^{5}\,M_{\odot}$), a twice as small clump
radius ($r_{clump}=0.5$pc), a three-times smaller core formation efficiency
per free-fall time ($\epsilon_{\rm ff,co}=0.1$) and an (unphysically) small
return fraction of the core mass ($f_{return}=0.3$). Results are presented
along with the fiducial model in the top panel of Fig. 1.

Models with a reduced return fraction $f_{return}$ (orange line) or a reduced
core formation efficiency per free-fall time $\epsilon_{\rm ff,co}$ (green
line) fail to replenish significantly the clump gas reservoir at
$t=t_{return}$, either because of too great a mass fraction remaining locked
inside cores (small $f_{return}$), or because the gas density has not
decreased enough to establish a strong gas density contrast between $t=0$ and
$t=t_{return}$ (small $\epsilon_{\rm ff,co}$).

Figure 1: Evolution with time of the clump gas density. Top panel: Results for
the fiducial model ($m_{clump}=10^{4}\,M_{\odot}$, $r_{clump}=1\,$pc,
$\epsilon_{\rm ff,co}=0.30$, $t_{return}=1$ Myr and $f_{return}=0.85$: solid
red line), for a ten-times more massive clump ($m_{clump}=10^{5}\,M_{\odot}$,
dash-dotted blue line), for a twice as compact clump ($r_{clump}=0.5$ pc:
dotted magenta line), for a three-times lower core formation efficiency per
free-fall time ($\epsilon_{\rm ff,co}=0.1$: dashed green line), and for a
$\sim 3$-times lower return fraction of the mass in cores ($f_{return}=0.3$:
dash-dotted orange line). Bottom panel: Results when varying the return time
of the core mass: $t_{return}/\tau_{\rm ff,0}=0.3,1.0,3.0,6.0,10.0$ (see key),
all other parameters being identical to the top panel. The fiducial model (red
line) is the same as in the top panel and corresponds to $t_{return}/\tau_{\rm
ff,0}=6.0$. Note the different extents of the $y$-axis in both panels.

In contrast, a higher initial gas density (obtained through either a larger
clump mass $m_{clump}$ or a smaller clump radius $r_{clump}$) promotes the
formation of successive star formation episodes. This is because the mass
fraction of a clump which gets turned into cores depends on time in units of
the free-fall time, $t/\tau_{\rm ff,0}$. The physical return time of cores
being fixed ($t_{return}=1$ Myr), the core return time $t_{return}/\tau_{\rm
ff,0}$ is longer for denser clumps. This builds a stronger contrast between
the gas densities at $t=0$ and $t=t_{return}$ since the clump gas density is
given a longer time-span $t_{return}/\tau_{\rm ff,0}$ to decrease. In turn,
this allows the core gas released at $t\gtrsim t_{return}$ to make a deeper
impact.

Figure 2: Same as Fig. 2 but for the core formation rates normalized to their
initial value

The importance of the return time of failed cores is highlighted in the bottom
panel of Fig. 1, where the fiducial model ($\tau_{\rm ff,0}=0.17$ Myr and
$t_{return}=1$ Myr $\equiv 6\,\tau_{\rm ff,0}$) is compared to models with
$t_{return}/\tau_{\rm ff,0}=0.3,1.0,3.0,10.0$. Clearly, longer return times
$t_{return}/\tau_{\rm ff,0}$ stimulate more core formation episodes.

Figure 2 shows the corresponding core formation rates normalized to their
initial value (same model parameters and color-coding as in Fig. 1). The
variations of the core formation rate are wider than those of the gas density,
since the core formation rate scales as the gas density to the power of 1.5
(Eq. 1). Fig. 2 highlights the ability of our model to create very diverse
core formation histories. This provides us with a hand to understand the
strong variance observed in cluster multiple stellar populations, an aspect we
shall discuss in Sec. 7.

The top panel of Fig. 2 may prompt us to conclude that dense and massive
clumps are more likely to generate multiple core formation episodes. This
conclusion is erroneous, however, and we will return to this important aspect
also in Sec. 7. Here, we stress that the model parameters shaping a core
formation history are the dimensionless parameters, namely, the core formation
efficiency per free-fall time $\epsilon_{\rm ff,co}$, the return fraction
$f_{return}$ and return time $t_{return}/\tau_{\rm ff,0}$ of the core mass. We
now develop a tool to quantify their impact.

###  2.2 Criterion for a second core formation episode

The rise of a second core formation episode requires the gas density, in the
aftermath of its steady decrease since core formation onset, to make an upturn
and to start increasing at $t=t_{return}$. This is only doable if non-
collapsing cores return to the clump gas a mass greater than that ”pumped” to
form new cores. Building on Eq. 1, the following equation expresses the clump
gas density at $t=t_{return}+dt$, $\rho_{\rm gas}(t_{return}+dt)$, as a
function of its counterpart at $t=t_{return}$, $\rho_{\rm gas}(t_{return})$:

$\displaystyle\rho_{\rm gas}(t_{return}+dt)$ $\displaystyle=\rho_{\rm
gas}(t_{return})$ (5) $\displaystyle-\sqrt{\frac{32G}{3\pi}}\epsilon_{\rm
ff,co}\rho_{\rm gas}(t_{return})^{3/2}\cdot dt$
$\displaystyle+f_{return}\sqrt{\frac{32G}{3\pi}}\epsilon_{\rm
ff,co}\rho_{0}^{3/2}dt\;.$

The second term on the right-hand side corresponds to the clump gas used by
ongoing core formation. The third term accounts for the cores formed over the
time-span $[0,{\rm d}t]$ that return to the clump gas over the time-interval
$[t_{return},t_{return}+{\rm d}t]$. This third term depends on the clump
initial gas density $\rho_{0}$, for this is the gas density at $t=0$ which
dictates the formation rate of the cores being dispersed at $t\gtrsim
t_{return}$. Following a few developments, we obtain the slope of the
evolution with time of the gas density at $t=t_{return}$:

$\displaystyle\left[\frac{d\rho_{\rm gas}(t)}{dt}\right]_{t=t_{return}}$
$\displaystyle=\frac{\rho_{\rm gas}(t_{return}+dt)-\rho_{\rm
gas}(t_{return})}{dt}$ (6) $\displaystyle=\frac{CoFR(t=0)}{V}\cdot
S(\epsilon_{\rm ff,co},f_{return},t_{return}/\tau_{\rm ff,0})$

with

$S(\epsilon_{\rm ff,co},f_{return},t_{return}/\tau_{\rm
ff,0})=\left[f_{return}-\left(1+\frac{1}{2}\epsilon_{\rm
ff,co}\frac{t_{return}}{\tau_{\rm ff,0}}\right)^{-3}\right]\;.$ (7)

In Eq. 6, $V$ is the clump volume and $CoFR(t=0)$ is the clump initial core
formation rate:

$\frac{CoFR(t=0)}{V}=\left[\frac{d\rho_{\rm
core}(t)}{dt}\right]_{t=0}=\sqrt{\frac{32G}{3\pi}}\epsilon_{\rm
ff,co}\rho_{0}^{3/2}\;.$ (8)

The function $S(\epsilon_{\rm ff,co},f_{return},t_{return}/\tau_{\rm ff,0})$,
given by Eq.7, determines whether the slope of $\rho_{gas}(t)$ turns positive
at $t=t_{return}$, or not. For instance, if the cores return no gas
($f_{return}=0$), then $S<0$ and $\left[\frac{d\rho_{\rm
gas}(t)}{dt}\right]_{t=t_{return}}<0$. That is, the gas density keeps
decreasing and there is no second core formation episode. Conversely, if
$S(\epsilon_{\rm ff,co},f_{return},t_{return}/\tau_{\rm ff,0})>0$,
$t=t_{return}$ corresponds to the onset of a second core formation episode.

Figure 3: The $S$-function in dependence of the core formation efficiency per
free-fall time $\epsilon_{\rm ff,co}$, and of the return fraction $f_{return}$
and return time $t_{return}/\tau_{\rm ff,0}$ of the mass in cores (Eq. 7). It
is color-coded according to the right-hand-side bar. When $S>0$ (red-to-yellow
area), the clump gas density at $t=t_{return}$ increases, possibly yielding a
second core formation episode. When $S<0$ (red-to-black area), the clump gas
density at $t=t_{return}$ keeps decreasing, making impossible a second core
formation episode. The white symbols depict the locus of points for which
$S=0$.

Figure 3 maps $S(f_{return},t_{return}/\tau_{\rm ff,0})$ for $\epsilon_{\rm
ff,co}=0.1$ and $\epsilon_{\rm ff,co}=0.5$. Its value is color-coded according
to the right-hand-side bar, and the white symbols part the region for which
$S>0$ from that with $S<0$. When $S>0$ (red-to-yellow area), the clump gas
gains, at $t\gtrsim t_{return}$, more mass from failed cores formed at
$t\gtrsim 0$ than it loses to ongoing core formation. When $S<0$ (red-to-black
area), the formation of new cores consumes more clump gas than what failed
cores provide. The higher the core formation efficiency per free-fall time,
the wider the area of the $(f_{return},t_{return}/\tau_{\rm ff,0})$ parameter
space for which $S>0$. For a given core formation efficiency per free-fall
time, a long return time of the core gas $t_{return}/\tau_{\rm ff,0}$ and a
large return mass fraction $f_{return}$ favor a gas density increase at
$t=t_{return}$ and, therefore, an increase of the core formation rate. Such
conclusions have already been suggested by Fig. 1.

We note that since, at core formation onset,

$\left[\frac{d\rho_{\rm gas}(t)}{dt}\right]_{t=0}=-\left[\frac{d\rho_{\rm
core}(t)}{dt}\right]_{t=0}=-\frac{CoFR(t=0)}{V}\,,$ (9)

Eq. 6 can be rewritten as:

$\displaystyle\left[\frac{d\rho_{\rm gas}(t)}{dt}\right]_{t=t_{return}}=$
$\displaystyle-\left[\frac{d\rho_{\rm gas}(t)}{dt}\right]_{t=0}$ (10)
$\displaystyle\times S(\epsilon_{\rm ff,co},f_{return},t_{return}/\tau_{\rm
ff,0})\;.$

As can be seen, $S$ compares how fast the gas density varies at $t=t_{return}$
and $t=0$. Since the upper limit on $S$ is $1$ (when $f_{return}=1$ and
$t_{return}/\tau_{\rm ff,0}\rightarrow\infty$), Eq. 10 expresses that, in
absolute value, gas density variations at $t\gtrsim t_{return}$ cannot be
faster than those at $t\gtrsim 0$.

The model parameters promoting a second core formation episode are a high core
formation efficiency per free-fall time $\epsilon_{\rm ff,co}$, a large mass
fraction $f_{return}$ of core gas not ending into stars, and a long return
time of the core gas to the clump gas $t_{return}/\tau_{\rm ff,0}$. A second
core formation episode can be compared to a bouncing ball. A ball bounces back
if (1) it is dropped from a high height, and (2) its encounter with the ground
is as elastic as possible, such that it does not lose much of the kinetic
energy accumulated through its initial fall. The height equivalent in our
model is $\epsilon_{\rm ff,co}t_{return}/\tau_{\rm ff,0}$, which determines
which fraction of the clump mass has been turned into cores at $t=t_{return}$:
the higher $\epsilon_{\rm ff,co}$, the longer $t_{return}/\tau_{\rm ff,0}$,
the greater the initial gas mass fraction turned into cores, and the more the
gas density has dropped compared to its initial value. As for our $f_{return}$
return fraction, its bouncing ball counterpart is the near-conservation of the
ball kinetic energy: the higher $f_{return}$, the higher the core gas mass
released at $t=t_{return}$, and the higher the gas density and core formation
rate - like the ball - bounce back.

The high core formation efficiencies per free-fall time required for multiple
core formation episodes (e.g. $\epsilon_{\rm ff,co}=0.30$) are to be
contrasted with the star formation efficiencies per free-fall time measured in
the Galactic disk. The latter are of order one per cent (Krumholz & Tan, 2007;
Pokhrel et al., 2021). Given that individual star-forming cores have a star
formation efficiency of $\simeq 1/3$ (Matzner & McKee, 2000; Alves et al.,
2007), the formation efficiency per free-fall time of prestellar cores should
not exceed a few per cent. Multiple core formation episodes therefore require
non-collapsing cores to dominate the core population. This is also indicated
by our high return fraction ($f_{return}=0.85$). We now review the evidence
for significant populations of failed cores in molecular gas clouds, as found
by observations and simulations.

##  3 Evidence for failed cores in molecular clouds

Belloche et al. (2011a) map the nearby molecular cloud Chamaeleon I in dust
continuum emission at 870$\micron$. They find that 17% of its starless cores
have a mass higher than the critical Bonnor-Ebert mass111The Bonnor-Ebert mass
is the critical mass beyond which isothermal cores cannot be in hydrostatic
equilibrium and collapse (Bonnor, 1956; Tan et al., 2014).. That is, only a
minority of the Chamaeleon I cores are gravitationally unstable, and they tend
to be the most massive ones. Additionally, the mass function of the starless
cores is over-populated at its low-mass end, while being consistent with the
Salpeter (1955) slope at the high-mass end. Belloche et al. (2011a) therefore
conclude that most starless cores in Chamaeleon I are not prestellar. Belloche
et al. (2011b) further suggest that some of these low-mass gravitationally-
stable cores may never reach their critical mass and ”get dispersed” in the
future. Their interpretation is supported by a number of numerical
simulations.

According to Clark & Bonnell (2005), most of the clumps (their terminology for
cores) produced by supersonic turbulence are unbound. In one of their
simulations, unbound clumps outnumber bound ones by two orders of magnitude.
Unbound clumps/cores tend to be the low-mass ones, while prestellar cores are
the most massive ones (since they have to reach the Jeans mass in the first
place). A consequence of this is the absence of a one-to-one mapping from
cores to stars and, therefore, from the core mass function to the stellar
initial mass function (Klessen, 2001; Clark & Bonnell, 2005).

The nearby Pipe molecular cloud presents us with yet another population of
failed cores. 159 cores with a mean volume density $\simeq 7\cdot
10^{3}\,cm^{-3}$ are detected by Lada et al. (2008). A comparison of their
escape velocity with their - mostly thermal - gas velocity dispersion show
most of them to be unbound (Lada et al., 2008, their fig. 4). The transition
between the bound and unbound regimes occurs at a core mass of
$2-3\,M_{\odot}$, with 140 cores being less massive than $3\,M_{\odot}$.
Unless they experience a decrease of their internal pressure, or an increase
of their external pressure (or a combination of both), such cores are unlikely
to be prestellar. Lada et al. (2008) stress that unbound cores should persist
for one or more sound crossing-times before they disperse in their
surroundings (see also Keto et al., 2006). The radius of the unbound cores is
of order 0.075 pc (Lada et al., 2008, their Fig. 6) and the sound speed in an
isothermal gas at a temperature of $10$K is $c_{s}=0.2\,km\cdot s^{-1}$. The
sound crossing time of the unbound cores is therefore
$\tau_{core,cross}=2\cdot 0.075\,{\rm pc}/(0.2\,{\rm pc}\cdot{\rm
Myr}^{-1})\simeq 0.75$ Myr. Unbound cores in the Pipe cloud should thus
persist as coherent entities for about 1 Myr, possibly more (Lada et al.,
2008).

Based on MHD simulations, Vázquez-Semadeni et al. (2005) too find that some
cores/clumps form and dissolve later on. Considering a turbulent isothermal
molecular cloud with a mean number density $n_{0}=500\,cm^{-3}$, they infer a
life-expectancy of 0.5 to 3.0 Myr for its non-collapsing cores, similar to
what Lada et al. (2008) estimate for the Pipe cloud cores. Note that the
simulated cloud and the Pipe cloud have mean densities within an order of
magnitude of each other (the mean density of the Pipe cloud is $\simeq
100\,cm^{-3}$, Kainulainen et al., 2014) .

The life-expectancy of failed cores equates with the return time $t_{return}$
of our preliminary model and, therefore, with the periodicity with which a
clump unfolds its successive star formation episodes. Both Vázquez-Semadeni et
al. (2005) and Lada et al. (2008) estimate this life-expectancy to be of order
1 Myr. Interestingly, this is similar to the time-span separating the three
stellar populations of the ONC (Beccari et al., 2017, their Fig. 6), and Fig.
5 of Schneider et al. (2018) reveals a similar periodicity for the R136
cluster in 30 Dor. We shall come back to this in Sec. 6.

##  4 Failed- versus Star-Forming Cores

The preliminary model of Section 2 assumes that star-forming and failed cores
release their gas on the same time-span $t_{return}$. Yet, both types of core
evolve on different time-scales. First, Vázquez-Semadeni et al. (2005) stress
that the re-expansion of failed cores is slower than their initial
compression, as a result of self-gravity. Second, their MHD simulations show
that failed cores are often less dense than prestellar ones, implying a longer
free-fall time. Therefore, for at least two reasons (self-gravity and longer
free-fall time), failed cores evolve more slowly than star-forming ones. This
will help us meet one of the criterions unveiled in Sec. 2, namely, that
failed cores should have long return times $t_{return}/\tau_{\rm ff,0}$ to
stimulate a second core formation episode. In this section, we assign
different time-scales and different return fractions to failed and star-
forming cores. In the case of star-forming cores, their envelope is the mass
they return to their surroundings (Matzner & McKee, 2000).

We assume that collapsing/star-forming cores have a mean density two orders of
magnitude higher than that of their host clump. Their free-fall time
$\tau_{ff,colco}$ is thus ten times shorter than the clump’s:
$\tau_{ff,colco}=0.1\tau_{\rm ff,0}$. Their lifetime ranges from one to ten
times their free-fall time, depending on the magnetic field (e.g. Vázquez-
Semadeni et al., 2005; Galván-Madrid et al., 2007; Ward-Thompson et al., 2007,
and references therein). We therefore assume that the time-span between the
formation of prestellar cores and the release of their envelope is 3 times
their free-fall time. This defines the return time of their envelope to the
clump gas as

$t_{ret,colco}=3\tau_{ff,colco}=0.3\tau_{\rm ff,0}\;.$ (11)

Models yielding multiple core and star formation episodes hinge on failed
cores outnumbering star-forming ones (see Fig. 3 and below). The exact values
of the density contrast between star-forming cores and their host clump and of
$t_{ret,colco}/\tau_{\rm ff,0}$ are thus minor issues.

As for the return time of the failed cores, namely, the time-span on which
they re-expand into their surroundings, we fix it to 1 Myr (see Sec. 3):

$t_{ret,fco}=1\,{\rm Myr}\;.$ (12)

We will lift this assumption in Sec. 6, a single physical return time
$t_{ret,fco}=1$ Myr for failed cores being an oversimplification, especially
when considering very different host clump densities $\rho_{0}$.

Let us now quantify the core mass fractions (i.e. return fractions) associated
to the return times defined above. We introduce the star formation efficiency
per free-fall time $\epsilon_{\rm ff,st}$ in our model and derive the return
fractions as a function of $\epsilon_{\rm ff,co}$ and $\epsilon_{\rm ff,st}$.
The mass $dm_{co}(t)$ in cores formed over a time interval $[t,t+dt]$ consists
of three components: (1,2) stars and envelopes formed and released,
respectively, by star-forming cores, and (3) failed cores. The following set
of equations relates the mass fraction of failed cores, $F_{fco}$, to the mass
fractions of envelopes and stars from star-forming cores, $F_{env}$ and
$F_{*}$, respectively:

$1=F_{fco}+F_{*}+F_{env}\;.$ (13) $F_{env}=2F_{*}\;.$ (14)
$\frac{\epsilon_{\rm ff,st}}{\epsilon_{\rm ff,co}}=F_{*}$ (15)

For star-forming cores, we assume that the mass of their envelope is twice the
mass of their stars (Eq. 14), i.e. the star formation efficiency of a star-
forming core is one-third (Matzner & McKee, 2000; Alves et al., 2007). Eq. 15
quantifies which fraction $F_{*}$ of the total mass in cores eventually ends
into stars. We define it as the ratio of the star formation efficiency per
free-fall time $\epsilon_{\rm ff,st}$ of the host clump to its core formation
efficiency per free-fall time $\epsilon_{\rm ff,co}$. Note that $\epsilon_{\rm
ff,co}$ accounts for all cores, both failed and star-forming ones. For an
assumed pair ($\epsilon_{\rm ff,co}$, $\epsilon_{\rm ff,st}$), the above
system of equations yields the respective contributions of stars, star-forming
core envelopes, and failed cores to the core mass $dm_{co}(t)$ formed over the
time-span $[t,t+dt]$. Note that in Sec. 2.2, $\epsilon_{\rm ff,co}$ and
$f_{return}$ were independent of each other. With the above set of equations
and a fixed $\epsilon_{\rm ff,st}$, $\epsilon_{\rm ff,co}$ and $F_{fco}$ are
now tied to each other, an increase of $\epsilon_{\rm ff,co}$ implying an
increase of $F_{fco}$.

Table 1: Mass fractions of the stars ($F_{*}$) and envelopes ($F_{env}$) from star-forming cores, and of the failed cores ($F_{fco}$), for given $\epsilon_{\rm ff,co}$ and $\epsilon_{\rm ff,st}$.  $\epsilon_{\rm ff,co}$ |  $F_{*}$ |  $F_{env}$ |  $F_{fco}$
---|---|---|---
$\epsilon_{\rm ff,st}=0.01$
0.030 | 0.33 | 0.67 | 0.00
0.100 | 0.10 | 0.20 | 0.70
0.300 | 0.03 | 0.07 | 0.90
0.500 | 0.02 | 0.04 | 0.94
$\epsilon_{\rm ff,st}=0.026$
0.078 | 0.333 | 0.667 | 0.000
0.300 | 0.087 | 0.173 | 0.740
0.500 | 0.052 | 0.104 | 0.844

Table 1 provides examples for various core formation efficiencies per free-
fall time $\epsilon_{\rm ff,co}$ and for $\epsilon_{\rm ff,st}=0.01$ (Krumholz
& Tan (2007); see also Fig. 4 in Parmentier & Pasquali (2020)) and
$\epsilon_{\rm ff,st}=0.026$ (Pokhrel et al., 2021). When $\epsilon_{\rm
ff,st}=0.01$ and $\epsilon_{\rm ff,co}=0.03$ (or when $\epsilon_{\rm
ff,st}=0.026$ and $\epsilon_{\rm ff,co}=0.078$), failed cores make no
contribution to the core mass budget (since $F_{*}=\epsilon_{\rm
ff,st}/\epsilon_{\rm ff,co}=1/3$ implies $F_{env}=2/3$ and $F_{fco}=0$). In
the absence of failed cores, the only gas returned to the intra-cluster medium
is that of star-forming core envelopes, which, per clump free-fall time,
amounts to a few per cent of the clump gas (i.e. twice the star formation
efficiency per free-fall time). This is an extremely small fraction compared
to the $\gtrsim 90$ % of the clump mass which, per free-fall time, remain in
its gas reservoir. Unsurprisingly, a model with no failed cores fails to
generate multiple peaks in the core formation history (see below). As the core
formation efficiency per free-fall time $\epsilon_{\rm ff,co}$ increases -
while retaining the same star formation efficiency per free-fall time
$\epsilon_{\rm ff,st}$ \- the contribution of star-forming cores
$F_{env}+F_{*}$ dwindles, while that of failed cores $F_{fco}$ rises (see
Table 1).

Figure 4: Evolution with time of the core formation rate (dash-dotted lines)
and stellar mass (solid lines) of a clump with an initial gas mass
$m_{clump}=10^{4}\,M_{\odot}$, a radius $r_{clump}=1$ pc and a star formation
efficiency per free-fall time $\epsilon_{\rm ff,st}=0.01$. With higher core
formation efficiencies per free-fall time $\epsilon_{\rm ff,co}$, failed cores
contribute more to the core mass budget (see Table 1), and multiple peaks get
carved in the core formation history

The four cases of Table 1 with $\epsilon_{\rm ff,st}=0.01$ are illustrated in
Fig. 4. The mass and radius of the model clump are
$m_{clump}=10^{4}\,M_{\odot}$ and $r_{clump}=1$ pc, respectively, yielding an
initial free-fall time $\tau_{\rm ff,0}=0.167$ Myr for the clump gas. The
return times of the core gas to the intra-cluster medium are as given by Eqs
11 ($t_{ret,colco}=0.3\tau_{\rm ff,0}=0.05$ Myr) and 12 ($t_{ret,fco}=1\,{\rm
Myr}=6\tau_{\rm ff,0}$). Dash-dotted and solid lines show the time evolutions
of the core formation rate and of the stellar mass, respectively. Note that
stellar masses start rising at $t=t_{ret,colco}$, that is, when the first
stars emerge from their gaseous cocoon.

A higher core formation efficiency per free-fall time yields a higher core
formation rate initially, but also a faster decrease of this one since the
greater gas fraction trapped in cores reduces the clump gas mass accordingly.
An increasingly steeper gas density gradient inside clumps would produce a
similar effect (Fig. 4 in Parmentier, 2019). At $t\gtrsim t_{ret,fco}=1$ Myr,
i.e. after one return time of the failed cores, the core formation rate may
make an upturn, a rise borne by the gas supplied by the failed cores formed at
$t\gtrsim 0$ Myr. The reasons as to why larger $\epsilon_{\rm ff,co}$ drive
more star formation episodes is two-fold: (1) as in Sec. 2.2, a larger
$\epsilon_{\rm ff,co}$ implies a greater decrease of the gas density after a
given time-span $t/\tau_{\rm ff,0}$ (that is, the bouncing ball is dropped
from a greater height), and (2) larger $\epsilon_{\rm ff,co}$ with a fixed
star formation efficiency per free-fall time $\epsilon_{\rm ff,st}$ implies
greater return fractions from failed cores (see Table 1; that is, the ball
retains more of its kinetic energy while bouncing). These two effects can be
visualized in Fig. 3 as (1) a widening of the region for which $S>0$ (i.e. the
white symbols and the red-to-yellow area are brought downwards for increasing
$\epsilon_{\rm ff,co}$), and (2) an increase of $f_{return}$ (i.e. moving to
the right of the horizontal axis).

That a higher core formation efficiency per free-fall time is more likely to
generate successive core formation episodes comes at a cost, however: higher
$\epsilon_{\rm ff,co}$ yield lower stellar masses at any given time $t$ (see
solid lines in Fig. 4). This reduced star formation efficiency reflects the
larger mass in failed cores produced for higher $\epsilon_{\rm ff,co}$, and
the correspondingly lower gas masses available for further core formation.
This effect is not tremendous as increasing $\epsilon_{\rm ff,co}$ by a factor
of $17$ (i.e. $0.50/0.03$) reduces stellar masses and star formation
efficiencies by a factor $\simeq 3$.

##  5 Peak number, gas density and core formation efficiency

The color-magnitude diagram of the ONC presents 3 well-defined pre-main
sequences, possibly suggesting three stellar populations with an age
separation of about 1 Myr (Beccari et al., 2017). That is, the ONC would have
experienced three successive core and star formation episodes. Building on the
model refined in the previous section, we now explore which core formation
efficiency per free-fall time $\epsilon_{\rm ff,co}$ and which initial gas
density $\rho_{0}$ are required to carve peaks in the star formation history
of a clump with an imposed periodicity of 1 Myr. The results will be used in
the next section to constrain the properties of the gaseous progenitor of the
ONC based on its observed star formation history.

Figure 5: Logarithmic stellar age distributions for 3 core formation
efficiencies per free-fall time (from top to bottom panel: $\epsilon_{\rm
ff,co}=0.1,0.3,0.5$), and 4 initial gas volume densities $\rho_{0}$ (see key
in middle panel). The adopted densities are representative of cluster-forming
regions of the Galactic disk. Model parameters are provided in Table 2.

We test three core formation efficiencies per free-fall time, namely,
$\epsilon_{\rm ff,co}=0.1,0.3,0.5$. The star formation efficiency per free-
fall time is fixed, and we adopt $\epsilon_{\rm ff,st}=0.01$. The
corresponding core mass fractions converted into stars ($F_{*}$) and released
to the clump gas as failed cores ($F_{fco}$) are given in Table 1. As for the
clump initial gas density, we test values ranging from $\rho_{0}=10^{2}$ to
$10^{5}\,M_{\odot}\cdot pc^{-3}$ in logarithmic steps of 1. The lower bound is
of order the current mean density of the ONC, namely, $4500\,M_{\odot}$ of
stars and gas enclosed within a radius of 2 pc (Hillenbrand & Hartmann, 1998).
The upper bound corresponds to the density of the ”ultra-high density
material” mapped at pc-scale by Hacar et al. (2020) in the active star-forming
cloud OMC-1, located in the vicinity of the ONC. Densities in the range
$10^{4}{\rm-}10^{5}\,M_{\odot}\cdot pc^{-3}$ are also those characterizing the
central regions of starburst clusters in the Galactic disk (e.g. Westerlund 1
and NGC 3603; see top panel of Fig. 10 in Parmentier et al., 2014). To vary
the clump density, we vary its mass $m_{clump}$ while keeping its radius
$r_{clump}$ constant. We assume $r_{clump}=2$ pc, which is approximately the
current ONC radius (Hillenbrand & Hartmann, 1998; Da Rio et al., 2017). The
corresponding clump masses are given in Table 2. All models are assigned a
star formation duration of 4 Myr. We note that stellar feedback is not taken
into account, and that stellar-feedback driven gas expulsion may terminate
star formation at an earlier time than assumed here (e.g. Rahner et al., 2017;
Kroupa et al., 2018). We impose a physical return time $t_{ret,fco}=1$ Myr for
the failed cores, as dictated by the time-span separating the ONC three
stellar populations (see Fig. 6 in Beccari et al., 2017). Combining the clump
density with the physical return time of the failed cores provides their
return time in units of the clump free-fall time $t_{ret,fco}/\tau_{\rm ff,0}$
(see Table 2). The return time of failed cores will be discussed further in
Sec. 6. With $t_{ret,fco}=1$ Myr and a star formation duration of 4 Myr, we
expect our model star formation histories to present from 0 to 4 peaks,
depending on the core formation efficiency per free-fall time $\epsilon_{\rm
ff,co}$ and failed core return time $t_{ret,fco}/\tau_{\rm ff,0}$.

Table 2: Input parameters $\epsilon_{\rm ff,co}$ and $\rho_{0}$ for the models of Fig. 5. $r_{clump}=2$ pc, $\epsilon_{\rm ff,st}=0.01$ and $t_{ret,fco}=1$ Myr are imposed. $m_{st,4Myr}$ is the stellar mass at $t=4$ Myr Model |  $\epsilon_{\rm ff,co}$ |  $\rho_{0}$ | $t_{ret,fco}$ |  $m_{clump}$ |  $m_{st,4Myr}$
---|---|---|---|---|---
ID |  |  $[M_{\odot}\cdot pc^{-3}]$ | [$\tau_{\rm ff,0}$] |  $[M_{\odot}]$ |  $[M_{\odot}]$
A1 | 0.10 | $10^{2}$ |  1.2 | $3.3\cdot 10^{3}$ | $1.3\cdot 10^{2}$
A2 | 0.10 | $10^{3}$ |  3.9 | $3.3\cdot 10^{4}$ | $3.5\cdot 10^{3}$
A3 | 0.10 | $10^{4}$ | 12.4 | $3.3\cdot 10^{5}$ | $6.9\cdot 10^{4}$
A4 | 0.10 | $10^{5}$ | 39.1 | $3.3\cdot 10^{6}$ | $1.0\cdot 10^{6}$
B1 | 0.30 | $10^{2}$ |  1.2 | $3.3\cdot 10^{3}$ | $1.1\cdot 10^{2}$
B2 | 0.30 | $10^{3}$ |  3.9 | $3.3\cdot 10^{4}$ | $2.2\cdot 10^{3}$
B3 | 0.30 | $10^{4}$ | 12.4 | $3.3\cdot 10^{5}$ | $3.4\cdot 10^{4}$
B4 | 0.30 | $10^{5}$ | 39.1 | $3.3\cdot 10^{6}$ | $0.4\cdot 10^{6}$
C1 | 0.50 | $10^{2}$ |  1.2 | $3.3\cdot 10^{3}$ | $0.9\cdot 10^{2}$
C2 | 0.50 | $10^{3}$ |  3.9 | $3.3\cdot 10^{4}$ | $1.6\cdot 10^{3}$
C3 | 0.50 | $10^{4}$ | 12.4 | $3.3\cdot 10^{5}$ | $2.3\cdot 10^{4}$
C4 | 0.50 | $10^{5}$ | 39.1 | $3.3\cdot 10^{6}$ | $0.3\cdot 10^{6}$

The logarithmic stellar age distributions are presented in Fig. 5. The
downward black arrows at an age of 5 Myr indicate the onset of core formation.
Note in that respect the slight shift with the onset of the stellar age
distribution, especially noticeable for the lowest-density models (i.e.
$\rho_{0}=100\,M_{\odot}\cdot pc^{-3}$; red lines). This small shift
corresponds to the time-delay $t_{ret,colco}=0.3\tau_{\rm ff,0}$ required by
star-forming cores to shed their envelope in the intra-cluster medium (Eq.
11).

Figure 5 reveals that the presence of one peak per Myr for $\geq 3$Myr
requires a clump density of at least $\rho_{0}\simeq 10^{3}\,M_{\odot}\cdot
pc^{-3}$ when the core formation efficiency per free-fall time is
$\epsilon_{\rm ff,co}=0.5$ (dashed green line in bottom panel). Lower core
formation efficiencies per free-fall time must be ”compensated” with higher
initial gas densities (e.g. $\rho_{0}\simeq 10^{5}\,M_{\odot}\cdot pc^{-3}$
when $\epsilon_{\rm ff,co}=0.1$, dashed magenta line in top panel). We note
that models with a density similar to that of the present-day ONC (i.e.
$\rho_{0}\simeq 10^{2}\,M_{\odot}\cdot pc^{-3}$) fail to carve multiple peaks
in their star formation history, irrespective of the adopted core formation
efficiency per free-fall time (solid red lines in all panels).

How can we disentangle between the models with
$\rho_{0}=10^{3}\,M_{\odot}\cdot pc^{-3}$ and $\epsilon_{\rm ff,co}=0.5$ and
$\rho_{0}=10^{5}\,M_{\odot}\cdot pc^{-3}$ and $\epsilon_{\rm ff,co}=0.1$
(models C2 and A4 in Table 2)? The initial gas mass $m_{clump}$ and the
stellar mass $m_{st}$ at $t=4$ Myr constitute powerful tools here (see last
two columns of Table 2). The C2 model ($\rho_{0}=10^{3}\,M_{\odot}\cdot
pc^{-3}$ and $\epsilon_{\rm ff,co}=0.5$) has formed $1.6\cdot
10^{3}\,M_{\odot}$ in stars over 4 Myr (or $1.4\cdot 10^{3}\,M_{\odot}$ over
3.3 Myr, which is the ONC star formation duration we will infer in the next
section). This mass agrees much better with the current stellar content of the
ONC ($1800\,M_{\odot}$; Hillenbrand & Hartmann, 1998) than the stellar masses
predicted with $\rho_{0}=10^{5}\,M_{\odot}\cdot pc^{-3}$. Given the adopted
radius $r_{clump}=2$ pc, an initial gas density
$\rho_{0}=10^{5}\,M_{\odot}\cdot pc^{-3}$ implies an initial gas mass of
several million solar masses. This is 5 orders of magnitude higher than the
mass of ”ultra-high density material” detected by Hacar et al. (2020) in OMC-1
($\simeq 30\,M_{\odot}$ at $\simeq 10^{5}\,M_{\odot}\cdot pc^{-3}$). The
comparison of the observed gas and star masses to their model counterparts
discards therefore the model A4. In other words, considering both the shape
and normalisation of a star formation history allows one to disentangle
whether the most appropriate model combines a high core formation efficiency
per free-fall time to a low-density clump, or the other way round.

In the case of the ONC, we discard models with
$\rho_{0}=10^{5}\,M_{\odot}\cdot pc^{-3}$ not just because of their too large
gas and stellar masses, but also because the imposed physical return time
$t_{ret,fco}=1\,$Myr equates with $t_{ret,fco}=40\tau_{\rm ff,0}$ ($\tau_{\rm
ff,0}=0.026$ Myr). Our model assumptions (e.g. constancy of the clump radius
$r_{clump}$) are unlikely to remain valid that long.

##  6 The inter-peak age-span as a density probe

A further observable which helps shed light on the initial gas density of a
cluster-forming clump is the age span separating successive peaks of its
stellar age distribution. Since high-density systems evolve faster than low-
density ones, we intuitively expect the frequency of the peaks to be higher
when the clump density is larger. Yet, we have so far assumed that the return
time of failed cores $t_{ret,fco}$ is independent of other model parameters.
We now lift our assumption of a fixed return time $t_{ret,fco}=1\,$Myr for the
failed cores, and we relate $t_{ret,fco}$ to the clump free-fall time
$\tau_{\rm ff,0}$, as already done for $t_{ret,colco}$, the time needed for
star-forming cores to shed their envelope (Eq. 11).

###  6.1 The return time of the failed cores

Vázquez-Semadeni et al. (2005) infer an analytical estimate of how long it
takes for a re-expanding core to double its radius, that is, to drop its mean
density by almost an order of magnitude, when the external pressure exerted by
the surrounding gas is negligible. We adopt this expansion time as a proxy to
the failed core return time $t_{ret,fco}$. The left panel of Fig. 1 in
Vázquez-Semadeni et al. (2005) shows that the return time $t_{ret,fco}$ is in
the following range

$2\frac{R_{eq}}{c_{s}}\lesssim t_{ret,fco}\lesssim 6\frac{R_{eq}}{c_{s}},$
(16)

with $R_{eq}$ the core equilibrium radius and $c_{s}$ the core gas sound
speed. The exact value depends on how far from equilibrium the core is
following its initial compression. Since:

$\frac{R_{eq}}{c_{s}}=\sqrt{\frac{\pi}{G\rho_{co}}}$ (17)

for a core density $\rho_{co}$ (e.g. Binney & Tremaine, 1987), we can
write222It is worth noting here that Vázquez-Semadeni et al. (2005) adopt
$R_{eq}/c_{s}$ as their definition of the free-fall time, which is therefore
longer than ours by a factor $\sqrt{32/3}$.:

$t_{ret,fco}=\top\sqrt{\frac{\pi}{G\rho_{fco}}}=\top\sqrt{\frac{32}{3}}\sqrt{\frac{3\pi}{32G\rho_{fco}}}=3.3\top\tau_{ff,fco}$
(18)

with $2\lesssim\top\lesssim 6$.

Now let us relate the mean density $\rho_{fco}$ and free-fall time
$\tau_{ff,fco}$ of the failed cores to those of the host clump. According to
Vázquez-Semadeni et al. (2005), core-to-cloud mean density ratios of $\sim
30\rm{-}100$ yield gravitationally unstable cores when these are modeled as
Bonnor-Ebert spheres. For our failed cores, we therefore adopt a weaker
contrast $\mathcal{C}=10$ between their density $\rho_{fco}$ and the host
clump density $\rho_{0}$:

$\mathcal{C}=\frac{\rho_{fco}}{\rho_{0}}=10\;.$ (19)

The free-fall time of failed cores then obeys $\tau_{ff,fco}=\tau_{\rm
ff,0}/\sqrt{\mathcal{C}}$, and we can modify Eq. 18 as

$t_{ret,fco}=3.3\top\frac{\tau_{\rm ff,0}}{\sqrt{\mathcal{C}}}\;.$ (20)

With $\mathcal{C}\simeq 10$, the return time/re-expanding time of the failed
cores becomes:

$t_{ret,fco}\simeq\top{\tau_{\rm ff,0}}$ (21)

with $2\lesssim\top\lesssim 6$.

###  6.2 The ONC as a test-bed

In Sec. 5, we have made a preliminary estimate of the initial gas density and
free-fall time of the ONC progenitor, based on the 1 Myr star formation
periodicity imposed by the observations. We have found:
$\rho_{clump}^{ONC}\simeq 10^{3}\,M_{\odot}\cdot pc^{-3}$, hence
$\tau_{ff,0}^{ONC}\simeq 0.26$ Myr. Does the corresponding $\top$-value make
sense? Taking the 1-Myr time-span separating the successive stellar
populations of the ONC as a proxy of the failed-core return time, we obtain
$\top\simeq t_{ret,fco}/\tau_{\rm ff,0}\simeq 1\,{\rm Myr}/0.26\,{\rm
Myr}\simeq 3.9$, which falls in the range of viable $\top$-values quoted
above. We conclude that the age-span between successive peaks of a stellar age
distribution looks as a promising probe of the initial gas density. A
prerequisite for this is that the return time of failed cores in units of the
clump free-fall time, namely, the $\top$ parameter (Eq. 21), is known. This
topic obviously deserves more scrutiny, observations as well as simulations.

Figure 6: Top panel: Logarithmic age distribution functions of the three
populations of pre-main sequence stars identified by Beccari et al. (2017) in
the ONC (blue, green and red histograms; see key). The sum of the three
distributions is given by the black line with triangles. Middle panel:
Logarithmic stellar age distribution function for the ONC-I model (parameters
given in Table 3). The vertical blue, green and red dashed lines indicate the
mean ages of the three populations as given in Kroupa et al. (2018). Bottom
panel: Model age distribution function of the middle panel convolved by a
Gaussian of standard deviation $\sigma_{age}$ accounting for stellar age
uncertainties. Results are presented for $\sigma_{age}=0.1,0.2$ and $0.3$ Myr
(see key). The solid black line with triangles is the ONC total age
distribution from the top panel normalized so as to present the same old-peak
amplitude as the model with $\sigma_{age}=0.2$ Myr

Figure 6 compares the observed age distribution of ONC pre-main sequence stars
with a model with initial gas density $\rho_{0}=10^{3}\,M_{\odot}\cdot
pc^{-3}$ and core formation efficiency per free-fall time $\epsilon_{\rm
ff,co}=0.50$. The top panel shows the age distributions of the three stellar
populations identified by Beccari et al. (2017), using the same color-coding
as in their Fig. 6. Note that old ages are to the left and young ages (i.e.
most recent star formation) are to the right. The black line with triangles is
the total age distribution. Following the formation of the oldest population,
the overall age distribution declines, highlighting a globally decreasing star
formation rate as found for other star-forming regions (e.g. Chamaeleon I,
Luhman, 2007; Belloche et al., 2011b) and as predicted by Parmentier et al.
(2014) (see also Preibisch, 2012). We note that the three populations hardly
show up when considering the total age distribution, an effect partly driven
by stellar age uncertainties (see below).

The middle panel presents the model age distribution for $\epsilon_{\rm
ff,co}=0.5$, $\rho_{0}=10^{3}\,M_{\odot}\cdot pc^{-3}$ ($m_{clump}=3.35\cdot
10^{4}\,M_{\odot}$ and $r_{clump}=2\,$pc), and an age of 3.3 Myr for the onset
of core formation. Model parameters are quoted in Table 3 under ”ONC-I”. The
model age distribution presents two well-distinct peaks, prolonged by a tail
of lower star formation activity. The amplitude of the second peak is lower
than that of the first/oldest one, as required to reproduce an overall
decreasing star formation history. The mean ages of the old, intermediate and
young populations as given by Beccari et al. (2017, their Table 1) are,
respectively, 2.87, 1.88, and 1.24 Myr, with the individual stellar ages taken
from Da Rio et al. (2016). More recent estimates of the three population mean
ages (based on a different set of pre-main-sequence isochrones) are mentioned
in Kroupa et al. (2018, their section 2.1) as 3.0, 1.8 and 0.8 Myr. These
three ages are highlighted as the blue, green, and red vertical dashed lines,
respectively. Note the shift of the youngest population towards even younger
ages (i.e. shift between the red histogram in the top panel and the dashed red
line in the middle panel). The most recent stellar population of the ONC would
thus correspond to the trail of low-level star formation activity in the model
of Fig. 6.

Figure 7: Same as bottom panel of Fig. 6 but for the ONC-II model:
$m_{clump}=10^{4}\,M_{\odot}$, $r_{clump}=1.32$ pc$,\epsilon_{\rm
ff,st}=0.026$ and $\epsilon_{\rm ff,co}=0.30$ (see also Table 3)

The bottom panel shows the impact of stellar age uncertainties. For this, we
convolve the model output with a Gaussian function of standard deviation
$\sigma_{age}$, with $\sigma_{age}$ the uncertainty of individual stellar
ages. The mean difference between the age estimates of the 3 populations as
given by Beccari et al. (2017) and Kroupa et al. (2018) amounts to 0.2 Myr. We
thus adopt $\sigma_{age}=0.2$ Myr. At an age of 1 Myr (3 Myr), this equates
with a logarithmic age uncertainty $\sigma_{log_{10}age}=0.09$ (0.03), in
agreement with the uncertainties quoted in Table 5 of Da Rio et al. (2016).
The convolved age distribution is given in the bottom panel of Fig. 6 as the
solid magenta line. We also show the result for a larger ($\sigma_{age}=0.3$
Myr) or a smaller ($\sigma_{age}=0.1$ Myr) uncertainty (dash-dotted and dashed
purple lines). At this stage, we refrain from a comprehensive browsing of the
parameter space and from quantifying the fit quality with e.g. a
$\chi^{2}$-test because of: (1) the uncertainties of the observed stellar age
distribution at young ages, and (2) the absence of an initial gas density
gradient in the current model. How a density gradient will qualitatively
affect our results is discussed later in this section and in Sec. 7.4.

The relative amplitudes of the peaks of an observed star age distribution are
of limited help to probe the history of a star-forming region. This is because
the peaks are unequally impacted by stellar age uncertainties. For instance,
in our model, the first/oldest peak is the sharpest, thus the one most
curtailed by stellar age uncertainties. A more robust indicator of the
relative strengths of the successive star formation episodes is the stellar
mass formed by each of them. The three populations identified by Beccari et
al. (2017) contain, respectively, 111, 63 and 16 pre-main-sequence
stars333Beccari et al. (2017) mention a census of 24 stars for the ONC
youngest population. However, counting the number of objects in the red
histogram of their Fig. 6 gives 16, and we use this latter number.. This is
equivalent to 58 %, 33 % and 9 % of the ONC stellar mass if the mean stellar
mass is the same for all three populations, and if the observed samples are
representative of their respective populations. The model-predicted fractions
are more even, however, with 40 %, 31.5 % et 28.5 % of the total stellar mass
in the first, second and third million year, respectively. Adding an initial
density gradient to the model clump will contribute to alleviating this (see
Sec. 7.4). The total stellar mass produced by the model is $1400\,M_{\odot}$,
in rough agreement with that observed in the ONC ($1800\,M_{\odot}$,
Hillenbrand & Hartmann, 1998).

Table 3: The ONC models of Sec. 6.2. Common parameters are: $\rho_{0}=10^{3}\,M_{\odot}\cdot pc^{-3}$, $t_{ret,fco}=1\,{\rm Myr}\equiv 3.9\,\tau_{\rm ff,0}$, and an age of 3.3 Myr for core formation onset. Also given are the stellar mass and star formation efficiency at $t=3$ Myr Model |  $m_{clump}$ |  $\epsilon_{\rm ff,co}$ | $\epsilon_{\rm ff,st}$ |  $m_{st,3Myr}$ | $SFE_{gl,3Myr}$
---|---|---|---|---|---
ID |  [$10^{4}M_{\odot}$] |  |  | [$M_{\odot}$] |
ONC-I | $3.3$ | 0.5 | 0.010 | $1400.$ | 0.04
ONC-II | $1.0$ | 0.3 | 0.026 | $1500.$ | 0.15

The gas mass required by the ONC-I model may still seem prohibitive
($m_{clump}=3.35\cdot 10^{4}\,M_{\odot}$). It is higher by a factor $\sim 3$
than any gas mass estimate proposed so far for the ONC progenitor (e.g.
$1.2\cdot 10^{4}\,M_{\odot}$; Kroupa et al., 2001). One parameter which can
still be adjusted to reduce the initial gas mass and whose variations have not
been explored yet is the star formation efficiency per free-fall time. By
increasing $\epsilon_{\rm ff,st}$, one can form a similar stellar mass out of
a smaller amount of gas. Pokhrel et al. (2021) infer, as an average for 12
nearby molecular clouds, $\epsilon_{\rm ff,st}=0.026$. A model with
$\epsilon_{\rm ff,st}=0.026$ and $\epsilon_{\rm ff,co}=0.30$ is shown in Fig.
7 (model ONC-II in Table 3). The clump mass and radius are, respectively,
$m_{clump}=10^{4}\,M_{\odot}$ and $r_{clump}=1.32$ pc (the radius has been
decreased so as to retain the same density $\rho_{0}$ as for the ONC-I model).
The contributions in stars, star-forming core envelopes and failed cores to
the core mass budget are provided in Table 1. The stellar mass formed after 3
Myr is $m_{st}\simeq 1500\,M_{\odot}$, similar to that of the ONC-I model. The
reason as to why a similar stellar mass forms out of a three-times smaller gas
mass is two-fold here. Not only is the star formation efficiency per free-fall
time $\epsilon_{\rm ff,st}$ higher ($\epsilon_{\rm ff,st}=0.026$ vs.
$\epsilon_{\rm ff,st}=0.01$), the core formation efficiency per free-fall time
$\epsilon_{\rm ff,co}$ is also smaller ($\epsilon_{\rm ff,st}=0.30$ vs.
$\epsilon_{\rm ff,st}=0.50$), which has been shown in Fig. 4 to increase the
clump star formation efficiency. Note that this latter effect comes at the
expense of reducing the trough between the old and intermediate-age
populations (compare the magenta and brown solid lines in bottom panel of Fig.
6 and Fig. 7, respectively).

Having used the observed stellar age distribution and stellar mass of the ONC
to constrain the mass, density and core formation efficiency per free-fall
time of its gaseous progenitor, it is tempting to compare our results with
earlier studies. Kroupa et al. (2001) and Farias et al. (2020) study the
kinematics of the ONC stars, including those unbound by gas expulsion and
stellar interactions. The ONC model of Kroupa et al. (2001) builds on an
initial gas mass of $1.1\cdot 10^{4}\,M_{\odot}$ with a half-mass radius of
$r_{0.5}=0.45$ pc, and a star formation efficiency of 1/3. It gives rise to a
bound cluster which retains one third of its initial stellar content in the
aftermath of gas expulsion and survives for about 150 Myr. Our ONC-II model
has a similar initial gas mass, but is more extended, the half-mass radius of
our uniform density clump being $1.32/2^{1/3}\simeq 1.05$ pc. A gas density
gradient will make our model more compact, however. For instance, the same gas
mass ($m_{clump}=10^{4}\,M_{\odot}$) distributed inside the same radius
($r_{clump}=1.32$ pc) but according to a singular isothermal sphere yields an
initial half-mass radius of $r_{clump}/2=0.66$ pc for the gas. The stars will
be even more concentrated (Parmentier & Pfalzner, 2013). Another difference
between ONC-II and Kroupa et al. (2001) is the global star formation
efficiency $SFE_{gl}$. Kroupa et al. (2001) assume $SFE_{gl}=0.33$ at gas
expulsion, while our model predicts $SFE_{gl}=m_{st}/m_{clump}=0.15$ at $t=3$
Myr. Does our low star formation efficiency imply that the ONC-II model
cluster is unable to survive gas expulsion? Not once our model will include a
gas density gradient. Stars will then form more centrally-concentrated than
the embedding gas (Parmentier & Pfalzner, 2013), and the locally higher star
formation efficiency achieved by its inner regions will help the cluster
retain part of its stars despite a globally low star formation efficiency (see
also Adams, 2000). Using $N$-body simulations building on the model of
Parmentier & Pfalzner (2013), Shukirgaliyev et al. (2017) conclude that the
star formation efficiency threshold for cluster survival is $SFE_{gl}=0.15$,
with instantaneous gas expulsion and an external tidal field (taken as that of
the Solar neighborhood). Not only will a gas density gradient help the cluster
survive gas expulsion, it will also drive more efficient star formation
(Parmentier, 2019). This will raise $SFE_{gl}$ above the threshold predicted
by Shukirgaliyev et al. (2017), and reduce the initial gas mass $m_{clump}$
accordingly (since the ONC stellar mass is imposed).

Farias et al. (2020) use their gradual cluster formation model (Farias et al.,
2019) to model the velocity distribution function of the ONC stars, including
run-aways and walk-aways. They conclude that slow star formation
($\epsilon_{\rm ff,st}=0.01{\rm-}0.03$, rather than $\epsilon_{\rm ff,st}\geq
0.1$) and a progenitor clump surface density $\Sigma_{clump}\simeq
7000\,M_{\odot}\cdot pc^{-2}$ (rather than $\Sigma_{clump}\simeq
700\,M_{\odot}\cdot pc^{-2}$) provide the best match to the observations
(their Fig. 10). The surface density of our ONC-II clump
$\Sigma_{clump}=10^{4}/(\pi 1.32^{2})\simeq 1800\,M_{\odot}\cdot pc^{-2}$ is 4
times lower than what they advocate. Both models are not entirely comparable,
however. While our star formation histories are intrinsically time-varying,
Farias et al. (2019) assume a constant star formation rate fixed by the
initial gas mass and free-fall time (their Eq. 7).

Although, in principle, the comparison of different modeling approaches should
help cast more light on the formation and evolution of star clusters, this is
no straightforward task, and the variety of model assumptions calls for a
careful comparison.

###  6.3 Do observed single-population clusters hide multiple star formation
episodes ?

Figure 8: Intrinsic (solid lines) and measured (dashed lines) stellar age
distributions for a cluster-forming clump with
$\rho_{0}=10^{4}\,M_{\odot}\cdot pc^{-3}$. The blue and magenta lines
correspond to $\epsilon_{\rm ff,co}=0.08$ and $\epsilon_{\rm ff,co}=0.5$,
respectively (see key). The red long-dashed line is the measured age
distribution for $\epsilon_{\rm ff,co}=0.5$ raised to the same amplitude as
for $\epsilon_{\rm ff,co}=0.08$: both measured age distributions are
undistinguishable in shape

The frequency at which cluster-forming clumps pile up successive star
formation episodes now depends on the clump density $\rho_{0}$ (Eqs 20 and
21). An interesting consequence is that the star formation histories of
starburst clusters could have close-by peaks that, in the observed age
distribution, get blended due to the short age-span separating them. To
illustrate this possibility, let us consider two clumps with the same mass
$m_{clump}=3.3\cdot 10^{5}\,M_{\odot}$, radius $r_{clump}=2$ pc (i.e. clump
density $\rho_{0}=10^{4}\,M_{\odot}\cdot pc^{-3}$) and star formation
efficiency per free-fall time $\epsilon_{\rm ff,st}=0.026$, but distinct core
formation efficiencies per free-fall time. With $\epsilon_{\rm ff,co}=0.08$,
the first clump form no failed cores ($F_{fco}=0.0$; see Table 1 for
$\epsilon_{\rm ff,st}=0.026$). In contrast, with $\epsilon_{\rm ff,co}=0.5$,
the second clump has a failed core mass fraction $F_{fco}=0.84$ (see Table 1).
Their respective intrinsic stellar age distributions (Fig. 8) reflect this
difference: a single spike when $\epsilon_{\rm ff,co}=0.08$, but 2-to-3 peaks
for $\epsilon_{\rm ff,co}=0.50$. Note also the lower amplitude of the stellar
age distribution for $\epsilon_{\rm ff,co}=0.50$, which reflects the lower
star formation efficiency induced by higher $\epsilon_{\rm ff,co}$, as already
noticed in Fig. 4.

The measured age distributions, that is, the intrinsic distributions convolved
with a Gaussian error, tell a different story, however. We adopt
$\sigma_{age}=0.2$ Myr, which is slightly smaller than what Kudryavtseva et
al. (2012) inferred for the Galactic starburst cluster Westerlund 1
($\sigma_{age}=0.25$ Myr). Given the very short clump free-fall time
$\tau_{\rm ff,0}=0.08$ Myr, the age span between two successive peaks is
$\simeq 0.3$ Myr ($\top=3.9$), and the multi-peak aspect shown by the
$\epsilon_{\rm ff,co}=0.50$ model is erased once stellar age uncertainties are
accounted for (see magenta short-dashed line in Fig. 8). Figure 8 also shows
as the red long-dashed line the measured age distribution for $\epsilon_{\rm
ff,co}=0.50$ normalized to the same maximum as for $\epsilon_{\rm
ff,co}=0.08$. The shapes of both measured age distributions are hardly
distinguishable, despite containing fundamentally different star formation
histories.

Starburst clusters in the Galactic disk - like Westerlund 1 - are observed to
be single-burst clusters (Kudryavtseva et al., 2012). The exercise above shows
that the actual underlying star formation history may be quite different,
however. We will revisit this issue in the future, once an initial gas density
gradient will have been integrated into our model. The reader may indeed ask
as to why the full-width at half-maximum of our convolved age distributions in
Fig. 8 is as large as $\simeq 1.5$ Myr, that is, significantly larger than the
observed full-width at half-maximum of $0.4$ Myr inferred by Kudryavtseva et
al. (2012) for Westerlund 1. This effect stems from the absence of an initial
gas density gradient. A density gradient $\rho_{gas}\propto r^{-2}$ will yield
a full-width at half-maximum reduced by a factor 3-to-4 compared to a uniform
density clump (see Fig. 4 in Parmentier, 2019), thereby bringing model and
observed full-widths at half-maximum in agreement.

##  7 Discussion

###  7.1 Diversity in star-formation histories

Our model is able to generate very diverse star formation histories, from a
star formation rate steadily decreasing with time to ’jagged’ star formation
histories. Additionally, the time-span between two successive peaks can vary
from clump to clump, since it is a function of the initial gas density (Eq.
20).

Figure 5 seems to suggest that more numerous star formation episodes arise in
denser and/or more massive clumps (in each panel, compare red solid and
magenta dashed lines). This conclusion is erroneous, however. In each panel of
Fig. 5, the varying parameter is the clump mass, while the clump radius and
the physical return time of the failed cores are kept constant. Consequences
are (1) a modification of the clump density $\rho_{0}$, hence of the clump
free-fall time $\tau_{\rm ff,0}$ and, therefore, (2) a modification of the
return time of the failed cores in units of the clump free-fall time
$t_{ret,fco}/\tau_{\rm ff,0}$. The seemingly greater ability of more massive
clumps to carve more peaks in their star formation history in Fig. 5 stems
from their failed core return time $t_{ret,fco}/\tau_{\rm ff,0}$ becoming
longer by virtue of an increasingly shorter free-fall time $\tau_{\rm ff,0}$
coupled to a constant physical return time $t_{ret,fco}$. A longer return time
of the failed cores $t_{ret,fco}/\tau_{\rm ff,0}$ allows a greater decrease of
the gas density by the time failed cores disperse, yielding thereby a greater
relative increase of the gas density and of the core formation rate at
$t\gtrsim t_{ret,fco}$. In other words, to give the first core formation
episode more time to subside increases the likelihood of the second core
formation episode to spike (see also Sec. 2.2 and Fig. 3). Contrary to what a
superficial glimpse at Fig. 5 may suggest, the initial gas mass has no
influence on the shape of the core and star formation histories.

Figure 9: Evolution of the core formation rate normalized with the initial
free-fall time $\tau_{\rm ff,0}$ and initial mass $m_{clump}$ of the clump
gas, as a function of time after core formation onset in units of the clump
free-fall time. Each line with given color and type corresponds to four model
clumps with different masses $m_{clump}$ and densities $\rho_{0}$ (see text
for details).

To illustrate this even more convincingly, let us consider a clump of mass
$m_{clump}=3.3\cdot 10^{5}\,M_{\odot}$ and radius $r_{clump}=2$ pc (density
$\rho_{0}=10^{4}\,M_{\odot}\cdot pc^{-3}$). We also consider three clumps with
a 10-times higher mass and either (i) the same radius ($m_{clump}=3.3\cdot
10^{6}\,M_{\odot}$, $r_{clump}=2$ pc, $\rho_{0}=10^{5}\,M_{\odot}\cdot
pc^{-3}$), or (ii) the same volume density ($m_{clump}=3.3\cdot
10^{6}\,M_{\odot}$, $r_{clump}=4.3$ pc, $\rho_{0}=10^{4}\,M_{\odot}\cdot
pc^{-3}$), or (iii) the same surface density ($m_{clump}=3.3\cdot
10^{6}\,M_{\odot}$, $r_{clump}=6.3$ pc, $\rho_{0}=0.3\cdot
10^{4}\,M_{\odot}\cdot pc^{-3}$).

The mass $m_{clump}$ and the density $\rho_{0}$ define the dimensional
parameters of the model. For each pair ($m_{clump}$, $\rho_{0}$), we vary the
dimensionless model parameters, namely: the core formation efficiency per
free-fall time $\epsilon_{\rm ff,co}$, the return time of the failed cores in
units of the clump free-fall time $t_{ret,fco}/\tau_{\rm ff,0}$, and the star
formation efficiency per free-fall time $\epsilon_{\rm ff,st}$. We test the
following values: $\epsilon_{\rm ff,co}=(0.3,0.5)$, $t_{ret,fco}/\tau_{\rm
ff,0}=(3.3,10.0)$, and $\epsilon_{\rm ff,st}=(0.0026,0.026)$. Twenty-four core
formation histories are presented in Fig. 9 (top panel: $\epsilon_{\rm
ff,co}=0.30$; bottom panel: $\epsilon_{\rm ff,co}=0.50$). Time on the $x$-axis
is in units of the clump free-fall time, and the core formation rate on the
$y$-axis is normalized by the initial mass and free-fall time of the gas.
Models with a given set of dimensionless parameters (i.e. lines with given
color and dashed type in a given panel) display identical normalized core
formation histories. For instance, the 4 red solid lines, each depicting a
($m_{clump}$, $\rho_{0}$) combination, coincide with each other. The
normalized core formation history $\tau_{\rm ff,0}CoFR/m_{clump}(t/\tau_{\rm
ff,0})$ is therefore independent of the clump mass $m_{clump}$ and density
$\rho_{0}$. It is shaped by the parameter triple $(\epsilon_{\rm
ff,co},t_{ret,fco}/\tau_{\rm ff,0},\epsilon_{\rm ff,st})$444For the sake of
completeness, we should also mention the return time of the envelopes of star-
forming cores $t_{ret,colco}/\tau_{\rm ff,0}$. However, given that star-
forming cores represent a minor fraction of the total core mass in models
generating successive star formation episodes, this fourth dimensionless
parameter is here of minor importance. , while the clump mass and density
(equivalently, the mass and free-fall time) act as scaling factors. A larger
clump mass raises the core formation history without modifying its overall
shape. A shorter free-fall time (i.e. denser clump) also provides a scaled-up
version of the core formation history, while simultaneously increasing its
”frequency”.

For the parameters tested, Fig. 9 shows that we can expect up to 5-6 star
formation episodes (for $\epsilon_{\rm ff,co}=0.50$ and $t_{ret,fco}/\tau_{\rm
ff,0}=10$). The choice of $t_{ret,fco}/\tau_{\rm ff,0}=10$ for the return time
of the failed cores may seem extreme. It requests $\top\simeq 10$, which goes
beyond the range of values advocated in Sec. 6.1 (see Eqs 20-21). Yet,
$t_{ret,fco}/\tau_{\rm ff,0}=10$ may be relevant for environments where the
external pressure is high enough to stabilize the cores for a longer time-span
than when it is negligible (Lada et al., 2008).

Although the numbers of star formation episodes in Fig. 9 are comparable to
the numbers of stellar populations identified in old globular clusters (often
$\geq$3-4 subpopulations; e.g. Carretta, 2015; Bastian & Lardo, 2018), we note
that the relation between the two is not that straightforward (see Sec.
7.2)555An even greater number of stellar populations - 16 - have been
extracted from the chromosome map of $\omega$ Cen, the most massive Galactic
globular cluster (Milone et al., 2017). $\omega$ Cen, however, is likely the
nucleus of a now disrupted dwarf galaxy (e.g. Hilker & Richtler, 2000) and its
early evolution cannot thus be explained with a cluster-based only scenario. .

The conclusion that the clump mass does not shape the core formation history
seems at odd with what we know about old globular clusters, namely, that more
massive globulars tend to present stronger multiple stellar population
patterns. For instance, the mass of Galactic globular clusters correlates with
their maximum helium abundance variation (Milone et al., 2018), with their
maximum nitrogen abundance variation (Lagioia et al., 2019), and with the
number fraction of stars enriched in hot hydrogen burning elements (Milone et
al., 2017). Here, it is crucial to stress that the multiple stellar
populations (generations if they fit a time sequence) of old globular clusters
are identified based on their chemical abundances. In contrast, the successive
star formation episodes of our models are identified based on the star
formation history and the failed core return time $t_{ret,fco}$.

The detailed star formation histories of old globular clusters have remained
unknown so far, and our model does not compound the star formation history of
a cluster with its chemical evolution. Yet, could our approach help shed light
on as to why more massive globular clusters present more prominent multiple-
stellar-population patterns?

###  7.2 Potential link with multiple stellar populations in old globular
clusters

Patterns of multiple stellar populations vary greatly among Galactic old
globular clusters. Specifically, the morphology of their chromosome
map666Chromosome maps trace multiple stellar populations based on their star
$UV$-photometry. They build on the property that ultraviolet colors are very
sensitive to star-to-star helium and light-element abundances (Milone et al.,
2015) is highly diverse - with no one cluster being identical to another. For
instance, the branch of the so-called second generation is sometimes more
clumpy (e.g. NGC2808), sometimes more continuous (e.g. NGC6652; see Fig. 4 in
Milone et al., 2017). These maps also reveal varying numbers of globular
cluster subpopulations. Similar conclusions can be drawn from the sodium-
oxygen anticorrelation (see e.g. Fig. 1 in Charbonnel, 2016). It is tempting
to connect this aspect to our model star formation histories, whose aspect is
sometimes ’wavy’, with well-defined peaks and troughs (when $\epsilon_{\rm
ff,co}$ and $t_{retfco}/\tau_{\rm ff,0}$ are large), sometimes smoothly
decreasing with time. But first, we have to cope with the apparent
contradiction noted in Sec. 7.1, that high-mass globular clusters present the
most striking patterns of multiple stellar populations and of enrichment in
hot hydrogen burning elements, while the shape of our model star formation
histories are independent of the progenitor clump mass (Fig. 9). Can we get
out of this predicament?

Yes, if we assume that the stellar mass built by the clump must reach a given
threshold $m_{th}$ to trigger the formation of the polluters providing the hot
hydrogen burning elements. This could be the case if the polluters are
supermassive stars, namely, stars with masses a few $10^{3}$ \- a few
$10^{4}\,M_{\odot}$ (Denissenkov & Hartwick, 2014). Their predicted yields are
in excellent agreement with the light-element abundance patterns observed in
globular clusters (e.g. the sodium-oxygen anti-correlation: compare Fig. 1 in
Charbonnel (2016) with the top-right panel of Fig. 1 in Denissenkov & Hartwick
(2014)). In addition, supermassive stars could pollute their environment
during a few million years (Denissenkov & Hartwick, 2014), similar to the
time-span for cluster development. Their formation mechanism would be via
stellar collisions, which requires either a high cluster density (e.g.
Portegies Zwart et al., 2004), or a high number of stars (at least a few
$10^{6}$ stars; Gieles et al., 2018). Here, we build on the scenario that
their formation requires a minimum number of stars, hence a minimum stellar
mass $m_{th}$.

Let us consider three gas clumps (1,2,3) of mass $m_{1}<m_{2}<m_{3}$. The
least massive clump is also less massive than the threshold mass $m_{th}$.
With $m_{1}<m_{th}$, clump 1 is unable to form a supermassive star, and the
corresponding cluster does not display the light-element abundance variations
typical of globular clusters. Yet, if the core formation efficiency
$\epsilon_{\rm ff,co}$ and return time $t_{ff,fco}/\tau_{\rm ff,0}$ of clump 1
are large enough, its star formation history does display successive star
formation episodes. This may explain why open clusters do not present light-
element abundance variations (MacLean et al., 2015), despite the fact that
some (future) open clusters present multiple star formation episodes (e.g. the
ONC, Beccari et al., 2017).

We now consider clump 2 and clump 3 as massive enough to form a supermassive
star, i.e. $m_{th}<<m_{2}<m_{3}$. We also assume that their mass is the only
parameter differentiating them (i.e. their free-fall times are equal, as well
as their dimensionless parameters). With $SFE_{2,th}$ and $SFE_{3,th}$ the
star formation efficiencies of clump 2 and clump 3 at the time they reach the
stellar mass threshold $m_{th}$, we have $SFE_{3,th}<SFE_{2,th}$ since
$m_{th}=SFE_{3,th}m_{3}=SFE_{2,th}m_{2}$. Because both clumps have identical
free-fall times and dimensionless parameters, the growths with time of their
respective star formation efficiencies are alike, $SFE_{2}(t)=SFE_{3}(t)$.
Therefore, the inequality $SFE_{3,th}<SFE_{2,th}$ implies that clump 3 reaches
the threshold $m_{th}$ at an earlier time than clump 2. This shift in time
brings forward the formation of a supermassive star in clump 3, thereby
providing a longer pollution time-span. As a result, the more massive clump
should form a larger fraction of polluted stars, and achieve more extreme
light-element abundance variations, as observed (Milone et al., 2017, 2018;
Lagioia et al., 2019).

As long as the supermassive star does not ”operate”, newly-formed stars
constitute the - chemically-defined - first stellar generation. We emphasize
that the end of this generation coincides with the onset of the supermassive
star pollution, not with the end of the first star formation episode of our
model. If the intra-cluster medium pollution starts during the third star
formation episode, the first, second and beginning of third star formation
episodes define the first stellar generation. Conversely, one could consider
the case of a protoglobular clump with $\epsilon_{\rm ff,co}<0.1$, but massive
enough to form a supermassive star. Its star formation history would be
peakless (small core formation efficiency per free-fall time, hence small
return fraction $f_{ret,fco}$). Yet, the chromosome map would present two
groups of stars corresponding to the stars born before and after the onset of
the supermassive star pollution. In other words, the shape of the star
formation history and the chemical evolution are decoupled from each other.
While the shape of the star formation history depends on the model
dimensionless parameters, the chemical evolution depends on a dimensional
parameter (in the present case, the threshold stellar mass $m_{th}$). This is
the very reason as to why, all through this paper, we have refrained from
using the term ”stellar generations” to describe the successive star formation
episodes of our model, this term being often used in studies of globular
cluster multiple stellar populations. A consistent comparison between the
number of star formation episodes predicted by our model and the number of
multiple stellar generations in old globular clusters calls for a comparison
between the return time of the failed cores with the chemical enrichment time-
scale, which is beyond the scope of the present study.

Finally, we stress that the early evolution of globular clusters may see
successive pollution episodes, driven by different polluters operating on
different time-scales (e.g. D’Antona et al., 2016). Despite its low-mass
($2.2\cdot 10^{4}\,M_{\odot}$, Baumgardt et al., 2019), the globular cluster
NGC6535 hosts two chemically-defined stellar populations (Milone et al.,
2017). How could such a low-mass globular cluster contain stars bearing the
imprint of hot hydrogen burning ? Either the pollution of protoglobular clumps
in hot hydrogen burning products is not caused by supermassive stars only. Or
the current mass of NGC6535 is a tiny residual of a much larger initial mass.
In either case, it is interesting to note that NGC6535 hosts two chemically-
defined populations, while the much less massive ONC presents three star
formation episodes. One more evidence that peaks in star formation histories
are not direct matches to globular cluster multiple stellar generations.

###  7.3 Caveats and Limitations

Our idealized model does obviously not capture all aspects of the star
formation process. In this section we assess its limitations and how more
realistic conditions could affect its outcomes. We start by discussing our
assumption of an isolated clump in near-equilibrium (i.e. constant clump mass
and radius).

The constancy of the clump radius constitutes a model central hypothesis since
the gas density and core formation rate depend sensitively on it. Under what
conditions can a cluster-forming clump remain in near-equilibrium once star
formation has started? A key ingredient here is the feedback from forming
stars, which stir their surrounding gas. Protostellar winds given off by low-
mass stars were already suggested by Norman & Silk (1980) to be prime drivers
of turbulence within molecular clouds. On the pc-scale of cluster-forming
clumps, Nakamura & Li (2007) reach the same conclusion by means of
hydrodynamical simulations. They show that even in case of fast turbulence
decay, turbulence gets replenished from inside the clump due to the gas
motions generated by protostellar outflows. This happens because protostellar
outflows are collimated (i.e bipolar outflows; see Bally, 2016, for a review),
allowing them to transport energy and momentum to greater distances from their
sources than if they were spherically-symmetric (Nakamura & Li, 2007). This
”protostellar outflow-driven turbulence” allows the model clump of Nakamura &
Li (2007) to remain close to dynamical equilibrium provided that the star
formation efficiency per free-fall time amounts to a few per cent. There is in
fact an interplay between the rates of turbulence decay and star formation.
Matzner (2001) shows that a near-equilibrium configuration arises when
turbulence decays slowly (over 10 free-fall times). When turbulence decays
fast (over one free-fall time), however, continuous star formation is required
for protostellar outflows to offset turbulence decay. If star formation
happens discretely (i.e. one star at a time due to a low star formation rate),
turbulence decay cannot be offset and the system collapses777An example of a
cluster-forming clump whose shrinking is unimpeded by protostellar outflows is
provided in Huff & Stahler (2006).. The star formation efficiency per free-
fall time of our models ($\epsilon_{\rm ff,st}=1$-$2.6$ %) is close to the
limit for turbulence replenishment inferred by Nakamura & Li (2007) and,
therefore, our assumption of a constant clump radius appears reasonable.
Additionally, we note that even if our model had to include the contraction of
a clump, this contraction is more likely to happen when the star formation
rate reaches a low point, that is, one return time after a peak in the star
formation history. Such a clump contraction would be timely to stimulate even
more the rise of the next star formation episode by increasing the clump gas
density.

One could now ask what the fate of these protostellar outflows is. Do they
leave the cluster-forming clump (e.g. Farias et al., 2019), thereby decreasing
the gas mass available to star formation, or do they remain trapped in the
clump (our model). That is, how valid is the assumption of an isolated clump?
This is a complex issue since it involves (1) the outflow speed and the clump
escape velocity (hence the clump mass, radius, density profile, and the
location of the source star inside the clump), and (2) the outflow morphology
(which often consists of an inner high-velocity jet, superimposed to a slower
wider-angle component, Bally, 2016). The escaping outflows can also sweep part
of the clump gas on their way out, thereby increasing the gas mass fraction
lost from the clump. Federrath et al. (2014) find that more than a quarter of
the clump mass can be ejected this way. On the one hand, the star formation
rate would subside more quickly than predicted by our model, as a result of
the corresponding gas mass decrease. On the other hand, additional star
formation may be triggered since jets and outflows also drive local turbulent
gas compressions (Federrath et al., 2014). Should the formation of failed
cores be enhanced too, the multiple stellar population patterns predicted by
our model will be strengthened. We also note that, when the leftover star-
forming gas is completely expelled out of protoclusters due to massive star
feedback, star formation may be triggered in the ensuing shell of expanding
gas. This will add to the phenomenon of star cluster multiple stellar
populations if the stellar component of the shell collapses back (Parmentier,
2004; Rahner et al., 2018).

Another way of modifying the mass enclosed within the adopted clump radius is
by gas inflows. Myers (2009) find that all nine young stellar groups within
300 pc of the Sun radiate multiple filaments of parsec length, designing
thereby a hub-filament structure. Gas flowing inwards along these filaments
can fuel the cluster star formation process. This may be the case for the
Serpens South embedded cluster (Kirk et al., 2013) and for the ONC (Hacar et
al., 2017a). In the later case, Hacar et al. (2017a) estimate that the
integral-shaped filament currently supplies the ONC with gas at a rate of
$\dot{M}_{acc}\simeq 385\,M_{\odot}\cdot Myr^{-1}$. This accretion rate cannot
impact significantly the clump gas reservoir of our ONC-II model whose initial
gas mass is $10^{4}\,M_{\odot}$. This conclusion, however, would be different
if the gas accretion rate was higher in the past (Kroupa et al., 2018). To
assess the importance of filamentary accretion onto forming clusters cannot be
performed without knowledge of the gas accretion rate history. In the case of
our model, the effect would be to fill in the troughs in the time evolution of
the gas mass and, therefore, to reduce the amplitude of the oscillatory
patterns seen in some of our simulations, which we now turn to.

Our model star formation histories can present clear-cut oscillations, which
happens when failed cores dominate the core mass budget. These oscillations
stem from the adoption, for each model clump, of a unique return time for all
failed cores, i.e. $t_{ret,fco}=3.3\top\tau_{\rm ff,0}/\sqrt{\mathcal{C}}$
(with $3.3\top$ the return time of the failed cores in unit of their free-fall
time, and $\mathcal{C}$ the core-to-clump mean density ratio; see Eqs 20, 19,
18). The amplitude of the oscillations as predicted by our model therefore
constitutes an upper limit: in case of spatial variations in the core return
time, the gas replenishment from simultaneously-formed cores is spread more
evenly in time, smoothing thereby model oscillations. The effect would be
especially significant for cores formed at a peak of core formation activity,
since these cores bear the subsequent core formation history ”upheaval”.
Although we have assumed that the clump gas remains uniform in density (albeit
declining with time) through the clump evolution, it is not the case in
practice. Hydrodynamical simulations of turbulent self-gravitating clouds
produce filaments (e.g. Federrath, 2016; Smith et al., 2016), implying that
the immediate environment of cores varies from one core to another, an effect
also reinforced by spatial variations of the ram pressure arising from
protostellar outflows (Norman & Silk, 1980). With re-expanding cores confined
by different external pressures, the $\top$ parameter and, therefore, the
return time $t_{ret,fco}$ must vary from core to core. The mean density
contrast $\mathcal{C}$ of failed cores can also vary (although remaining
markedly smaller than that of star-forming ones). More realistic simulations
may thus show that the formation of clearly distinct star formation episodes
is the exception rather than the rule, even when failed cores dominate the
core mass budget by far. Interestingly, inspection of figures 3-7 in Milone et
al. (2017) suggests that globular clusters whose enriched population draws a
smooth distribution of points in the chromosome map outnumber those with
clearly clumpy distributions (although photometric uncertainties contribute to
this as well).

Another aspect that our present model does not account for is that of core
coalescence (e.g. Norman & Silk, 1980; Dib et al., 2010). The corresponding
mass increase may allow some of them to reach their Jeans mass, thereby
turning failed cores into star-forming ones. This too would reduce the
amplitude of the star formation history oscillations.

Coalescence and return-time variations of failed cores contribute to
explaining why the strong oscillatory patterns observed in some of our models
do not show up in hydrodynamical simulations888At least we are not aware of
such simulations. Unfortunately, such publications do not always present their
star or core formation history.: these are characterized by many more degrees
of freedom than our model. But we also stress that such simulations need to
cover at least 4-6 free-fall times of the modeled star-forming region for two
subsequent oscillations to be uncovered.

Finally, a note about magnetic fields. The higher their initial strength in
the cluster-forming clump, the lower the core formation efficiency per free-
fall time (Dib et al., 2010). Magnetic field differences among cluster-forming
clumps may thus contribute to whether or not subsequent star formation
episodes emerge in clusters.

###  7.4 Impact of an initial gas density gradient

In this contribution, we have restricted our attention to uniform density
clumps. Yet, observed cluster-forming clumps are centrally-concentrated (e.g.
Mueller et al., 2002). Compared to a top-hat profile, a gas density gradient
strengthens the star formation activity of a clump by locating a greater
fraction of its gas mass in its central regions, where the gas density is
higher and the free-fall time shorter. To quantify this effect, Parmentier
(2019) introduces the magnification factor $\zeta$, defined as the ratio
between the star formation rate of a centrally-concentrated clump and that of
its top-hat equivalent999The top-hat equivalent is a clump of uniform gas
density containing the same gas mass as the centrally-concentrated clump,
enclosed within the same radius. (her equations 8-9; see also Tan et al.,
2006). For instance, a gas density profile with a logarithmic slope of $-2$
(Mueller et al., 2002) and a central plateau of extent $10^{-2}$ the clump
radius yields $\zeta\simeq 2.5$ (see top panel in Fig. 1 of Parmentier &
Pasquali, 2020). The higher the gas density contrast between the clump center
and the clump edge, the higher the magnification factor is, and the more
invigorated star formation activity is.

In the framework of the present model, how would a density gradient affect our
results? The star formation rate of a centrally-concentrated clump is given by
Eq. 8 in Parmentier (2019), where $\epsilon_{\rm ff,int}$ is the intrinsic
star formation efficiency per free-fall time, namely, the efficiency of the
clump top-hat equivalent. The same equation can be rewritten as a core
formation rate, by substituting the core formation efficiency per free-fall
time, $\epsilon_{\rm ff,co}$, for the star formation efficiency per free-fall
time $\epsilon_{\rm ff,int}$:

$CoFR_{clump}=\zeta\epsilon_{\rm ff,co}\frac{m_{gas}}{\tau_{\rm ff}}\;.$ (22)

In this equation, the magnification factor $\zeta$ can be associated either to
the intrinsic core formation efficiency per free-fall time $\epsilon_{\rm
ff,co}$

$CoFR_{clump}=(\zeta\epsilon_{\rm ff,co})\frac{m_{gas}}{\tau_{\rm ff}}\,,$
(23)

or to the clump gas mass

$CoFR_{clump}\propto\epsilon_{\rm
ff,co}\left(\frac{\zeta^{2/3}m_{gas}}{r_{clump}}\right)^{3/2}\;.$ (24)

In the first case, the magnification factor raises the efficiency per free-
fall time, implying that an initial gas density gradient strengthens the
emergence of successive star formation episodes. In the second case, a given
core formation rate can be achieved with a gas mass reduced by a factor
$\zeta^{2/3}$. Equivalently, the clump star formation efficiency $SFE_{gl}$
increases by a factor $\zeta^{2/3}$. For instance, the star formation
efficiency of the ONC-II model ($SFE_{gl}=0.15$, Sec. 6.2, Table 3) will be
increased by a factor $1.6{\rm-}2.1$ if $\zeta=2{\rm-}3$.

We also expect a gas density gradient to modify the distribution of the
stellar mass among the successive star formation episodes. The ”weight” of the
first/oldest star formation episode will increase at the expense of the
subsequent ones, since the gas density contrast from clump center to clump
edge plummets with time, causing in turn a decrease of the magnification
factor (see Figs 2, 5 and 6 in Parmentier, 2019).

We will investigate these aspects in future papers.

##  8 Conclusions

We have proposed a novel approach to understand the multiple stellar
populations of star clusters. Our model hinges on the presence of failed cores
in star-forming regions. Failed cores are density peaks of the interstellar
medium which do not end collapsing and forming stars. Instead, they re-expand
and disperse back in the inter-core gas. Although their existence has been
known for more than two decades, both in observations (Taylor et al., 1996;
Lada et al., 2008; Belloche et al., 2011a, b) and simulations (Clark &
Bonnell, 2005; Vázquez-Semadeni et al., 2005; Gómez et al., 2007), they have
not received the same scrutiny as their star-forming counterparts. The
formation of failed cores necessarily depletes the gas reservoir of a star-
forming region. We therefore propose that their formation prompt the gas
density, core and star formation rates to dwindle. Their subsequent
dispersion, however, allows the density of the surrounding gas to re-increase,
thereby rejuvenating the core and star formation rates. When failed cores
dominate the total core mass, this controlling mechanism of the available gas
mass and density yields peaks and troughs in the time evolution of the gas
density, hence of the core and star formation histories. The time-span on
which failed cores disperse - which we coin their return time - determines the
age-span separating the successive peaks of a star formation history.

Multiple peaks in the stellar age distribution of cluster-forming clumps (e.g.
the ONC; Beccari et al., 2017) imply that a star formation episode must at
some point recede, to allow the next one to rise. To implement the scenario
described above, we therefore build on the volume-density driven cluster-
formation model of Parmentier & Pfalzner (2013), which predicts a declining
star formation activity in cluster-forming clumps under the assumptions of
constant clump mass and size.

We start in Sec. 2 with a preliminary model in which star-forming and failed
cores release their gas on the same time-scale. We find that three
dimensionless parameters determine whether a core formation history unfolds at
least two peaks: the core formation efficiency per free-fall time
$\epsilon_{\rm ff,co}$, the fraction $f_{return}$ of the mass in cores that
returns to the clump gas, and the associated time-span $t_{return}/\tau_{\rm
ff,0}$, this one being expressed in units of the clump free-fall time
$\tau_{\rm ff,0}$. High core formation efficiencies per free-fall time
$\epsilon_{\rm ff,co}$, high return fractions $f_{return}$, and long return
times $t_{return}/\tau_{\rm ff,0}$ favor the rise of a second core formation
episode (see red-to-yellow areas in Fig. 3).

Following a review of the evidence for failed cores in observed and simulated
star-forming regions in Sec. 3, we start to augment our model in Sec. 4.
First, we introduce the star formation efficiency per free-fall time
$\epsilon_{\rm ff,st}$, whose combination with the core formation efficiency
per free-fall time $\epsilon_{\rm ff,co}$ defines the respective fractions of
the core mass turned into stars ($F_{*}$), given off to the clump gas as star-
forming core envelopes ($F_{env}$), and as failed/rebounding cores ($F_{fco}$)
(see Eqs 13-15, and Table 1). Secondly, star-forming cores lose their envelope
on a shorter time-span than failed cores disperse. We thus lift our assumption
of identical return times for the failed cores and the star-forming core
envelopes, and we assign them different time-scales (Eqs 11-12). We assume
that failed cores disperse on a time-span $t_{ret,fco}=1$ Myr, which is the
estimated life-expectancy of failed cores in the Pipe molecular cloud (Lada et
al., 2008), and the core re-expansion time in the magnetohydrodynamical
simulations of Vázquez-Semadeni et al. (2005). This is also the age-span
separating the successive populations of pre-main sequence stars in the ONC
(Beccari et al., 2017),

Next, in Sec. 5, we generate stellar age distributions for various clump
masses and core formation efficiencies per free-fall time (Fig. 5). All model
clumps are given the same star formation efficiency per free-fall time
$\epsilon_{\rm ff,st}$, radius $r_{clump}$, and physical return time
$t_{ret,fco}=1$ Myr. Based on the number of peaks being carved in the star
formation history with the imposed 1-Myr-periodicity and the corresponding
stellar mass, we estimate that the ONC gaseous progenitor had a mean gas
density $\rho_{0}=10^{3}\,M_{\odot}\cdot pc^{-3}$ and a core formation
efficiency per free-fall time $\epsilon_{\rm ff,co}\geq 0.3$. It bears
repeating that this is the long return time of the failed cores
$t_{ret,fco}/\tau_{\rm ff,0}$ which drives the ’jagged’ structure of the star
formation history of the most massive clumps, not the clump mass itself.

In Sec. 6, we lift the assumption of a constant physical return time
$t_{ret,fco}$ for the failed cores. Using the analytical estimate of Vázquez-
Semadeni et al. (2005, see left panel of their Fig. 1), we infer that the re-
expanding time of failed cores amounts to a few free-fall times of their host
clump (Eq. 21, with the free-fall time defined as in Eq. 2). We also propose
two models for the star formation history of the ONC (Figs 6-7 and Table 3),
and emphasize that the star formation history of starburst clusters may
contain several star formation episodes, these remaining concealed (i.e.
blended) due to their inter-peak age-span being comparable to stellar age
uncertainties.

In Sec. 7, we discuss various aspects of our model. First, we emphasize the
great diversity in the core and star formation histories we have generated.
The shapes of the core and star formation histories depend on the
dimensionless model parameters ($t_{ret,fco}/\tau_{\rm ff,0}$, $\epsilon_{\rm
ff,co}$, $\epsilon_{\rm ff,st}$). That is, the evolution of the normalized
core formation rate $\tau_{\rm ff,0}CoFR/m_{clump}$ as a function of time
$t/\tau_{\rm ff,0}$ depends exclusively on ($t_{ret,fco}/\tau_{\rm ff,0}$,
$\epsilon_{\rm ff,co}$, $\epsilon_{\rm ff,st}$) (Fig. 9). The dimensional
parameters ($m_{clump}$, $\rho_{0}$) (i.e. clump mass and density) act as
normalizing factors. Second, we discuss the potential relation between our
model and the multiple stellar populations/generations ubiquitously observed
in Galactic old globular clusters. That more massive globular clusters present
stronger patterns of multiple stellar populations seems, at first sight, to
contradict our finding of mass-independent star formation histories. Our star
formation episodes, however, are defined based on the failed core return time,
while the multiple stellar generations of old globular clusters are mapped
based on their light-element abundances (e.g. chromosome map, Milone et al.
(2015), and sodium-oxygen anti-correlation, Carretta et al. (2010)). We
qualitatively show that our model can produce massive clusters with a greater
fraction of enriched stars if the formation of the stars polluting the cluster
with hot hydrogen burning elements requires a threshold stellar mass. Third,
we discuss the caveats and limitations of our model, emphasizing that more
realistic conditions will smooth the strong oscillatory patterns carved in
some of the predicted star formation histories. Finally, we discuss how the
inclusion of a density gradient in our model clumps will impact our present
results. We will chart in greater depth these additional aspects, along with
that of multiple stellar populations, in future papers.

GP acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation – Project-ID 138713538 – SFB 881 (”The Milky Way System”,
subproject B05). This research has made use of NASA ’s Astrophysics Data
System. We thank the anonymous referee for a constructive report.

## References

  * Adams (2000) Adams, F. C. 2000, ApJ, 542, 964. doi:10.1086/317052
  * Alves et al. (2007) Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17. doi:10.1051/0004-6361:20066389
  * Anderson et al. (2009) Anderson, J., Piotto, G., King, I. R., et al. 2009, ApJ, 697, L58. doi:10.1088/0004-637X/697/1/L58.
  * André et al. (2014) André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, Protostars and Planets VI, 27. doi:10.2458/azu_uapress_9780816531240-ch002.
  * Ballesteros-Paredes et al. (2020) Ballesteros-Paredes, J., André, P., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 76. doi:10.1007/s11214-020-00698-3.
  * Bally (2016) Bally, J. 2016, ARA&A, 54, 491. doi:10.1146/annurev-astro-081915-023341.
  * Bastian & Lardo (2018) Bastian, N. & Lardo, C. 2018, ARA&A, 56, 83. doi:10.1146/annurev-astro-081817-051839
  * Baumgardt et al. (2019) Baumgardt, H., Hilker, M., Sollima, A., et al. 2019, MNRAS, 482, 5138. doi:10.1093/mnras/sty2997
  * Beccari et al. (2017) Beccari, G., Petr-Gotzens, M. G., Boffin, H. M. J., et al. 2017, A&A, 604, A22. doi:10.1051/0004-6361/201730432
  * Belloche et al. (2011a) Belloche, A., Schuller, F., Parise, B., et al. 2011, A&A, 527, A145. doi:10.1051/0004-6361/201015733
  * Belloche et al. (2011b) Belloche, A., Parise, B., Schuller, F., et al. 2011, A&A, 535, A2. doi:10.1051/0004-6361/201117276
  * Binney & Tremaine (1987) Binney, J. & Tremaine, S. 1987, Princeton, N.J. : Princeton University Press, c1987., p291
  * Bonnor (1956) Bonnor, W. B. 1956, MNRAS, 116, 351. doi:10.1093/mnras/116.3.351
  * Carretta et al. (2006) Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2006, A&A, 450, 523. doi:10.1051/0004-6361:20054369.
  * Carretta et al. (2010) Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2010, A&A, 516, A55. doi:10.1051/0004-6361/200913451.
  * Carretta (2015) Carretta, E. 2015, ApJ, 810, 148. doi:10.1088/0004-637X/810/2/148
  * Clark & Bonnell (2005) Clark, P. C. & Bonnell, I. A. 2005, MNRAS, 361, 2. doi:10.1111/j.1365-2966.2005.09105.x.
  * Charbonnel (2016) Charbonnel, C. 2016, EAS Publications Series, 80-81, 177. doi:10.1051/eas/1680006
  * Cordoni et al. (2018) Cordoni, G., Milone, A. P., Marino, A. F., et al. 2018, ApJ, 869, 139. doi:10.3847/1538-4357/aaedc1
  * Da Rio et al. (2016) Da Rio, N., Tan, J. C., Covey, K. R., et al. 2016, ApJ, 818, 59. doi:10.3847/0004-637X/818/1/59.
  * Da Rio et al. (2017) Da Rio, N., Tan, J. C., Covey, K. R., et al. 2017, ApJ, 845, 105. doi:10.3847/1538-4357/aa7a5b.
  * D’Antona et al. (2002) D’Antona, F., Caloi, V., Montalbán, J., et al. 2002, A&A, 395, 69. doi:10.1051/0004-6361:20021220.
  * D’Antona et al. (2016) D’Antona, F., Vesperini, E., D’Ercole, A., et al. 2016, MNRAS, 458, 2122. doi:10.1093/mnras/stw387
  * Denissenkov & Hartwick (2014) Denissenkov, P. A. & Hartwick, F. D. A. 2014, MNRAS, 437, L21. doi:10.1093/mnrasl/slt133
  * D’Ercole et al. (2008) D’Ercole, A., Vesperini, E., D’Antona, F., et al. 2008, MNRAS, 391, 825. doi:10.1111/j.1365-2966.2008.13915.x
  * Decressin et al. (2007) Decressin, T., Charbonnel, C., & Meynet, G. 2007, A&A, 475, 859. doi:10.1051/0004-6361:20078425.
  * D’Ercole et al. (2008) D’Ercole, A., Vesperini, E., D’Antona, F., et al. 2008, MNRAS, 391, 825. doi:10.1111/j.1365-2966.2008.13915.x.
  * Dib et al. (2010) Dib, S., Hennebelle, P., Pineda, J. E., et al. 2010, ApJ, 723, 425. doi:10.1088/0004-637X/723/1/425
  * Dupree et al. (2017) Dupree, A. K., Dotter, A., Johnson, C. I., et al. 2017, ApJ, 846, L1. doi:10.3847/2041-8213/aa85dd
  * Galván-Madrid et al. (2007) Galván-Madrid, R., Vázquez-Semadeni, E., Kim, J., et al. 2007, ApJ, 670, 480. doi:10.1086/522081
  * Gómez et al. (2007) Gómez, G. C., Vázquez-Semadeni, E., Shadmehri, M., et al. 2007, ApJ, 669, 1042. doi:10.1086/521620.
  * Farias et al. (2020) Farias, J. P., Tan, J. C., & Eyer, L. 2020, ApJ, 900, 14. doi:10.3847/1538-4357/aba699
  * Farias et al. (2019) Farias, J. P., Tan, J. C., & Chatterjee, S. 2019, MNRAS, 483, 4999. doi:10.1093/mnras/sty3470
  * Federrath et al. (2014) Federrath, C., Schrön, M., Banerjee, R., et al. 2014, ApJ, 790, 128. doi:10.1088/0004-637X/790/2/128.
  * Federrath (2016) Federrath, C. 2016, MNRAS, 457, 375. doi:10.1093/mnras/stv2880.
  * Gieles et al. (2018) Gieles, M., Charbonnel, C., Krause, M. G. H., et al. 2018, MNRAS, 478, 2461. doi:10.1093/mnras/sty1059.
  * Hacar et al. (2017a) Hacar, A., Tafalla, M., & Alves, J. 2017, A&A, 606, A123. doi:10.1051/0004-6361/201630348
  * Hacar et al. (2020) Hacar, A., Hogerheijde, M. R., Harsono, D., et al. 2020, A&A, 644, A133. doi:10.1051/0004-6361/202038138.
  * Hatchell et al. (2005) Hatchell, J., Richer, J. S., Fuller, G. A., et al. 2005, A&A, 440, 151. doi:10.1051/0004-6361:20041836.
  * Hilker & Richtler (2000) Hilker, M. & Richtler, T. 2000, A&A, 362, 895.
  * Hillenbrand & Hartmann (1998) Hillenbrand, L. A. & Hartmann, L. W. 1998, ApJ, 492, 540. doi:10.1086/305076.
  * Huff & Stahler (2006) Huff, E. M. & Stahler, S. W. 2006, ApJ, 644, 355. doi:10.1086/503357.
  * Kainulainen et al. (2014) Kainulainen, J., Federrath, C., & Henning, T. 2014, Science, 344, 183. doi:10.1126/science.1248724
  * Keto et al. (2006) Keto, E., Broderick, A. E., Lada, C. J., et al. 2006, ApJ, 652, 1366. doi:10.1086/508251
  * Kirk et al. (2013) Kirk, H., Myers, P. C., Bourke, T. L., et al. 2013, ApJ, 766, 115. doi:10.1088/0004-637X/766/2/115
  * Klessen (2001) Klessen, R. S. 2001, ApJ, 556, 837. doi:10.1086/321626.
  * Kraft (1994) Kraft, R. P. 1994, PASP, 106, 553. doi:10.1086/133416.
  * Kroupa et al. (2001) Kroupa, P., Aarseth, S., & Hurley, J. 2001, MNRAS, 321, 699. doi:10.1046/j.1365-8711.2001.04050.x
  * Kroupa et al. (2018) Kroupa, P., Jeřábková, T., Dinnbier, F., et al. 2018, A&A, 612, A74. doi:10.1051/0004-6361/201732151.
  * Krumholz & Tan (2007) Krumholz, M. R. & Tan, J. C. 2007, ApJ, 654, 304. doi:10.1086/509101
  * Kudryavtseva et al. (2012) Kudryavtseva, N., Brandner, W., Gennaro, M., et al. 2012, ApJ, 750, L44. doi:10.1088/2041-8205/750/2/L44
  * Lada et al. (2008) Lada, C. J., Muench, A. A., Rathborne, J., et al. 2008, ApJ, 672, 410. doi:10.1086/523837.
  * Lagioia et al. (2019) Lagioia, E. P., Milone, A. P., Marino, A. F., et al. 2019, AJ, 158, 202. doi:10.3847/1538-3881/ab45f2.
  * Langer et al. (1993) Langer, G. E., Hoffman, R., & Sneden, C. 1993, PASP, 105, 301. doi:10.1086/133147
  * Luhman (2007) Luhman, K. L. 2007, ApJS, 173, 104. doi:10.1086/520114
  * Matzner & McKee (2000) Matzner, C. D. & McKee, C. F. 2000, ApJ, 545, 364. doi:10.1086/317785
  * Matzner (2001) Matzner, C. D. 2001, From Darkness to Light: Origin and Evolution of Young Stellar Clusters, 243, 757.
  * MacLean et al. (2015) MacLean, B. T., De Silva, G. M., & Lattanzio, J. 2015, MNRAS, 446, 3556. doi:10.1093/mnras/stu2348
  * Marino et al. (2018) Marino, A. F., Milone, A. P., Casagrande, L., et al. 2018, ApJ, 863, L33. doi:10.3847/2041-8213/aad868.
  * Milone et al. (2008) Milone, A. P., Bedin, L. R., Piotto, G., et al. 2008, ApJ, 673, 241. doi:10.1086/524188.
  * Milone et al. (2015) Milone, A. P., Marino, A. F., Piotto, G., et al. 2015, ApJ, 808, 51. doi:10.1088/0004-637X/808/1/51.
  * Milone et al. (2017) Milone, A. P., Piotto, G., Renzini, A., et al. 2017, MNRAS, 464, 3636. doi:10.1093/mnras/stw2531.
  * Milone et al. (2018) Milone, A. P., Marino, A. F., Renzini, A., et al. 2018, MNRAS, 481, 5098. doi:10.1093/mnras/sty2573
  * Mueller et al. (2002) Mueller, K. E., Shirley, Y. L., Evans, N. J., et al. 2002, ApJS, 143, 469. doi:10.1086/342881
  * Myers (2009) Myers, P. C. 2009, ApJ, 700, 1609. doi:10.1088/0004-637X/700/2/1609.
  * Nakamura & Li (2007) Nakamura, F. & Li, Z.-Y. 2007, ApJ, 662, 395. doi:10.1086/517515.
  * Norman & Silk (1980) Norman, C. & Silk, J. 1980, ApJ, 238, 158. doi:10.1086/157969
  * Oliveira et al. (2020) Oliveira, R. A. P., Souza, S. O., Kerber, L. O., et al. 2020, ApJ, 891, 37. doi:10.3847/1538-4357/ab6f76
  * Parmentier (2004) Parmentier, G. 2004, MNRAS, 351, 585. doi:10.1111/j.1365-2966.2004.07804.x.
  * Parmentier & Pfalzner (2013) Parmentier, G. & Pfalzner, S. 2013, A&A, 549, A132. doi:10.1051/0004-6361/201219648.
  * Parmentier et al. (2014) Parmentier, G., Pfalzner, S., & Grebel, E. K. 2014, ApJ, 791, 132. doi:10.1088/0004-637X/791/2/132.
  * Parmentier (2019) Parmentier, G. 2019, ApJ, 887, 179. doi:10.3847/1538-4357/ab53d6
  * Parmentier & Pasquali (2020) Parmentier, G. & Pasquali, A. 2020, ApJ, 903, 56. doi:10.3847/1538-4357/abb8d3.
  * Piotto et al. (2007) Piotto, G., Bedin, L. R., Anderson, J., et al. 2007, ApJ, 661, L53. doi:10.1086/518503.
  * Pokhrel et al. (2021) Pokhrel, R., Gutermuth, R. A., Krumholz, M. R., et al. 2021, arXiv:2104.04551
  * Portegies Zwart et al. (2004) Portegies Zwart, S. F., Baumgardt, H., Hut, P., et al. 2004, Nature, 428, 724. doi:10.1038/nature02448
  * Prantzos & Charbonnel (2006) Prantzos, N. & Charbonnel, C. 2006, A&A, 458, 135. doi:10.1051/0004-6361:20065374.
  * Preibisch (2012) Preibisch, T. 2012, Research in Astronomy and Astrophysics, 12, 1. doi:10.1088/1674-4527/12/1/001.
  * Rahner et al. (2017) Rahner, D., Pellegrini, E. W., Glover, S. C. O., et al. 2017, MNRAS, 470, 4453. doi:10.1093/mnras/stx1532.
  * Rahner et al. (2018) Rahner, D., Pellegrini, E. W., Glover, S. C. O., et al. 2018, MNRAS, 473, L11. doi:10.1093/mnrasl/slx149.
  * Schneider et al. (2018) Schneider, F. R. N., Ramírez-Agudelo, O. H., Tramper, F., et al. 2018, A&A, 618, A73. doi:10.1051/0004-6361/201833433
  * Shukirgaliyev et al. (2017) Shukirgaliyev, B., Parmentier, G., Berczik, P., et al. 2017, A&A, 605, A119. doi:10.1051/0004-6361/201730607
  * Smith et al. (2016) Smith, R. J., Glover, S. C. O., Klessen, R. S., et al. 2016, MNRAS, 455, 3640. doi:10.1093/mnras/stv2559
  * Tan et al. (2014) Tan, J. C., Beltrán, M. T., Caselli, P., et al. 2014, Protostars and Planets VI, 149. doi:10.2458/azu_uapress_9780816531240-ch007
  * Tan et al. (2006) Tan, J. C., Krumholz, M. R., & McKee, C. F. 2006, ApJ, 641, L121. doi:10.1086/504150.
  * Taylor et al. (1996) Taylor, S. D., Morata, O., & Williams, D. A. 1996, A&A, 313, 269
  * Vázquez-Semadeni et al. (2005) Vázquez-Semadeni, E., Kim, J., Shadmehri, M., et al. 2005, ApJ, 618, 344. doi:10.1086/425951.
  * Ward-Thompson et al. (2007) Ward-Thompson, D., André, P., Crutcher, R., et al. 2007, Protostars and Planets V, 33