Star Cluster Formation in Cosmological Simulations. II. Effects of Star Formation Efficiency and Stellar Feedback

The full description of the simulation setup is presented in Paper I. Here we first briefly summarize the common aspects and then describe the improvements to the model. Most of the simulations in this paper are new and have distinct features from those discussed in Paper I.

The simulations are run with the Eulerian gas dynamics and N-body adaptive refinement tree (ART) code (Kravtsov et al. 1997; Kravtsov 1999, 2003; Rudd et al. 2008), which includes several key physical ingredients, such as three-dimensional radiative transfer (Gnedin & Abel 2001; Gnedin 2014), a non-equilibrium chemical network of hydrogen and helium, non-equilibrium cooling and photoionization heating, phenomenological molecular hydrogen formation and destruction (Gnedin & Kravtsov 2011), and a subgrid-scale turbulence model (Semenov et al. 2016). We run the cosmological simulations from the initial condition that contains a main halo with a total mass ${M}_{200}\approx {10}^{12}\,{M}_{\odot }$ at z = 0 in a periodic box of 4 comoving Mpc in size. All simulations start on a 128³ root grid, which sets the dark matter particle mass ${m}_{\mathrm{DM}}=1.05\times {10}^{6}\,{M}_{\odot }$. The high dynamic range of spatial resolution is achieved by adaptive mesh refinement, where additional cell levels are added to the root grid. We apply a Lagrangian refinement criterion for both dark matter and gas components so that the mass of all gas cells varies only within a narrow range at all times. In addition, we adopt a Jeans refinement criterion with which cells larger than twice the local Jeans length will be refined.

In Paper I, we developed a novel algorithm for the formation of star clusters in cosmological simulations: continuous cluster formation (CCF). This algorithm is different from most of the star formation methods implemented in current mesh- and particle-based cosmological simulations, where star particles are spawned instantaneously by a Poisson process with their masses set beforehand. In CCF, each star particle represents a single star cluster formed at a local density peak of the molecular gas. Cluster particles grow their masses continuously by accreting ambient material at every local timestep (∼1000 yr), meaning that the star formation process is resolved with high time resolution. The growth of a cluster particle is terminated by its own energy and momentum feedback. Thus, the final particle mass is set self-consistently and represents the mass of a single star cluster formed within the natal GMC.

In order to avoid dependence on the time-variable physical size of a cell, L_cell, the cluster particle is allowed to grow its mass via gas accretion within a spherical region of fixed physical size. We interpret this sphere as the dense part of a GMC that would form a bound stellar system in freefall collapse. The optimal value of the sphere size was investigated in Paper I, ${R}_{\mathrm{GMC}}=5\,$ pc, which is similar to the sizes of observed massive cluster-forming clouds or clumps (e.g., Urquhart et al. 2014). We will henceforth refer to this star-forming sphere as the "GMC."

At the highest refinement level, the GMC sphere fully includes the peak-density ("central") cell and partially overlaps with its 26 neighbor cells from a 3 × 3 × 3 cube configuration. Accessing more than these immediate neighbors is computationally prohibitive in the ART code. The growth rate of a given cluster depends on the H₂ density of the overlapping cells,

$$\begin{eqnarray}&&\dot{M}=\displaystyle \sum _{\mathrm{cell}}{f}_{\mathrm{GMC}}\,{V}_{\mathrm{cell}}\,{\dot{\rho }}_{* }=\displaystyle \frac{{\epsilon }_{\mathrm{ff}}}{{\tau }_{\mathrm{ff}}}\,\displaystyle \sum _{\mathrm{cell}}{f}_{\mathrm{GMC}}\,{V}_{\mathrm{cell}}\,{f}_{{{\rm{H}}}_{2}}\,{\rho }_{\mathrm{gas}},\end{eqnarray} \tag{ 1 }$$

where V_cell is the volume of each neighbor cell, f_GMC is the fraction of V_cell that overlaps with the GMC, ${f}_{{{\rm{H}}}_{2}}$ is the mass fraction of hydrogen in molecular phase, ρ_gas is the total gas density, and ${\epsilon }_{\mathrm{ff}}$ is the star formation efficiency per freefall time. The freefall time is defined as

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}\equiv {\left(\displaystyle \frac{32G\rho }{3\pi }\right)}^{-1/2}\approx 1.6\,\mathrm{Myr}{\left(\displaystyle \frac{{n}_{{\rm{H}}}}{{10}^{3}{\mathrm{cm}}^{-3}}\right)}^{-1/2},\end{eqnarray} \tag{ 2 }$$

where n_H is the total number density of hydrogen of the central cell. The mass growth of the cluster particle is calculated and added to the particle mass at each local timestep, typically Δt ∼ 100 yr. The mass accumulation history of each cluster is thus temporally resolved with many thousands of steps.

Cluster formation is allowed only in cells above the number density ${n}_{{\rm{H}},\mathrm{th}}={10}^{3}\,$ cm⁻³. This is close to the observational estimate by Kainulainen et al. (2014) of the H₂ density threshold for star formation in nearby GMCs of ∼5 × 10³ cm⁻³. In addition, to avoid the creation of very small and numerous clusters, new particles are created only if their expected final mass is above threshold M_th. The expected mass is calculated as the initial rate $\dot{M}$ times the maximum allowed formation time, ${\tau }_{\max }=15$ Myr. Since the actual duration of cluster formation is significantly shorter in the new runs, we discuss and revise M_th in Section 2.2.4.

2.2.1. Gas Cell Refinement

2.2.2. Redshift-independent Star Formation Efficiency

The mass accretion rate of a cluster is given by Equation (1), which contains the parameter ${\epsilon }_{\mathrm{ff}}$. However, the meaning of ${\epsilon }_{\mathrm{ff}}$ is quantitatively different from what is commonly used in galaxy formation simulations. In the traditional prescription, stellar particles are given the mass calculated as the star formation rate density ${\dot{\rho }}_{* }$ times the volume of the cell containing the particle. The cell size is fixed in comoving coordinates but expands in proper physical coordinates, so that for the same gas density (typically, near the threshold for star formation), the particle mass increases with time. In contrast, in our model, each star cluster grows its mass over several million yr by accreting material within a cloud of fixed physical radius, ${R}_{\mathrm{GMC}}$. Our ${\epsilon }_{\mathrm{ff}}$ is applied to a fixed volume and does not vary with cosmic time. This is an important improvement.

We can relate our value of ${\epsilon }_{\mathrm{ff}}$ to that on the scale of a cell via the differences in volume. Consider the case when the diameter of the GMC sphere is smaller than a star-forming cell, and therefore the GMC is completely embedded in the cell. The smaller volume allowed to participate in star formation translates into the smaller efficiency for the whole cell,

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{ff},\mathrm{cell}}=\displaystyle \frac{{V}_{\mathrm{GMC}}}{{V}_{\mathrm{cell}}}\,{\epsilon }_{\mathrm{ff}}=\displaystyle \frac{4\pi }{3}{\left(\displaystyle \frac{{R}_{\mathrm{GMC}}}{{L}_{\mathrm{cell}}}\right)}^{3}\,{\epsilon }_{\mathrm{ff}}.\end{eqnarray} \tag{ 3 }$$

For example, the seemingly high value ${\epsilon }_{\mathrm{ff}}=50 \%$ in our CCF model is equivalent to ${\epsilon }_{\mathrm{ff},\mathrm{cell}}\approx 1 \%$ for a cell of 30 pc, a typical size for current highest-resolution cosmological simulations.

We allow cluster formation at the three finest refinement levels, which is consistent with our adopted star formation density threshold: n_crit = 10³ cm⁻³. For perfectly Lagrangian refinement, the average number density of hydrogen plus helium atoms for cells at level l and redshift z in our simulation box is

$$\begin{eqnarray}&&n(z)\approx 900\,{\mathrm{cm}}^{-3}\times {8}^{l-9}{\left(\displaystyle \frac{1+z}{3}\right)}^{3}.\end{eqnarray} \tag{ 4 }$$

We can see that at z = 2, cells at level 9 (the third-finest level at that time) have an average density close to ${n}_{\mathrm{crit}}$. Cells create new clusters as soon as they exceed the density threshold, which usually happens already at the third-finest level. The physical cell size at that level is in the range L_cell = 12–24 pc, which is larger than $2{R}_{\mathrm{GMC}}$. Thus, the corresponding cell efficiency initially is ${\epsilon }_{\mathrm{ff},\mathrm{cell}}=(0.038\mbox{--}0.30){\epsilon }_{\mathrm{ff}}$.

This is an important point for the comparison with the low efficiency per freefall time inferred in Galactic star-forming regions, ${\epsilon }_{\mathrm{ff}}\sim 1 \%$ (Kennicutt 1998; Krumholz et al. 2012; Vutisalchavakul et al. 2016). In our model, clusters begin forming with a similarly low efficiency. As the gas continues to collapse in freefall, the density increases and cells are refined to the highest level. Then the effective ${\epsilon }_{\mathrm{ff},\mathrm{cell}}$ increases as well, to match ${\epsilon }_{\mathrm{ff}}$, in agreement with the results of Murray (2011) and Lee et al. (2016).

