Determining Star Formation Thresholds from Observations

TL;DR

This study uses detailed simulations to evaluate observational techniques for detecting star formation thresholds, finding that some methods can indicate their presence but cannot reliably determine their exact values due to resolution limitations.

Contribution

The paper demonstrates through simulations that current observational diagnostics are limited in accurately identifying star formation thresholds in molecular clouds.

Findings

Projected star formation efficiency plots can detect thresholds.

Typical survey resolutions are insufficient to measure threshold values.

Common diagnostics may falsely suggest thresholds without physical basis.

Abstract

Most gas in giant molecular clouds is relatively low-density and forms star inefficiently, converting only a small fraction of its mass to stars per dynamical time. However, star formation models generally predict the existence of a threshold density above which the process is efficient and most mass collapses to stars on a dynamical timescale. A number of authors have proposed observational techniques to search for a threshold density above which star formation is efficient, but it is unclear which of these techniques, if any, are reliable. In this paper we use detailed simulations of turbulent, magnetised star-forming clouds, including stellar radiation and outflow feedback, to investigate whether it is possible to recover star formation thresholds using current observational techniques. Using mock observations of the simulations at realistic resolutions, we show that plots of projected star formation efficiency per free-fall time $\epsilon_{\rm ff}$ can detect the presence of a threshold, but that the resolutions typical of current dust emission or absorption surveys are insufficient to determine its value. In contrast, proposed alternative diagnostics based on a change in the slope of the gas surface density versus star formation rate surface density (Kennicutt-Schmidt relation) or on the correlation between young stellar object counts and gas mass as a function of density are ineffective at detecting thresholds even when they are present. The signatures in these diagnostics sometimes taken as indicative of a threshold in observations, which we generally reproduce in our mock observations, do not prove to correspond to real physical features in the 3D gas distribution.

Tables

Table 1. Table 1: Column 1: Simulation name. Columns 2-4: type of turbulence driving, 3D velocity dispersion, rms Mach number. Columns 5-6: ratio of thermal to magnetic pressure ( β 𝛽 \beta ) and Alfvén Mach number. Column 7: Dimensionless mass to flux ratio ( Equation 1 ). Columns 8-10: length of the simulation box, mean density, the sink particle threshold density. Columns 11-12: Time at which the simulation snapshot is taken and total number of stars present in the simulation snapshot used. Column 13: SFR per mean global freefall time ( ϵ ff subscript italic-ϵ ff \epsilon_{\rm ff} ).

Name	Turb. ($b$)	${\sigma }_v$	$ℳ$	$\beta$	${ℳ}_A$	${\mu }_{\Phi }$	$L$	${\rho }_0$	${\rho }_{\mathrm{sink}}$	Time	$N_{\mathrm{stars}}$	${ϵ}_{\mathrm{ff}}$
			[km/s]				[pc]	[g cm-3]	[g cm-3]	[Myr]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)
GTBJR	Mix ($b$=0.4)	1.0	5.0	0.33	2.0	3.8	2	3.28 $\times$ ${10}^{-21}$	8.3 $\times$ 10-17	4.254	12	0.031
C18Decay	None	1.254	6.6	0.046	1.0	1.56	0.65	4.46 $\times$ ${10}^{-20}$	1.16 $\times$ 10-15	1.91	55	0.12
C18Drive	Sol. ($b$=0.33)	1.254	6.6	0.046	1.0	1.56	0.65	4.46 $\times$ ${10}^{-20}$	1.16 $\times$ 10-15	1.842	63	0.048

Equations

$${\mu }_{\Phi }=\frac{M}{M_{\Phi }}=2\pi \sqrt{G}\left(\frac{M}{\Phi }\right),$$

$$t_{\mathrm{ff},2D}=\sqrt{\frac{3\pi }{32G{\rho }_{2D}}},$$

$${ϵ}_{\mathrm{ff},3D}=\frac{t_{\mathrm{ff}}({\rho }^{\prime })}{M_{\mathrm{gas},>}}\times \mathrm{SFR},$$

$${ϵ}_{\mathrm{ff},2D}=\frac{t_{\mathrm{ff},2D}}{M_{\mathrm{gas}}}\times \text{mbox}SFR,$$

$${\Sigma }_{\mathrm{th}}\sim {\rho }_{\mathrm{th}}{\lambda }_{\mathrm{J}}({\rho }_{\mathrm{th}})⟹{\rho }_{\mathrm{th}}\sim \frac{{\Sigma }_{\mathrm{th}}^{2}G}{\pi c_s^2},$$

$${ϵ}_{\mathrm{ff},\mathrm{obs}}=\left(\frac{0.5{\mathrm{M}}_{\odot }{N}_{\mathrm{YSO}}}{M_{\mathrm{gas}}}\right)\left(\frac{t_{\mathrm{ff},2D}}{t_{\mathrm{YSO}}}\right)$$

