On the Observed Diversity of Star Formation Efficiencies in Giant Molecular Clouds

TL;DR

This study uses hydrodynamical simulations to investigate the wide diversity of star formation efficiencies in GMCs, attributing it mainly to the unique evolutionary paths of individual clouds rather than global properties.

Contribution

It demonstrates that the observed scatter in star formation efficiency is due to diverse cloud evolutionary histories, not global cloud properties, using realistic simulations.

Findings

Efficiencies follow a wide distribution consistent with observations.

No strong correlation between efficiency and global cloud properties.

Unique evolutionary paths lead to diversity in efficiencies.

Abstract

Observations find a median star formation efficiency per free-fall time in Milky Way Giant Molecular Clouds (GMCs) on the order of $\epsilon_{\rm ff}\sim 1\%$ with dispersions of $\sim0.5\,{\rm dex}$. The origin of this scatter in $\epsilon_{\rm ff}$ is still debated and difficult to reproduce with analytical models. We track the formation, evolution and destruction of GMCs in a hydrodynamical simulation of a Milky Way-like galaxy and by deriving cloud properties in an observationally motivated way, measure the distribution of star formation efficiencies which are in excellent agreement with observations. We find no significant link between $\epsilon_{\rm ff}$ and any measured global property of GMCs (e.g. gas mass, velocity dispersion). Instead, a wide range of efficiencies exist in the entire parameter space. From the cloud evolutionary tracks, we find that each cloud follow a \emph{unique} evolutionary path which gives rise to wide diversity in all properties. We argue that it is this diversity in cloud properties, above all else, that results in the dispersion of $\epsilon_{\rm ff}$.

Figures

Equations

$${\dot{\rho }}_{\star }={ϵ}_{\mathrm{ff},\mathrm{SF}}{f}_{{\mathrm{H}}_2}\frac{{\rho }_{\mathrm{g}}}{t_{\mathrm{ff}}},$$

$${ϵ}_{\mathrm{ff}}=\frac{t_{\mathrm{ff}}}{t_{\star ,\mathrm{y}}}\frac{M_{\star ,\mathrm{y}}}{(M_{\mathrm{GMC}}+M_{\star ,\mathrm{y}})}$$

$${ϵ}_{\mathrm{ff},\lim }=\frac{t_{\mathrm{ff}}}{t_{\star ,\mathrm{y}}}\frac{1}{(1+M_{\mathrm{GMC}}/300{\mathrm{M}}_{\odot })},$$

$$\frac{d{M}_{\mathrm{GMC}}}{dt}=-{ϵ}_{\mathrm{ff},0}{\left(\frac{t_{\mathrm{age}}}{t_{\mathrm{ff},\mathrm{fixed}}}\right)}^{\delta }\frac{M_{\mathrm{GMC}}(t_{\mathrm{age}})}{t_{\mathrm{ff},\mathrm{fixed}}}-\alpha M_{\star ,\mathrm{tot}}(t_{\mathrm{age}})+\gamma ,$$

$$\frac{d{M}_{\star ,\mathrm{tot}}}{dt}={ϵ}_{\mathrm{ff},0}{\left(\frac{t_{\mathrm{age}}}{t_{\mathrm{ff},\mathrm{fixed}}}\right)}^{\delta }\frac{M_{\mathrm{GMC}}(t_{\mathrm{age}})}{t_{\mathrm{ff},\mathrm{fixed}}},$$

$$M_{\star ,\mathrm{y}}(t_{\mathrm{age}})=M_{\star ,\mathrm{tot}}(t_{\mathrm{age}})-M_{\star ,\mathrm{tot}}(t_{\mathrm{age}}-t_{\star ,\mathrm{y}}),$$

$${ϵ}_{\mathrm{ff},L16}(t_{\mathrm{age}})=\frac{t_{\mathrm{ff},\mathrm{fixed}}}{t_{\star ,\mathrm{y}}}\frac{M_{\star ,\mathrm{y}}(t_{\mathrm{age}})}{M_{\mathrm{GMC}}(t_{\mathrm{age}})+M_{\star ,\mathrm{tot}}(t_{\mathrm{age}})},$$

Full text

1 Introduction

It is within Giant Molecular Clouds (GMCs) that galaxies form the vast majority of their stars (Myers et al., 1986; Shu et al., 1987; Scoville & Good, 1989; McKee & Ostriker, 2007). In local spiral and dwarf galaxies, star formation on galactic scales is known to be a slow process (e.g. Bigiel et al., 2008), with gas depletion time scales on the order of billions of years. This inefficiency is also found on scales of individual GMCs (Krumholz & Tan, 2007), with a median star formation efficiencies per cloud free-fall time $\epsilon_{\rm ff}\sim 1\%$ (Myers et al., 1986; Murray, 2011; Krumholz et al., 2012). While most observations of $\epsilon_{\rm ff}$ in GMCs find very similar median value, different surveys find spreads in $\epsilon_{\rm ff}$ of different sizes. Using the median absolute deviation to robustly estimate the standard deviations in $\epsilon_{\rm ff}$ the observational data presented in Evans et al. (2014); Heyer et al. (2016); Lee et al. (2016); Vutisalchavakul et al. (2016); Ochsendorf et al. (2017) and Utomo et al. (2018) yields values ranging from $\sim 0.21$ to $\sim 0.83\,{\rm dex}$. In their recent review Krumholz et al. (2018) compiled data from 13 papers over the last decade (see their Fig. 10) and discuss the impact of the method used to estimate $\epsilon_{\rm ff}$ on measurements. The exact size and distribution of the spread in $\epsilon_{\rm ff}$ is therefore still a debated topic, particularly as the origins of this diversity is not yet understood.