As in every numerical simulation, at the highest refinement level, we are not properly resolving gas collapse and therefore underestimate the density and overestimate the freefall time. This numerical effect also requires adopting a larger value of ${\epsilon }_{\mathrm{ff}}$ to maintain the correct star formation rate (SFR).

These algorithmic and numerical issues make it difficult to directly compare our parameter ${\epsilon }_{\mathrm{ff}}$ with the observed efficiency on ∼30 pc scales of GMCs. We will discuss instead the distribution of the integral efficiency, which is the fraction of gas mass turned into stars, in Figure 8 and Section 4.8.

2.2.3. Restriction on Creation of New Clusters Near Existing Active Clusters

In Paper I, cluster particles were seeded in cold dense cells that contain local density peaks. We follow similar procedure here but with additional algorithmic improvements and constraints.

As we discussed in Section 2.2.1, the physical size of gas cells at the finest refinement level is always kept in the range 3–6 pc. This means that clusters in the smallest cells can accrete gas from all 27 cells that overlap with the GMC sphere (see Figure 1). To be consistent with our interpretation that the peak-density cell represents the star-forming part of a single GMC, the neighboring cells should not produce new clusters that would compete for gas supply with an existing cluster. Therefore, we prohibit a cell from creating a new particle if it already has an actively growing cluster in a neighbor cell.

Zoom In Zoom Out Reset image size

This restriction makes sense only for the finest refinement level cells. If cluster particles are formed in coarser cells, which enclose a large fraction of volume of the GMC sphere, then there is no competition with neighboring cells. In this case, we allow new clusters to be created.

We also do not apply this restriction to the finest-level cells that are separated by more than one cell size from the central cell. These are the diagonal cells in the cube. Thus, the only cells that are prohibited from creating new clusters are the six cells "face-touching" the central cell.

To describe the above criteria for cluster seeding more quantitatively, we define an overlap fraction, f_over, of the volume of a given "face-touching" neighbor that is occupied by the central GMC sphere. If a given cell satisfies the star formation criterion but has any "face-touching" neighbor with a significant overlap, f_over > 20%, that already contains an actively growing cluster, the above cell is not allowed to create a new particle. If the overlap fraction is small, f_over < 20%, we do not apply this restriction.

All of these checks are meant to minimize the impact of the neighbor restriction as much as is reasonable. However, we show in Section 3 that it still has a significant impact on the shape of the CIMF, relative to the results in Paper I.

2.2.4. No Removal of Low-mass Clusters

2.2.5. No Molecular Fraction Threshold for Continuing Cluster Growth

2.2.6. Gap Elimination

2.2.7. Mass-loss Rate Due to Stellar Evolution

In Paper I, a stellar particle lost mass due to passive stellar evolution and stellar winds at the rate given by

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dt}}=\displaystyle \frac{\eta }{t+{\tau }_{\mathrm{loss}}}\,{M}_{t=0},\end{eqnarray} \tag{ 5 }$$

where η is the typical fraction of mass loss and ${\tau }_{\mathrm{loss}}$ is the characteristic timescale, both of which depend on the IMF of stars. For a Kroupa (2001) IMF, we used η = 0.046 and ${\tau }_{\mathrm{loss}}=2.76\times {10}^{5}$ yr. The time evolution of the cumulative fraction of stellar mass loss is then obtained by integrating Equation (5) over time:

$$\begin{eqnarray}&&{f}_{\mathrm{loss}}(t)=\displaystyle \frac{{M}_{\mathrm{loss}}}{{M}_{t=0}}={\int }_{0}^{t}\displaystyle \frac{\eta }{t+{\tau }_{\mathrm{loss}}}{dt}=\eta \mathrm{ln}\left(\displaystyle \frac{t}{{\tau }_{\mathrm{loss}}}+1\right).\end{eqnarray} \tag{ 6 }$$

Figure 2 shows this mass-loss history, together with the results from the detailed stellar population synthesis model FSPS⁶ for the same IMF (Conroy et al. 2009; Conroy & Gunn 2010). We find that the previous implementation overestimates the mass-loss rate over the whole stellar lifetime. The overestimate is more prominent for younger stellar populations. For example, after 3 Myr, about 10% of the stellar mass is lost in Paper I modeling but less than 1% in FSPS modeling.

Zoom In Zoom Out Reset image size

Here we introduce a new mass-loss model that fits the time evolution of the cumulative fraction in FSPS with a second-order polynomial. We write ${f}_{\mathrm{loss}}(t)={{ax}}^{2}+{bx}+c$, where $x={\mathrm{log}}_{10}(t/\mathrm{yr})$, and obtain the best-fit parameters: a = −0.010, b = 0.288, and c = −1.42. The mass-loss rate at a given age t is obtained by differentiating f_loss(t) with respect to t:

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dt}}=\displaystyle \frac{b+2a{\mathrm{log}}_{10}(t)}{\mathrm{ln}10}\,\displaystyle \frac{{M}_{t=0}}{t}.\end{eqnarray} \tag{ 7 }$$

Note that this expression works for stars in the age range from 2.75 Myr to 13.7 Gyr. For ages younger than 2.75 Myr, we assume no appreciable mass loss.

2.2.8. Changes to Stellar Feedback

As in Paper I, the feedback subgrid model consists of mass, momentum, and energy injections from stellar winds, radiation pressure, and SN explosions. The main update in this paper is in the implementation of SN remnant (SNR) feedback. In the new SNR model, we estimate the partition of the thermal, kinetic, and turbulence energies of SNRs using a parameterization model calibrated by the high-resolution hydrodynamic simulations in an inhomogeneous turbulent medium by Martizzi et al. (2015). The energy and momentum input from SNRs depends on the ambient density and spatial resolution, as described in detail in Semenov et al. (2016).

One caveat of the SNR model is that it is calibrated by the simulations of SN explosion in isolation, rather than in a more complex star-forming environment. In reality, clusters produce a large number of massive stars that undergo SN explosion over several Myr. Gentry et al. (2017) found that such clustering of SNe can enhance momentum feedback by an order of magnitude relative to that delivered by an isolated SN. Another important motivation of boosting the injected momentum is to compensate for the momentum loss due to advection errors as the SN shell moves across the simulation grid (e.g., Semenov et al. 2017). Therefore, in our new runs, we boost the momentum feedback of the Martizzi et al. (2015) SNR model by a factor ${f}_{\mathrm{boost}}=3$, 5, or 10, as listed in Table 1.

Table 1.  Model Runs

Name	zf	${\epsilon }_{\mathrm{ff}}$	Feedback	${L}_{\mathrm{cell}}$ (pc) at z = 2
SFE10	1.7	0.1	earlya+5*SNRb	5
SFE50	1.5	0.5	early+5*SNR	5
SFE50-SNR3	2.0	0.5	early+3*SNR	5
SFE100	1.5	1.0	early+5*SNR	5
SFE200	1.5	2.0	early+5*SNR	5
SFEturb	2.7	Variable	early+5*SNR	5
fid-PaperI	3.2	0.1	early+old SNc	20
SNR10-old	2.0	0.1	early+10*SNR	10

Notes.

^a"early": early feedback schemes including stellar wind and radiative pressure. ^b"SNR": resolution-dependent SNR feedback scheme in Martizzi et al. (2015); the number before "SNR" is the boosting factor. ^c"old SN": previous SN thermal and momentum feedback with the amount calibrated by Agertz & Kravtsov (2015).

Download table as:  ASCIITypeset image

2.2.9. Initial Bound Fraction

In Paper I, we assumed that the whole accreted mass of cluster particles during their active growth is gravitationally self-bound when the particles emerge from their natal clouds. That is, the stellar particle mass is the bound cluster mass. This assumption does not take into account the complex dynamical evolution of star clusters in the early phase, when the boundedness is affected by the hierarchical structure of the ISM and gas expulsion due to stellar feedback. Recent observations and numerical simulations of turbulent clouds suggest that the fraction of mass in a given star-forming complex that remains bound after a few Myr depends strongly on the star formation efficiency on the scale of GMCs (Goodwin 1997; Geyer & Burkert 2001; Goodwin & Bastian 2006; Smith et al. 2011; Kruijssen et al. 2012).

In this paper we introduce the initial bound fraction, f_i, and assign it to each stellar particle. We adopt a linear dependence of f_i on the local star formation efficiency,

$$\begin{eqnarray}&&{f}_{i}\equiv \min \left(\displaystyle \frac{{\epsilon }_{\mathrm{int}}}{{\epsilon }_{\mathrm{core}}},1\right),\end{eqnarray} \tag{ 8 }$$

where _core = 0.5 is the correction factor for mass loss due to protostellar outflows, suggested by Kruijssen (2012), and _int is the integral star formation efficiency, defined as the ratio between the full mass of the active stellar particle (here we write it explicitly as M_f but later omit the subscript "f") and the maximum baryon mass of the GMC throughout the whole course of cluster accretion:

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{int}}\equiv \displaystyle \frac{{M}_{{\rm{f}}}}{{\max }_{\mathrm{time}}({M}_{\mathrm{star}}+{M}_{\mathrm{gas}})}.\end{eqnarray} \tag{ 9 }$$