$${\Sigma }_{\mathrm{SFR}}=\frac{0.5{\mathrm{M}}_{\odot }{\mathrm{N}}_{\mathrm{YSO}}}{t_{\mathrm{YSO}}A},$$

Full text

††pubyear: 2019††pagerange: Determining Star Formation Thresholds from Observations–References

1 Introduction

Understanding the physical factors behind the formation of stars from interstellar gas is key to developing a predictive theory of star formation and understanding the evolution of galaxies. Star formation is known to occur in filamentary structures (Goldsmith et al., 2008; André et al., 2014) in molecular clouds (Wong & Blitz, 2002; Kennicutt et al., 2007; Blanc et al., 2009; Krumholz, 2014) and it can be characterized by a quantity known as the star formation efficiency per free-fall time $\epsilon_{\rm ff}$ (Krumholz & McKee, 2005), which measures the fraction of the gas that is converted to stars per free-fall time. On molecular cloud scales, the average value of $\epsilon_{\rm ff}$ is known to be small, $\approx$ 0.01111As discussed in the Krumholz et al. (2019) review, the amount of spread about this average is a subject of current debate; Krumholz et al. argue that the weight of evidence favours a relatively small spread of $\approx 0.3$ dex, but some authors argue for larger spreads of $\gtrsim 1$ dex. (e.g., Krumholz & Tan 2007; Krumholz et al. 2012; Federrath 2013a; Evans et al. 2014; Salim et al. 2015; Vutisalchavakul et al. 2016; Heyer et al. 2016; Leroy et al. 2017; Sharda et al. 2018; see Krumholz et al. 2019 for a recent review) This means that on giant molecular cloud (GMC) or molecular cloud length scales ($\sim 100$ pc down to $1$ pc), star formation is very inefficient. However, as we move towards smaller length scales ($\sim$ 0.1 pc down to AU scales) tracing gas at densities exceeding $\sim 10^{7}$ cm*-3*, eventually there must be some density or size scale beyond which most of the mass in gas would wind up in a star $\sim 1$ dynamical time later. Therefore, there should exist a point after which $\epsilon_{\rm ff}$ does not remain small anymore and approaches unity.

Nearly every physical model of star formation predicts the existence of a threshold of this type. Knowing its value would tell us a great deal about how star formation works and what regulates it. For example, one could imagine that the dense clumps traced by HCN emission are gravitationally bound structures, whereas gas in GMCs is largely unbound, and this is why star formation is inefficient on GMC scales (e.g., Heiderman et al., 2010; Lada et al., 2012). If that explanation were correct, one would expect to find low $\epsilon_{\rm ff}$ on the scales of GMCs, but high $\epsilon_{\rm ff}$ in regions traced by HCN. Such an observation would be powerful evidence that a change in boundedness is what is regulating star formation. Alternately, a number of authors have proposed models in which star formation is regulated by supersonic turbulence (Krumholz & McKee, 2005; Federrath & Klessen, 2012, 2013; Padoan & Nordlund, 2011; Hennebelle & Chabrier, 2011; Hopkins, 2012, 2013; Federrath, 2015). A generic feature of such models is the existence of a characteristic density scale (or a range of them in some cases) at which gas becomes bound; this too represents a predicted threshold at which we would expect a change in $\epsilon_{\rm ff}$. Similar arguments can be used to derive critical densities, column densities, or length scales from magnetic-regulation models of star formation (Shu et al., 1987; McKee, 1989; Basu & Ciolek, 2004).

In order to detect a threshold in observations, we must know what it looks like, and for this purpose it is helpful to envision the process of star formation in a Lagrangian sense: as a conveyor belt moving mass from low to high to stellar densities. The value of $\epsilon_{\rm ff}$ characterises the speed of this flow for any particular fluid parcel: $\epsilon_{\rm ff}\ll 1$ at a given density $\rho$ means either that fluid parcels of density $\rho$ require many free-fall times $t_{\rm ff}\approx 1/\sqrt{G\rho}$ to substantially increase in density, or that fluid parcels that are increasing in density are nearly balanced out by those decreasing in density, so the net flow of mass to higher density is small. On the other hand, values of $\epsilon_{\rm ff}\approx 1$ mean that the typical fluid parcel requires only a time $\approx t_{\rm ff}$ to go up in density. Since the mass flow from low to high density must be continuous, at least in a time-averaged sense, the value of $\epsilon_{\rm ff}$ determines how much mass is “stuck” at a given density along the conveyor belt: low $\epsilon_{\rm ff}$ corresponds to places where the conveyor belt moves slowly and thus large amounts of mass build up, while high $\epsilon_{\rm ff}$ corresponds to places where the conveyor belt moves rapidly and there is relatively little mass. Thus the observational signature of a threshold is as illustrated in the cartoon in Figure 1: if one selects a density such that $\epsilon_{\rm ff}\ll 1$, then there is a large amount of mass stuck at densities from $n$ to $2n$, and contours drawn at these two densities will enclose a great deal of mass. If one selects a higher density $n$ such that $\epsilon_{\rm ff}\approx 1$, there is little mass between $n$ and $2n$. The threshold $n_{\rm thresh}$ is the density where the contours switch from being widely-spaced to narrowly-spaced in mass.