As the dispersion in $\epsilon_{\rm ff}$ is found independently of the method of observation, it is likely physical. Analytical models of star formation in supersonic turbulent flows (e.g. Padoan & Nordlund, 2011; Hennebelle & Chabrier, 2011; Federrath & Klessen, 2012) have successfully explained the low mean $\epsilon_{\rm ff}$ in GMCs and provided insight into how it can scale with global cloud properties, such as density and virial parameter. Krumholz & McKee (2005), for example, postulate that GMCs are turbulent gas structures which are characterised by a log-normal distribution (determined by the Mach number) and only regions within the cloud with density above some threshold are able to form stars. In their model the threshold density is determined by both the virial parameter of the cloud and its Mach number. Lee et al. (2016) (henceforth L16, see also Ochsendorf et al., 2017) have been critical of such models for failing to predict a sufficiently large dispersion in $\epsilon_{\rm ff}$. Furthermore, observations indicate a decreasing $\epsilon_{\rm ff}$ with increasing cloud mass ($M_{\rm GMC}$). Ochsendorf et al. (2017) argued that this is a result of massive clouds having a diffuse, non-star forming outer envelopes. This observed $M_{\rm GMC}$-$\epsilon_{\rm ff}$ relationship presents an additional constraint on any model attempting to explain the distribution in $\epsilon_{\rm ff}$. However it is worth noting that the $M_{\rm tot}$-$\epsilon_{\rm ff}$ relation could be, at least in part, the result of limitations in current observations to detect low mass clouds.

Feldmann & Gnedin (2011) developed a toy model where the growth of stellar mass in a GMC is determined only by the mass and free-fall time of available gas and an intrinsic efficiency per free-fall time ($\epsilon_{\rm ff,0}$), meanwhile the change in gas mass results from a combination of mass converted into stars, gas removed by feedback and accretion of new gas from the surrounding environment. Using this model they found that a $\sim 2$ dex spread in the measured $\epsilon_{\rm ff}$ can be explained by cloud evolution by adopting a constant input $\epsilon_{\rm ff}$. L16 explored the same effect but found that a fixed input $\epsilon_{\rm ff}$ fails to reproduce their observed distribution of $\epsilon_{\rm ff}$, with too few clouds predicted at high ($\gtrsim 10\%$) and low ($\lesssim 0.01\%$) efficiencies at any given time. By assuming that all GMCs evolve in a similar fashion, but observed at a different stage of evolution and by allowing for a time dependent $\epsilon_{\rm ff}$, this problem was mitigated (see their figure 8). While each GMC observed is at a different stage of its life, it is currently unclear if all clouds follow the same evolution, and how (or if) the galactic environment plays a role.

A number of authors have studied the dispersion in $\epsilon_{\rm ff}$ using simulations of individual GMCs as well as global disc simulations: the former looking for an explanation in the internal properties of clouds, while the latter allows for the impact of environment to be studied. Semenov et al. (2016) studied the impact of an explicit treatment of small scale gas turbulence, using simulations of entire Milky Way-like galaxies, on the parsec-scale star formation efficiency. From the local properties of the gas, they used the simple parametrisation of Padoan et al. (2012), calibrated on magneto-hydrodynamical simulations of star formation in supersonic turbulence, to compute $\epsilon_{\rm ff}$. They found that their simulations produced values of $0.01\%\lesssim\epsilon_{\rm ff}\lesssim 10\%$. While this is an encouraging result, we will in this work demonstrate that the mapping between the $\epsilon_{\rm ff}$ computed from local properties and what is actually derived from observations is complex and depends on the star formation history of the cloud, not just its instantaneous properties.

Recently Grudić et al. (2018) carried out 17 magnetohydrodynamic simulations of isolated GMCs of varying mass, radii and feedback models but identical surface density. They found that the spread in $\epsilon_{\rm ff}$ seen in observations is similar to the spread in $\epsilon_{\rm ff}$ measured throughout a cloud’s lifetime. However, as pointed out by the authors, their GMCs have “fairly artificial” initial conditions and lack the effects of the larger galactic environment in which GMCs are found. As shown by observations (e.g. Rosolowsky et al., 2003; Heyer et al., 2009), GMCs have a large range of properties that can not be captured by 17 overlapping models. Therefore full galactic (disc) simulations, which produce self-consistent GMCs with a range of properties, are needed to better model and investigate the entire evolution of such clouds.

In this study we will go beyond previous work by using parsec resolution simulations of entire disc galaxies and investigate the emerging GMC star formation efficiencies and how they evolve. This allows us to explore whether such a diversity in $\epsilon_{\rm ff}$ and the observed mass–$\epsilon_{\rm ff}$ relation is an artefact of observational methods or a physical result. The large number of GMCs found in our simulations allows us to look for correlations between different properties of a GMC and its $\epsilon_{\rm ff}$, thus determining if a single property is responsible for the observed scatter in $\epsilon_{\rm ff}$. Furthermore, taking advantage of the high temporal resolution of our simulations, we explore how GMCs evolve over their lifetime in a number of different properties and how this contributes to the diversity in $\epsilon_{\rm ff}$. Finally, combining our simulated data with analytical models we determine whether a single model is able to explain our simulations or observations.

This paper is organised as the following: in §2 we summarise our simulations and methods for identifying and tracking GMCs, in §3 we present the measured values of $\epsilon_{\rm ff}$ and how it relates to global GMC properties, in §4 we explore the source of the dispersion in $\epsilon_{\rm ff}$ and finally we present our conclusions in §5.

2 Method

2.1 Simulations

We make use of the two Milky Way-like galactic disc simulations in Grisdale et al. (2017), henceforth G17. The simulations are identical, apart from one being run with stellar feedback (our fiducial simulation) and one without. The simulations account for a dark matter halo, stellar and gaseous disc and a bulge. The initial conditions of both simulations are identical to the AGORA disc initial conditions described in Kim et al. (2016). They were run using the hydro+$N$-body, Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier, 2002). A cell is refined if it reaches a threshold mass of $9300{\,{\rm M}}_{\odot}$ and the minimum allowed cell size is $\Delta x\sim 4.6{\,\rm pc}$.

The adopted cooling, feedback and star formation models are outlined in G17 and Grisdale et al. (2018), hence forth G18 (see also Agertz et al., 2013; Agertz et al., 2015). Briefly, the feedback model accounts for the injection of momentum, energy, mass loss and enrichment from stellar winds, supernovae (II and Ia) and radiation pressure from young stars. Star formation occurs on a cell-by-cell basis according to the star formation law:

[TABLE]