Note that the linear relation between local star formation efficiency and initial bound fraction is based on the analysis of the star formation simulations of Bonnell et al. (2008). Since those simulations only cover high-efficiency cases with f_i > 0.4, it is unclear whether such linearity is valid for low-efficiency GMCs. Recent simulations of isolated GMCs (Goodwin 2009; Dale et al. 2014; Gavagnin et al. 2017) suggest that the bound fraction depends not only on ${\epsilon }_{\mathrm{int}}$ but also on the dynamical state of the young stars before gas expulsion by feedback. Including this effect in future work could improve the model.

It is important to emphasize that the GMC mass in our simulations varies on a timescale of less than a Myr, as the central part is converted into stars and outside gas flows in, in essentially freefall. The sum of stellar and gas mass usually increases over several Myr before declining as the feedback of young stars disperses the remaining gas. This situation is very different from simulations of star formation in isolated clouds, where the total baryon mass is fixed as the initial cloud mass.

Bound cluster mass in our model is then defined as the stellar particle mass at the end of the active growth period times the initial bound fraction: $M={f}_{i}\,{M}_{{\rm{f}}}$. The fraction of stars remaining bound to the cluster continues to evolve due to dynamical evaporation and tidal stripping. We calculate this process with a subgrid model and will discuss the evolution of the cluster mass function in a follow-up paper. We attach the initial bound fraction variable to every cluster particle and assume that the unbound part remains near the cluster, so that both the bound and unbound parts move together under the same gravitational potential.

2.2.10. Alternative Definition of Cluster Formation Timescale

Although the actual mass growth history of individual star clusters is difficult to determine in observations, recent theoretical models and simulations suggest a dynamic, time-variable process. Self-gravitating collapse of the dense parts of GMCs is thought to accelerate star formation until stellar feedback, in the form of stellar winds and radiative pressure/heating, pushes the accreting gas out. The whole process happens quickly, as indicated by the observed age spread of stars in embedded clusters within 3–4 Myr, only a couple freefall times of their natal clouds; see Section 2.4.1.

One way to quantify the age spread is to estimate the timescale for the formation of a given fraction of stars within a cluster, e.g., the period when the cluster mass grows from 10% to 90% of its final mass. However, since the final mass of the clusters in our model is not known until their growth is terminated by the feedback, obtaining the 10%–90% timescale requires tracking the detailed growth history of all clusters, which is computationally prohibitive. Therefore, in Paper I, the average duration of cluster formation was calculated as the mass-weighted timescale,

$$\begin{eqnarray}&&{\tau }_{\mathrm{ave}}\equiv \displaystyle \frac{{\int }_{0}^{{\tau }_{\max }}t\,\dot{M}(t)\,{dt}}{{\int }_{0}^{{\tau }_{\max }}\dot{M}(t)\,{dt}},\end{eqnarray} \tag{ 10 }$$

where $\dot{M}(t)$ is the cluster SFR at time t. This definition best describes steady mass accretion. For example, in the case of constant SFR, ${\tau }_{\mathrm{ave}}={\tau }_{\max }/2$.

One caveat of this definition is that, in the case of increasing and then decreasing SFR, ${\tau }_{\mathrm{ave}}$ actually records the epoch when $\dot{M}$ reaches its peak, rather than the duration of the process. Here we introduce a new quantity that better describes the width of such SFH. The age spread is defined as the ratio between the full particle mass M_f and the mass-weighted SFR over the whole growth history $\langle \dot{M}\rangle$:

$$\begin{eqnarray}&&{\tau }_{\mathrm{spread}}\equiv \displaystyle \frac{{M}_{{\rm{f}}}}{\langle \dot{M}\rangle }=\displaystyle \frac{{M}_{{\rm{f}}}}{{\int }_{0}^{{\tau }_{\max }}{\dot{M}}^{2}{dt}/{M}_{{\rm{f}}}}.\end{eqnarray} \tag{ 11 }$$

For a power-law mass accretion history, $\dot{M}\propto {t}^{\alpha }$, $\alpha \ne -1$, both ${\tau }_{\mathrm{ave}}$ and ${\tau }_{\mathrm{spread}}$ can be explicitly evaluated as

$$\begin{eqnarray}&&{\tau }_{\mathrm{ave}}=\displaystyle \frac{\alpha +1}{\alpha +2}{\tau }_{\max }\quad \mathrm{and}\quad {\tau }_{\mathrm{spread}}=\displaystyle \frac{2\alpha +1}{{(\alpha +1)}^{2}}{\tau }_{\max }.\end{eqnarray} \tag{ 12 }$$

However, in the case when $\dot{M}$ exhibits a peak, e.g., a Gaussian function $\dot{M}\propto \exp (-{(t-{t}_{0})}^{2}/2{\sigma }^{2})$,

$$\begin{eqnarray}{\tau }_{\mathrm{ave}} & = & {t}_{0}+\sqrt{\displaystyle \frac{2}{\pi }}\sigma \\ & & \times \displaystyle \frac{\exp {(-{t}_{0}^{2}/2{\sigma }^{2})-\exp (-({t}_{0}-{\tau }_{\max })}^{2}/2{\sigma }^{2})}{\mathrm{erf}({t}_{0}/\sqrt{2}\sigma )+\mathrm{erf}(({\tau }_{\max }-{t}_{0})/\sqrt{2}\sigma )},\end{eqnarray} \tag{ 13 }$$

while the new definition gives

$$\begin{eqnarray}&&{\tau }_{\mathrm{spread}}=\sqrt{\pi }\sigma \displaystyle \frac{{[\mathrm{erf}({t}_{0}/\sqrt{2}\sigma )+\mathrm{erf}(({\tau }_{\max }-{t}_{0})/\sqrt{2}\sigma )]}^{2}}{\mathrm{erf}({t}_{0}/\sigma )+\mathrm{erf}(({\tau }_{\max }-{t}_{0})/\sigma )}.\end{eqnarray} \tag{ 14 }$$

The relationship between the intrinsic width of the mass accretion history σ and the derived ages ${\tau }_{\mathrm{ave}}$ and ${\tau }_{\mathrm{spread}}$ is shown in Figure 3.

Zoom In Zoom Out Reset image size

When $\sigma \to 0$, the SFH reduces to a δ function located at $t={t}_{0}$. In this case, ${\tau }_{\mathrm{ave}}\to {t}_{0}$, while ${\tau }_{\mathrm{spread}}\to 2\sqrt{\pi }\sigma$. This extreme case illustrates the problem with our previous definition. Rather than recording the age spread of stars, ${\tau }_{\mathrm{ave}}$ actually reflects the epoch of the peak of star formation. Instead, the new definition follows the true width of the age distribution.

As σ increases, the SFH becomes flatter. At the other extreme, if $\sigma \gt {\tau }_{\max }$, the SFR can be considered constant. The above equations then reduce to ${\tau }_{\mathrm{ave}}\to {t}_{0}/2$ and ${\tau }_{\mathrm{spread}}\to {t}_{0}$.

In this paper, we present new runs with the above improvements. All start with the same initial conditions as in Paper I. Their key physical parameters are listed in Table 1. The number after "SFE" in their names corresponds to the local ${\epsilon }_{\mathrm{ff}}$ in percent. One exception is the "SFEturb" run, with variable turbulence-dependent ${\epsilon }_{\mathrm{ff}}$ based on the parameterization of ${\epsilon }_{\mathrm{ff}}$ as a function of cell turbulence energy (Padoan et al. 2012). The evolution of the turbulence energy is calculated by a subgrid-scale model from Schmidt et al. (2014), implemented in the ART code by Semenov et al. (2016). The standard value of the SNR momentum boost factor is ${f}_{\mathrm{boost}}=5$. In run SFE50-SNR3, we adopt ${f}_{\mathrm{boost}}=3$ to test the dependence of global SFR on this factor.

In addition to these new runs with full updates, we include for comparison the fiducial run "old-fid" from Paper I and the "SNR10-old" run with ${f}_{\mathrm{boost}}=10$ and some intermediate degree of update. These latter runs illustrate progression in the development of our algorithm.

(a) fid-PaperI.—The fiducial run in Paper I with ${\epsilon }_{\mathrm{ff}}=10 \%$. As shown in Figure 5, this run overestimates the expected SFR of the main galaxy by more than one order of magnitude at z = 4–8. Due to the high SFR, a large number of clusters more massive than 10⁵$\,{M}_{\odot }$ are formed at high z. The inability to suppress SFR at high z is caused by the weak feedback that is used in this run. Another sign of the ineffectiveness of the feedback comes from the mass-weighted cluster formation timescale ${\tau }_{\mathrm{ave}}$. In Paper I, we showed that the median value of ${\tau }_{\mathrm{ave}}$ for all clusters in this run is about 3 Myr, consistent with the observational constraint described in Section 2.4.1. Detailed analysis suggests that this median value is dominated by a large number of less massive cluster particles. For clusters more massive than ${10}^{5}\,{M}_{\odot }$, however, ${\tau }_{\mathrm{ave}}$ is much longer, sometimes even longer than ${\tau }_{\max }/2=7.5\,$ Myr, suggesting that the stellar feedback cannot terminate the gas accretion process for massive clusters within the allowed accretion timescale that is assigned in the simulation.