While this manifestation is relatively straightforward to investigate in a simulation where we have access to the full 3D density field, for observed clouds where we have access to quantities only in projection, it is obviously not possible to search for a threshold as illustrated in Figure 1. Instead, one must use some observable proxy. Onishi et al. (1998), Johnstone et al. (2004), and Hatchell et al. (2005) were among the first to search for a threshold using observations. Onishi et al. (1998) studied the relations between cores and young stellar objects in Taurus, and found that C18O cores that contained either compact dense centres traced by H13CO or cold Infrared Astronomical Satellite (IRAS) sources have high column densities, and therefore concluded that gravitational collapse and subsequent mass accretion occurs once a core’s column density exceeds $\sim 8\times 10^{21}$ cm*-2* ($A_{V}\sim 7-9$ mag). Johnstone et al. (2004) studied substructures in the Ophiuchus clouds using submillimeter continuum maps and compared them with visual extinction maps. These authors note the absence of any structures for $A_{V}\leq 7$ mag, and argue that this indicates that substructures only form where $A_{V}\geq 15$ mag.

Heiderman et al. (2010) and Lada et al. (2010) lend support to the idea that there is a threshold based on data for a larger sample of local molecular clouds. They argue that the correlation between the star formation rate and gas mass becomes increasingly tight as one considers the gas mass at higher and higher column densities. Könyves et al. (2015), André (2015), and André (2017) find that the objects they identify as cores in Herschel maps (which are in practice defined by contours at a certain signal to noise threshold) are found almost exclusively on background material characterized by an extinction of $A_{V}\sim 7-9$ mag. They argue that this extinction level is characteristic of the densities at which gravitational collapse of filamentary structures takes place. Lada et al. (2012) propose that the existence of a density threshold explains why external galaxies show a near-linear correlation between star formation rate and mass of gas traced by HCN line emission, but a superlinear correlation between star formation and CO line emission. In their model, HCN traces gas above the threshold, while CO traces gas below it.

However, a number of authors have questioned these conclusions. García-Burillo et al. (2012), Usero et al. (2015), and Bigiel et al. (2016) show that the correlation between HCN emission and star formation in extragalactic systems is not in fact linear, calling into question the idea that HCN traces gas that is especially closely linked to star formation (though see Shimajiri et al. 2017 for an opposing view). Elmegreen (2018), following up on earlier arguments by Krumholz & Thompson (2007), argues that the near-linearity of the HCN-star formation correlation is an observational selection effect rather than a physical threshold. For Galactic measurements, Burkert & Hartmann (2013) argue that the correlations observed by Heiderman et al. (2010) and Lada et al. (2010) do not require a particular density threshold, and variation in the dependence of $\Sigma_{\rm SFR}$ on $\Sigma_{\rm gas}$ is instead due to the increasing importance of gravity at higher densities. Consistent with this picture, Gutermuth et al. (2011) show that the density of young stellar objects increases as a smooth powerlaw with gas surface density, again suggesting a smooth rise in star formation rate in dense gas without any special threshold density value. Clark & Glover (2014) find that the clouds in their simulations can still form stars at cloud-averaged densities which are lower than the “threshold" value of $\sim 120\ \mathrm{M_{\odot}\ pc^{-2}}$ proposed by Heiderman et al. (2010) and Lada et al. (2010) and suggest that their threshold for star formation is more likely a consequence of the star formation process, rather than a prerequiste for star formation. In other words, regions where star formation is more active and there are more YSOs present also tend to be regions where a great deal of gas has collapsed to high surface density, but the latter is not a direct cause of the former.

In this paper, we aim to first confirm that there is a threshold density above which $\epsilon_{\rm ff}$ approaches unity in detailed simulations of star formation, and then to investigate whether we can recover this value from observations using a variety of proposed diagnostic methods. We answer these questions by analyzing simulations that include some of the most detailed descriptions of the different physical processes relevant for star formation, and mimic conditions prevalent in nearby molecular clouds. We use these simulations to create mock observations, and test whether various proposed methods for detecting thresholds in observations can in fact recover results that match what we obtain using the full 3D spatial and temporal information to which we have access for the simulations.

The rest of the paper is structured in the following manner. We begin by describing the simulations and the methodology used to create mock observations in Section 2. In Section 3, we first search for star formation thresholds using the full simulation information (in particular the gas volume density), and then assess the capability of various observational analysis methods to recover such a threshold using the projected information to which observers have access. We summarize our findings and the main conclusion of our study in Section 4.

2 Simulations and Analysis Methods