where $f_{\rm H_{2}}$ is the local mass fraction of molecular hydrogen (H2), $\rho_{\rm g}$ is the gas density, $t_{\rm ff}=\sqrt{3\pi/32G\rho_{\rm g}}$ is the local free-fall time and $\epsilon_{\rm ff,SF}$ is the local star formation efficiency per free-fall time of gas in the cell. For all star forming cells $\epsilon_{\rm ff,SF}$ is set to $10\%$ in the simulation with feedback and $1\%$ in the simulation without. As shown in G17, these choices lead to comparable galactic star formation histories in the two simulations. All star particles form with an initial mass of $300{\,{\rm M}}_{\odot}$.

G17 demonstrated that the simulation with stellar feedback give rise to a supersonically turbulent ISM, with a density and velocity structure in close agreement with local spiral galaxies. Furthermore, the resulting GMC population has masses, sizes, velocity dispersions and scaling relations (‘Larson’s relations’) closely matching that of the Milky Way (e.g. Heyer et al., 2009), as shown in G18. This makes our simulations a suitable platform for investigating the evolution and star formation properties of GMCs.

2.2 Cloud identification and analysis

We identify GMCs in two separate ways: 1) in projection (2D) using the CLUMPFIND algorithm (Williams et al., 1994) as implemented in the clump finding identification and analysis package CUPID111part of the Starlink Project (see Manset & Forshay, 2014; Starlink, 2015, for details) and 2) in 3D using the on-the-fly clump finding module PHEW (Parallel HiErarchical Watershed, Bleuler et al., 2015) built into RAMSES. The methods yield similar distributions of GMC properties, albeit with the 3D method giving slightly larger masses and sizes. The adopted cloud finding parameters and the resulting GMC properties are discussed in detail in G18. We adopt both methods when comparing simulated GMC star formation efficiencies in §3.1 but restrict all other analysis in this work to the 3D approach. This allows for higher time resolution, hence allowing for cloud tracking as well as better statistics (due to the number of clouds identified).

Throughout this work, data from the 2D clump finder were obtained from simulation snapshots between $t=150-450\,{\rm Myr}$, separated by $\Delta t=25\,{\rm Myr}$. For the 3D clump finder, data was obtained at $t=325-380\,{\rm Myr}$, with a temporal spacing of, on average, $\Delta t\sim 25,000$ years (see G17 for more details). Clouds that lie within the central kiloparsec of the galaxy are removed from the analysis. In total, during the period of analysis we identify 8,201 (3,434) GMCs using 2D clump finding in the simulation with(out) feedback and 655,499 (212,056) clouds with the 3D clump finder.

To accurately compare simulations to observations (e.g. L16), we use the estimator

[TABLE]

for each GMC, where $M_{\rm\star,y}$ is the mass of stars with the GMC which has an age less than $t_{\rm\star,y}$, $M_{\rm GMC}$ is the (molecular) gas mass of the GMC, $t_{\rm ff}=\sqrt{3\pi/32G\rho_{\rm GMC}}$ is the mean free fall time across the GMC and $\rho_{\rm GMC}$ is $M_{\rm GMC}$ divided by the GMC’s volume (i.e. its mean density). Stellar masses are calculated by considering only stars that overlap with gas belonging to a GMC222Other methods for matching stars to GMCs were explored and found to have little impact on the results presented in this work., either in projection (2D method) or in 3D.

To allow for a comparison to the stellar clusters detected by free–free emission (e.g. L16), we adopt $t_{\rm\star,y}=4\,{\rm Myr}$. We emphasise that $\epsilon_{\rm ff}$ should not be confused with $\epsilon_{\rm ff,SF}$ in Eq. 1. The former is the efficiency per free fall time averaged over the whole GMC ($\sim 10-70{\,\rm pc}$), while the latter is the efficiency at which an individual computational cell converts gas into stars ($\sim 4.6{\,\rm pc}$).

2.3 Cloud tracking

All GMC quantities are followed over time by employing the cloud tracking algorithm described in Tasker & Tan (2009). Briefly, for clouds found using the 3D333We only apply the tracking routine to clouds identified with the 3D clump finding method, as the time resolution in the 2D method is insufficient for accurate tracking., we use the position ($\mathbf{x}$) and velocities ($\mathbf{v}$), of each GMC at a time $t_{0}$, we adopt a linear approximation, where the change in position vector over one tracking step is $\Delta\mathbf{x}=\mathbf{v}\Delta t$, to predict where the cloud should be at the next cloud finder output, $t_{1}=t_{0}+\Delta t$. Next, separations ($\mathcal{S}$) between a cloud’s predicted position at $t=t_{1}$ and the positions of all clouds found at this time are calculated. Cloud identities are matched for clouds with the smallest $\mathcal{S}$ and satisfying either $\mathcal{S}\leq 2R_{{\rm GMC},t_{0}}$ or $\mathcal{S}\leq\langle R_{{\rm GMC},t_{1}}\rangle$, where $R_{{\rm GMC},t_{0}}$ is the radius of the cloud at $t_{0}$ and $\langle R_{{\rm GMC},t_{1}}\rangle$ is the mean cloud radius at $t_{1}$. In the case of multiple clouds from $t_{0}$ being linked to the same cloud at $t_{1}$, the cloud at $t_{1}$ inherits the identity of the most massive cloud from $t_{0}$, while the other progenitor cloud(s) are considered to have merged and are not tracked further.

To ensure that only clouds with a complete life-cycle are considered we exclude those formed in (or before) the first snapshot of our analysis, or those destroyed after the last.

Furthermore, it is important to note that because we are detecting clouds based on a fixed density threshold ($100{\,\rm{cm^{-3}}}$, see G18) it is possible for a cloud to drop below the detection limit but remain a coherent structure and then, at later times, pass back above the threshold. In such situations, our methods would register the cloud as having been destroyed and a new cloud forming. A “new” cloud of this type may be detected with a significant young stellar mass ($M_{\rm\star,y}$), hence giving the appearance of beginning its life with a high $\epsilon_{\rm ff}$. To mitigate this, we only considered GMCs with an initial $M_{\rm\star,y}\leq 1500{\,{\rm M}}_{\odot}$, equivalent to five (or less) star particles. An alternative method, which would better reflect the complex cycle of gas ending up in GMCs, would be to employ tracer particles to track the gas of each GMC (e.g. Semenov et al., 2016), or to identify stellar clusters and their associated molecular gas. We leave this for a future investigation.