(b) SFE10-old.—Same as fid-PaperI, except (1) a refinement strategy that keeps the ratio of ${R}_{\mathrm{GMC}}/{L}_{\mathrm{cell}}$ roughly constant but one level coarser than that described in Section 2.2.1 and (2) a new SNR feedback prescription described in Section 2.2.8 with a momentum boosting factor ${f}_{\mathrm{boost}}=5$. We found that the SFH of this run follows the abundance matching results until z ≈ 2.5, when a major merger happens in the main galaxy. This merger brings so much cold gas into the inner regions of the galaxy that stellar feedback is unable to disperse it. An intense starburst occurs, and a prominent stellar spheroid forms at the galactic center. As we discuss in the next section, adding one additional refinement level prevents this starburst and leads to a reasonable SFH. Another result is that, although the SFR of this run is much smaller than that of fid-PaperI due to the adoption of the new feedback model, there are still many clusters, especially massive ones, that have very long formation timescales. This finding inspired a detailed investigation of the origin of the long timescale and the implementation of gap elimination algorithms, which we describe in Sections 2.2.5 and 2.2.6.

(c) SNR10-old.—Same as SFE10-old but with a larger boosting factor ${f}_{\mathrm{boost}}=10$. The goal of this run is to probe the upper limit of ${f}_{\mathrm{boost}}$. As shown in Figure 5, the SFH is always under the abundance matching result. This underestimate of SFR sets an upper limit on the reasonable choice of ${f}_{\mathrm{boost}}$.

For comparison, we have also run a simulation with the standard ${f}_{\mathrm{boost}}=5$ but boosted the momentum feedback from radiative pressure by a factor of 10. This boost produced no apparent effect on SFH, unlike the SN momentum boost.

2.4.1. Observed Age Spread of Young Clusters

To use the predicted formation duration as an additional test of the models, in Table 2, we compile recent measurements of the spread of the relative ages of stars in several young star clusters. Although obtaining accurate age measurements is still challenging, current observations suggest that the age spread should be less than 6 Myr for various star-forming regions in different galaxies. This upper limit also agrees with small-scale hydrodynamic simulations of cluster formation, which suggest that the star formation process should proceed on a timescale comparable to the freefall time (e.g., Hartmann et al. 2012; Grudić et al. 2018). This relatively short timescale provides a strong constraint on the implementation of star formation and stellar feedback.

Table 2.  Compilation of the Observed Properties of Nearby Young Star Clusters

Name	Age (Myr)	Age Spread (Myr)	Mass ($\,{M}_{\odot }$)	References
Orion Nebula Cluster	4–5	1–3	2000	Da Rio et al. (2010b), Jeffries et al. (2011)
R136	4–5	3	4.5e5	Massey & Hunter (1998)
LH95	4	2.8–4.4	?	Da Rio et al. (2010a)
NGC 1569A	$\lt 5$	Free of dust extinction	1e6	Maoz et al. (2001)
Westerlund 1	5	0.4 or $\lt 1$	6.3e4	Negueruela et al. (2010), Kudryavtseva et al. (2012)
NGC 4103	?	2–4	?	Forbes (1996)
NGC 3603 YC	2	0.1	?	Kudryavtseva et al. (2012)
NGC 3603 HII	1	$\lt 1$	1.9e4	Pang et al. (2013)
W3 Main	?	2–3	?	Bik et al. (2012)
Antennae clusters	?	$\lt 6,{A}_{V}=1$ mag	?	Whitmore & Zhang (2002)
NGC 4449 clusters	?	$\lt 5,{A}_{V}=0.5\mbox{--}1.5$ mag	0.5–5e4	Reines et al. (2008)
M83 clusters	?	$\lt 4$	?	Hollyhead et al. (2015)

Download table as:  ASCIITypeset image

2.4.2. SFR versus Maximum Cluster Mass in Nearby Galaxies

In Figure 5, we show the SFH of the main halo for several runs with different values of the star formation and feedback parameters. The SFH is calculated from all cluster particles within 0.5 R_vir of the main galaxy from the last available snapshot of each run, where R_vir is the virial radius of the dark matter halo. The rate of star formation is averaged over 100 Myr around a given epoch in order to smooth out stochasticity.

Zoom In Zoom Out Reset image size

In general, simulations with ${f}_{\mathrm{boost}}=5$ agree well with the abundance matching results, in the sense of both SFR and cumulative stellar mass growth. This suggests that the strength of stellar feedback used in these runs is appropriate to reproduce the inefficient star formation at high redshift. Changing the intensity of feedback has a measurable effect on the SFH. The run with ${f}_{\mathrm{boost}}=3$ systematically overestimates SFR at all times, while the run with ${f}_{\mathrm{boost}}=10$ underestimates SFR by an order of magnitude. The resulting SFHs of the two runs are well beyond the 1σ confidence intervals of the abundance matching result. Therefore, the choice of ${f}_{\mathrm{boost}}$ in our simulations is bracketed between 3 and 10.

On the other hand, varying ${\epsilon }_{\mathrm{ff}}$ from ∼0.01 to 2.0 has almost no systematic effect on SFH. This insensitivity to the value of ${\epsilon }_{\mathrm{ff}}$ has already been found in several recent cosmological simulations with the implementation of strong stellar feedback and high spatial resolution (Agertz et al. 2013; Hopkins et al. 2013; Agertz & Kravtsov 2015). The feedback-controlled low-efficiency star formation activity can be interpreted as the short lifetime of star-forming regions, which then requires a large number of star-forming cycles to convert all the gas into stars (Semenov et al. 2017).

However, the lack of sensitivity of SFH to the local star formation efficiency does not mean that one can assign arbitrary values to ${\epsilon }_{\mathrm{ff}}$ in galaxy formation simulations. As we show in the rest of this section, ${\epsilon }_{\mathrm{ff}}$ has dramatic effects on the properties of individual star-forming regions and, in turn, the properties of young star clusters formed within.

Above, we showed that with ${f}_{\mathrm{boost}}=5$, our simulations can reproduce the SFH of Milky Way–sized galaxies, suggesting that the galaxy-wide star formation activity is well suppressed with our feedback implementation. Another commonly used global representation of star formation is the KSR, which is the relationship between the gas surface density and the surface density of the SFR.

To derive the KSR for molecular gas in our simulations, we split the main galaxy disk into several 1 × 1 kpc square areas and calculate ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ and ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ within each area. Here ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ is estimated by using clusters younger than 20 Myr. Varying this timescale in the range of 15–50 Myr does not affect the results. We considered different spatial smoothing scales and found that using smaller scales between 100 and 500 pc produces a larger scatter of the KSR but with similar median values to the 1 kpc case. We collect measurements from z ≈ 9 to the last available snapshot in each run and calculate the median value of ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ in a given range of ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$. Figure 6 shows the KSR in six simulations with different ${\epsilon }_{\mathrm{ff}}$.

Zoom In Zoom Out Reset image size

When expressed through molecular H₂ gas, the observed KSR in nearby spiral galaxies at z = 0 is consistent with linear and does not vary significantly with metallicity or galaxy type (Bigiel et al. 2008, 2011). In our simulations, we do not find a systematic redshift evolution of the KSR in various snapshots from z ≈ 9–1.5. Our simulations produce an approximately linear relation but with systematically higher ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ than observed at z = 0 by about a factor of 2–10. We also find that the simulations with higher ${\epsilon }_{\mathrm{ff}}$ tend to have somewhat higher ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ for a given ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$, but this trend is overwhelmed by the very large scatter.

The reason for this overestimate is not clear yet. It may be caused by the averaging scheme. The morphology of our simulated galaxies is irregular (see Figure 4), unlike the axisymmetric disks of low-redshift galaxies. Big gaps in the gas density distribution, caused by strong feedback, make the estimate of ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ less reliable. Interestingly, the FIRE simulations with spatial resolution comparable to ours show similar overestimate of Σ_SFR, pushing the gas depletion timescale into the range between 100 Myr and 1 Gyr. This depletion timescale, although shorter than the measurements from the local universe, is, however, not inconsistent with high-z results (Sharon et al. 2013; Tacconi et al. 2013; Rawle et al. 2014; Hodge et al. 2015; Sharda et al. 2018).

Figure 7 shows the initial bound fraction predicted for our model clusters. We find a strong positive correlation between f_i and stellar particle mass for all runs with a wide range of model parameters. This trend is not sensitive to either the formation epoch or host galaxy mass, suggesting that f_i is mostly independent of the global galactic environment and reflects only local properties of individual star-forming regions.