In our study we use the simulations described in Federrath (2015) and Cunningham et al. (2018) (hereafter F15 and C18, respectively). We choose these simulations since they include detailed treatments of gravity, turbulence, magnetic fields, mechanical jet/outflow and radiation feedback and obtain star formation rates (F15) and IMFs (C18) that are among the closest matches to observations to date. We briefly highlight their main features strictly relevant to this work and refer the reader to F15 and C18 for more details. We also summarize some key simulation parameters in Table 1.

2.1 F15 Simulations

F15 use the AMR code (Berger & Colella, 1989) FLASH (Fryxell et al., 2000; Dubey et al., 2008) to solve the compressible MHD equations. While the simulations described in F15 include physical processes in steps of increasing complexity, we only use the one with the most complete set of physical processes, which includes self-gravity, turbulence, magnetic fields, and jet/outflow feedback. The simulation that we use here is not directly from F15, but is constructed with the same initial and boundary conditions as all the simulations in F15, is identical to the most complex simulation in F15 labelled GvsTMJ and additionally includes radiation feedback based on the implementation by Federrath et al. (2017) who use the model described in Offner et al. (2009) for protostellar evolution. This simulation is labelled GTBJR in Onus et al. (2018) and we adopt the same name for the purposes of this work. The turbulence driving in the simulation excites a natural mixture of solenoidal and compressible modes, corresponding to a turbulence driving parameter $b=0.4$ (Federrath et al., 2010a). We refer the reader to F15 for more details on the implementation of sink particle formation, turbulence, magnetic fields, jets/outflows (Federrath et al., 2014) and to Federrath et al. (2017) for the implementation of radiation feedback. Fragmentation, star formation and accretion are modelled with the sink particle technique by Federrath et al. (2010b).

The simulation box has a total cloud mass $M=388\ \mathrm{M}_{\odot}$ and has a box length of $L=2$ pc, with a mean density of $\rho_{0}=3.28\times 10^{-21}$ g cm*-3*, corresponding to a global free-fall time of $t_{\rm ff}=1.16$ Myr. The velocity dispersion is $\sigma_{v}=1$ km s*-1*, and the sound speed is $c_{s}=0.2$ km s*-1* at the initial temperature $T=10$ K, which results in an rms Mach number of $\mathcal{M}=5$. The simulation starts with an initial uniform magnetic field of $B=10$ ${\rm\mu}$G. These simulations can be characterized by the magnetic field strength parameter $\mu_{\Phi}$, constant throughout for the whole cloud (albeit varying locally within the cloud) since mass and magnetic flux are conserved for the entire simulation box, defined as

[TABLE]

where $M_{\Phi}=\Phi/2\pi\sqrt{G}$ is the magnetic critical mass (Mouschovias & Spitzer 1976; the mass below which a cloud cannot collapse and above which collapse cannot be prevented by magnetic fields alone), and $\Phi$ is the total magnetic flux threading the cloud. The value of $\mu_{\Phi}$ for this simulation is 3.8. The resulting virial ratio is $\alpha_{\rm vir}=1.0$ and the plasma beta is $\beta=0.33$ (corresponding to an Alfvén Mach number of $\mathcal{M}_{\rm A}=2.0$). These physical conditions are chosen to mimic those found in nearby, low-mass, star-forming regions such as Perseus or Taurus.

Sink particles in the simulation are formed dynamically when a local region undergoes gravitational collapse. Once the gas density in a cell exceeds a density of $\rho_{\rm sink}=\pi c_{s}^{2}/G\lambda_{\rm J}^{2}$, a control volume of radius $r_{\rm sink}=\lambda_{\rm J}/2$ is formed around it and it is checked whether all the gas in that volume is Jeans-unstable, is gravitationally bound and is collapsing towards the central cell. If all these additional checks are passed, a sink particle is formed in the central cell. Performing these additional checks suppresses spurious sink formation in transient shocks (Federrath et al., 2010b).

2.2 C18 Simulations

C18 use the ORION2 AMR code (Li et al., 2012) to solve the equations of ideal MHD along with treatments of coupled self-gravity (Truelove et al., 1998; Klein et al., 1999), and radiation transfer (Krumholz et al., 2007). The simulations also include feedback due to protostellar outflows following the procedure in Cunningham et al. (2011); stars form following the sink particle algorithm of Krumholz et al. (2004), and protostellar evolution uses the model described in Offner et al. (2009). The gas is initially evolved under the action of a turbulent driving force for two crossing times, $t_{\rm cross}=L/v_{\rm rms}=0.51$ Myr, after which self gravity is switched on. The evolution of the system thereafter is categorized into two cases, one where turbulence is allowed to decay and the other where a constant rate of energy is injected to balance the rate of turbulent decay (Mac Low, 1999). Although C18 carry out simulations with a wide range of $\mu_{\Phi}$ values, we focus on the simulations with $\mu_{\Phi}$ = 1.56 because this value is comparable to observed mass to flux ratios in nearby molecular clouds, and because these simulations run long enough to produce enough protostars to yield meaningful statistics. These simulations are also inefficient in forming stars for the driven turbulence case. We label the two cases C18Decay and C18Drive (based on whether turbulence is allowed to decay or being driven) for the purposes of this work; these correspond to the runs in Rows 1 and 3 in Table 1 of C18. For more details of the physics included and the limitations of the simulations, we refer the reader to C18.