Finally, clouds with lifetimes shorter than a million years are also discarded. After the tracking is complete, and the above criteria are applied, 1879 unique GMCs evolutionary tracks remain, which we focus on in §3.2.

3 Results

3.1 Distribution of star formation efficiencies

3.1.1 The GMC Mass - Star Formation Efficiency Relation

We begin our analysis by calculating $\epsilon_{\rm ff,2D}$ (Eq. 2) for the clouds identified in projection, as outlined in §2.3. In the left hand panel of Fig. 1 we show how $\epsilon_{\rm ff,2D}$ varies with total cloud mass ($M_{\rm tot}=M_{\rm GMC}+M_{\rm\star,y}$, as defined in L16) in our fiducial simulation (i.e. including feedback) and compare these to the observational data in L16. The simulated GMC population has a wide range of $\epsilon_{\rm ff,2D}$ values which agree well with observations, see Fig.1.

Both simulated and observed $\epsilon_{\rm ff,2D}$ decrease with increasing $M_{\rm tot}$ (see also Ochsendorf et al., 2017), raising the question as to whether this is due to a physical process, or a result of the cloud identification method. In our simulation $M_{\rm\star,y}=N\cdot M_{\star}$, where $M_{\star}=300{\,{\rm M}}_{\odot}$444Particles in the feedback simulation lose mass, however during the first $4\,{\rm Myr}$ of their evolution they only lose a maximum of $10\%$ of their initial mass. is the mass resolution of the star particles and $N$ is the number of star particles in the GMC. This defines a lower limit of the estimated $\epsilon_{\rm ff,2D}$ in our simulation,

[TABLE]

below which our simulation cannot sample star formation. This ‘sampling limit’ is shown in Fig. 1 for $t_{\rm ff}=t_{\rm\star,y}$, illustrating how this introduces a bias in how $\epsilon_{\rm ff,2D}$ relates to $M_{\rm tot}$. Star clusters identified via free–free emission, as done by L16, have a similar bias; such identification requires the presence of UV-emitting massive stars to ionise the surrounding ISM, hence setting a lower limit to the detectable star formation efficiency for all cloud masses, which scales in a similar fashion as Eq. 3 (see also Murray, 2011; Grudić et al., 2018). This is likely a contributing factor as to why the simulation and observations agree on the low $\epsilon_{\rm ff,2D}$ end of the distribution. Indeed, with the exception of a single GMC,555The exception, found with $\epsilon_{\rm ff}=1.8\times 10^{-4}$ and $M_{\rm tot}=2.8\times 10^{4}{\,{\rm M}}_{\odot}$, only has $16{\,{\rm M}}_{\odot}$ in stellar mass. the observed GMCs shown in Fig. 1 (from L16) have $M_{\rm\star,y}\gtrsim 100{\,{\rm M}}_{\odot}$, close to the mass resolution of the simulation.

3.1.2 Dispersion of Star Formation Efficiencies

To quantify the spread of $\epsilon_{\rm ff,2D}$ we show the normalised distribution, the median value of $\epsilon_{\rm ff,2D}$ ($\langle\epsilon_{\rm ff,2D}\rangle_{\rm med}$) and the standard deviation of $\log{\epsilon_{\rm ff,2D}}$, $\sigma_{\log{\epsilon_{\rm ff,2D}}}$666The standard deviation ($\sigma_{\log{\epsilon_{\rm ff,2D}}}$) is robustly estimated via the median absolute deviation (MAD): $\sigma_{\log{\epsilon_{\rm ff,2D}}}=1/0.6745$ MAD (Müller, 2000; Romeo & Fathi, 2016)., for both the simulation and L16’s observations in the middle panel of Fig. 1.

The shape of the $\epsilon_{\rm ff,2D}$ distribution in our simulation is in excellent agreement with observations. We find a median star formation efficiency per free-fall time of $\sim 2~{}\%$, also matching observations. The simulated $\sigma_{\log{\epsilon_{\rm ff,2D}}}$ is not as large a value as L16’s data ($0.72\,{\rm dex}$ compared to $0.83\,{\rm dex}$) but is within the range of values reported in the literature (see §1). This is likely due to the simulated population of GMCs being better sampled at low masses compared to observation (see also G18), where the spread in $\epsilon_{\rm ff,2D}$ is smaller, as well the observed population having a couple of extremely inefficiently star forming clouds ($\epsilon_{\rm ff}\sim 10^{-4}$). From Fig. 2, which shows how $\langle\epsilon_{\rm ff,2D}\rangle_{\rm med}$ and $\sigma_{\log{\epsilon_{\rm ff,2D}}}$ vary with $M_{\rm tot}$, we find that $\langle\epsilon_{\rm ff,2D}\rangle_{\rm med}$ in the simulation is compatible with observations at almost all cloud masses: low mass clouds ($\sim$ few $\times 10^{4}{\,{\rm M}}_{\odot}$) reaching almost $\epsilon_{\rm ff}\sim 10\%$, whereas clouds in excess of $10^{6}{\,{\rm M}}_{\odot}$ have $\epsilon_{\rm ff}\lesssim 1\%$, on average. This result can be summarised as $\epsilon_{\rm ff}\propto M_{\rm tot}^{\beta}$ with $\beta\sim-0.34$ and $-0.36$ for our simulation and L16’s observations respectively. Finally, a qualitative agreement is found for $\sigma_{\log{\epsilon_{\rm ff,2D}}}-M_{\rm tot}$, with an increasing scatter with increasing cloud mass.