Zoom In Zoom Out Reset image size

Instead, we find that the initial bound fraction varies strongly with ${\epsilon }_{\mathrm{ff}}$. At $M={10}^{5}\,{M}_{\odot }$, f_i can reach above 50% for ${\epsilon }_{\mathrm{ff}}\geqslant 0.5$ but is limited to only 1%–10% for runs with ${\epsilon }_{\mathrm{ff}}$ < 0.5. This difference in normalization is due to the corresponding scaling of the integral star formation efficiency because of our assumption ${f}_{i}\propto {\epsilon }_{\mathrm{int}}$ (Equation (8)).

The dependence of ${\epsilon }_{\mathrm{int}}$ on stellar particle mass has already been noticed in recent GMC-scale simulations. For example, Grudić et al. (2018, hereafter G17) ran a series of MHD simulations of isolated turbulent clouds with different initial gas surface density. They found a tight relation between ${\epsilon }_{\mathrm{int}}$ and ${{\rm{\Sigma }}}_{\mathrm{gas}}$ and parameterized it by the relation

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{int}}={\left(\displaystyle \frac{1}{{\epsilon }_{\max }}+\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{crit}}}{{{\rm{\Sigma }}}_{\mathrm{gas}}}\right)}^{-1},\end{eqnarray} \tag{ 15 }$$

with best-fit parameters ${\epsilon }_{\max }=0.77$ and ${{\rm{\Sigma }}}_{\mathrm{crit}}\,=2800\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$. The critical surface density ${{\rm{\Sigma }}}_{\mathrm{crit}}$ above which the efficiency rises significantly corresponds to the cloud mass ${M}_{\mathrm{crit}}\approx 2.2\times {10}^{5}\,{M}_{\odot }{(R/5\mathrm{pc})}^{2}$. Using our definition of the initial bound fraction, the above equation can be rewritten as a relationship between f_i and scaled particle mass $m\equiv M/({\epsilon }_{\max }{M}_{\mathrm{crit}})$:

$$\begin{eqnarray}&&{f}_{i}=\min (\sqrt{{m}^{2}+4m}-m,1).\end{eqnarray} \tag{ 16 }$$

We emphasize that the radius R used in G17 is not the same as the radius of our GMC sphere, ${R}_{\mathrm{GMC}}$. In their simulations of isolated clouds, both the mass and size of the cloud are fixed by the initial condition, while our star-forming GMC has an open boundary through which gas flows in and out, controlled by gravity and stellar feedback. Our ${R}_{\mathrm{GMC}}$ is a lower limit to R, since the cluster can accrete more distant gas, e.g., from all 27 neighbor cells. We can estimate a radius of the corresponding accretion region by taking it to be a sphere of the same volume as the 27 gas cells. Since the physical size of our cells at the finest refinement level varies between 3 and 6 pc (see Section 2.2.1), we take the average length to be with 4.5 pc. This gives the effective radius R ≈ 8.4 pc. We use this value to compare the G17 scaling with ours.

The right panel of Figure 7 shows a good agreement between the relation given by Equation (16) and our runs with ${\epsilon }_{\mathrm{ff}}\geqslant 1$. The lower efficiency runs cannot reach the relatively high values of f_i predicted by G17. We fit the relation between stellar particle mass and initial bound fraction as ${f}_{i}\propto {M}^{a}$ and find the power-law slope for our clusters in the range a = 0.43–0.51. It is consistent with the G17 result a ≈ 0.5 for clusters in the mass range $M\lt {10}^{5}\,{M}_{\odot }$.

In the SFEturb run, f_i has a nonmonotonic downturn at $M\sim {10}^{5}\,{M}_{\odot }$. This downturn is caused by the relatively small number of massive clusters in this run, so that the median value of f_i is dominated by some rare cluster formation events. Indeed, these massive clusters are formed in a small satellite galaxy that contains only one compact GMC with mass ∼10⁸$\,{M}_{\odot }$ at z ∼ 4. The low ${\epsilon }_{\mathrm{ff}}$ in the SFEturb run leads to slow gas consumption, which allows the cloud to persist for a long time. Because of the high density of material in the vicinity of the GMC, young clusters acquire large velocity dispersion and leave the GMC. New clusters then appear in their stead and draw gas from the same cloud. This large cloud mass in the denominator of Equation (9) then leads to lower integrated star formation efficiency ${\epsilon }_{\mathrm{int}}$, which causes the downturn of f_i. This effect is less visible for the higher-efficiency runs, where GMCs are consumed and disrupted quicker.

The initial bound fraction in our models depends somewhat on the SNR boosting factor ${f}_{\mathrm{boost}}$. The run with ${f}_{\mathrm{boost}}=3$ shows systematically higher f_i than the corresponding ${f}_{\mathrm{boost}}=5$ run, although only by ∼20%. The dependence of f_i on ${f}_{\mathrm{boost}}$ can be understood as the balance between gravity and feedback: with lower ${f}_{\mathrm{boost}}$, a cluster in one GMC of a given mass requires more stellar mass to reach the same amount of SN momentum to fully disperse the gas. Therefore, a cluster in the SFE50-SNR3 run reaches a higher final mass for the same GMC mass, which implies higher values of ${\epsilon }_{\mathrm{int}}$ and f_i.

There are several consequences of introducing the initial bound fraction to determine the initial cluster mass. First, due to the mass dependence of f_i, the shape and normalization of the CIMF is different from the stellar particle mass function. Second, the sensitivity of f_i to the choice of ${\epsilon }_{\mathrm{ff}}$ leads to a different integrated cluster formation efficiency and maximum cluster mass for runs with different ${\epsilon }_{\mathrm{ff}}$. These observables can be used to constrain ${\epsilon }_{\mathrm{ff}}$, as we discuss below.

Because of the time-dependent accretion of gas onto GMCs and dispersal by stellar feedback, the integral efficiency ${\epsilon }_{\mathrm{int}}$ of cluster formation deviates from the adopted efficiency per freefall time ${\epsilon }_{\mathrm{ff}}$ in a given run.

Figure 8 shows the particle mass-weighted (${\epsilon }_{\mathrm{int}}M$) distributions of ${\epsilon }_{\mathrm{int}}$/${\epsilon }_{\mathrm{ff}}$ in all of our main runs. In the case of the SFEturb run, the instantaneous efficiency ${\epsilon }_{\mathrm{ff}}$ in a given cell depends on the local properties of supersonic turbulence. We calculate the mass-weighted mean efficiency $\langle {\epsilon }_{\mathrm{ff}}\rangle \approx 0.03$ for this run and use it in the leftmost panel, which shows ${\epsilon }_{\mathrm{int}}/\langle {\epsilon }_{\mathrm{ff}}\rangle$.

Zoom In Zoom Out Reset image size

There exists a weak trend that the median value of the distribution of ${\epsilon }_{\mathrm{int}}$/${\epsilon }_{\mathrm{ff}}$ decreases with ${\epsilon }_{\mathrm{ff}}$. For the SFEturb, SFE10, SFE50, SFE100, and SFE200 runs, the median values of ${\epsilon }_{\mathrm{int}}$/${\epsilon }_{\mathrm{ff}}$ are 0.52, 0.44, 0.31, 0.16, and 0.15, which correspond to the values of ${\epsilon }_{\mathrm{int}}=0.016$, 0.044, 0.155, 0.16, and 0.3, respectively. Larger ${\epsilon }_{\mathrm{ff}}$ leads to faster initial growth of a given cluster, which produces more intense stellar feedback. Stronger feedback removes more material from the GMC and leads to larger deviation of ${\epsilon }_{\mathrm{int}}$ from ${\epsilon }_{\mathrm{ff}}$.

Similar to Paper I, all distributions show a large scatter of ${\epsilon }_{\mathrm{int}}$. The 10th–90th percentile range around the median value is ≈0.7 dex for all of the runs, except the SFEturb run, which has a larger range of 1.1 dex. The equivalent 1σ scatter is about 0.25 dex for the fixed ${\epsilon }_{\mathrm{ff}}$ runs and 0.43 dex for the SFEturb run. The larger scatter in the latter run is caused by the additional variation of instantaneous ${\epsilon }_{\mathrm{ff}}$ over the course of cluster growth because of the variation of local turbulence.

The CIMF is one of the key properties of young star clusters. Here we examine the CIMF of model clusters in the main halo, summed over the central galaxy and satellite galaxies. In Paper I, we showed that the Schechter function provides a good description of the shape of the cluster mass function. Here we also fit all CIMFs with a Schechter function using the maximum-likelihood method. Because the shape of the CIMF for small clusters now deviates more from a single power law, we restrict the fit only to clusters more massive than a certain minimum mass, above which the CIMF is best described by the Schechter form. The best-fit slopes, cutoff masses, and choice of minimum mass at different epochs for all runs are listed in Table 3. The typical 1σ uncertainty for the Schechter fit is around 0.1–0.15 for α and around 0.06–0.13 dex for M_cut.

Table 3.  CIMF Best-fit Parameters