The initial temperature and sound speed are the same as for the F15 simulations. The turbulent driving force is purely sinusoidal (b$\sim$0.33) and is scaled to maintain an rms Mach number of $\mathcal{M}=6.6$, leading to a virial parameter $\alpha_{\rm vir}=1.05$. The simulation box has a total cloud mass of $M$ = 185 $\mathrm{M}_{\odot}$ and has a box length of $L=0.65$ pc, with a mean density of $\rho_{0}=4.46\times 10^{-20}$ g cm*-3* corresponding to a global free-fall time of $t_{\rm ff}=0.315$ Myr. The value $\mu_{\Phi}=1.56$ corresponds to a magnetic field of $B=81.3$ ${\rm\mu}$G and an Alfvén Mach number of $\mathcal{M}_{A}$ = 1.4.

Sink particles in the simulation are formed only on the finest AMR level when the gas becomes dense enough to exceed a local Jeans number of $J=(G\rho\Delta x^{2}/\pi c_{s}^{2})^{1/2}>1/4$ or equivalently $\rho>\rho_{\rm sink}=\pi c_{s}^{2}/16G\Delta x^{2}$, where $\Delta x$ is the cell width at the finest AMR level.

2.3 Creating Mock Observations

To create mock observations of column density we choose snapshots from GTBJR, C18Decay and C18Drive at times corresponding to 4.254, 1.91 and 1.842 Myr respectively and make projection maps along each of the three cardinal axes. To compare these maps to observations of nearby molecular clouds, we smooth them with a Gaussian kernel with a FWHM corresponding to typical resolutions in observations using dust-based tracers; we focus on dust rather than molecular lines because dust measurements generally provide the most accurate estimates of total column density on small scales (Goodman et al., 2009). There are a wide range of distances and resolutions found in observational studies available in the literature. The physical resolution depends on target distance, observation wavelength, instrument, and technique. For this reason we consider two representative cases: (i) a resolution of 3.0 arcmin (typical of NIR extinction measurements; Juvela & Montillaud 2016) for a cloud at a distance of 140 pc (typical of Taurus; Torres et al. 2007), corresponding to an absolute resolution of 0.070 pc and (ii) a resolution of 36.9 arcsec (typical for Herschel observations at 500 $\mu$m; André et al. 2010) at a distance of 260 pc (typical of Aquila; Straižys et al. 2003; Könyves et al. 2015), corresponding to an absolute resolution of $\approx$ 0.046 pc. We shall denote these two resolutions as Res1 and Res2 respectively.

We draw contours on the resulting smoothed image at specific surface density values (which we also refer to as “contour levels”), $\Sigma_{\rm gas}$, evenly spaced throughout the entire range available for a given projection map along each axis. Since the contour shapes can vary and some contours can stretch over the entire length of the simulation box, we use only the ones that form closed curves. In cases where there are multiple closed contours for a given contour level, we combine the areas of the different unique contours and treat them as a single entity. We then project the positions of all the stars (sink particles) formed in the simulations onto these maps and count the number of Young Stellar Objects (YSOs), i.e., stars having an age less than 0.5 Myr (Class 0/I), enclosed within each contour (denoted by $N_{\rm YSO}$). We compute the following quantities for each contour level on both the smoothed and un-smoothed maps: (i) the enclosed area ($A$), (ii) the enclosed gas mass ($M_{\rm gas}$), (iii) the total mass of the YSOs ($M_{*}$), and (iv) the free-fall time

[TABLE]

where $\rho_{\rm 2D}=3\sqrt{\pi}M_{\rm gas}/4A^{1.5}$ is the mean density that would be estimated by an observer under the assumption that the line-of-sight size of the region is comparable to the size projected on the plane of the sky (e.g., Krumholz et al., 2012). We show an example of the method in Figure 2.

3 Results

Having generated our simulated observations, we now investigate how well one can diagnose star formation thresholds from them using various techniques. We begin in Section 3.1 by investigating the physical threshold in the simulation using the full 3D and time-dependent information to which we have access. In the remainder of this section we investigate various methods for trying to recover this information from 2D projected data and mock observations.

3.1 The true threshold: variation of $\mathbf{\epsilon_{\rm ff}}$ with $\mathbf{\rho_{\rm gas}}$

As discussed in Section 1, the dimensionless quantity $\epsilon_{\rm ff}$ is the fraction of an object’s gas mass that is transformed into stars in one free-fall time at the object’s mean density. Since denser regions invariably form stars more quickly than more diffuse ones (indeed, it would be very surprising if they did not), any claim of any type of volume or column density “threshold” for star formation below which star formation is suppressed, must reveal itself as a significant change in $\epsilon_{\rm ff}$ across the threshold.