It is important to note that Fig. 2 and the values of $\beta$ given above do not account for the sampling limit (see §3.1.1). As shown by the L16’s data in Fig. 1, clouds below this sampling limit do exist and therefore need to be accounted for when calculating true value of $\beta$. For example, if we ignore clouds in our feedback simulation that sit on or close to the sampling limit we find $\beta$ decreases, i.e. the relationship become steeper. The steeping in due to the fact that sampling limit preferentially impacts lower mass ($<10^{5}{\,{\rm M}}_{\odot}$ ) clouds, indeed $\langle\epsilon_{\rm ff,2D}\rangle_{\rm med}$ at these masses is only $\sim 1\sigma_{\log{\epsilon_{\rm ff,2D}}}$ from the sampling limit. Given that $\beta$ will be biased in both simulation and observation, determining the impact of the sampling limit is of significant importance. One way to assess the how the sampling limit biases measurements of $\beta$ would be to compare data with a simple Bayesian model. Such a comparison is beyond the scope of this study and we leave it for future work.

We calculate the efficiency per free-fall time ($\epsilon_{\rm ff}$) for GMCs identified using the 3D clump finder and repeat the above analysis, see right most panel of Fig. 1. We find similar median efficiencies as before ($\langle\epsilon_{\rm ff}\rangle_{\rm med}\sim 1\%$), but the spread in $\epsilon_{\rm ff}$, which still covers several orders of magnitude, is smaller: $\sigma_{\log{\epsilon_{\rm ff}}}=0.37$.

In summary, the measured star formation efficiency per free-fall time in a cloud tells us little, if anything at all, about the input efficiency on smaller scales, i.e. $\langle\epsilon_{\rm ff}\rangle_{\rm med}\neq\epsilon_{\rm ff,SF}$. The former depends not only on the turbulent density substructure (as discussed below) and global properties of the cloud but also on its evolution and therefore the rate of star formation, that results in the young stellar population observed at any instance.

3.1.3 Role of Feedback on $\epsilon_{\rm ff}$

Fig. 1 and 2 include analysis from the simulation without feedback. We find very similar trends as before, i.e. $\epsilon_{\rm ff,2D}\propto M_{\rm tot}^{\beta}$ with $\beta\sim-0.69$ and a smaller scatter in 3D compared to 2D (i.e. $\sigma_{\log{\epsilon_{\rm ff,2D}}}=1.1\rightarrow\sigma_{\log{\epsilon_{\rm ff}}}=0.2$). As with simulation with feedback we emphasise that the value of $\beta$ given above does not account for the sampling limit and therefore may not accurately represent the correlation between $\epsilon_{\rm ff,2D}$ and $M_{\rm tot}$. Given that the simulation without feedback has a propensity to produce unphysical (see G18), massive ($M_{\rm GMC}>10^{7}{\,{\rm M}}_{\odot}$), low $\epsilon_{\rm ff,2D}(\lesssim 10^{-5})$ clouds, we expect that clouds will be below the sampling limit at all masses, but particularly at high mass (see left panel of Fig. 1) which makes predicting the impact of the sampling limit on $\beta$ much more difficult in this case.

The primary difference between the GMCs from the simulation with feedback and the simulation without is the latter has a population of long lived massive ($\gtrsim 10^{7}{\,{\rm M}}_{\odot}$) and very inefficiently star forming ($\epsilon_{\rm ff}\lesssim 10^{-3}$) clouds. Interestingly, we find that in 3D the simulation with feedback yields GMCs able to reach higher $\epsilon_{\rm ff}$ values than clouds in the simulation without, which is a result of stellar feedback removing gas from the clouds. From these results we conclude that, while stellar feedback plays a role in determining the shape of the distribution of measured $\epsilon_{\rm ff,2D}$ and $\epsilon_{\rm ff}$, it is not the source of the dispersion in either.

Given the reasonable match between GMCs in the simulation without feedback and the observations, it might seem just as reasonable to use this simulation as the simulation with feedback in further analysis. However, as shown in both G17 and G18, this simulation fails to produce a realistic neutral ISM and distribution of GMC properties, therefore we focus all further analysis on the simulation with feedback.

Having established that star formation efficiencies in the simulated GMC population closely match observations, we next aim to quantify why this is the case.

3.1.4 Role of Cloud Properties on $\epsilon_{\rm ff}$

From this point forward we focus our analysis to clouds identified using the 3D clump finder. Fig. 3 shows the relation between $\epsilon_{\rm ff}$ of GMCs and gas mass ($M_{\rm GMC}$), galactocentric radius ($R_{\rm gal,2D}$), velocity dispersion ($\sigma_{\rm v,GMC}$), size ($R_{\rm GMC}$), gas density ($\rho_{\rm GMC}$) and virial parameter ($\alpha_{\rm vir,GMC}$) in the simulation with feedback (we refer the reader to G18 for discussion on how these quantities are calculated).

No strong correlation is found between any cloud property and $\epsilon_{\rm ff}$; clouds with similar physical parameters show a great diversity of $\epsilon_{\rm ff}$. Not surprisingly, there is a trend for dense GMCs to have higher $\epsilon_{\rm ff}$ and likewise, as $t_{\rm ff}\propto\rho^{-0.5}$, a trend for clouds with short free-fall times to be efficient at forming stars. We find a signature of clustering in the different $\epsilon_{\rm ff}$-spaces, e.g. in the $\epsilon_{\rm ff}-M_{\rm GMC}$ space, clouds cluster around $[\epsilon_{\rm ff},M_{\rm GMC}]\sim[-1.25,10^{6}{\,{\rm M}}_{\odot}]$ and $\sim[-2.3,10^{5.4}{\,{\rm M}}_{\odot}]$, suggesting there are some preferential values that clouds are drawn to.

Given that $\epsilon_{\rm ff}$ in the former of these two regions is centred on $\epsilon_{\rm ff}\sim 5\%$ and therefore within a factor of 2 from $\epsilon_{\rm ff,SF}$, this region could be an artefact of the star formation model employed in the simulations. To test this we reproduce Fig. 3 using clouds identified in the simulation without feedback (Fig. 4 shows the $\epsilon_{\rm ff}$-$\rho_{\rm GMC}$ panel). While the vast majority ($\sim 70\%$) of clouds are highly inefficient at forming stars at higher density, i.e. $\rho_{\rm GMC}\gtrsim 10{\,{\rm M}}_{\odot}{\,\rm pc}^{-3}$, there is a tendency for $\epsilon_{\rm ff}\rightarrow\epsilon_{\rm ff,SF}$. By directly comparing the $\epsilon_{\rm ff}$-$\rho_{\rm GMC}$ panel in Fig. 3 to Fig. 4 we find that general shape of the two distributions is very similar, but that latter is missing most of the $\epsilon_{\rm ff}>10^{-1.5}$ clouds. We therefore conclude that the high $\rho_{\rm GMC}$-high $\epsilon_{\rm ff}$ region is at least partially a result of the star formation prescription employed in our simulations. In future work we explore in detail, to what the degree the star formation prescription drives clouds into this region and if other factors play a role.