Runs	${M}_{\min }/\,{M}_{\odot }$	α	${M}_{\mathrm{cut}}/\,{M}_{\odot }$
	Full-particle
SFE10	6 × 104	3.4	$1.3\times {10}^{5}$
SFE50	6 × 104	3.3	$1.6\times {10}^{5}$
SFE50-SNR3	4 × 104	3.7	$4.5\times {10}^{5}$
SFE100	4 × 104	3.1	$4.5\times {10}^{5}$
SFE200	4 × 104	3.6	NAa
	Full-cluster
SFE10	5 × 103	2.6	$3.6\times {10}^{4}$
SFE50	2 × 104	2.5	$3.1\times {10}^{5}$
SFE50-SNR3	1 × 104	2.9	$1.3\times {10}^{5}$
SFE100	2 × 104	2.5	$2.3\times {10}^{5}$
SFE200	4 × 104	3.3	NAa
	$z\approx 2.0$-cluster
SFE10	3 × 103	2.6	$2.4\times {10}^{4}$
SFE50	8 × 103	2.2	$4.0\times {10}^{4}$
SFE50-SNR3	1 × 104	2.0	$1.3\times {10}^{4}$
SFE100	1 × 104	1.7	$3.5\times {10}^{4}$
SFE200	6 × 103	1.2	$2.5\times {10}^{4}$
	$z\approx 5.3$-cluster
SFE10	$2.5\times {10}^{3}$	1.1	$1.6\times {10}^{4}$
SFE50	6 × 103	1.3	$1.3\times {10}^{5}$
SFE50-SNR3	6 × 103	1.8	$1.7\times {10}^{5}$
SFE100	5 × 103	1.4	$2.5\times {10}^{5}$
SFE200	3 × 104	1.5	$8.5\times {10}^{5}$

Note.

^aThe cutoff mass cannot be obtained, since the CIMF of all clusters in the SFE200 run is consistent with a pure power law.

Download table as:  ASCIITypeset image

In Figure 9, we show the combined CIMF of all clusters formed at all times up to the last available snapshot of a given run. Here we distinguish between the stellar particle mass (M) and the cluster mass that takes into account the initial bound fraction (${f}_{i}M$). Because f_i itself depends positively on particle mass, the power-law slope of the CIMF is in general shallower than that of the particle mass function. Among different runs, we find that the high-mass end of the CIMF is also strongly affected by the initial bound fraction. Simulations with lower ${\epsilon }_{\mathrm{ff}}$ tend to have lower f_i for a given cluster mass, thus modifying the mass function more strongly at all masses. We notice that the CIMFs here extend only to $\sim {10}^{6}\,{M}_{\odot }$ for ${\epsilon }_{\mathrm{ff}}\geqslant 0.5$ runs, rather than $\sim {10}^{7}\,{M}_{\odot }$ as in Paper I. This is mainly due to the lower SFR caused by stronger feedback implementation in the current simulations. The maximum mass is reduced to $\sim {10}^{5}\,{M}_{\odot }$ for the SFEturb and SFE10 runs because of even smaller values of f_i.

Zoom In Zoom Out Reset image size

To compare more directly with observations of young clusters formed in the same star formation episode, in Figures 10 and 11, we show the CIMF of clusters younger than 100 Myr at different epochs, z ≈ 2 and 5.3, respectively. At z ≈ 2, the main galaxy has not experienced any major mergers for more than 500 Myr. The CIMF of young clusters reflects the properties of the ISM with quiescent star formation. We find a wide range of power-law slopes for different runs, from α ≈ 2.6 for the SFE10 run to α ≈ 1.2 for the SFE200 run. There exists a systematic trend that higher ${\epsilon }_{\mathrm{ff}}$ leads to smaller α. A shallower CIMF slope means that, for a given galactic SFR, clusters with higher mass are more likely to be created. That is the reason why the runs with higher ${\epsilon }_{\mathrm{ff}}$ tend to have higher-mass tails in the CIMF.

Zoom In Zoom Out Reset image size

Figure 11 shows the CIMF of young clusters around z ≈ 5.3, when the main galaxy experiences a major merger. The overall normalization of the mass function is lower than that at z ≈ 2, because the galaxy at z ≈ 5.3 is much smaller and contains less gas. However, the maximum cluster mass at this epoch is roughly one order of magnitude higher than that at z ≈ 2. This is because the power-law slopes of the CIMF, α ≈ 1.1–1.8, are much shallower than those in the non-merger case. The cutoff masses of these CIMFs are all above ${10}^{5}\,{M}_{\odot }$, which is several times larger than those at z ≈ 2. We already saw such an enhancement of the formation of massive clusters in gas-rich galaxy mergers in the previous runs presented in Paper I. This effect is even more pronounced in the new runs with stronger feedback.

Zoom In Zoom Out Reset image size

In Paper I, we concluded that cluster formation is an environmentally dependent process where the high-mass end of the CIMF varies strongly with the star formation activity of the host galaxy. Both recent observations (Adamo et al. 2015; Johnson et al. 2016, 2017) and theoretical work (Kruijssen 2015) suggest that SFR surface density (${{\rm{\Sigma }}}_{\mathrm{SFR}}$), rather than SFR itself, better represents the intensity of star formation and physical conditions of the galactic ISM. Therefore, here we explore the effects of ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ on such cluster properties as the fraction of clustered star formation, Γ, and maximum cluster mass, M_max.

Analogous to constructing the KSR in Section 3.2, we split the disk of the main galaxy into the 1 × 1 kpc square grid. For each square area, we calculate the surface density of the SFR averaged over 20 Myr. The fraction of clustered star formation is defined as the ratio between the mass in bound clusters and the total stellar mass formed within the same time interval, 20 Myr:

$$\begin{eqnarray}&&{\rm{\Gamma }}\equiv \displaystyle \frac{\sum {f}_{i}M}{\sum M}.\end{eqnarray} \tag{ 17 }$$

We also test a larger averaging timescale from 20 to 50 Myr and find that the results are not sensitive to the choice of the timescale.

This Γ is different from the quantity we used in Paper I in two ways. First, we estimated Γ in Paper I by splitting the galactic disk into concentric circular bins. However, due to the stronger feedback implemented in this paper, no well-defined gas disk is present in the main galaxy. Instead, the disks are clumpy and asymmetric, often showing features of strong outflows. Therefore, it is not possible to find a definitive center and cylindrical-symmetric axis to construct circular bins. Second, in Paper I, we defined Γ simply as the fraction of clusters more massive than ${10}^{4}\,{M}_{\odot }$, not all of which are necessarily gravitationally bound. In this paper, with the introduction of the initial bound fraction, we have a more physical way of estimating the clustered fraction as it is defined in observations.

Figure 12 shows a positive correlation between Γ and ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ for all runs. We find that Γ changes by about one order of magnitude over the range ${{\rm{\Sigma }}}_{\mathrm{SFR}}={10}^{-3}-1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}$. At ${{\rm{\Sigma }}}_{\mathrm{SFR}}\gt 1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, Γ saturates at the maximum value set by the initial bound fraction of the representative massive clusters. A similar saturation of Γ at high ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ is also found in the analytical model of Kruijssen et al. (2012); however, the physical origin of this saturation is different. In Kruijssen et al. (2012), the saturation is caused by the "cruel cradle effect," in which young clusters formed within less-dense regions are destroyed by the tidal interaction with other star-forming regions during their embedded phase, while in our model, the saturation is caused by the mass-dependent initial bound fraction.

Zoom In Zoom Out Reset image size

It is clear that Γ also depends strongly on ${\epsilon }_{\mathrm{ff}}$. Run SFE200, with ${\epsilon }_{\mathrm{ff}}=2.0$, shows a very high Γ up to 60%, while runs SFE10 and SFEturb do not have any star-forming regions with Γ > 10%. This sensitivity makes Γ an excellent observable to constrain the choice of ${\epsilon }_{\mathrm{ff}}$.

Figure 13 shows the relationship between ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ and the maximum mass of clusters formed within a 100 Myr interval. We find a clear positive correlation going roughly as

$$\begin{eqnarray*}&&{M}_{\max }\propto {{\rm{\Sigma }}}_{\mathrm{SFR}}^{2/3}.\end{eqnarray*}$$

The power-law slope is similar to the best-fit value, ≈0.7, for the observations of maximum cluster mass described in Section 2.4.2. The similarity of the slope among all runs reveals a robust result that the high-mass end of the CIMF depends uniquely on the intensity of star formation activity. The normalization of this relation also scales monotonically with ${\epsilon }_{\mathrm{ff}}$, because of the ${\epsilon }_{\mathrm{ff}}$ dependence of the initial bound fraction.

Zoom In Zoom Out Reset image size

Note that the runs with ${\epsilon }_{\mathrm{ff}}\leqslant 0.1$ cannot produce clusters more massive than ${10}^{5}\,{M}_{\odot }$ at all considered epochs, z > 1.5. This indicates that no clusters in these runs would survive dynamical disruption to become globular clusters at the present time. This is another evidence against such low values of ${\epsilon }_{\mathrm{ff}}$.