At a minimum the simulations must have a threshold because they include sink particles – by construction gas that meets the conditions for sink particle creation $\epsilon_{\rm ff}\sim 1$. A necessary but not sufficient condition for sink particles to form is that the density exceed a threshold value $\rho_{\rm sink}$ (or number density $n_{\rm sink}$), which is $2.13\times 10^{7}$ cm*-3* for GTBJR and $2.9\times 10^{8}$ cm*-3* for C18Decay and C18Drive. However, a transition from low efficiency to free-fall collapse could also occur at densities considerably below these values. To determine if this is the case, in Figure 3 we plot $\epsilon_{\rm ff}$ as a function of $\rho_{\rm gas}$ using

[TABLE]

where SFR is the true, time-averaged star formation rate in the simulations, $M_{\mathrm{gas,>}}$ is the mass of gas above a certain density value chosen as the threshold ($\rho_{\rm gas,th}$), and $t_{\mathrm{{ff}(\rho^{\prime})}}$ is the free-fall time evaluated at $\rho^{\prime}$, which is the mean density for all the gas with $\rho>\rho_{\rm gas,th}$. This method is equivalent to 3D contouring (the 3D version of the method used in Section 2) with the assumption that all the stars lie in the densest regions.

Figure 3 shows a few interesting features. $\epsilon_{\rm ff,3D}$ is roughly constant to within a factor of $\sim 3$ from $n\approx 10^{4}-10^{6.5}$ cm*-3* for all the three simulations. (The GTBJR simulation also shows a slight bump in $\epsilon_{\rm ff}$ at $n\sim 10^{5}$ cm*-3*, but this appears to be a transient feature; similar features appear and disappear at a range of intermediate $n$ values in other F15 simulations (not shown), but do not appear at a consistent density or time, and no similar features appear in the C18 simulations.) Note that, while $\epsilon_{\rm ff}$ is constant from $n\approx 10^{4}-10^{6.5}$ cm*-3*, it begins to rise for even lower densities, particularly in the C18 simulations. We disregard this effect because it is an artifact of the periodic box used in the simulations. For density thresholds approaching the mean simulation box density ($\approx 10^{4}$ cm*-3* for C18, $\approx 10^{3}$ cm*-3* for GTBJR), further decreases in density do not add any more mass or star formation inside the contour. It is this effect that drives the rise in $\epsilon_{\rm ff}$ at the lowest densities.

Despite these artifacts at low density, Figure 3 provides clear evidence of a physical threshold: $\epsilon_{\rm ff,3D}$ is roughly constant and $\lesssim 0.1$ over the range from $n\approx 10^{3}-10^{6.5}$ cm*-3*, then rises sharply at $n\approx 10^{6.5}$ cm*-3*. This is clearly well below the threshold density at which sink particles form,222Note that, for GTBJR, due to the large number of additional checks (see e.g. Federrath et al. (2010b)) imposed before sink particles form, the mean density at which sink particles form in practice is considerably higher than $\rho_{\rm sink}$. indicating that this rise is not due to the sink particle algorithm. Instead, Figure 3 provides clear evidence that star formation in the simulations transitions from inefficient, $\epsilon_{\rm ff}<0.1$, to efficient, $\epsilon_{\rm ff}\sim 1$, at $n\approx 10^{6.5}$ cm*-3*. Thus, we confirm that there exists a true volume density threshold for star formation in the simulations. We therefore next investigate how well various observational methods can recover it.

3.2 Thresholds in projection: $\mathrm{\epsilon_{\rm ff}}$ vs. $\mathbf{\Sigma_{\rm gas}}$

In this section, we consider a first possible method to recover the volume density threshold using essentially the same method proposed by Krumholz et al. (2012), Federrath (2013b), Salim et al. (2015), and Sharda et al. (2018), which simply represents a projected version of the volumetric analysis used above.

We first estimate $\epsilon_{\rm ff}$ (denoted $\epsilon_{\rm ff,2D}$ here, to distinguish the projected, 2D from the true, 3D value) for each contour level in the un-smoothed column density maps for each simulation (see e.g. top panel in Figure 2) as

[TABLE]

where $t_{\rm ff,2D}$ is the free-fall time estimated from projected quantities, per Equation 2. Figure 4 shows the variation of $\epsilon_{\rm ff}$ as a function of contour level $\Sigma_{\rm gas}$ for this un-smoothed image. The curves are qualitatively similar to those shown in Figure 3 and show the “rising $\epsilon_{\rm ff}$” feature at the highest densities. Therefore, evidence for a density threshold can be recovered from such a plot. To estimate the threshold volume density $\rho_{\rm th}$ associated with the column density $\Sigma_{\rm th}$ where $\epsilon_{\rm ff,2D}$ begins to rise, we can hypothesise that the characteristic line of sight depth should be of order the Jeans length, i.e.,

[TABLE]