Next we explore how cloud properties evolve and the role this plays in establishing the wide range of observed star formation efficiencies.

3.2 Cloud evolution

3.2.1 Individual clouds

Fig. 3 includes evolutionary tracks of eight randomly selected GMCs, with their positions shown every $0.25\,{\rm Myr}$. Each cloud has a unique path through the seven parameter spaces explored. For example, most clouds tend to become more gravitationally bound over their lifetime ($\alpha_{\rm vir,GMC}$ decreases) while the crimson cloud becomes less bound ($\alpha_{\rm vir,GMC}$ increases).

GMCs are not confined to a single area of parameter space but can move from one region to another and in a variety of ways. Comparing the black and grey clouds with either of the blue clouds shows that some clouds explore only a small fraction of a given parameter space while others might explore a large portion. The one possible exception to this is a cloud’s progression in $\epsilon_{\rm ff}$–$R_{\rm gal,2D}$ space, where we see that clouds are “born” and “die” at approximately the same galactic radius ($R_{\rm gal,2D}$), with very little variations over the clouds lifetime. This is due to the short cloud evolution timescale ($\sim 10\,{\rm Myr}$) compared to the galactic dynamical timescale ($\sim 100-200\,{\rm yr}$). Cloud evolution in the $\epsilon_{\rm ff}$–$\sigma_{\rm v,GMC}$ space is similar to $\epsilon_{\rm ff}$–$R_{\rm gal,2D}$, i.e. evolution in $\epsilon_{\rm ff}$ occurs while $\sigma_{\rm v,GMC}$ remains largely unchanged.