The upper panels of Figure 14 show the cumulative distributions of cluster formation timescales ${\tau }_{\mathrm{ave}}$. We find a clear trend that the higher the ${\epsilon }_{\mathrm{ff}}$, the shorter the timescale, especially for low-mass clusters ($M\lt {10}^{5}\,{M}_{\odot }$). This means that the growth history of these clusters is dominated by the gas accretion, which is in turn controlled by ${\epsilon }_{\mathrm{ff}}$. This result, obtained after many algorithmic updates described in Section 2, is still consistent with that in Paper I.

Zoom In Zoom Out Reset image size

More quantitatively, we split the whole cluster sample into low- ($\lt {10}^{5}\,{M}_{\odot }$) and high- ($\gt {10}^{5}\,{M}_{\odot }$) mass clusters and calculate the 10%, 50%, and 90% percentiles of the distribution of ${\tau }_{\mathrm{ave}}$. The results for both low- and high-mass clusters are listed in Table 4.

Table 4.  Cluster Formation Timescales in Myr

Runs	$M\lt {10}^{5}\,{M}_{\odot }$	$M\gt {10}^{5}\,{M}_{\odot }$
	10%–50%–90%	10%–50%–90%
	${\tau }_{\mathrm{ave}}$
SFEturb	0.25–1.07–3.13	0.43–2.07–8.62
SFE10	0.14–0.57–1.83	0.96–2.12–3.85
SFE50	0.07–0.26–0.84	0.52–1.44–2.53
SFE50-SFR3	0.09–0.26–0.80	0.45–1.26–2.35
SFE100	0.06–0.19–0.60	0.32–1.43–2.42
SFE200	0.04–0.15–0.46	0.55–1.59–2.29
	${\tau }_{\mathrm{spread}}$
SFEturb	0.09–0.29–1.13	0.84–1.99–5.76
SFE10	0.11–0.28–0.64	0.84–1.61–2.67
SFE50	0.08–0.17–0.35	0.23–0.51–1.20
SFE50-SFR3	0.09–0.19–0.40	0.48–0.91–1.78
SFE100	0.08–0.15–0.30	0.28–0.53–0.90
SFE200	0.07–0.13–0.25	0.28–0.49–0.91

Download table as:  ASCIITypeset image

In general, the timescales for low-mass clusters are very short. The median values of ${\tau }_{\mathrm{ave}}$ are all below ∼1 Myr. To better understand the changes across different runs, instead of just comparing the median values, we quantitatively evaluate the differences of the probability density distributions of ${\tau }_{\mathrm{ave}}$. We introduce a "shift factor" ${f}_{\mathrm{best}}$ that would make cumulative distributions of samples s₁ and s₂ most similar to each other. Specifically, we iteratively solve for the factor ${f}_{\mathrm{best}}$ that maximizes the p-value of the Kolmogorov–Smirnov test between s₁ and ${f}_{\mathrm{best}}\,{s}_{2}$. We find that the ${\tau }_{\mathrm{ave}}$ distributions in the SFE200, SFE100, and SFE50 runs are shifted by 0.27, 0.34, and 0.45, respectively, relative to the SFE10 run. This suggests a strong anticorrelation between ${\epsilon }_{\mathrm{ff}}$ and ${\tau }_{\mathrm{ave}}$, with a scaling relation that can be best described as ${\tau }_{\mathrm{ave}}\propto {\epsilon }_{\mathrm{ff}}^{-0.45}$. This is interesting because, although the instantaneous gas accretion rate of clusters is linearly correlated with ${\epsilon }_{\mathrm{ff}}$ in Equation (1), the correlation between ${\epsilon }_{\mathrm{ff}}$ and ${\tau }_{\mathrm{ave}}$ is nonlinear.

On the other hand, the comparison of the distributions of ${\tau }_{\mathrm{ave}}$ for the SFE50 and SFE50-SNR3 runs shows that the boosting factor of SN feedback ${f}_{\mathrm{boost}}$ has a negligible effect. This is expected, since ${f}_{\mathrm{boost}}$ only controls the intensity of the momentum feedback from SNe that explode after 3 Myr, which is longer than the typical formation timescale. Feedback from SNe marks the very end of the formation of low-mass clusters.

For massive clusters ($M\gt {10}^{5}\,{M}_{\odot }$), however, the distributions of the cluster formation timescale are very similar for runs with ${\epsilon }_{\mathrm{ff}}\geqslant 0.5$. They all have a median timescale around 1.5 Myr and truncate at ${\tau }_{\mathrm{ave}}\approx 3\,\mathrm{Myr}$, the time when the first SNe explode. This suggests that the growth history of massive clusters is determined by both gas accretion and stellar feedback. However, an immediate termination of gas accretion by SNe only happens with high ${\epsilon }_{\mathrm{ff}}$. Slower star formation cannot accumulate enough massive stars and SNe to disperse the natal GMCs within the first 3 Myr. This can be seen from the long tail of ${\tau }_{\mathrm{ave}}$ in the SFE10 and SFEturb runs. Although the median value is only ≈2–3 Myr, there are many clusters with ${\tau }_{\mathrm{ave}}\gt 3\,\mathrm{Myr}$ in the SFE10 run and ${\tau }_{\mathrm{ave}}\gt 5\,\mathrm{Myr}$ in the SFEturb run. Such long timescales are inconsistent with the observed age spread of young star clusters in nearby star-forming regions, as discussed in Section 2.4.1.

The alternative definition of the formation timescale, ${\tau }_{\mathrm{spread}}$, shows even smaller values than ${\tau }_{\mathrm{ave}}$. The median age spread ranges from 0.1 to 0.3 Myr for low-mass clusters and from 0.5 to 2.0 Myr for high-mass ones (see Table 4). However, the cumulative distribution of ${\tau }_{\mathrm{spread}}$ for the SFEturb run still has a long tail toward a large age spread, which is again inconsistent with the observations. In general, ${\tau }_{\mathrm{spread}}$ correlates with ${\epsilon }_{\mathrm{ff}}$ but less strongly than ${\tau }_{\mathrm{ave}}$.

In Section 3.1, we found that the SFH of the main galaxy in runs with ${f}_{\mathrm{boost}}=5$ is consistent with the abundance matching results. This SFH is very sensitive to the choice of ${f}_{\mathrm{boost}}$, which gives a strong constraint on its value: $3\lt {f}_{\mathrm{boost}}\lt 10$.

On the other hand, the global SFH depends little on the value of ${\epsilon }_{\mathrm{ff}}$ over the whole two orders of magnitude range. The inefficiency of star formation on large scales is not set by hand using small values of instantaneous ${\epsilon }_{\mathrm{ff}}$. Rather, as proposed recently by Semenov et al. (2017), the inefficiency may come from multiple star formation cycles on small scales, where stellar feedback processes quickly destroy individual star-forming regions. As long as the feedback is strong enough to trigger these star formation–feedback cycles, the global properties of galaxies and the inefficient star formation activity can be reproduced, and the outcome is not sensitive to the detailed implementation. This phenomenon is also found in high-resolution cosmological simulations by Hopkins et al. (2014).

In contrast to the good match of the SFH, the normalization of the KSR in all of our runs is overestimated by a factor of 2–20. In the analysis, we used the 20 Myr averaging timescale to calculate ${{\rm{\Sigma }}}_{\mathrm{SFR}}$. We chose this short timescale because we found that the structure of the ISM in our simulations is vulnerable to the feedback from young clusters, so that even the kpc-scale gas distribution changes very rapidly. Also, since we used the high density threshold for cluster formation, ${n}_{\mathrm{th}}={10}^{3}\,{\mathrm{cm}}^{-3}$, the star-forming regions are concentrated in the densest parts of the galaxy. Therefore, the $1\times 1\,\mathrm{kpc}$ grid used to derive the KSR sometimes contains only a couple of star-forming complexes. This situation is very different from the observations of low-redshift galaxies, where comparable areas contain many small star-forming regions, and both ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ and ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ are averaged over different phases of star formation.

Analogous to the instantaneous SFR in Equation (1), higher ${\epsilon }_{\mathrm{ff}}$ runs tend to have higher ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ for a given ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$. However, the slope of the molecular KSR is nearly linear and consistent with the observed linear relation, despite the nonlinear scaling ${\dot{\rho }}_{* }\propto {\rho }_{\mathrm{gas}}^{3/2}$ used in Equation (1).

In Section 3.5, we showed that the power-law slope of the CIMF is sensitive to the choice of ${\epsilon }_{\mathrm{ff}}$. Observations of the CIMF in nearby galaxies reveal a Schechter-like function with a power-law slope close to −2 (Portegies Zwart et al. 2010). The SFE10 run shows the slope as steep as ≈2.5, while the SFE200 run has a very shallow slope smaller than 1.5. Therefore, under the current simulation setup, we find that the range ${\epsilon }_{\mathrm{ff}}=0.5\mbox{--}1.0$ matches the observations best.