where $\lambda_{\rm J}(\rho_{\mathrm{th}})$ is the Jeans length evaluated at $\rho_{\mathrm{th}}$. The top axis in Figure 4 shows the values of $n_{\rm th}$ evaluated using Equation 5. We find that $\rho_{\mathrm{th}}$ estimated from Equation 5 and Figure 4 is reasonably close to the actual values found in Figure 3, i.e., $\epsilon_{\rm ff,2D}$ begins to rise at a projected column density $\approx$ $2\times 10^{3}\ \mathrm{M}_{\odot}$ pc*-2*, and this corresponds to an estimated volume density $n_{\rm th}\sim 10^{6.5}$ cm*-3*, not far off the 3D value one would infer directly from Figure 3.

Next we investigate whether this “rising $\epsilon_{\rm ff}$” feature can be recovered from the mock observations (see Section 2) we use in this work. As is common in observational work, we use the areas $A$, gas mass $M_{\rm gas}$, YSO count $N_{\rm YSO}$, free-fall time $t_{\rm ff}$ of each contour to determine $\epsilon_{\rm ff}$ as

[TABLE]

where $t_{\rm ff,2D}$ is again computed from Equation 2, but this time for the smoothed maps, and $t_{\rm YSO}=0.5$ Myr, our nominal estimate of the class 0/I lifetime. Figure 5 shows the $\epsilon_{\rm ff}$ versus $\Sigma_{\rm gas}$ for the mock observations at the two different resolutions, Res1 and Res2. We do not probe extremely low surface density values in the mock observations since such contour levels do not possess any closed curves (contours) in the projection map. Therefore, there may be variation in the estimated value of $\epsilon_{\rm ff}$ for surface density values below the first closed contour. We see that the rising $\epsilon_{\rm ff}$ feature can still be distinguished, and thus even in the mock observations we can detect clear evidence for a threshold using this method; however, the value of the inferred threshold column density, and thus the volume density that one would infer using Equation 5, depends strongly on resolution for the two resolutions we have explored. For the lower resolution case Res1, $\epsilon_{\rm ff,obs}$ appears to begin rising at columns of several hundred M⊙ pc*-2*, whereas in the un-smoothed map the rise is at several times $10^{3}\ \mathrm{M}_{\odot}$ pc*-2*. The corresponding volume density threshold one would infer from Equation 5 is thus close to two orders of magnitude too small. The error is much smaller for the Res2 case. We conclude that plots of $\epsilon_{\rm ff,obs}$ versus $\Sigma_{\rm gas}$ are a reliable method for detecting the existence of thresholds, but that inferred threshold values can be strongly affected by beam-smearing, which washes out the small, high-density structures one would need to measure in order to infer the threshold column density reliably.

3.3 Thresholds from the Kennicutt-Schmidt relation: $\Sigma_{\rm SFR}$** vs. $\Sigma_{\rm gas}$**

We next examine the method for detecting column density thresholds proposed by Heiderman et al. (2010, hereafter H10). H10 propose searching for a threshold for star formation by plotting a Kennicutt-Schmidt (KS) relation, i.e., by plotting $\Sigma_{\rm SFR}$ versus $\Sigma_{\rm gas}$, for individual star-forming clumps. They argue that this relationship turns over from superlinear to near-linear at a surface density $\Sigma_{\rm th}\approx 129\ \mathrm{\mathrm{M}_{\odot}}$ pc*-2*, and that this provides evidence for a threshold star formation surface density. To test this method on our simulations, we compute the observationally-estimated star formation rate per unit area for the different column density contour levels as

[TABLE]

where our method of estimating the star formation rate per unit area is identical to that used by H10. We similarly compute the mean surface density of the gas inside each contour, $\Sigma_{\rm gas}=M_{\rm gas}/A$.

In Figure 6, we plot the relationship between $\Sigma_{\rm gas}$ and $\Sigma_{\rm SFR}$ from our simulations for the lower resolution Res1 alongside the data from H10. Changing the resolution to Res2 does not result in a qualitative change in the features of the plot. While all three simulations, GTBJR, C18Decay and C18Drive extend significantly above the observations at higher column densities, GTBJR substantially overlaps the observations at lower column densities. The discrepancy at higher columns is not all that surprising, since H10’s data represent measurements only for whole “clouds” (defined as the outermost detectable contours given their sensitivity), not sub-regions within them as do our high column density lines.

The main point to take from Figure 6 is that it provides no evidence for a volume or column density threshold for star formation in our simulations, despite the fact that there is one, as is apparent from examining Figures 3 and 4. To the extent that the observed data shown in Figure 6 provide evidence for any change in star formation behaviour above a column density threshold, as opposed to being purely an observational artefact (Krumholz et al., 2012), this threshold does not appear to be related to a change in star formation efficiency.

3.4 Thresholds from the star to gas ratio: $N_{\rm YSO}/M_{\rm gas}$ vs. $\Sigma_{\rm gas}$