A visual inspection of the evolution of these eight clouds (see additional material777Additional material can be found at https://www.youtube.com/playlist?list=PLO1VCxfFwCjDefBvKCF7-CptcZLn-gGMT) reveals that the environment of a GMC is as important as its internal processes. For example, the red cloud is situated in a particularly dense spiral arm which feeds the cloud with gas, allowing $M_{\rm GMC}$ to increase by almost an order of magnitude during the first $0.5\,{\rm Myr}$ of its life. In contrast the grey cloud forms in a much lower density environment, which is quickly disrupted by shear from galactic rotation. Furthermore only one of the eight clouds (dark blue) is clearly destroyed by supernovae, with the other seven instead appearing to be destroyed by shear, demonstrating that environment plays a role in the evolution of a cloud throughout its life.

3.2.2 General Trends in Evolution

To aid in teasing out general trends in cloud evolution, we create 2D histograms of $M_{\rm GMC},\,M_{\rm\star,y}$ and $\epsilon_{\rm ff}$ as functions of GMC age ($t_{\rm age}$) for all 1,879 evolutionary tracks (see §2.3) in Fig. 5. There is a wide range of different evolutionary paths taken by GMCs and that most clouds only live for $3$–$4\,{\rm Myr}$. Furthermore, this demonstrates that the eight randomly selected GMCs overlaid in both Fig. 3 and Fig. 5 are not the only GMCs with unique evolutionary paths.

The gas mass of a GMC when it is first detected tends to be the maximum gas mass ($M_{\rm GMC,max}$) that the cloud reaches. However, some GMCs reach their $M_{\rm GMC,max}$ at $t_{\rm age}\neq 0$ and have therefore gained mass through cloud-cloud collisions and accretion. Despite this, the general trend is for GMCs to lose mass as they evolve.

Interestingly, Kawamura et al. (2009) infer from observations that $M_{\rm GMC}(t_{\rm age}=0)\neq M_{\rm GMC,max}$ and that over a period of $\sim 30\,{\rm Myr}$ that $M_{\rm GMC}$ increases by a up to a factor of three, even after star formation begins (see their §4.2 and Table 4). This discrepancy could be due to any number of factors such as: how we account for cloud-cloud interactions, clouds dropping below our detection threshold (see §2.3), differences in how we define and detect clouds or inaccuracies in the model used to infer evolution of GMCs from observations. We leave further exploration of this difference to future work.

An obvious assumption would be that the lost gas mass is converted into stars. This is at least partially true, as as the young stellar mass ($M_{\rm\star,y}$) tends to increase during the first few million years of a clouds life. At $t_{\rm age}\gtrsim 4\,{\rm Myr}$, $M_{\rm\star,y}$ decreases with increasing age, yet clouds continue to lose significant fractions of their gas mass. Therefore, the lost gas mass is removed from the GMC by other means (e.g. feedback and shear).

In general we find that $\epsilon_{\rm ff}$ tends to increase during the first $4$–$6\,{\rm Myr}$ of a GMC’s life after this point we see that $\epsilon_{\rm ff}$ tends to either plateau or decrease. This corresponds to the age at which star particles will experience their first supernova event and therefore is a strong indicator of feedback limiting $\epsilon_{\rm ff}$.

To summarise, at any given $t_{\rm age}$ it is possible for clouds to have a large variety in $\epsilon_{\rm ff}$, $M_{\rm GMC}$ and $M_{\rm\star,y}$ and the exact value of given property at a given age is unique to each cloud.

4 Discussion

4.1 Cloud Conformity or Diversity?

As shown by the observational data present in L16 and included in Fig. 1 and discussed in §1, there is a significant dispersion in the values of $\epsilon_{\rm ff}$ for GMCs in the Milky Way. Feldmann & Gnedin (2011) put forward a model that is able to produce a wide spread in measured $\epsilon_{\rm ff}$ values ($1\lesssim\epsilon_{\rm ff}\lesssim 100\%$, see their Fig. 2) by adopting a fixed universal efficiency per free-fall time $\epsilon_{\rm ff,0}$ (analogous to $\epsilon_{\rm ff,SF}$ used in our simulations) and allowing the GMC to evolve with time. L16 combined Feldmann & Gnedin (2011)’s model with the star formation prescription given in Krumholz & McKee (2005) to allow for a time-dependant $\epsilon_{\rm ff,0}$ which results in a pair of coupled ordinary differential equations:

[TABLE]

and

[TABLE]

where $\gamma$ is the rate of gas accretion on to the GMC and $\alpha$ is a parametrisation of the disruption of GMCs due to feedback. Having solved the above equations the young stellar mass

[TABLE]

and the model equivalent of $\epsilon_{\rm ff}$

[TABLE]

can be calculated at each at every $t_{\rm age}$. In this model, the evolution of $M_{\rm GMC}/M_{\rm GMC,max},\,M_{\rm\star,y}/M_{\rm GMC,max}$ and $\epsilon_{\rm ff,\,L16}$ is the same for all clouds, for a given choice of $\delta$, $\epsilon_{\rm ff,0}$ and fixed free-fall time ($t_{\rm ff,\,fixed}$). In essence, this model requires conformity in the evolution of all GMCs and the dispersion in $\epsilon_{\rm ff}$ is produced by observing a population of GMCs, with each cloud at a different stage in its evolution.

We adopt the values of $\alpha,\,\gamma$ and $t_{\rm ff,\,fixed}$ given in L16 ($3.5,\,0$ and $6.7$ respectively) and reproduce their model for $\delta=0,1$ and $2$. L16 employed $\epsilon_{\rm ff,0}=0.014$ to ensure that all models produce clouds with lifetimes of $\sim 20\,{\rm Myr}$ and argued that $\epsilon_{\rm ff,0}\propto t_{\rm age}^{\delta}$ with $\delta=2$ was required to match observations. We adopt several different values for $\epsilon_{\rm ff,0}$. Firstly, we use $0.27,\,0.52$ and $0.47$ for the three values of $\delta$ respectively as these values ensure that all models convert $10\%$ of their gas mass to stars by the time the cloud is destroyed ($M_{\rm GMC}/M_{\rm GMC,max}\leq 0.01$). This ensures that all models have the same initial and final conditions (i.e. enforces conformity between models). Thus allowing us to determine if the diversity in our simulated GMCs can be explained by conformity to a single evolutionary path. The second set of values, $0.091,\,0.019$ and $0.006$, produce cloud lifetimes of $10,\,15$ and $20\,{\rm Myr}$ respectively for $\delta=2$ (i.e. assuming that GMCs will have different evolutionary paths). This set of models will allow for a determination on whether simply allowing for different GMC lifetimes is enough to explain the diversity in clouds efficiencies.

Testing this model against the evolutionary tracks of all GMCs in our simulation, i.e the red and blue lines in Fig. 5, shows that due to the large spread in our data, the model (independent of $\delta$ and $\epsilon_{\rm ff,0}$) overlaps with our simulated GMCs in each parameter space. However, that to reproduce the diversity seen in the simulated GMCs would require a significant number of models, each with a different values for $\delta$ and $\epsilon_{\rm ff,0}$: therefore a model which produces diversity is required.

In their recent work, Grudić et al. (2018) carried out isolated GMC simulations for three different mass clouds ($2\times 10^{4,\,5,\,6}{\,{\rm M}}_{\odot}$). They found that by observing a population of clouds, all with the same mass, at random points during their lifetime they could reasonably reproduce L16’s observed distribution in $\epsilon_{\rm ff}$. Furthermore, this was found independently of cloud mass (see their Fig. 4). Their simulations and conclusion support the models presented in L16: i.e. all clouds follow (nearly) identical evolutions and the spread in $\epsilon_{\rm ff}$ is a result of observing a population of different aged clouds.

Fig. 6 shows how $\epsilon_{\rm ff}$ evolves for clouds with gas mass at “birth”888Defined as the first time a cloud is detected by the clump finder ($M_{\rm birth}$). The clouds with the largest $M_{\rm birth}$ tend to reach higher $\epsilon_{\rm ff}$ and have a (marginally) higher $\langle\epsilon_{\rm ff}\rangle_{\rm med}$. The most noteworthy result from Fig. 6 is that birth mass plays only a small role in determining the initial value of $\epsilon_{\rm ff}$ as shown by $\epsilon_{\rm ff}(t_{\rm age}\leq 0.8\,{\rm Myr})$ having at least a 2 dex spread in each mass bin. Furthermore, while the dispersion in $\epsilon_{\rm ff}$ decrease as clouds age, it never reaches zero. This implies that knowing the birth mass of a GMC is not enough to predict the efficiency at which it converts gas into stars, and thus its star formation history. Carrying out a similar experiment using $R_{\rm birth}$ (not shown), yields nearly identical results: clouds are born with a range of different $\epsilon_{\rm ff}$ independent of their initial size.

Unlike the model presented in L16 and the isolated GMCs simulations of Grudić et al. (2018), the GMCs in this work are simulated in a (realistic999As shown by the analysis of the simulation in G17 and G18.) galactic environment which can heavily influence their evolution. In §3.2 we found that different GMCs were effected by their environment in different ways, i.e. some experience mergers and others are sheared apart, etc. Indeed recent observations of NGC 2276 have found that galactic-scale tidal forces and ram pressure has lead to large variations in molecular content of the galaxy, resulting in some regions with variations in the depletion time scale (the ratio of the molecular gas mass to the star formation rate) as large as several orders of magnitude when measured on scales of $\sim 450{\,\rm pc}$ (Tomic̆ić et al., 2018). It is therefore likely that measurements of $\epsilon_{\rm ff}$ on cloud scales in such a galaxy also find large variations. We therefore argue that the initial, intrinsic properties of GMCs ($M_{\rm birth}$ and $R_{\rm birth}$) are not sufficient to set the initial value and evolution of $\epsilon_{\rm ff}$: other factors such as the galactic environment (e.g shear) must be taken into account. It is the combination of a wide range of possible cloud properties and the environment in which GMCs live that naturally give rise to the observed and simulated spread in $\epsilon_{\rm ff}$.

Additionally, we note that Fig. 5 and 6 show that in our simulation $\epsilon_{\rm ff}$ does not follow a smooth, systematic evolution, instead it is able to both increase and decrease as clouds age. This is a direct contradiction of the prediction made by the L16 model, which predicts a continually increase in $\epsilon_{\rm ff,\,L16}$ with $t_{\rm age}$. The evolution seen in $\epsilon_{\rm ff}$ for our simulated clouds, as discussed above, result from combination of time dependant factors (e.g. mergers, gas accretion, feedback, shear and galactic tides) that are difficult to model as constant parameter in any given model and thus further evidence that the environment of clouds need to be accounted for when exploring their evolution.

Finally, we note that the analytical model we have adopted from L16 is not the only model. For example, the model present in Vázquez-Semadeni et al. (2018) predicts different star formation histories (and hence instantaneous $\epsilon_{\rm ff}$) for clouds of different mass. Their results and conclusions support this work and our conclusion that it is the evolution history of a cloud needs to be know to able to determine its $\epsilon_{\rm ff}$ at any given age. Therefore any model (numerical or analytical) must be able to capture the full range of physical processes that occur within a GMC and its interactions with its environment.

4.2 Limitations of the Simulations

The instantaneous H2 fraction calculated at run-time to determine the star formation rate of a computational-cell (see §2.1) is not stored or advected through the simulation. As a result we have to determine the molecular content within the simulation in post processing. For simplicity we chose to adapt a density threshold of $\rho_{\rm mol}=100{\,\rm{cm^{-3}}}$, with all gas above this value considered to be molecular. This limits the maximum value of $t_{\rm ff}$, which in turn acts as a limiting factor in determining $\epsilon_{\rm ff}$ from the simulations. If larger values of $t_{\rm ff}$ could be reached smaller values of $\epsilon_{\rm ff}$ maybe detected. A simple solution would be to rerun the simulations but including a treatment of the chemistry and thus allowing the molecular fraction of the gas to be self-consistently determined by the simulation, which could then be used to identify GMCs. However given that the current simulation is able to not only reproduce the median value of $\epsilon_{\rm ff}$ but also the size and distribution of the spread in values (and does so using a universal efficiency on the scale of computational cells, i.e. $\epsilon_{\rm ff,SF}$) we leave re-simulation for future work.

Isolated GMC simulations are able to completely resolve the internal structure of the GMC but at the cost of the galactic environment (e.g. see Padoan et al., 2016; Grudić et al., 2018). The simulations used throughout this work have such an environment but they are limited in spatial resolution (i.e. $\Delta x\sim 4.6{\,\rm pc}$, see §2.1). This resolution results in GMCs being made up of several computational cells and thus stars form and inject feedback into specific regions with the clouds. This allows one generation of stars to alter the gas structure within a GMC and even remove gas, thus determining where the next generation of stars form and how a GMC evolves. Despite the limited resolution, our simulation is able to accurately reproduce the galaxy wide gas probability distribution function (PDF), the range and distribution of cloud properties, including $\epsilon_{\rm ff}$ (see G17 and G18). We therefore argue that accurately fully resolving the internal structure of GMCs is not as important as accurately reproducing the large ($\geq 100{\,\rm pc}$) scale galactic environment.

5 Conclusions

In this work we explore the efficiency of Giant Molecular Clouds (GMCs) at forming stars in hydrodynamical simulations of Milky Way-like galaxies. The primary goal of this work is to explain the observed spread in the star formation efficiency per free fall time ($\epsilon_{\rm ff}$). To this end we calculate $\epsilon_{\rm ff}$ for each GMC found within two simulations, one with stellar feedback and one without. Using a tracking algorithm we follow the evolution of $\epsilon_{\rm ff}$ (and other properties) of each cloud throughout its lifetime. Our key results are:

Galactic disc simulations where star formation is determined by a Schmidt star formation law (applied to molecular gas) are able to produce the observed spread in the measured values of $\epsilon_{\rm ff}$ for GMCs. A large spread in values is found independently of the presence of stellar feedback, however the inclusion of feedback in the simulation prevents highly inefficient ($\epsilon_{\rm ff}<10^{-4}$) massive ($M_{\rm tot}>10^{7}{\,{\rm M}}_{\odot}$) clouds from forming. Stellar feedback is not the main source of the dispersion in $\epsilon_{\rm ff}$. 2. 2.

No single GMC property determines the $\epsilon_{\rm ff}$ of a cloud. Comparing seven key properties (gas mass, free-fall time, galactic radius, velocity dispersion, radius, density and virial parameter) of GMCs with $\epsilon_{\rm ff}$ shows no significant correlation. Instead we find that a cloud with a given value in any of the above properties is able to have a wide range of values in $\epsilon_{\rm ff}$. 3. 3.

Each GMC evolves in a unique way, determined by both its initial properties and its environment. It is therefore not possible to describe the evolution of all clouds by a single analytical model neglecting environmental effects. Furthermore, the evolution of a particular property for a given cloud is not smooth or uniform: a cloud is able to explore a wide range of values during its lifetime. This leads to a natural spread in properties and in particular the value of $\epsilon_{\rm ff}$. 4. 4.

The evolution of $\epsilon_{\rm ff}$ throughout a GMC’s lifetime does not follow a systematic increase, contrary to predictions of simple analytical models. Instead the measured value of $\epsilon_{\rm ff}$ for a cloud is driven by a number of time dependant factors, including stellar feedback and galactic environment, which can cause both increases or decreases as the clouds ages. This allows for a variety of different star formation histories.

In future work we will explore the processes that drive cloud evolution and how this leads to diversity in cloud properties.

acknowledgments

We thank the anonymous referee for their valuable and insightful comments. KG acknowledge support from the Science and Technology Facilities Council (grant ST/N002717/1), as part of the UK E-ELT Programme at the University of Oxford. KG also thanks and acknowledges support from New College, University of Oxford via the Balzan Fellowship. OA acknowledges support from the Swedish Research Council (grant 2014- 5791). OA and FR acknowledge support from the Knut and Alice Wallenberg Foundation. The research of JD and AS is supported by Adrian Beecroft and the STFC. This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure.