{'fraction_non_alphanumeric': 0.03626763491227584, 'fraction_numerical': 0.02148063970730896, 'mean_word_length': 4.250600305610129, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The X-factor in Galaxies: II. The molecular hydrogen -- star formation relation

There is ample observational evidence that the star formation rate (SFR) surface density, Sigma_SFR, is closely correlated with the surface density of molecular hydrogen, Sigma_H2. This empirical relation holds both for galaxy-wide averages and for individual>=kpc sized patches of the interstellar medium (ISM), but appears to degrade substantially at a sub-kpc scale. Identifying the physical mechanisms that determine the scale-dependent properties of the observed Sigma_H2-Sigma_SFR relation remains a challenge from a theoretical perspective. To address this question, we analyze the slope and scatter of the Sigma_H2-Sigma_SFR relation using a set of cosmological, galaxy formation simulations with a peak resolution of ~100 pc. These simulations include a chemical network for molecular hydrogen, a model for the CO emission, and a simple, stochastic prescription for star formation that operates on ~100 pc scales. Specifically, star formation is modeled as a Poisson process in which the average SFR is directly proportional to the present mass of H2. The predictions of our numerical model are in good agreement with the observed Kennicutt-Schmidt and Sigma_H2-Sigma_SFR relations. We show that observations based on CO emission are ill suited to reliably measure the slope of the latter relation at low (<20 M_sun pc^-2) H2 surface densities on sub-kpc scales. Our models also predict that the inferred Sigma_H2-Sigma_SFR relation steepens at high H2 surface densities as a result of the surface density dependence of the CO/H2 conversion factor. Finally, we show that on sub-kpc scales most of the scatter in the relation is a consequence of discreteness effects in the star formation process. In contrast, variations of the CO/H2 conversion factor are responsible for most of the scatter measured on super-kpc scales.


INTRODUCTION
The formation of stars is far from being a well understood and solved problem. The complex interplay between star formation, the formation of molecular clouds, and the physical processes operating within the interstellar medium (ISM), such as turbulence, magnetic fields, self-gravity, or feedback from the stellar population, poses many formidable challenges to the proper interpretation and modeling of the star formation process (e.g., McKee & Ostriker 2007). Fortunately, observations have revealed a number of empirical relations which provide some guidance for the development of a theoretical model. This includes the well-studied Kennicutt-Schmidt relation, a correlation between the surface density of the neutral ISM and Σ SFR (Schmidt 1959;Kennicutt 1989Kennicutt , 1998. The realization that star formation is more closely correlated with molecular gas than neutral gas (Talbot 1980;Wong & Blitz 2002, but cf. Kennicutt 1989 led to considerable interest in the study of the Σ H2 − Σ SFR relation (e.g., Heyer et al. 2004;Gao & Solomon 2004;Greve et al. 2005;Krumholz & Thompson 2007;Kennicutt et al. 2007;Bigiel et al. 2008;Robertson & Kravtsov 2008;Gnedin et al. 2009;Krumholz et al. 2009;Gnedin & Kravtsov 20101 Bigiel et al. 2011;Schruba et al. 2011;Rahman et al. 2011Rahman et al. , 2012. The properties of the Σ H2 − Σ SFR relation are tied to the small scale (scales of molecular clouds and below) connection between star formation and the supply in form of molecular gas. Therefore, measuring the slope and scatter 5 of the Σ H2 − Σ SFR relation allows to probe potentially relevant physical mechanisms that are involved in the star formation process. This requires, of course, that biases in the observational tracers used to derive Σ H2 and Σ SFR are properly accounted for. The main goal of this paper is therefore to study and isolate the importance of both various physical mechanisms and observational biases for the slope and scatter of the Σ H2 − Σ SFR relation.
Our theoretical predictions are based on high resolution, cosmological galaxy formation simulations equipped with subgrid models to estimate H 2 abundances and SFRs. We study observational biases by post-processing our simulations with a model for the J = 1 → 0 line emission of CO (Feldmann et al. 2012; paper I from now on). CO line emission is the most commonly used tracer to measure H 2 abundances (e.g., Wilson et al. 1970;Scoville & Sanders 1987;Brown & Vanden Bout 1991;Young et al. 1995;Regan et al. 2001;Helfer et al. 2003;Solomon & Vanden Bout 2005;Tacconi et al. 2006;Riechers et al. 2006;Kuno et al. 2007;Tacconi et al. 2010;Genzel et al. 2010;Daddi et al. 2010a,b;Ivison et al. 2011).
Detecting H 2 in emission is difficult due to the high energy gap between the ground state and the lowest excited levels of H 2 and the lack of a permanent electric dipole moment of the H 2 molecule. Our CO emission model allows us to determine the CO/H 2 conversion factor (also called the X-factor or X CO ) which is formally defined as i.e., as the ratio between H 2 column density, N H2 , and the CO velocity integrated intensity, W CO , of the J = 1 → 0 rotational transition. As shown in paper I, the X-factor depends on the conditions of the ISM, in particular on the dust-to-gas ratio and the gas surface density. Hence, a crucial question that we address in this paper is whether the use of a constant value for X CO , as done in many observational studies, introduces significant biases in the estimates of the slope and scatter of the Σ H2 −Σ SFR relation. Similarly, we use stellar population synthesis modeling to compute the FUV and Hα fluxes of our model galaxies in order to assess whether the use of these star formation tracers leads to observational biases. From a theoretical perspective, the slope of the Σ H2 − Σ SFR relation provides us with a means to distinguish between (and potentially rule out) different models of star formation. In general, the slopes of the Σ H2 − Σ SFR relation and that of the underlying relation between the spatial densities of SFRs and H 2 (Guibert et al. 1978;Talbot 1980), the ρ H2 −ρ * relation, could be different due to variations in the gas scale height (Talbot 1971) or changes in the density distribution . However, if the ρ H2 −ρ * relation is linear, so is the Σ H2 − Σ SFR relation (Schaye & Dalla Vecchia 2008;. We note that a linear relation has been interpreted as evidence that star formation occurs in dense clumps within molecular clouds and that the number of such clumps scales with the total amount of molecular gas (Bigiel et al. 2008). An alternative possibility that has been suggested is that the timescale for star formation is controlled by small scale processes, e.g., stellar feedback, ambipolar diffusion, etc. which do not vary much from one location to the next, and hence ensure a steady (when averaged over sufficiently large spatial and temporal scales) conversion of H 2 into stars (Wong & Blitz 2002). In contrast, a slope of ∼ 1.3 can be explained by models that are based on a constant star formation efficiency per free-fall time (Krumholz & McKee 2005;Krumholz et al. 2009). Star formation driven by cloudcloud collisions could lead to even steeper slopes (Tasker 2011).
While the slope of the Σ H2 − Σ SFR relation has been the subject of a number of observational studies, there is, unfortunately, no clear consensus yet on its exact value. While many studies find slopes that are close to unity (∼ 0.8 − 1.3, e.g., Wong & Blitz 2002;Komugi et al. 2005;Bigiel et al. 2008;Blanc et al. 2009;Bigiel et al. 2011;Schruba et al. 2011;Rahman et al. 2011Rahman et al. , 2012, a number of works prefer a steeper slope, e.g., ∼ 1.4 (Heyer et al. 2004;Kennicutt et al. 2007), ∼ 1.6 (Thilker et al. 2007), ∼ 1.2 − 1.9 (Liu et al. 2011). Various factors may contribute to these differences, in-cluding choices in the fitting methodology, how star formation maps are corrected for dust extinction and diffuse emission, the resolution of the observations, or the range in surface densities over which the fit is done. In particular, assumptions about the amount of diffuse emission present in the star formation maps and its subsequent treatment in the data analysis have a strong impact on the derived slope Liu et al. 2011;Rahman et al. 2011). However, if the fitting is restricted to regions of sufficiently high Σ H2 , then the impact of any diffuse emission is minimized and the Σ H2 −Σ SFR relation is found to be close to linear (Rahman et al. 2012).
Besides the slope, the normalization of the Σ H2 − Σ SFR relation is an important parameter that needs to be matched by the theoretical models.
Fortunately, studies that find an approximatively linear slope of the Σ H2 − Σ SFR relation also find little variation of the normalization (Σ H2 /Σ SFR ∼ few Gyr) with galaxy mass or other global galaxy properties (e.g., Braine et al. 2001;Wong & Blitz 2002;Bigiel et al. 2008;Genzel et al. 2010;Bigiel et al. 2011;Bolatto et al. 2011). A notable exception is the claim that the normalization depends on the specific star formation rate (Saintonge et al. 2011), although there is a potential worry that this result is, to some extent at least, an artifact of the large scatter in the Σ H2 − Σ SFR relation and the use of correlated quantities (Rahman et al. 2012). Driven by these observational findings we adopt in this paper a (stochastic) star formation model that is based on a linear ρ H2 −ρ * relation with a constant gas depletion time of a few Gyr.
The scatter in the Σ H2 −Σ SFR relation has received significantly less attention compared to the slope or the normalization. The reason is probably that the scatter depends even more on the specific details of the respective observational survey and the data analysis procedures. The finding that the scatter increases with increasing spatial resolution of a survey is sometimes used to argue that star formation scaling relations "break down" on small scales (Momose et al. 2010;Onodera et al. 2010;Schruba et al. 2010). However, we will argue in this paper that the proper interpretation is that on small scales the intrinsically stochastic nature of the scaling relations becomes simply more evident.
A first step to understand and quantify the different contributors to the scatter was made by . Here, we extend our previous analysis in several ways. We consider the scatter that results from spatial X-factor fluctuations. We also analyze the scatter related to star formation tracers (Hα and F U V ). Finally, we highlight the importance of stochasticity in the star formation process and propose a simple model of star formation that allows a classification of the various mechanisms that contribute to the scatter in the Σ H2 − Σ SFR relation.
The outline of the paper is as follows. In §2 we present the detail of our numerical approach, including the set-up of the simulations ( §2.1), the star formation model ( §2.2), details of sub-grid physics ( §2.3), and the modeling of the H 2 and star formation tracers ( §2.4). We then show in §3.1 that our simulations are able to reproduce the observed Kennicutt-Schmidt and Σ H2 − Σ SFR relations. Next, in §3.2, we discuss observational claims of a nonlinear ρ H2 −ρ * relation based on slope measurements of Note. -The first five columns provide the name, box size, numerical resolution, and the adopted cosmological parameters of each zoom-in simulation. The penultimate column states whether the simulation is run with gas-to-dust ratios and interstellar radiation fields fixed to the corresponding values in the solar neighborhood. The parameter ∆t SF in the last column refers to the average time between individual star formation events, see §2.2. The fiducial SF model uses ∆t SF = 10 7 yr. the Σ H2 − Σ SFR relation at low Σ H2 . Then, in §3.3 we analyze the importance of the X-factor for high redshift galaxies with high gas surface densities. We address the origin of the scatter in the Σ H2 − Σ SFR relation in §3.4. In §4 we discuss our star formation model in the context of various observational constraints. Finally, in §5, we summarize our results and conclusions.

Simulations
All simulations have been run with the Eulerian hydrodynamics + N-body code ART (Kravtsov et al. 1997(Kravtsov et al. , 2002 that uses an adaptive mesh refinement (AMR) technique to increase the resolution selectively in specified regions of interest. We also use the standard method of embedding these regions in layers of subsequently lower dark matter resolution to further reduce the computational cost, but still capture the impact of large scale tidal fields correctly (Katz 1991;Bertschinger 2001).
Simulation MW-fid focusses its computational resources on a Lagrangian region that encloses five virial radii of a MW-sized halo at z = 0 (total mass ∼ 10 12 M ⊙ ) in a 6 Mpc h −1 box. The simulation is started from cosmological initial conditions with the parameters given in Table 1. It is run fully self-consistently down to redshift z = 4. At this point "fixed ISM conditions" are imposed on the simulation as described in §2.3 and it is continued for additional 600 Myr before it is analyzed. By then the high-resolution Lagrangian region harbors a large disk galaxy with a virial mass of ∼ 4.2 × 10 11 M ⊙ and several less massive galaxies. More details on the setup of the MW-fid simulation can be found elsewhere paper I).
The simulations MW-dt5, MW-dt6, and MW-sl2 are spawned from the z = 4 snapshot of MW-fid and continued for additional 600 Myr with fixed ISM conditions. The first two simulations are run with reduces values of ∆t SF , the average time scale between individual star formation events, see §2.2, but are otherwise identical to MW-fid. The run MW-sl2 uses a modified star formation model, see §3.2 for details.
Simulation HZ-csm refines on seven regions within a 25 Mpc h −1 box. By z = 0, these regions have collapsed to halos in the 10 11 −10 13 M ⊙ mass range. The simulation is started at z = 100 and is continued down to z = 1.8 fully self-consistently. The setup of the HZ-csm simulation is discussed in detail in Zemp et al. (2012). Simulation HZ-fid is spawned from the z = 2 snapshot of HZ-csm.
It is continued with fixed ISM conditions for additional 200 Myr (down to z = 1.8) to allow the gas to react to the changes in the ISM properties and to reach its new equilibrium state.
The main properties of our simulations are summarized in Table 1.

Star formation model
Since individual resolution elements in our simulations correspond to (at best) GMC scales, it is clear that we are not able to follow in any realistic detail the formation and evolution of individual bound star clusters, let along individual stars. Hence, our approach is to marginalize over the complexities involved in star formation by describing star formation on a statistical level. Such an approach will have its limitations, of course, and we are not aiming at reproducing observations perfectly, but rather try to isolate and explain some of the trends in the observational data.
There are several lines of evidence suggesting that a statistical description of star formation is a plausible ansatz, at least on the level of star clusters and in relatively quiescently star-forming galaxies. One is the observation that the scaling between Σ SFR and Σ H2 is roughly linear and holds over a wide range of surface densities and galaxy metallicities (Bigiel et al. 2008;Genzel et al. 2010), although it may break down in strongly outof-equilibrium situations such a major galaxy mergers (Teyssier et al. 2010;Bournaud et al. 2011). Another one is provided by the similarity of the shape of the mass function of young star clusters among non-starbursting galaxies (Portegies Zwart et al. 2010), which indicates that the conversion of interstellar gas into star clusters proceeds on average in the same way in different galaxies with different global properties. Also, the existence of a correlation between the total number of clusters and the luminosity of the most luminous cluster in a galaxy can be understood in a purely statistical formation scenario of star clusters (Larsen 2002).
Our model assumes that H 2 mass is a good tracer of the SFR (Pelupessy et al. 2006;Robertson & Kravtsov 2008;Gnedin et al. 2009;Papadopoulos & Pelupessy 2010;Kuhlen et al. 2012) and can be summarized by the following three statements: (1) star formation is a stochastic process and the number of individual star formation events 6 in a given time interval is described by a homogeneous Poisson process, (2) the average stellar mass formed per unit time is proportional to the present H 2 mass, and (3) the factor of proportionality, τ −1 dep , is a constant. Here, average refers to the ensemble average over independent patches of the ISM of the same size and with the same H 2 masses, which may, however, differ in other properties. In other words, we marginalize over most of the GMC parameters which may be potentially relevant for star formation and keep only the explicit dependence of the SFR on H 2 mass.
We give a detailed description of the model in the appendix. Here, we only summarize the properties that we will need later in the paper.
The discreteness of star formation leads to scatter in the Σ H2 − Σ SFR relation because the SFR that is realized in a given region and over a given period differs from the ensemble average SFR tied to the present H 2 mass. The scatter increases (up to a point) with increasing average time between individual star formation events, ∆t SF , but decreases with increasing life-time of a given star formation tracer ∆t * . In fact, the scatter depends on these quantities only in form of the dimensionless ratio ∆t * /∆t SF , the average number of individual star formation events during ∆t * .
Our model contains three parameters: the average H 2 depletion time scale, τ dep , the spatial scale on which our model operates, l SF , and ∆t SF . We fix τ dep to 2.9 Gyr in order to fit the observed normalization of the Σ SFR -Σ H2 relation, see §3.1. We stress that τ dep is not the H 2 gas depletion time of individual star forming molecular clouds, but an average depletion time that includes the non star forming molecular gas. Our model is assumed to operate on scales of l SF = 60 − 100 pc scales, which is the peak resolution of our numerical simulations, see Table 1. Our default value for ∆t SF is 10 Myr. This choice is based on the following considerations.
First, the crossing time (or the free-fall time) of a molecular region should be a lower limit on the time period between two bursts of star formation in that region. This time t cross ∼ L/σ can be estimated from the Larson relation between cloud size L and cloud velocity dispersion σ (Larson 1981;Solomon et al. 1987). It scales as ∼ 1 Myr (L/pc) 0.5 and, hence, t cross ∼ 8 − 10 Myr for regions of size l SF = 60 − 100 pc. Since the largest molecular cloud complexes in the Milky Way have radii of ∼ 100 pc (Dame et al. 1986), this line of reasoning cannot be used to constrain ∆t SF on larger spatial scales. Despite having a similar numerical value on ∼ 100 pc scales, ∆t SF and t cross are quite different physical quantities. In fact, t cross , a measure of the duration of an individual star formation event, increases with scale, while ∆t SF , the period between successive star formation events, decreases with scale due to the smoothing effect of spatial averaging on Poisson noise.
Secondly, we have checked that with the choice ∆t SF = 10 Myr individual star formation events in the simulations cover a broad range of masses up to 10 6 M ⊙ . The embedded star clusters that form during a star formation event of total mass 10 6 M ⊙ will have masses up to, but not exceeding, ∼ 10 5 M ⊙ for a reasonable choice of the slope of the embedded cluster mass function (∼ 2.3−2.4). This compares well with the fact that, at least in the is as the simultaneous formation of a number of embedded star clusters in a given region in space (∼ 60−100 pc in our simulations).
Finally, ∆t SF ∼ 10 Myr is consistent with the observed relation between SFR and maximum embedded star cluster mass (Weidner et al. 2004). We note that these considerations do not rule out a moderately larger value, e.g, ∆t SF ∼ 20 Myr.
2.3. Further sub-grid modeling All simulations include a photo-chemical network, metal enrichment from supernova (type Ia and type II), but no thermal energy injection, optically thin radiative cooling by hydrogen (including H 2 ), helium, and metal lines, and 3D radiative transfer of UV radiation. The details of the implementation can be found in Gnedin et al. (2009) and . Here, we give a brief recount.
The photo-chemical network in ART follows the formation and destruction of the five major atomic and ionic species of hydrogen and helium. The formation of H 2 on dust grains and the destruction of H 2 via photodissociation in the Lyman-Werner bands are taken into account. The transfer of ionizing and non-ionizing UV radiation from stellar sources is computed in the OTVET approximation (Gnedin & Abel 2001). Radiative stellar feedback is important for the heating and cooling balance of the gas and for the abundance of H 2 and CO. Unlike the re-ionized intergalactic medium, the dense interstellar medium within galaxies may well be opaque to ionizing photons of all but the nearest stars and radiative transfer effects cannot be neglected.
Most of our simulations are run with "fixed ISM conditions" (see . First of all this means that the dust-to-gas ratio D MW , which is normally assumed to scale linearly with the local gas metallicity, is kept fixed at a value that corresponds to 7 Z ⊙ , and we will denote this as D MW = 1. The gas-to-dust ratio is a crucial parameter that enters the formation rates and the dust shielding of molecular hydrogen and also our CO emission model. The metallicities of the selfconsistently enriched gas are still used to compute the cooling rates that enter the hydrodynamical solver. This is done in order to avoid numerical artifacts such as sudden increases in the gas accretion rates resulting from excess, non-equilibrium cooling. Furthermore, in the fixed ISM runs, the normalization of the radiation field at 1000Å is fixed to J MW = 10 6 photons cm −2 s −1 sr −1 eV −1 , a value typical for the solar neighborhood in the Milky Way (Draine 1978;Mathis et al. 1983). We use the notation U MW = 1, where U MW is the intensity of the radiation field at 1000 A in units of J MW . We stress that only the normalization of the radiation field computed with the OTVET solver is fixed. The shape of the radiation spectrum is not modified.
2.4. Postprocessing: CO, FUV and Hα emission 7 In this paper Z ⊙ refers to the metallicity of the solar neighborhood. Specifically, Z ⊙ = 0.02 or 12 + log 10 (O/H) = 8.92, which is somewhat larger than the metallicity of the Sun according to recent estimates (Allende Prieto et al. 2001;Asplund et al. 2004Asplund et al. , 2009  -Relation between the molecular gas (neutral gas) surface density Σ H 2 (Σ HI+H 2 ) and the surface density of the star formation rate Σ SFR as predicted by our fiducial SF and CO emission model for Milky-Way like ISM properties (D MW = 1, U MW = 1). The left panel shows the Σ H 2 vs Σ SFR relation, the right panel the classical Kennicutt-Schmidt relation (Σ HI+H 2 vs Σ SFR ). The spatial resolution is chosen to match as closely as possible the resolution of the observational studies (Bigiel et al. 2008. The gas surface densities include a factor 1.36 that accounts for the presence of Helium. The red dot-dashed line shows the median Σ SFR (time averaged over the past 20 Myr) for a given "true" H 2 or total gas surface density (as computed in the simulation). The blue solid line shows instead the median of Σ SFR as a function of the "inferred" H 2 or total gas surface density, i.e., a gas density in which Σ H 2 is inferred from the CO intensity, as predicted by our model, using the galactic X-factor. The magenta hashed region indicates the typical scatter in Σ SFR (25 and 75-th percentiles in the left panel, 16 and 84-th percentiles in the right panel, ). The black dashed line (best fit) and the contour lines (containing 90% and 50% of the data) in the left panel are observational measurement of the Σ H 2 -Σ SFR relation by Bigiel et al. (2011) on kpc scales. The contour line in the right panel shows the distribution of sub-kpc sized patches in the sample of nearby galaxies by Bigiel et al. (2008), see their Fig. 7. The figure demonstrates that our modeling of the star formation and the CO emission is consistent with the observed relations between Σ SFR and the molecular and neutral gas surface density, respectively.
The J = 1 → 0 12 C 16 O emission is computed as described in paper I. In short, a sub-grid model for the CO emission is constructed based on the results of a suite of small scale magneto-hydrodynamical ISM simulations (Glover & Mac Low 2011). This model contains two free parameters.
One is the CO brightness temperature that is related to both the temperature of the CO emitting gas and the temperature of the cosmic microwave background (CMB). In §3.1, §3.2, and §3.4 we study the Σ H2 -Σ SFR in the local Universe and we adopt a gas temperature of 10 K (a typical temperature of molecular clouds in the Milky Way). The corresponding brightness temperature is 6.65 K. In §3.3 we predict the CO emission for galaxies at z ∼ 2. We assume that the increase in the CMB temperature at those redshifts is compensated by an increase in the gas temperature from 10 K to 14.5 K, such that the brightness temperature remains approximatively constant. Such a moderate increase in gas temperature is consistent with the detailed modeling of the gas temperature in non-starbursting high redshift galaxies using photon-dissociation regions codes coupled with a semi-analytic galaxy evolution model (Lagos et al. 2012).
The other free parameter is the scaling of the CO line width (either a constant line width or a virial scaling). By default we show results for the virial line width scaling, which is probably the more realistic of the two (see Fig. 3 of paper I), but if relevant we will point out the changes that result from the assumption of a constant line width. The sub-grid model is applied to the highest refined resolution elements in the simulation. The contributions from these individual, ∼ 60 − 100 pc sized resolution elements are then combined in the optically thin limit to derive the CO emission from larger regions.
The FUV and Hα emission of each stellar particle is computed with Starburst-99 (Leitherer et al. 2010) assuming solar metallicity, high mass loss Geneva tracks. We checked that switching to, e.g., Padua tracks does not affect any of our results in a significant way. The Hα luminosities are directly taken from the Starburst-99 output, while the broad-band FUV luminosities are computed using the Galex FUV transmission curve (Morrissey et al. 2005) and the UV spectra provided by Starburst-99.

Actual and inferred star formation relations
It has been shown (Gnedin et al. 2009;Krumholz et al. 2009;Gnedin & Kravtsov 2010) that an H 2 -based star formation prescription is able to reproduce the relation between the surface densities of neutral gas and SFRs (the Kennicutt-Schmidt relation) both for galaxies in the local Universe and for galaxies at high redshift. In fact, the connection between star formation and molecular hydrogen provides a simple physical interpretation for the drop in the SFR surface densities at low gas surface densities (Robertson & Kravtsov 2008). In this  Fig.1, but for a SF model in whichρ * ∝ ρ 2 H 2 . The left panel uses surface densities measured on kpc scales, while the right panel shows the corresponding results for spatial averaging scales of 250 pc. This SF model leads to a steepening of the inferred Σ H 2 − Σ SFR relation at sufficiently high H 2 surface densities ( 20 M ⊙ pc −2 ). However, this steepening is suppressed for Σ H 2 < 20 M ⊙ pc −2 , especially when measured at sub-kpc spatial resolution, as a result of the increase of X CO with decreasing Σ H 2 at low surface densities (paper I). It will therefore be challenging to extract an unbiased estimate for the slope of the relation between SFR and gas density from CO observations of MW-like galaxies at low surface densities.
picture the drop is a manifestation of a transition between a neutral and a molecular hydrogen phase. Specifically, above a characteristic gas surface density the gas is shielded from the interstellar radiation field by a sufficiently large dust optical depth and the molecular phase prevails. The characteristic surface density depends on the dust-to-gas ratio, the strength of the interstellar radiation field, and the density structure of the ISM, see, e.g, .
However, this interpretation of the Kennicutt-Schmidt relation neglects a potentially relevant detail. While theoretical models can predict H 2 surface densities directly, observations often rely on CO observations to infer H 2 surface densities. A crucial test of our understanding of the Kennicutt-Schmidt relation is therefore whether this agreement still holds even if we take the effects of the CO/H 2 conversion factor into account. In other words, we need to compare theory and observations on an equal footing, i.e., using H 2 surface densities that are derived from CO emission in both cases. Our star formation model contains the gas depletion time τ dep as a free parameter. This parameter is the conversion factor between H 2 mass and the ensemble average SFR, see §2.2, but, since we assume τ dep =const, it also corresponds to the normalization of the Σ H2 − Σ SFR relation. We find that a depletion time τ dep of 2.9 Gyr leads to an inferred (CO based) depletion time of ∼ 2.35 Gyr (including Helium, Bigiel et al. 2011) and, thus, to excellent agreement between predictions and observations, see Fig. 1a. We note that the proper choice of τ dep is largely degenerate with the value of the X-factor. Specifically, our CO emission model prefers a median X CO value for a Milky Way like ISM that that is ∼ 65% larger than the value X CO,MW = 2 × 10 20 cm −2 K −1 km −1 s used by Bigiel et al. (2011), explaining the need for a somewhat larger H 2 depletion time. Fig. 1a shows the inferred (CO based) Σ H2 − Σ SFR relation based on our simulation MW-fid. It is an approximatively linear relation over two orders of magnitude in H 2 surface density. This is not an entirely obvious result since, at least on sufficiently small scales, the CO/H 2 conversion factor is strongly dependent on surface density (paper I). However, it turns out that the spatial averaging over a large set of regions with different X CO values erases most of this surface density dependence on kpc and larger scales. Hence, the actual (linear) slope of the Σ H2 − Σ SFR relation is recovered (but see §3.3).
We show our predictions for the Kennicutt-Schmidt relation in Fig. 1b, finding good agreement between our theoretical predictions and observations of the inferred (CO based) Kennicutt-Schmidt relation. We thus conclude that a star formation model that depends linearly on the abundance of H 2 is consistent with both the observed Kennicutt-Schmidt relation and the observed Σ H2 − Σ SFR relation of normal star-forming galaxies in the local Universe.
3.2. The slope of the Σ H2 − Σ SFR relation The slope of the Σ H2 −Σ SFR relation has been a subject of significant debate over the last years. While a number of studies find a slope of about unity, other observational works favor steeper slopes. It has been argued that these discrepancies arise from observational obstacles such as diffuse emission in the infrared, or uncertainties in the dust absorption that make it challenging to obtain reliable SFR estimates, especially in regions of low star formation activity and gas surface density (Liu et al. 2011;Rahman et al. 2011;Leroy et al. 2012).
A varying CO/H 2 conversion factor is another of those potential obstacles, but has been largely neglected since there is no general observational handle on this quantity in different galaxies. The aim of this section is to test the implications for the observed Σ H2 − Σ SFR relation if thė ρ * − ρ H2 relation is highly non-linear and the X-factor is taken into account. In the following we assume that any observational uncertainties related to the estimates of SFR (e.g., diffuse infrared emission, dust absorption, etc.) are under control and corrected for. We focus instead on the role of the CO/H 2 conversion factor.
For this test we use simulation MW-sl2, which uses the same general approach to star formation as presented in §2.2, but with a density dependent gas depletion time.
The results of this test are shown in Fig. 2. While the actual Σ H2 − Σ SFR relation is very steep with a slope n ∼ 2, the inferred (based on CO observations) Σ H2 − Σ SFR relation flattens significantly at Σ H2 20 M ⊙ pc −2 . This effect is already pronounced on kpc scales, but becomes very strong on scales of 250 pc. In other words, we predict that the H 2 gas depletion time derived from CO emission would appear to decrease, not increase, with decreasing surface density at Σ H2 20 M ⊙ pc −2 on sufficiently small scales. This effect is a consequence of the anti-correlation between X CO and Σ H2 at low gas surface densities, caused by the presence of large amounts of CO-dark molecular gas (Wolfire et al. 2010;Krumholz et al. 2011;Narayanan et al. 2011b;paper I;Shetty et al. 2011;Narayanan et al. 2012).
For 20 M ⊙ pc −2 Σ H2 100 M ⊙ pc −2 , on the other hand, the X-factor plays only a small role and the steep slope (n ∼ 2) is recovered. At such gas surface densities the diffuse emission in the star formation maps is also expected to be less of an obstacle (Rahman et al. 2012). Hence, the good news is that a robust measurement of the slope of the Σ H2 − Σ SFR relation based on CO data will be possible in areas of moderately high gas and SFR surface density. At low Σ H2 , however, and especially on sub-kpc scales, an accurate determination of the slope on the basis of CO observations will be difficult.
3.3. Comparison with z ∼ 2 galaxies Studies of the Σ H2 − Σ SFR relation that are based on samples of galaxies from both the local Universe and from high redshift indicate that the slope is close to (but not quite) linear over a many orders of magnitude of Σ H2 . In particular, if one selects only non-interacting galaxies, the slope of the Σ H2 − Σ SFR relation is approximatively ∼ 1.1 − 1.2. This, however, appears to be driven by a steepening at high surface densities, because the slope measured on intermediate H 2 surface densities is much closer to unity (Bigiel et al. 2008).
In our analysis we focus on non-interacting galaxies, that presumably form stars in an approximatively steady state. We do not include interacting or merging galaxies in our discussion, as such objects can be scattered off the Σ H2 − Σ SFR relation Daddi et al. 2010b).
While there are several ways to explain a non-linear slope of the Σ H2 − Σ SFR relation, we demonstrate below that our star formation and CO emission model coupled with an actually linear Σ H2 −Σ SFR relation does predict a slightly non-linear slope of the inferred (CO-based) Σ H2 − Σ SFR relation at high Σ H2 , consistent with observations. In essence, we argue that the observed super-linear slope can be understood as an X-factor effect. We note that the following predictions rest on the assumption that there is no significant trend of the CO line ratios or the CO brightness temperature with gas surface density. Our predictions are specifically for the J = 1 → 0 rotational transition line of CO, while galaxies at z ∼ 1 − 2 are typically observed in CO line emission resulting from J = 2 → 1 or J = 3 → 2 transitions. Systematic trends in the emission line ratios could thus modify the slope of the Σ H2 − Σ SFR relation. A systematic variation of the brightness temperature with surface density would have a similar effect. So far there is no clear observational evidence that either of these assumptions is violated in steady-state, normally star forming galaxies at z ∼ 1 − 2.
Our predictions for the inferred (CO based) Σ H2 −Σ SFR relation at z ∼ 2 are shown in Fig. 3a. Galaxies with inferred molecular gas densities Σ H2 100 M ⊙ pc −2 are shifted off the actually linear Σ H2 − Σ SFR relation. The origin of this offset is the small, but systematic, variation of the X-factor with surface density.
As discussed in paper I, at high Σ H2 the CO emission from a small ISM patch (∼ 20 − 100 pc) ceases to scale linearly with the gas column density (and thus the H 2 surface density) due to the increased optical thickness at the line center of the CO emission line. This effect is somewhat, but not fully, compensated by an increase in the width of the emission line, so that overall there is a remaining increase of the X-factor with increasing H 2 surface density. For this reason, the use of a constant CO/H 2 conversion factor leads to systematic shifts in the inferred Σ H2 −Σ SFR relation. At intermediate H 2 surface densities (the precise range depends on the spatial scale) the conversion factor is essentially constant and, hence, in this case the measured slope is predicted to be very close to linear, as observed (Bigiel et al. 2008).
Star forming galaxies at higher redshifts typically have lower Z than galaxies of a similar mass in the local Universe, e.g., Maiolino et al. (2008). Our fully selfconsistent cosmological simulation HZ-csm predicts that metallicities are only ∼ 0.5 Z ⊙ at z ∼ 2. Hence, D MW ∼ 0.5 and one may wonder whether the combination of a high redshift, low metallicity sample and a low redshift, high metallicity sample could explain the observed superlinear relation. In order to test the importance of the metallicity dependence of CO/H 2 conversion factor we rerun the HZ-csm simulation with fixed Milky-Way like ISM conditions (simulation HZ-fid, see §2.1). Fig. 3b shows that galaxies lie along the observed, slightly super-linear Σ H2 − Σ SFR relation, even if they had D MW = 1. We therefore conclude that the superlinear slope is caused primarily by the scaling of X CO with H 2 surface density, and that metallicity and dustto-gas ratio variations have only a small effect.
We stress again that our CO emission model rests on a number of assumptions as pointed out above. If either of these assumptions were broken, an alternative explanation for the super-linearity of the Σ H2 − Σ SFR relation The Σ H 2 − Σ SFR relation as predicted by a fully cosmological, hydrodynamical simulation at z ∼ 2 (simulations HZ-csm, HZ-fid). Each panel shows the predictions of the "inferred" Σ H 2 − Σ SFR relation, i.e., the one where H 2 surface densities are derived from 12 CO emission maps using a galactic conversion factor X CO,MW = 2 × 10 20 K −1 cm −2 km −1 s. (Left) the interstellar radiation fields and metallicities as computed self-consistently within the simulation at z = 1.8. The H 2 weighted mean metallicity at this epoch is only ∼ 0.5 Z ⊙ but spans a broad range from (1-σ) 0.25 to 0.75 Z ⊙ . This results in a value of X CO that is above the galactic conversion factor and hence shifts the median of the inferred Σ H 2 − Σ SFR relation slightly toward lower surface densities when compared with the observations by Bigiel et al. (2011) for galaxies in the local Universe. The more interesting result is, however, that the slope of the Σ H 2 − Σ SFR relation appears to steepen, in particular at high Σ SFR . (Right) The simulation is restarted at a slightly earlier epoch and continued for ∼ 200 Myr down to z = 1.8, but this time with dust-to-gas ratios and UV radiation fields fixed to D MW = 1 and U MW = 1, respectively. The steepening of the Σ H 2 − Σ SFR relation remains visible and is therefore not a result of changes in the dust-to-gas ratios or interstellar radiation fields. Symbols and lines are as in Fig. 1b. In addition, the gray dashed line shows the fit to the observed Σ H 2 − Σ SFR relation based on a large sample of low and high-z galaxies by Genzel et al. (2010). The red triangles mark the individual positions of simulated galaxies with stellar masses exceeding 10 10 M ⊙ . All simulation predictions are based on our fiducial SF and CO model. The latter assumes a virial scaling of the CO line width. The red arrow at the top indicates the median shift in the inferred H 2 column density of the simulated galaxies if the CO line width would be fixed to a constant value of 3 km s −1 . The figure shows that galaxies with high gas or SFR surface densities appear to deviate from a linear Σ H 2 − Σ SFR relation, consistent with the observations of Genzel et al. (2010), despite the fact that the underlying relation between SFR and H 2 mass is perfectly linear. The super-linear slope (∼ 1.1 − 1.2) of the inferred Σ H 2 − Σ SFR relation is caused primarily by the increase of X CO with increasing Σ H 2 at high gas column densities. Our results apply only to galaxies that are in an equilibrium mode of star formation, not to starbursting galaxies. In the latter environments our CO model becomes unreliable and the galaxy-wide ratio between total molecular gas and the for SF relevant dense molecular gas (n > 10 4 cm −3 ) may change (Gao & Solomon 2004;Lada et al. 2012;Papadopoulos et al. 2012).
would be required.
Systematic variations of the CO line ratios or the CO brightness temperature with gas surface density could be responsible for a difference between the intrinsic slope of the Σ H2 − Σ SFR relation and the slope derived from CO observations (Narayanan et al. 2011a(Narayanan et al. , 2012. Alternatively, a variation of the CO line ratios or brightness temperature with redshift could lead to a systematic z dependence of the normalization of the observed Σ H2 − Σ SFR relation. This can produce an artificial trend with surface density if a galaxy sample is used in which Σ H2 (and Σ SFR ) strongly correlate with galaxy redshift. Finally, the Σ H2 − Σ SFR relation may actually get steeper at high column densities. For instance, it has been suggested that star formation becomes more efficient at high column densities because external pressure on molecular clouds shifts the balance between gravity and turbulent support (Krumholz & McKee 2005). Upcoming observations with the Atacama Large Millimeter Array will hopefully enable us to distinguish between our model and these alternatives (e.g., see Fu et al. 2012).
3.4. The scatter in the Σ H2 − Σ SFR relation There are many effects and processes that could, in principle, contribute to the scatter in the Σ H2 − Σ SFR relation. Clearly, scatter can arise from (1) uncertainties related to the method of estimating SFRs, (2) uncertainties related to the estimation of H 2 masses and surface densities, (3) a possible non-linearity of the star formation process, and (4), any systematic uncertainties in the observables that were not accounted for. We will focus in this paper on the scatter sources (1) and (2). The potential role of (3) is discussed in detail in . We do not attempt to model sources that fall under category (4), since they are not intrinsic but depend on the specifics of the observational survey.
SFR estimates: The stellar mass that was formed over some past time interval will, in general, not coincide with the SFR that is expected based on the present H 2 mass. In other words, as shown in , the use of a time-averaged SFR as an estimator of the ensemble-average SFR introduces scatter. The following (not necessarily distinct) mechanisms fall under this category: • discreteness of star formation: star formation occurs in individual star formation events, i.e., is clustered in time, • stochasticity of star formation: star formation relations on small scales hold only in an (ensemble) average sense.
• fluctuations in the H 2 abundance: H 2 densities and, thus, the ensemble average SFRs may fluctuate on short time scales (Glover & Mac Low 2007), while the SFRs derived from tracers are smoothed because of the inherent time averaging, • evolutionary processes: the conversion efficiency from gas to stars may change over the lifetime of molecular clouds (Murray 2011;cf. Feldmann & Gnedin 2011; see also Schruba et al. 2010;Onodera et al. 2010), • imperfect tracers: observational tracers of SFRs (e.g., Hα luminosities) may provide only approximate estimates of time-averaged SFRs because the proper conversion factor is not known exactly (e.g., it depends on the precise star formation history).
Observations indicate that most stars form in quantized units (embedded star clusters, Lada & Lada 2003) and hence that star formation is a discretized process. While this inevitably makes star formation a stochastic process (in the above sense), it does not mean that all stochasticity in star formation arises from discreteness effects. For instance, the SFR in a molecular cloud depends on more than just its molecular mass. Magnetic fields, cosmic ray density, or the virialization state of the cloud will play a role to some extent. Marginalizing over these additional control parameters will give rise to an apparent stochasticity of star formation. In the formalism of §2.2, the gas depletion time τ dep becomes a function of these additional parameters and its replacement with some appropriately averaged, constant depletion time leads to the appearance of stochasticity. For instance, a star formation efficiency that evolves over the lifetime of molecular clouds can be interpreted in this sense. Here, the age of the cloud is the additional control parameter that determines τ dep .
Our numerical models account for scatter caused by fluctuations in the H 2 abundance, the use of observational SFR tracers, the time discreteness of star formation, and any incidental stochasticity, but do not include the potential contributions from any additional control parameters. H 2 mass estimates: As mentioned above, uncertainties in the estimates of the H 2 surface densities contribute to scatter in the Σ H2 − Σ SFR relation. Hence, variations in the CO/H 2 conversion factor are a potential source of scatter for observational studies that are based on CO emission. As discussed in paper I, the X-factor can vary significantly, even for a fixed dust-to-gas ratio and H 2 column density, because Σ H2 depends on the product of H 2 mass fraction and total gas density, while the CO emission depends on the latter but not the former.

Systematic variations of observables:
We measure the scatter under the condition that the dust-to-gas ratio and the interstellar radiation field are kept fixed (in the sense of §2.3). This is an important point since the Xfactor depends strongly on the former quantity (but only weakly on the latter, see paper I). Hence, variations in the dust-to-gas ratio lead to systematic modulations of conversion factor and, if unaccounted for, to scatter in the Σ H2 − Σ SFR relation.
Often metallicity variations within a given galaxy are a strong function of galacto-centric radius (Searle 1971) and can be large from one galaxy to another. Hence, such systematic X-factor variation will appear as a galaxy-to-galaxy scatter . Since the amount of scatter that is created in this way depends on the particular selection function of the galaxy sample, we do not include galaxy-to-galaxy scatter in this analysis. Our predictions should therefore be compared with observations based on samples of galaxies with approximatively the same dust-to-gas ratio.
In Fig. 4a we plot the scatter of the Σ H2 −Σ SFR relation as a function of spatial averaging scale and separate the contributions that result from the use of time-averaging tracers of star formation and uncertainties in X CO , respectively. The scatter is computed from all regions with an (inferred) H 2 column density between 10 and 100 M ⊙ pc −2 , a CO velocity integrated intensity equal to or larger than 0.2 K km s −1 , and a minimum SFR surface density of 3×10 −4 M ⊙ yr −1 kpc −2 . These limits are chosen to roughly mimic typical values encountered in observational studies and it is clear that the exact numerical predictions will depend to some extent on these limits. In particular, the choice of the minimal Σ SFR is crucial, since the scatter is measured in log Σ H2 − log Σ SFR space and, if not removed, regions with very low star formation would contribute enormously to the scatter (even worse, regions with zero star formation would make the scatter formally infinite).
Our numerical modeling predicts that X-factor variations induce a scatter of the order of ∼ 0.1 − 0.2 dex. Fig. 4a shows that this scatter may be relevant on scales ∼ kpc and above, but on smaller scales the total scatter is primarily due to the use of time averaged SFRs, at least for our fiducial choice ∆t SF = 10 Myr. Again we stress that this analysis assumes that variations in the dust-to-gas ratio and CO brightness temperature are small or accounted for. Fig. 4a also shows that the scatter that results from the use of time-averaged SFRs is a strong function of scale. It reaches ∼ 0.5 dex at ∼ 100 pc, but is only ∼ 0.1 dex at kpc scales.
Furthermore, there is little difference between the use of FUV luminosity-based SFRs and the use of the actual time-averaged SFR over the last 20 Myr. Hence, a varying FUV-to-stellar mass conversion factor (see imperfect tracers above) contributes little to the scatter. This rules out the suggestion by Leroy et al. (2012) that the variation of the FUV luminosity over the lifetime of a single stellar population (SSP) dominates the scatter, at least on scales of ∼ 100 pc and above. In fact, if the luminosity-weighted SFRs differ little from the timeaveraged SFRs, e.g., if SFRs are constant, then the luminosity evolution of the tracer becomes completely irrelevant. This can be seen from the following simple analysis. In both panels the scatter is measured over the range 10 M ⊙ pc −2 < Σ H 2 < 100 M ⊙ pc −2 . Regions with Σ SFR < 3 × 10 −4 M ⊙ yr −1 kpc −2 or I CO < 0.2 K km s −1 are excluded from the analysis. (Left) The blue dot-dashed line shows the scatter that arises when Σ H 2 is know exactly, but SFR are inferred from FUV luminosities. The blue dotted lines marks the analogous result when SFR are derived from the stellar masses formed within the last 20 Myr. Scatter in X CO at fixed CO emission leads by itself to a scatter in the inferred Σ H 2 − Σ SFR relation of the order of 0.1-0.2 dex and is shown as the magenta hashed region. The lower and upper boundaries of this region correspond to the cases of virial scaling of the CO line width vs constant line width, respectively (see text). Fluctuations in X CO are not an important source of scatter on sub-kpc scales (at fixed dust-to-gas ratio and interstellar radiation field), but become increasingly relevant on scales of ∼ kpc and above. Finally, the black horizontally hashed region shows the combined scatter that takes into account both the scatter in X CO and the scatter associated with the estimations of SFRs. (Right) This panel shows how the scatter depends on the SF tracer (FUV, Hα, or simple time averaged SFR) and on assumptions about the stochasticity of the SF process (see legend). The critical parameter is the average time ∆t SF between SF events at a given site within the galaxy (see text). Specifically, the upper 4 lines show the scatter as derived from the various SF tracers for ∆t SF = 10 Myr (our fiducial value), while the two lines just below correspond to ∆t SF = 1 Myr. The gray hashed region shows the scatter for a run with ∆t SF = 0.1 Myr. This scatter results from the mismatch in time scales between SFRs that are averaged over the lifetime of a particular tracer (4 Myr -lower boundary; 20 Myr -upper boundary) and H 2 masses that are observed at a given instant. The figure shows that stochastic effects play a crucial role in determining the overall scatter in the Σ H 2 − Σ SFR relation. Furthermore, modulo X CO effects, the scatter decreases with scale l roughly as a power law ∝ l −α , with α ≈ 0.5 − 0.7, consistent with the findings and interpretation given by  The total spectral luminosity from a Lagrangian volume element at time t is given as Here, φ ν (t) is the spectral luminosity from a SSP of unit mass and age t. This can be rewritten as where are the total spectral energy emitted by an SSP of unit mass and the luminosity weighted SFR, respectively. Hence, if time-averaged and luminosity weighted SFRs trace each other closely, then the constant E ν is the perfect (i.e., scatter-free) conversion factor between tracer luminosity and time-averaged SFR. The validity of this statement does not depend on the form of φ ν . We show in Fig. 4b that the scatter in the Σ H2 − Σ SFR relation depends on both the star formation tracer and the time discreteness parameter ∆t SF . We find that Hα-based SFR estimates lead to more scatter in the Σ H2 −Σ SFR relation than the use of FUV flux as a tracer. Hence, the star formation tracer with the shorter lifetime (Hα) leads to larger scatter. Similarly, when we estimate the SFR based on the actual stellar mass formed within the past 4 Myr and the 20 Myr, we find that the use of a shorter averaging time leads to more scatter in the Σ H2 − Σ SFR relation.
The averaging timescales of 4 and 20 Myr correspond roughly to the luminosity weighted timescales of Hα and FUV emission . It is therefore not entirely surprising that the scatter predictions computed using these time-averaged SFRs are similar to the predictions that use Hα and FUV based tracers. It demonstrates that Hα and FUV based tracers can, to a good degree of approximation, be treated as a top-hat filter with a width of ∼ 4 Myr and ∼ 20 Myr, respectively. This correspondence will break, however, if the scales are small enough and the lifetimes of the particular tracer short enough such that luminosity weighted SFRs and time-averaged SFRs begin to differ substantially. For Hα based SFR estimates this appears to happen on scales of < 400 pc, while for FUV based tracers the effect is small even on spatial averaging scales of ∼ 100 pc. Fig. 4b also shows that the scatter depends on the average time between star formation events ∆t SF . As expected a shorter ∆t SF means that individual star formation events involve less stellar mass but occur at a higher rate, which reduces the scatter. The sharp drop in the scatter when ∆t SF is reduced to 1 Myr proves that much of the scatter in our fiducial ∆t SF = 10 Myr model arises from a single source, the discreteness of star formation. Hence, observational estimates of the scatter can be used to put tight constraints on the value of ∆t SF . In addition, a systematic observational study of the scatter in the Σ H2 − Σ SFR relation as function of ISM environment would allow to determine whether (and how) ∆t SF depends on ISM properties.
In order to assess how much of the total scatter is caused by time variations in the H 2 abundance (and not related to time discreteness of star formation), we also show the scatter for a run with ∆t SF = 0.1 Myr. This timescale is close to the smallest dynamical time step in our simulation and, hence, effectively eliminates any discreteness (beyond that dictated by the simulation time step) in the star formation model. We then measure the scatter in the true Σ H2 − Σ SFR relation using time-averaged SFRs (using 4 Myr and 20 Myr as averaging times). In this way we isolate the scatter that is caused by the observational actuality that SFRs are time-averaged quantities, but H 2 masses are observed at a particular instant. The scatter that results in this way is relatively small, see Fig. 4b. On super-kpc scales it is dominated by the scatter that is caused by X CO fluctuations and on sub-kpc scales by the scatter due to the discreteness of star formation. However, we point out that time variations in the H 2 abundance couple in a non-linear way to the Poisson noise of individual star formation events (see appendix). Hence, they contribute to the scatter caused by the discreteness of star formation and, hence, cannot be neglected.
In  we found that noise inserted by hand on small scales leads to scatter in the Σ H2 −Σ SFR relation that decreases roughly as power law ∝ l −α , where α ≈ 0.5, with increasing spatial averaging scale l. Fig. 4b shows that the scatter due to the stochastic nature of star formation follows this scaling approximatively. In particular, it is clearly less steep than a α = 1 scaling which would be the naive expectation if the gas were arranged in disk of uniform density. As discussed in  the scaling deviates from α = 1 because the density distribution of the ISM, determined by turbulence, is far from being uniform.
How does the scatter that is predicted for our fiducial choice ∆t SF = 10 Myr compare with observations? Such a comparison is not straightforward since observational estimates of the scatter depend on choices in the methodology. For instance, the treatment of diffuse emission does not only affect the inferred slope of the Σ H2 − Σ SFR relation, but also the scatter. Furthermore, the measured scatter depends on the surface density range of the fit. With these caveats in mind we will now compare our predictions to observational studies that give quantitive estimates of the scatter at a given scale. Rahman et al. (2011) infer a scatter of about 0.3-0.4 dex for SFRs based on Hα luminosities on ∼ 0.5 kpc scales. The scatter is lower (∼ 0.1 − 0.3 dex) when they use tracers with longer lifetimes (FUV + 24µm). Our model predicts a scatter of ∼ 0.4 dex for Hα based SFRs and a scatter of ∼ 0.2 dex if FUV luminosities are used to trace star formation on ∼ 0.5 kpc scales. This quantitative agreement is a further justification of the choice ∆t SF = 10 Myr. Verley et al. (2010) investigate how the scatter increases with increasing spatial resolution. Their Fig. 4 shows a scatter of ∼ 0.4 dex on 360 pc scales, ∼ 0.35 dex on 720 pc scales, and ∼ 0.3 dex on 1.4 kpc scales. This is a somewhat shallower scaling than our predictions (∼ 0.4 dex on 500 pc scales, ∼ 0.25 dex on kpc scales and 0.1-0.15 dex on super-kpc scales), see Fig. 4b. However, their scatter is computed by an iterative clipping method that would tend to underestimate the scatter if the scatter is large, i.e., on smaller scales. Schruba et al. (2010) compare the gas depletion times τ CO and τ H2 in apertures centered on peaks of CO and Hα emission, respectively. Unfortunately, they do not report the scatter in the Σ H2 − Σ SFR relation. However, we can compare the change of log 10 (τ CO − τ H2 ), a crude estimator of the scatter, with changing spatial averaging scale using their Fig. 3. This results in a scaling similar to the predictions given in our Fig. 4.
To summarize, we find that for galaxies with Milky-Way like ISM conditions most of the scatter in the Σ H2 − Σ SFR relation is a consequence of the time discreteness of star formation. Systematic variations in the SFR tracer conversion factors are only relevant for tracers with short lifetimes (e.g., Hα) and when observations are done on sufficiently small scales (< 400 pc). Variations of the H 2 /CO conversion factor can dominate the scatter on kpc scales and above, but are unimportant on much smaller scales. Finally, fluctuations in the H 2 abundance play a supporting role, enhancing the scatter caused by the discreteness of the star formation process.

DISCUSSION
How does the star formation model presented in this paper relate to alternative interpretations of the scatter in the Σ H2 − Σ SFR relation?
A commonly made suggestion is that the relation "breaks down" on small scales (Momose et al. 2010;Onodera et al. 2010;Schruba et al. 2010) because molecular clouds pass through evolutionary phases in which they transform from CO-bright, but star-less, clouds to star forming regions with little surrounding molecular gas. One of the difficulties with a picture in which most of the molecular gas in a galaxy goes through a well defined sequence of stages is the following. It does not explain why the H 2 depletion time of molecular clouds with embedded young stellar objects is only a few 100 Myr (Lada et al. 2010), while the depletion time on galactic scales is an order of magnitude larger . In fact, this observation, originally used as evidence to support long lifetimes of molecular clouds (Zuckerman & Evans 1974), implies that the majority of the molecular gas in the galaxy has to be in a non-star forming state (Elmegreen 2000), possibly in either non-star forming clouds that will be dispersed before star formation has a chance to begin, or in unbound molecular associations. The way our model addresses this problem is that the observed gas depletion time of a molecular and star forming region is smaller than the average depletion τ dep by a factor ∼ ∆t * /∆t SF (see appendix). Observations of individual molecular clouds often derive SFRs based on counts of young stellar objects (∆t * = 1 − 2 Myr) or by using tracers of ionizing radiation from massive stars (∆t * ∼ few Myr), e.g., Hα emission or free-free radio emission. Consequently, ∆t * < ∆t SF = 10 Myr and the observed gas depletion time in star forming regions is shorter than τ dep . This is balanced by the large (formally infinite) depletion time in molecular regions that are not currently star forming.
It has also been suggested that the scatter in the Σ H2 − Σ SFR relation is caused by variations in the ratio between molecular gas and dense ( 10 4 cm −3 ) molecular gas (Lada et al. 2012). Since our model does not explicitly follow dense gas, such variations could be included as an additional stochastic component in the star formation model, see §3.4. However, the tight correlation between dense gas and star formation in many different environments (Lada et al. 2012;Papadopoulos et al. 2012) is suggestive of an alternative way. In the context of our model the only thing required is to re-interpret the word "individual star formation event" as "individual dense gas formation event" and to assume that, once gas becomes very dense, star formation is inevitable and will proceed on the ∼ free-fall time of the respective dense gas clump. This modification of our model does not specify the physical mechanism for the sudden increase in gas density, but several plausible options exist, e.g., cloud collisions (Tasker & Tan 2009).
This re-interpretation of our star formation model accounts, by construction, for the observed tight correlation between dense gas and star formation rate. Furthermore, the ratio between dense (i.e., star forming) molecular gas and all molecular gas in a given region will vary depending on how many "dense gas formation events" have occurred in the region. In this model both the scatter in the Σ H2 − Σ SFR relation and the scatter in the mass ratio between dense and not-so-dense gas on small scales are a consequence of the discrete formation of dense gas clumps out of molecular gas.
Further potential contributors to the scatter include the incomplete sampling of the IMF and the drifting of stars out of their parent molecular clouds. These scatter sources are unlikely to be relevant given the relatively large spatial scales ( 100 pc) and SFR surface densities (> 10 −3 M ⊙ yr −1 kpc −2 ) in our study (see Onodera et al. 2010).
We conclude that in the context of our star formation model there is no "break down" of the Σ H2 −Σ SFR scaling relation. Instead, the proper interpretation is that the discrete and stochastic nature of star formation becomes evident as observations probe smaller and smaller scales.

SUMMARY AND CONCLUSIONS
In this paper we studied the slope and the scatter of the Σ H2 −Σ SFR relation using cosmological galaxy formation simulations coupled with models for star formation, H 2 chemistry, and CO emission. Our focus is especially on the role of the CO/H 2 conversion factor. We found that X-factor variations with surface density can result in significant biases of the measured slope. In particular, at high spatial resolution (few 100 pc or better) and sufficiently low surface densities (Σ H2 < 20 M ⊙ pc −2 ) the slope inferred from CO observations is shallower than the actual slope. In contrast, the inferred slope becomes steeper than the true slope if galaxies with high H 2 surface densities (above 100 M ⊙ pc −2 ) are included in the sample, providing an possible explanation for the slightly super-linear slope of the Σ H2 −Σ SFR relation seen at high gas surface densities (e.g., Genzel et al. 2010). Yet, we also showed that measurements at 500 pc resolution over a surface density range often studied in samples of nearby galaxies, 10 M ⊙ pc −2 < Σ H2 < 100 M ⊙ pc −2 , are essentially unbiased.
Variations in the X-factor contribute to the scatter in the Σ H2 − Σ SFR relation (of the order of ∼ 0.1 − 0.2 dex), dominating over many other scatter sources when the spatial resolution of the survey is ∼ kpc or larger. This even holds if there are no significant spatial variations of the dust-to-gas ratio or the interstellar radiation field. Such variations are expected in a heterogeneous sample of galaxies, leading to additional galaxy-to-galaxy scatter with an amount that depends on the properties of the particular galaxy sample (e.g., Schruba et al. 2011). On sub-kpc scales, however, spatial variations in the CO/H 2 conversion factor contribute little to the overall scatter (assuming a fixed dust-to-gas ratio). On such scales much of the scatter is a consequence of the fact that the measured, time-averaged SFRs differ from the SFRs that are expected based on the present amount of H 2 .
We demonstrated that the scatter in the Σ H2 − Σ SFR relation (on scales of ∼ 100 pc and larger) is primarily a consequence of the discreteness of the star formation process. The luminosity evolution of SFR tracers can become relevant for tracer with short lifetimes, e.g., Hα, on small scales (less then a few 100 pc). For FUV-based SFRs, however, the scatter does not differ significantly from the scatter based on a hypothetical tracer with a plain 20 Myr lifetime. The differences in the timescales between SFR tracers (at least a few Myr, up to 100 Myr) and that of H 2 masses (essentially instantaneous measurements) does lead to some scatter, but it is typically dominated by scatter that results from X-factor variations (on super-kpc scales) and by scatter due to the discreteness of star formation (on sub-kpc scales).
The predictions made in the paper suggest a number of observational tests that could be used to constrain the presented numerical models. Fortunately, most of the more obvious tests, e.g., looking for a change in slope of the Σ H2 −Σ SFR relation at very low and very high surface densities (Fig 2b, Fig 3), should be feasible with future ALMA observation. A clear test of the predictions of the Poisson star formation model should be possible with a systematic, observational study of how the scatter in the Σ H2 − Σ SFR relation scales with spatial scale and how it depends on the lifetimes of star formation tracers. Such a study would allow to constrain the discreteness of star formation and, more generally, would be a crucial guide for the development of the theoretical underpinnings of star formation in a galactic context. This work was supported in part by the DOE at Fermilab, by the NSF grant AST-0708154, by the NASA grant NNX-09AJ54G, and by the Kavli Institute for Cosmological Physics at the University of Chicago through the NSF grant PHY-0551142 and PHY-1125897 and an endowment from the Kavli Foundation. The simulations used in this work have been performed on the Joint Fermilab -KICP Supercomputing Cluster, supported by grants from Fermilab, Kavli Institute for Cosmological Physics, and the University of Chicago. This work made extensive use of the NASA Astrophysics Data System and arXiv.org preprint server.

APPENDIX
In this appendix we give a more formal definition of the Poisson star formation model introduced in §2.2 and describe its implementation in the ART code. The model assumes that the number of individual star formation events N ∆t in a given time interval ∆t is a Poisson random variable with a mean and variance λ = ∆t ∆tSF . Here, ∆t SF is the average time interval between two star formation events. The continuum, non-stochastic limit corresponds to ∆t SF → 0.
The ensemble average SFR density in a given region is proportional to the H 2 density, i.e., where the gas depletion time τ dep is assumed to be a constant. The ensemble average SFR in a given fixed volume l 3 SF is simply S = ρ * l 3 SF . We denote the time average of S over some interval ∆t as S ∆t . The actual SFR, S ∆t , that occurred in the volume l 3 SF during the time interval ∆t is a random variable and defined as which ensures that the ensemble average of S ∆t equals S ∆t . This equations shows that the stellar mass formed during an individual star formation event is simply given by S ∆t ∆t SF , i.e., it scales linearly with the average SFR. The time interval ∆t may be some small, arbitrarily chosen time step. We are typically interested in measuring the mean and variance of the SFR (or of its logarithm) over some physical time interval ∆t * ≫ ∆t, e.g., over the lifetimes of star formation tracers such as Hα or FUV emission. The actual SFR over ∆t * is given by with indices i running over the ∆t * /∆t time intervals of length ∆t. It is worthwhile to have a closer look at (3). First, the expression on the right hand side implies that S ∆t * does not depend on the time step ∆t provided ∆t is sufficiently short compared with the typical time over which S fluctuates. Moments of S ∆t * and log 10 S ∆t * depend on the timescales ∆t SF and ∆t * only via the ratio ∆t * /∆t SF . The equation also allows us to estimate the H 2 depletion time ρ H2 l 3 SF /S ∆t * that an observer would measure in a star forming region (N ∆t * = i N ∆t,i ≥ 1). If we assume that the H 2 density in the region remains constant during ∆t * , then the depletion time that the observer would infer is τ dep (∆t * /∆t SF )/N ∆t * , which is smaller than τ dep if ∆t * < ∆t SF .
Equation (3) further shows that S ∆t * depends in a non-linear way on both the time evolution of the H 2 mass in a given region and the number of individual star formation events. Hence, we expect that the scatter in the Σ H2 − Σ SFR relation that arises from the discreteness of star formation depends in a non-linear way on both the time variations in the H 2 mass and on Poisson shot noise in the number of individual star formation events. To be quantitative let us approximate the scatter in the Σ H2 − Σ SFR relation with σ log 10 S∆t * , the standard deviation of the logarithm of S ∆t * > 0. Let us further assume that S ∆t,i fluctuates over time ∆t * as a log-normal random variable, i.e., S ∆t,i ∝ e σX , where X is a Gaussian random variable with mean 0 and variance 1. The choice ∆t * = 20 Myr, ∆t SF = 10 Myr and σ = 1.4 results in σ log 10 S∆t * ≈ 0.60. With the same set of parameters but ∆t SF → 0 we find σ log 10 S∆t * ≈ 0.07. If we assume σ = 0, i.e., S ∆t,i constant, we obtain σ log 10 S∆t * ≈ 0.24. Hence, the variance of log 10 S ∆t * is not simply the sum of the variances caused by Poisson noise and H 2 fluctuations, respectively, but is determined to a large extent by their covariance. Using this simple test setup we also find that the scaling of σ log 10 S∆t * with ∆t * /∆t SF depends on the value of σ, although, for sufficiently large ∆t * /∆t SF , the scatter decreases with increasing ∆t * /∆t SF .
Our simulations adopt the Poisson star formation model in the following way. Every ∆t = 10 5 yr the code computes the ensemble average SFR S in each resolution element (l SF ∼ 100 pc, see Table 1) based on the present H 2 mass. Then, a random realization of N ∆t is drawn from a Poisson distribution with the mean ∆t/∆t SF and, if N ∆t > 0, the code creates a stellar particle of mass N ∆t S∆t SF .