We next consider the method for detecting thresholds, and the associated evidence for a threshold, presented by Lada et al. (2010, hereafter L10). L10 survey 11 nearby clouds and produce IR extinction maps. They count the number of YSOs within contours drawn on clouds and show that the relation $N_{\rm YSO}\approx M_{\rm gas}^{0.96\pm 0.11}$ is satisfied for these clouds when $M_{\rm gas}$ is taken to be the total mass of gas above an extinction threshold of 0.8 $A_{\mathrm{K\ mag}}$. L10 find that the cloud-to-cloud scatter in the ratio $M_{\rm gas}/N_{\rm YSO}$ (or equivalently $N_{\rm YSO}/M_{\rm gas}$) between the 11 local clouds surveyed is minimised if one counts mass and YSOs within a contour corresponding to an extinction $A_{\mathrm{K}}=0.8\ \mathrm{mag}$. They argue based on this that gas at extinction above $A_{\mathrm{K}}=0.8\ \mathrm{mag}$ constitutes the active material for star formation, so that $A_{\mathrm{K}}=0.8\ \mathrm{mag}$ constitutes a threshold.

In Figure 7, we show the $N_{\rm YSO}/M_{\rm gas}$ as a function of contour level $\Sigma_{\rm gas}$ for the lower resolution (Res1) mock observations from simulations alongside the data from L10; we use the Res1 rather than Res2 case because this is more typical of NIR extinction methods. We find that the mock data from the simulations lie in close agreement with the observed data, and that the scatter between the results for the three projection axes for a single simulation is comparable to the observed cloud-to-cloud scatter. However, examining Figure 7 it is clear that the simulation points do not display any feature through which one could recover the physical volume density threshold that we have shown is present. The simulations and observations both do show a rise in $N_{\rm YSO}/M_{\rm gas}$ at the highest contour levels. This could plausibly be attributed to a rise in star formation efficiency at the highest column densities, but we note that these features are (1) at extinction values far higher than the purported threshold of $0.8$ mag in K-band luminosity for which L10 argue, and (2) as a result of low resolution are at column densities that are significantly below the true threshold discussed in Section 3.1 and Section 3.2.

4 Summary and Conclusions

In this paper we use simulations of star formation including a wide range of physical processes to probe the existence of star formation thresholds, and we investigate whether it is possible to measure them using current observational techniques. We show that the simulations we analyse do possess a true local volume density threshold (different from the sink particle threshold) of $n_{\rm th}\sim 10^{6}-10^{7}$ cm*-3*; gas above this threshold appears to form stars efficiently, without substantial opposition from turbulence, magnetic fields or outflow feedback. The evidence for this is a drastic change in the star formation efficiency per free-fall time $\epsilon_{\rm ff}$ across the threshold, a feature that naturally defines a true star formation threshold in the simulations.

We then investigate whether it is possible to recover this information using only the 2D, static information that is available in observations. We create projection maps from the simulations, which we then blur to resolutions typical of mid-infrared dust emission or near-infrared dust extinction studies. We find that when we estimate $\epsilon_{\rm ff}$ from un-smoothed projection maps, we can still recover the characteristic rising $\epsilon_{\rm ff}$ feature at column density values of $\Sigma_{\rm gas,th}\sim 2\times 10^{3}\ \mathrm{M}_{\odot}\ \mathrm{pc}^{-2}$, and we propose a method to estimate an associated volume density that recovers the true, 3D threshold volume density reasonably well. For the mock observations at realistic resolution, we find that plots of projected $\epsilon_{\rm ff}$ remain a reliable indicator of the presence of a threshold, as shown by a rising $\epsilon_{\rm ff}$ feature, but that the resolution currently typical of such observations renders recovery of the true column or volume density a challenge. In contrast, methods based on comparing $\Sigma_{\rm SFR}$ and $\Sigma_{\rm gas}$, or on measuring the ratio of young stellar object number to gas mass, provide no evidence for a volume or column density threshold for star formation in our simulations, despite the fact that one exists. We therefore conclude that these methods are not reliable indicators for the presence of a threshold for star formation.

Acknowledgements

We would like to thank the anonymous referee whose useful comments have helped improve the manuscript. MRK and CF both acknowledge support from the Australian Research Council’s (ARC) Discovery Projects and Future Fellowship funding schemes, awards DP160100695 (MRK), DP170100603 (CF), FT180100375 (MRK), and FT180100495 (CF), and from the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). AJC’s work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. The simulations used in this work were made possible by grants of high-performance computing resources from the following: the National Center of Supercomputing Application through grant TGMCA00N020, under the Extreme Science and Engineering Discovery Environment, which is supported by National Science Foundation grant number OCI-1053575; the NASA High-End Computing Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center; the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391); the Partnership for Advanced Computing in Europe (PRACE grant pr89mu); the National Computational Infrastructure, which is supported by the Australian Government (grants ek9 and jh2); the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. The simulation software FLASH was in part developed by the DOE-supported Flash Center for Computational Science at the University of Chicago. LLNL-JRNL-766677