The normalization of the CIMF also depends strongly on ${\epsilon }_{\mathrm{ff}}$. Because of the positive correlation between ${\epsilon }_{\mathrm{ff}}$ and the initial bound fraction, runs with lower ${\epsilon }_{\mathrm{ff}}$ create clusters that have lower f_i for a given particle mass. Therefore, applying the initial bound fraction shifts the particle mass function to smaller masses and leads to the apparent lower normalization of the CIMF for lower ${\epsilon }_{\mathrm{ff}}$ runs.

This conclusion relies on the adopted linear model for the initial bound fraction, Equation (8), which is an extrapolation of a few simulation results available in the literature (Kruijssen 2012). If the true relation is not linear or has a large scatter, as indicated recently by Goodwin (2009), Dale et al. (2014), and Gavagnin et al. (2017), the derived CIMF would be affected. The overall galaxy dynamics is unaffected by the model for f_i because we assume that the bound and unbound components of a given cluster move together.

In Section 3.5, we find that there is a significant difference in CIMFs produced during major mergers and quiescent periods. In a merger, the slope of the CIMF becomes shallower, while the truncation mass rises much higher. This promotes the formation of most massive clusters, $M\gtrsim 3\times {10}^{5}\,{M}_{\odot }$. It is similar to the result we found in Paper I, despite the many updates to the algorithm in this paper.

The existence of more than 100 globular clusters in the Milky Way at present provides an additional constraint on the creation of most massive clusters at high redshift. As can be seen in Figures 9 and 13, clusters with masses larger than $\sim 3\times {10}^{5}\,{M}_{\odot }$, which are potential globular clusters progenitors, can only be created in the runs with ${\epsilon }_{\mathrm{ff}}\geqslant 0.5$.

In Paper I, we presented a positive correlation between ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ and the fraction of young clusters with mass above ${10}^{4}\,{M}_{\odot }$. Roughly speaking, this fraction is a reasonable proxy for the fraction of clustered star formation that we estimated in Section 3.3. This is because massive clusters tend to have higher f_i, and therefore, most of the bound cluster mass is contributed by the most massive clusters. In Paper I, the shape of the particle mass function was not affected strongly by the value of ${\epsilon }_{\mathrm{ff}}$, and similarly, the fraction of massive clusters was not sensitive to the choice of ${\epsilon }_{\mathrm{ff}}$. In the new runs, however, the strong correlation between ${\epsilon }_{\mathrm{ff}}$ and the initial bound fraction f_i leads to a positive correlation between ${\epsilon }_{\mathrm{ff}}$ and the integrated cluster fraction, Γ. Thus, Γ as a function of ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ provides us with a new diagnostic to constrain ${\epsilon }_{\mathrm{ff}}$.

The positive correlation between ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ and Γ found in Section 3.6 is caused by two effects. First, in high-${{\rm{\Sigma }}}_{\mathrm{SFR}}$ areas, the CIMF is likely to extend to higher masses, as demonstrated in Figure 13. Second, high-mass clusters typically have larger initial bound fractions, which leads to higher Γ. Therefore, this correlation comes about from a combination of the variation of CIMF in different environments and the mass-dependent initial bound fraction.

We find that the simulations with ${\epsilon }_{\mathrm{ff}}=0.5\mbox{--}1.0$ reproduce the observed values of Γ over a wide range of ${{\rm{\Sigma }}}_{\mathrm{SFR}}\,={10}^{-3}-1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}$. However, there are three data points from Adamo et al. (2015) that show very high ${\rm{\Gamma }}\sim 0.5$ at ${{\rm{\Sigma }}}_{\mathrm{SFR}}\sim 1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}$, which cannot be reached by either the SFE50 or SFE100 run. It is important to keep in mind that these three data points are from galaxy-integrated rather than spatially resolved measurements, which we did in this paper. Note also that the data compilation shown in Figure 12 combines measurements of different SFR tracers, cluster identification criteria, and averaging spatial and temporal scales. Future observations with consistent methodology and higher spatial resolution are required to place better constrains on the star formation models.

In Section 3.7, we find a positive correlation between the SFR surface density and maximum mass of young clusters. The normalization of this relation scales with ${\epsilon }_{\mathrm{ff}}$. We find that runs with ${\epsilon }_{\mathrm{ff}}\geqslant 0.5$ are consistent with the observations of clusters in nearby galaxies over a large range of ${{\rm{\Sigma }}}_{\mathrm{SFR}}={10}^{-4}-2\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}$. Due to the small size of observational samples, it is hard to constrain ${\epsilon }_{\mathrm{ff}}$ into a narrower range, but at least ${\epsilon }_{\mathrm{ff}}\leqslant 0.1$ is ruled out.

In Section 2.2.10, we discussed how the two definitions of the cluster formation timescale probe different mass accretion behavior of the clusters.

First, we check whether the mass accretion history can be described by a simple power law, $\dot{M}\propto {t}^{\alpha }$, as suggested by Murray & Chang (2015). Such power-law mass accretion gives ${\tau }_{\mathrm{ave}}/{\tau }_{\max }=(\alpha +1)/(\alpha +2)$ and ${\tau }_{\mathrm{spread}}/{\tau }_{\max }\,={(2\alpha +1)/(\alpha +1)}^{2}$. Therefore, if the accretion history of individual clusters can indeed be described by a power law, the indexes α derived from the ratios ${\tau }_{\mathrm{ave}}$/${\tau }_{\max }$ and ${\tau }_{\mathrm{spread}}$/${\tau }_{\max }$ should be consistent with each other. However, we find that majority of the clusters do not show consistent α from these two definitions, suggesting that most clusters experience accretion histories that cannot be described by a simple power law.

Then we explore the detailed output of the mass accretion rate at each local timestep (∼100 yr) for a fraction of massive clusters above ${10}^{5}\,{M}_{\odot }$. We find a general pattern that the accretion rate is relatively low at the beginning of cluster formation, when the gas density is near the threshold ${n}_{\mathrm{crit}}$. As the parent GMC collapses to higher density, the central gas cell where the cluster resides is refined to one or two higher levels, which increases the mass growth rate. We find that this collapse phase typically lasts for ∼2/3 of the total cluster formation duration. Such superlinear mass growth is qualitatively consistent with the theoretical prediction of the collapse of turbulence clouds (Murray & Chang 2015; Murray et al. 2017). As the mass accretion reaches its peak by either exhausting the star-forming gas or removing it by stellar feedback, the density quickly declines, and, once it is below the threshold ${n}_{\mathrm{crit}}={10}^{3}\,{\mathrm{cm}}^{-3}$, the cluster growth is shut down completely. This pattern clearly cannot be fitted by a single power law of time, and it is better described by a peaked distribution. In this case, ${\tau }_{\mathrm{ave}}$ and ${\tau }_{\mathrm{spread}}$ reflect distinct characteristics of the accretion of individual clusters. While ${\tau }_{\mathrm{ave}}$ indicates the timescale for a given cluster to reach the peak of its mass accretion, ${\tau }_{\mathrm{spread}}$ quantifies the width of that peak.

Similar to the results of Paper I, in the new runs, the cluster formation timescale depends strongly on ${\epsilon }_{\mathrm{ff}}$. This robust relationship could be used to constrain ${\epsilon }_{\mathrm{ff}}$. In practice, however, the clusters in most runs have formation timescales that are shorter than 3 Myr, so that they all are within the range of the observed age spread. One exception is the SFEturb run, in which ∼30% of massive clusters have ${\tau }_{\mathrm{ave}}\gt 6\,$ Myr, which is inconsistent with observations.

Here we summarize the constraints for the star formation and feedback parameters, ${\epsilon }_{\mathrm{ff}}$ and ${f}_{\mathrm{boost}}$, resulting from the different observables discussed above.

• 1.  
  Global SFH: $3\lt {f}_{\mathrm{boost}}\lt 10$.
• 2.  
  Slope of CIMF: ${\epsilon }_{\mathrm{ff}}=0.5\mbox{--}1.0$.
• 3.  
  Fraction of clustered star formation: ${\epsilon }_{\mathrm{ff}}=0.5\mbox{--}2.0$.
• 4.  
  Maximum cluster mass: ${\epsilon }_{\mathrm{ff}}\geqslant 0.1$.
• 5.  
  Cluster formation timescale: ${\epsilon }_{\mathrm{ff}}\geqslant 0.1$.

Based on the joint constraints, we find that ${f}_{\mathrm{boost}}\approx 5$ and ${\epsilon }_{\mathrm{ff}}=0.5\mbox{--}1.0$ are favored in the current framework of star formation and stellar feedback subgrid models and the current spatial and mass resolution.

We emphasize again that the ${\epsilon }_{\mathrm{ff}}$ used in this paper is different from that in other cosmological simulations: (1) it is only applied within a star-forming sphere of the fixed physical volume, rather than a time-variable cell volume; and (2) the mass growth rate of a given cluster varies significantly during the active accretion period due to changes of the local gas density, which leads to broad scatter of the integrated efficiency ${\epsilon }_{\mathrm{int}}$ at a fixed ${\epsilon }_{\mathrm{ff}}$.