The Ratio of CO to Total Gas Mass in High Redshift Galaxies

Walter et al. (20012) have recently identified the J=6-5, 5-4, and 2-1 CO rotational emission lines, and [C_{II}] fine-structure emission line from the star-forming interstellar medium in the high-redshift submillimeter source HDF 850.1, at z = 5.183. We employ large velocity gradient (LVG) modeling to analyze the spectra of this source assuming the [C_{II}] and CO emissions originate from (i) separate unvirialized regions, (ii) separate virialized regions, (iii) uniformly mixed unvirialized region, and (iv) uniformly mixed virialized regions. We present the best fit set of parameters, including for each case the ratio $\alpha$ between the total hydrogen/helium gas mass and the CO(1-0) line luminosity. We also present computations of the ratio of H_{2} mass to [C_{II}] line-luminosity for optically thin conditions, for a range of gas temperatures and densities, for direct conversion of [C_{II}] line-luminosities to"dark-H_{2}"masses. For HDF 850.1 we find that a model in which the CO and C^{+} are uniformly mixed in gas that is shielded from UV radiation, requires a cosmic-ray or X-ray ionization rate of $\zeta \approx$ 10^{-13} s^{-1}, plausibly consistent with the large star-formation rate ($\sim$ 10^{3} M$_{\odot}$ yr^{-1}) observed in this source. Enforcing the cosmological constraint posed by the abundance of dark matter halos in the standard $\Lambda$CDM cosmology and taking into account other possible contributions to the total gas mass, we find that three of these four models are less likely at the 2$\sigma$ level. We conclude that modeling HDF 850.1's ISM as a collection of unvirialized molecular clouds with distinct CO and C^{+} layers, for which $\alpha$ = 0.6 M$_{\odot}$ (K km s^{-1} pc^{2})^{-1} for the CO to H_{2} mass-to-luminosity ratio, (similar to the standard ULIRG value), is most consistent with the $\Lambda$CDM cosmology.


INTRODUCTION
Observations of high-redshift CO spectral line emissions have greatly increased our knowledge of galaxy assembly in the early Universe. At redshifts z∼2, this includes the discovery of turbulent star-forming disks with cold-gas mass fractions and star-formation rates significantly larger than in present day galaxies (Daddi et al. 2010;Genzel et al. 2012;Magnelli et al. 2012;Tacconi et al. 2010Tacconi et al. , 2012. The high star-formation rates are correlated with large gas masses and luminous CO emission lines (Kennicutt & Evans 2012) observable to very high redshifts. A prominent example is the luminous submillimeter and Hubble-Deep-Field source HDF 850.1, which is at a redshift of z = 5.183 as determined by recent detections of CO(6-5), CO(5-4), and CO(2-1) rotational line emissions, and also [CII] fine-structure emission in this source (Walter et al. 2012). The high redshift of HDF 850.1 offers the opportunity of setting cosmological constraints on the conversion factor from CO line luminosities to gas masses, via the implied dark matter masses and the expected cosmic volume density of halos of a given mass. Such analysis is the subject of our paper.
CO emitting molecular clouds, which provide the raw material for star formation, are usually assumed to have undergone complete conversion from atomic to molecular hydrogen. However, since H2 has strongly forbidden rotational transitions and requires high temperatures (∼500 K) to excite its rotational lines, it is a poor tracer of cold ( 100 K) molecular gas. Determining H2 gas masses in the interstellar medium (ISM) of galaxies has therefore relied on tracer molecules. In particular, 12 CO is the most commonly employed tracer of ISM clouds; aside from being the most abundant molecule after H2, CO has a weak dipole moment (µe = 0.11 Debye) and its rotational levels are thus excited and thermalized by collisions with H2 at relatively low molecular hydrogen densities (Solomon & Vanden Bout 2005).
The molecular hydrogen gas mass is often obtained from the CO luminosity by adopting a mass-to-luminosity conversion factor α = MH 2 /L ′ CO(1−0) between the H2 mass and the J =1-0 115 GHz CO rotational transition (Bolatto et al. 2013). The value for α has been empirically calibrated for the Milky Way Galaxy by three independent techniques: (i) correlation of optical extinction with CO column densities in interstellar dark clouds (Dickman 1978); (ii) correlation of gamma-ray flux with the CO line flux in the Galactic molecular ring (Bloemen et al. 1986;Strong et al. 1988);and (iii) observed relations between the virial mass and CO line luminosity for Galactic GMCs (Solomon et al. 1987). These methods have all arrived at the conclusion that the conversion factor in our Galaxy is fairly constant. The standard Galactic value is α = 4.6 M⊙/(K km −1 pc 2 ). Subsequent studies of CO emission from unvirialized regions in Ultra-Luminous Infrared Galaxies (ULIRGs) found a significantly smaller ratio. For such systems, α = 0.8 M⊙/(K km −1 pc 2 ) (Downes & Solomon 1998). These values have been adopted by many (Sanders et al. 1988;Tinney et al. 1990; Wang et al. 1991;Walter et al. 2012) to convert CO J=1-0 line observations to total molecular gas masses, but without consideration of the dependence of α on the average molecular gas conditions found in the sources being considered. Since all current observational studies of α leave its range and dependence on the average density, temperature, and kinetic state of the molecular gas still largely unexplored, its applicability to other systems in the local or distant Universe is less certain (Papadopoulos et al. 2012).
In this paper, we estimate the CO emitting gas masses in HDF 850.1 using the large-velocity-gradient (LVG) formalism to fit the observed emission line spectral energy distribution (SED) for a variety of model configurations, and for a comparison to the Galactic and ULIRG conversions. In §2 we outline the details of the LVG approach, including an overview of the escape probability method and a derivation of the gas mass from the line intensity of the modeled source. In §3, we present the best fit set of parameters that reproduce HDF 850.1's detected lines and calculate the corresponding molecular gas mass, assuming the CO and [CII] emission lines originate from (i) separate virialized regions, (ii) separate unvirialized regions, (iii)uniformly mixed virialized regions, and (iv) uniformly mixed unvirialized regions. The inferred gas masses enable us to set lower-limits on the dark-matter halo mass for HDF 850.1. In §4 we compare the estimated halo-masses to the number of such objects expected at high redshift in ΛCDM cosmology, and show that models with lower values of α are favored. We conclude with a discussion of our findings and their implications in §5.

LARGE VELOCITY GRADIENT MODEL
We start by describing our procedure for quantitatively analyzing the [CII] and CO emission lines detected at the position of HDF 850.1 using the large velocity gradient (LVG) approximation. We consider a multi-level system with population densities of the i th level given by ni. The equations of statistical equilibrium can then be written as : where l is the total number of levels included; since the set of l statistical equations is not independent, one equation may be replaced by the conservation equation where ntot is the number density of the given species in all levels. In our application, ntot = nCO. Following the notation of Poelman & Spaans (2005), Rij is given in terms of the Einstein coefficients, Aij and Bij , and the collisional excitation (i < j )and de-excitation rates (i > j ) Cij : where Jij is the mean radiation intensity corresponding to the transition from level i to j averaged over the local line profile function φ(νij). The total collisional rates Cij depend on the individual temperature-dependent rate coefficients and collision partners, usually H2 and H for CO rotational excitation. The difficulty in solving this problem is that the mean intensity at any location in the source is a function of the emission and varying excitation state of the gas all over the rest of the source, and is thus a nonlocal quantity. To obtain a general solution of the coupled sets of equations describing radiative transfer and statistical equilibrium, we adopt the approach developed by Sobolev (1960) and extended by Castor (1970) and Lucy (1971) and assume the existence of a large velocity gradient in dense clouds. This assumption is justified given the interstellar molecular line widths which range from a few up to a few tens of kilometers per second, far in excess of plausible thermal velocities in the clouds (Goldreich & Kwan 1974). They suggest that these observed velocity differences arise from large-scale, systematic, velocity gradients across the cloud, a hypothesis that lies in accord with the constraints provided by observation and theory.
In the limit that the thermal velocity in the cloud is much smaller than the velocity gradient across the radius of the cloud, the value of Jij at any point in the cloud, when integrated over the line profile, depends only upon the local value of the source function and upon the probability that a photon emitted at that point will escape from the cloud without further interaction. Thus Jij becomes a purely local quantity, given by: where Sij is the line source function, Detection of four lines tracing the star-forming interstellar medium in HDF 850.1 (Walter et al. 2012) assumed constant through the medium. In this expression gi and gj are the statistical weights of levels i and j respectively, βij is the "photon escape probability", and Bij (νij, TB) is the background radiation with temperature TB. In our models we set TB to the CMB temperature of 16.9 K at z = 5.183. We ignore contributions from warm dust (da Cunha et al. 2013). For a spherical homogenous collapsing cloud, the probability that a photon emitted in the transition from level i to level j escapes the cloud is given by where τij is the optical depth in the line, The equations of statistical equilibrium are therefore reduced to the simplified form and can be solved through an iterative process to give the fractional level populations ni/ntot (for a given choice of densities for the collision partners, usually H2 but also H, and kinetic temperature T kin ). Assuming the telescope beam contains a large number of these identical homogeneous collapsing clouds (e.g., Hailey-Dunsheath 2009), the corresponding emergent intensity of an emission line integrated along a line of sight is then simply where χ is the abundance ratio, χ ≡ ntot/n, where n is the total hydrogen gas volume density (cm −3 ) and N is the hydrogen column density (cm −2 ). Given a series of observed lines with frequency νij , one can identify a set of characterizing parameters that best reproduces the observed line ratios and intensity magnitudes; among these parameters are the cloud's kinetic temperature (T kin ), velocity gradient (dv/dr ), gas density n and collision partner (H and H2) gas fractions, abundance ratio (χ), and column density (N ). For a spherical geometry, this column density can then be further related to the molecular gas mass of the cloud in the following way where the factor µ = 1.36 takes into account the helium contribution to the molecular weight and R = DAθ/2 is the effective radius of the cloud, with DA being the angular diameter distance to the source and θ the beam size of the line observations. N is defined in terms of the column density obtained from the LVG calculation, N ′ = N (1 + z) 4 where the (1 + z) 4 multiplicative factor reflects the decrease in surface brightness of a source at redshift z in an expanding universe.
In the case where the emitting molecular clouds are gravitationally bound, applying the virial theorem to a homogenous spherical body yields the following constraint (Goldsmith 2001 where a = 7.77×10 −3 if n = nH 2 and a = 5.50×10 −3 if n = nH . The velocity gradient is inversely proportional to the dynamical time scale. In models where the clouds are assumed to be virialized, dv/dr and n are no longer independent input parameters of the model, but rather vary according to equation (11). For the optically thick 12 CO J=1-0 transition line (β ≈ 1/τ ), carrying out the LVG calculations with this additional virialization condition leads to a simple relation between the gas mass and the CO(1-0) line luminosity, (12) where the excitation temperature Texc ≈ T kin when the emission line is thermalized. (In this expression we assume complete conversion to H2 so that the gas mass is the H2 mass.) Empirically, the Galactic value of this mass-toluminosity ratio for virialized objects bound by gravitational  n H 2 = 10 cm − 3 n H 2 = 10 2 cm − 3 n H 2 = 10 3 cm − 3 n H 2 = 10 4 cm − 3 Figure 1. The gas mass to [CII] luminosity ratio for a C + to H 2 abundance ratio of χ C + =10 −4 as a function of the kinetic temperature T kin at a molecular hydrogen number density of n H 2 = 10, 10 2 , 10 3 , and 10 4 cm −3 . Calculations were made in the optically thin regime, assuming the fine-structure transition J=3/2→1/2 is due solely to spontaneous emission processes and collisions with ortho-and para-H 2 .

ANALYSIS OF HDF
1900.1 GHz (158 µm), one of the main cooling lines of the star-forming in stellar medium, has also been detected (Table 1). These lines have placed HDF 850.1 in a galaxy over density at z = 5.183, a redshift higher than those of most of the hundreds of submillimetre-bright galaxies identified thus far. Walter et al. (2012) used an LVG model to characterize the CO spectral energy distribution of HDF 850.1 and found that the observed CO line intensities could be fit with a molecular hydrogen density of 10 3.2 cm −3 , a velocity gradient of 1.2 km s −1 pc −1 , and a kinetic temperature of 45 K. Then, assuming α = 0.8 M⊙(K km s −1 pc 2 ) −1 as for ULIRGs, they used the 1-0 line luminosity inferred from their LVG computation to infer that MH 2 = 3.5×10 10 M⊙. However, as Papadopoulos et al. (2012) argue, adopting a uniform value of α for ULIRGs neglect its dependence on the density, temperature, and kinematic state of the gas; this may limit the applicability of computed conversion factors to other systems in the local or distant Universe.
Here, we broaden the LVG analysis carried out in Walter et al. 2012 and use the LVG-modeled column density to estimate the total gas mass of the source. We present several alternative models, each subject to a slightly different constraint. In particular, we first consider the case where the CO and [CII] lines originate from different regions of the molecular cloud. This picture is consistent with the standard structure of PDRs in which there is a layer of almost totally ionized carbon at the outer edge, intermediate regions where the carbon is atomic, and internal regions where the carbon is locked into CO (Sternberg & Dalgarno 1995;Tielens & Hollenbach 1985;Wolfire et al. 2010). In this pic-ture then, the hydrogen is fully molecular in the CO emitting regions. We then consider models in which the CO molecules and C + ions are uniformly mixed, such that the line emissions originate in gas at the same temperature and density. These models resemble conditions found in UV-opaque cosmic-ray dominated dark cores of interstellar molecular clouds where the chemistry is driven entirely by cosmic-ray ionization. For HDF850.1, the ionization rates could be significantly higher than in the Milky Way, leading to enhanced C + in the UV-shielded regions. In both instances, we perform our LVG computations assuming (a) virialized clouds for which the virialization condition (Equation [11]) has been imposed and (b) gravitationally-unbound clouds.
To carry out these computations, we use the Mark & Sternberg LVG radiative transfer code described in Davies et al. (2012). Energy levels, line frequencies and Einstein A coefficients are taken from the Cologne Database for Molecular Spectroscopy (CDMS).The excitation and deexcitation rates of the CO rotational levels that are induced by collisions with H2 are taken from Yang et al. (2010) while the C + collisional rate coefficients come from Flower & Launay (1977) and Launay & Roueff (1977).

Separate CO, C + Virialized Regions
We first consider a model in which the CO and [CII] emission lines detected at the position of HDF 850.1 originate in separate regions of the molecular gas cloud, regions which are not necessarily at the same temperature and number density. For self-gravitating clouds in virial equilibrium, the velocity gradient is no longer an independent input parameter of the LVG model, but varies with nH 2 according to equation (11). To find the unique solution that yields the two observed line ratios, I CO(6−5) /I CO(2−1) and I CO(6−5) /I CO(5−4) , we assume a canonical value of χ CO = 10 −4 for the relative CO to H2 abundance and vary the remaining two parameters, temperature and molecular hydrogen density, over a large volume of the parameter space. We find, under this virialization constraint, that the observed CO lines are best fit with a kinetic temperature of 70 K and a molecular hydrogen number density of 10 2.6 cm −3 (with a corresponding velocity gradient of ≈ 0.16 km s −1 pc −1 ). The column density that yields the correct line intensity magnitudes is 4.2×10 19 cm −2 , corresponding to a molecular hydrogen gas mass of MH 2 ≈ 2.13×10 11 M⊙. The H2 mass to CO luminosity conversion factor obtained in this model is α = 5.1 M⊙(K km s −1 pc 2 ) −1 , a value similar to the Galactic conversion factor observed for virialized molecular clouds in the Milky Way, suggesting that HDF 850.1 may have some properties in common with our Galaxy.
Reducing the relative CO to H2 abundance by a factor of two, to χ CO = 5×10 −5 , results in a best fit solution with a molecular hydrogen gas mass of MH 2 ≈ 2.68×10 11 M⊙, nearly 25% larger than the value obtained assuming χ CO = 10 −4 . Ranges on the fit parameter consistent with the observational uncertainties are listed in Table 2.

Separate CO, C + Unvirialized Regions
We then consider a model in which the CO and [CII] emission lines are assumed to originate from separate regions of gravitationally-unbound molecular clouds. Since unvirialized clouds generally demonstrate a higher degree of turbulence relative to their virialized counterparts, we expect the velocity gradient in this model to be greater than the velocity gradient obtained for the virialized model, (dv/dr) virialized ≈ 0.16 km s −1 pc −1 . We therefore fix the velocity gradient to be ten times the virialized value, (dv/dr) unvirialized = 1.6 km s −1 pc −1 , and, assuming a canonical value of χ CO = 10 −4 , find the solution that yields the two observed line ratios by varying T and nH 2 . We find, under these assumptions, that the CO SED is best fit with a molecular hydrogen density of 10 3 cm −3 and a kinetic temperature of 100 K. For this set of parameters, the beamaveraged H2 column density is NH 2 ≈ 1.0×10 19 cm −2 , giving an associated molecular gas mass of MH 2 ≈ 5.16×10 10 M⊙. This estimate of the gas mass is nearly 50% larger than that obtained by Walter et al. (2012) by applying the H2 massto-CO luminosity relation with the typically adopted conversion factor for ULIRGs, α = 0.8 M⊙ (K km s −1 pc 2 ) −1 . Given our inferred H2 gas mass and predicted CO(1-0) line luminosity from our LVG fit, we find that α = 1.2 M⊙ (K km s −1 pc 2 ) −1 , in this model. The increase in molecular hydrogen density and the reduction in inferred mass (relative to the values obtained in the previous model where the virialization condition was imposed) arise from our assumption that the velocity gradient in this model is greater than the velocity gradient obtained in the virialized case. For a fixed χ, the optical depth drops with increasing dv/dr (Equation [7]); since β, the probability of an emitted photon escaping, correspondingly increases, the radiation is less "trapped" and a higher density, nH 2 , is required to produce the observed CO excitation lines. Furthermore, since Ii,j ∝ βi,j M , a larger β implies that less mass is required to reproduce an observed set of line intensities. Therefore, assuming (dv/dr) unvirialized = 10 (dv/dr) virialized causes the optical depth, and consequently, the inferred mass, to drop by a factor of nearly 4 in this model.

Optically thin [CII]
In the two models above, where the CO and [CII] lines are assumed to be emitted from separate regions of the molecular clouds, the single detected ionized carbon line is insufficient in constraining the parameters of the LVG modeled [CII] region. We thus consider the optically thin regime of the [CII] line (β ≃ 1), such that where χ C + is the C + to H2 abundance ratio (fixed at a value of 10 −4 in these calculations) and xi represents the fraction of ionized carbon molecules in the i th energy level. Assuming that the fine-structure transition J=3/2→1/2 is due solely to spontaneous emission processes and collisions with orthoand para-H2, the equations of statistical equilibrium reduce to a simplified form and can be solved to obtain x (J =3/2) as a function of temperature and H2 number density. In Fig. 1, we have plotted the resulting mass-to-luminosity ratio, MH 2 /L [CII] , as a function of the kinetic temperature for several different values of nH 2 . Given the high [CII]/farinfrared luminosity ratio of L [CII] /LF IR = (1.7±0.5)×10 −3 in HDF 850.1 (Walter et al. 2012), it is reasonable to assume the ionized carbon is emitting efficiently and to thus consider the high temperature, large number density limit (T kin ≃ 500 K, nH 2 ≃ 10 4 cm −3 ). In this limit, the massto-luminosity ratio is ≈ 1.035 M⊙(K km −1 pc 2 ) −1 and the corresponding molecular gas mass of the C + region, using the detected line intensity of the ionized carbon line L [CII] = 1.1×10 10 L⊙, is found to be MH 2 ≈ 5.2×10 10 M⊙ .

Uniformly Mixed CO, C + Virialized Region
We next consider models in which the CO molecules and C + ions are mixed uniformly, such that the corresponding line emissions originate in gas at the same temperature, density, and velocity gradient. For these conditions, the chemistry is driven by cosmic-ray ionization and the density fractions ni/n for CO and C + depend on a single parameter, the ratio of the cloud density to the cosmic-ray ionization rate ζ (Boger & Sternberg 2005). For Galactic conditions with ζ ∼10 −16 s −1 , the C + to CO ratio is generally very small. However, in objects such as HDF 850.1 where the ionization rate may be much larger due to high star formation rates, this ratio may be enhanced significantly.  Figure 2. Dependence of the density ratios, n C + /n CO (solid line) and n H 2 /n H (dashed line), on the H 2 ionization rate ζ at a fixed solar metallicity Z' =1, for our best-fit LVG parameters T kin = 160 K and n H 2 = 10 3 cm −3 ). As ζ increases, the abundance of C + relative to CO in the cosmic-ray dominated dark cores of interstellar clouds grows while that of H 2 to atomic hydrogen decreases. The desired value n C + /n CO ≈ 13 is obtained for ζ ≈ 2.5×10 −14 s −1 , at which point n H 2 ≈ 0.4n H .
C + -C -CO interconversion in a purely ionization-driven chemical medium. For solar abundances of the heavy elements (Z ′ = 1), a fairly reasonable assumption given the high star formation rate observed in HDF 850.1, the cosmic-ray ionization rate required to achieve this high C + to CO ratio in a cloud with temperature T =160 K and density nH 2 =10 3 cm −3 is of order ζ ≃ 2.5×10 −14 s −1 (Figure 2). This is significantly enhanced compared to the Milky Way value, ζ ∼10 −16 s −1 , and is plausibly consistent with the fact that HDF 850.1 has a star-formation rate of 850 solar masses per year, a value which is larger than the measured Galactic SFR by a factor of 10 3 (Walter et al. 2012).
Furthermore, for such a high cosmic-ray ionization rate of this magnitude, the hydrogen is primarily atomic for the implied LVG gas density. We find that the ratio of molecular hydrogen to atomic hydrogen in the interstellar clouds is nH 2 /nH ≈ 0.4 (Figure 2). We thus replace H2 with H as the dominant collision partner. Given the uncertain CO-H rotational excitation rates (see Shepler et al. 2007), we assume that the rate coefficients are equal to the rates for collisions with ortho-H2 (Flower & Pineau Des For 2010). We find that the observed CO and [CII] lines together are best fit with a temperature of 160 K, an atomic hydrogen number density of nH = 10 3 cm −3 (with a corresponding velocity gradient of 0.17 km s −1 pc −1 ), and an abundance ratio (relative to H) of 9.3×10 −5 and 7×10 −6 for C + and CO respectively. The column density that yields the correct line intensity magnitudes is 8.5×10 19 cm −2 , corresponding to an atomic hydrogen gas mass of MH ≈ 2.14×10 11 M⊙. The molecular hydrogen mass is then MH 2 = 2(nH 2 /nH)MH ≈ 1.72×10 11 M⊙ and the total gas mass estimate is Mgas ≈ 3.86×10 11 M⊙. The corresponding conversion factor in this model is α = 9.8 M⊙(K km s −1 pc 2 ) −1 , where α is now defined as the total gas mass to CO luminosity ratio.
Reducing the fixed sum of abundance ratios (relative to H) of CO and C + by a factor of two, to χ CO + χ C + = 5×10 −5 , results in a best fit solution with a total gas mass of Mgas ≈ 4.57×10 11 M⊙, nearly 20% larger than the value obtained assuming χ CO + χ C + = 10 −4 .

Uniformly Mixed CO, C + Unvirialized Region
In the case where we assume gravitationally-unbound molecular clouds, we again find a best fit model with a C + to CO ratio of χ C + /χ CO ≈ 13, indicating that atomic hydrogen is the dominant collision partner in the LVG calculations. We thus calculate a three-dimensional grid of model CO and [CII] lines, varying T kin , nH , and the relative abundances χ [CII] and χ CO , with the constraints that χ CO + χ C + = 10 −4 and (dv/dr) unvirialized = 10(dv/dr) virialized = 1.7 km s −1 pc −1 . The observed set of CO and [CII] lines, assumed to have been emitted from the same region, are fit best with a temperature of 180 K, an atomic hydrogen number density of 10 3.4 cm −3 and an abundance ratio of 9.3×10 −5 and 7×10 −6 for C + and CO respectively. For this set of parameters, the beam-averaged H column density is well constrained to be NH ≈ 2.9×10 19 cm −2 . This corresponds to an atomic and molecular gas mass of MH ≈ 7.33×10 10 M⊙ and MH 2 ≈ 5.20× 10 10 M⊙ respectively, yielding a total gas mass estimate of Mgas ≈ 1.25×10 11 M⊙. The ratio of the a For each model parameter, the top row represents the unique, best fit value obtained for the specified model. The bottom row provides the range of parameter values that yield results consistent with the observed line intensity ratios within the error bars of the observed data points (Walter et al. 2012). b α here is defined as the ratio of total gas mass to CO(1-0) luminosity, α = Mgas/L' CO(1−0) . In the models where the CO and [CII] lines are assumed to be originating from separate regions, Mgas = M H 2 since estimates of the atomic hydrogen gas mass could not be obtained via the LVG calculations. In the models where the CO molecules and C + ions are assumed to be uniformly mixed, the total gas mass is the sum of the molecular and the atomic gas masses, Mgas = M H 2 +M H .
total gas mass to the CO luminosity in this model is α = 5.1 M⊙(K km s −1 pc 2 ) −1 .

COSMOLOGICAL CONSTRAINTS
Our inferred gas masses enable us to set cosmological constraints. For a particular set of cosmological parameters, the number density of dark matter halos of a given mass can be inferred from the halo mass function. The Sheth-Tormen mass function, which is based on an ellipsoidal collapse model, expresses the comoving number density of halos n per logarithm of halo mass M as, (1 + 1 (aν 2 ) p )e −aν 2 /2 (14) where a reasonably good fit to simulations can be obtained by setting A = 0.322, a = 0.707, and p = 0.3 (Sheth & Tormen 1999). Here, ρm is the mean mass density of the universe and ν = δcrit(z)/σ(M ) is the number of standard deviations away from zero that the critical collapse overdensity represents on mass scale M. Integrating this comoving number density over a halo mass range and volume element thus yields N, the expectation value of the total number of halos observed within solid angle A with mass greater than some M h and redshift larger than some z, where H0 is the Hubble constant, DA is the angular diameter distance, and Ωm and ΩΛ are the present-day density parameters of matter and vacuum, respectively. Under the assumption that the number of galaxies in the field of observation follows a Poisson distribution, the probability of observing at least one such object in the field is then P = 1 − F (0, N ) where F (0, N ) is the Poisson cumulative distribution function with a mean of N . Given the detection of HDF 850.1, we can say that out of the hundreds of submillimetre-bright galaxies identified so far, at least one has been detected in the Hubble Deep Field at a redshift z > 5 with a halo mass greater than or equal to the halo mass associated with this source. This observation, taken together with the theoretical number density predicted by the Sheth-Tormen mass function, implies that an atomic model that yields an expectation value N can be ruled out at a confidence level of where the solid angle covered by the original SCUBA field in which HDF850.1 was discovered is ≃ 9 arcmin 2 (Hughes et al. 1998) and a ΛCDM cosmology is assumed with H0 = 70 km s −1 Mpc −1 , ΩΛ = 0.73, and Ωm = 0.27 (Komatsu et al. 2011). M h,min , the minimum inferred halo mass for HDF 850.1, is related to the halo's minimum baryonic mass component, a quantity derived in §3 via the LVG technique, in the following way, where the baryonic and the total matter density parameters are Ω b = 0.05 and Ωm = 0.27 respectively. Each model's estimate of the minimum baryonic mass associated with HDF 850.1 therefore corresponds to an estimate of N(5,M h,min ) and respectively yields the certainty with which the model can be discarded. The confidence with which models can be ruled out on this basis is plotted as a function of the minimum baryonic mass estimated by the model (solid curve in Figure 3). To check the consistency of the four LVG models considered in §3 with these results, dashed lines representing the masses derived from each model are included in the plot (upper left panel). Assuming the CO and [CII] molecules are uniformly mixed in virialized clouds results in a baryonic mass M b,min = 3.86×10 11 M⊙. The probability of observing at least one such source, with a corresponding halo mass M h 2.1×10 12 M⊙, is ∼ 7×10 −2 ; this model can thus be ruled out at the 1.8σ level (solid). The model which postulate separate virialized regions (dashed) can be ruled out with relatively less certainty, at the 1σ level. On the other hand, modeling the CO and [CII] emission lines as originating from mixed (dotted) or separate (dash-dot) unvirialized regions, results in minimum baryonic masses which are consistent with the constraint posed by equation (15). We expect to find N ∼ 1.5 and 10 such halos, respectively, at a redshift z 5. The fact that the expected average number of observed sources for the latter model is much higher than the actual number of sources observed may be due to the incompleteness of the conducted survey and is therefore not grounds for ruling out this model.
The baryonic masses, obtained via the LVG method, represent conservative estimates of the total baryonic content associated with HDF 850.1. In particular, mass contributions from any ionized gas or stars residing in the galaxy are not taken into account, and in the models where a layered cloud structure is assumed, the atomic hydrogen mass component is left undetermined. Furthermore, ejection of baryons to the IGM through winds may result in a halo baryon mass fraction that is smaller than the cosmic ratio, Ω b /Ωm, used in this paper, resulting in a conservative estimate of the minimum halo mass. We therefore consider the effects on the predicted number of observed halos, N(M h,min ), if the minimum baryonic mass derived for each model is doubled, e.g. assuming a molecular-gas mass fraction of ∼ 1/2 (Tacconi et al. 2013). Using these more accurate estimates of HDF 850.1's baryonic mass, we find that the two models in which the virialization condition is enforced can be ruled out at the ∼ 2.7σ and 2σ level for the mixed and separate models, respectively (upper right panel).
For comparison, we also consider the gas masses implied by the CO(1-0) line luminosities predicted by the two models where distinct CO and C + layers are assumed (lower left panel). Adopting a CO-to-H2 conversion factor of α = 0.8 (K km s −1 pc 2 ) −1 , we find that enforcing the cosmological constraint posed by the abundance of dark matter halos does not rule out either of the two models in which separate CO and C + regions were assumed, even if the minimum baryonic masses are doubled to account for neglected contributions to the total gas mass (lower right panel). If a conversion factor of α = 4.6 (K km s −1 pc 2 ) −1 is used, both models can be ruled out at the ∼ 1σ level and increasing these obtained masses by a factor of two drives up the sigma levels to ∼ 1.8σ for both models.

SUMMARY
In this paper, we employed the LVG method to explore alternate model configurations for the CO and C + emission lines regions in the high-redshift source HDF 850.1. In particular, we considered emissions originating from (i) separate virialized regions, (ii) separate unvirialized regions, (iii) uniformly mixed virialized regions, and (iv) uniformly mixed unvirialized regions. For models (i) and (ii) where separate CO and C + regions were assumed, the kinetic temperature, T kin , and the molecular hydrogen density, nH 2 , were fit to reproduce the two observed line ratios, I CO(6−5) /I CO(2−1) and I CO(6−5) /I CO(5−4) , for a fixed canonical value of the CO abundance (relative to H2), χ CO = 10 −4 . The column density of molecular hydrogen, NH 2 , was then fit to yield the correct line intensity magnitudes and the molecular gas mass was derived for each respective model. In models (iii) and (iv) where the CO molecules and C + ions were assumed to be uniformly mixed with abundance ratios that satisfied the constraint, χ CO + χ C + = 10 −4 , we found that a relatively  for each model were doubled to account for neglected contributions to the total gas mass; models (a) and (c) can now be ruled out at the 2σ and 2.7σ levels respectively. Lower left panel : The vertical lines represent the mass values implied by the predicted CO(1-0) line luminosities from models (a) and (b) with a CO-to-H 2 conversion factor of α = 0.8 and 4.6 M ⊙ (K km s −1 pc 2 ) −1 . If these minimum baryonic mass values are then doubled (lower right panel ), both models can be ruled out at the ∼ 1.8σ level in the case where α = 4.6 M ⊙ (K km s −1 pc 2 ) −1 is adopted as the conversion factor.
high ionization rate of ζ ≃ 2.5×10 −14 s −1 is necessary to reproduce the set of observed line ratios, {I [CII] /I CO(2−1) , I [CII] /I CO(5−4) , I [CII] /I CO(6−5) }. Since the hydrogen in a cloud experiencing an ionization rate of this magnitude is primarily atomic, two additional parameters, nH and NH , were introduced and the set of LVG parameters, {T kin , nH 2 , nH , χ CO , χ C + , NH 2 , NH }, were fit and used to obtain both a molecular and an atomic hydrogen gas mass for each model. The gas masses derived by employing the LVG technique thus represent conservative estimates of the minimum baryonic mass associated with HDF 850.1.
These estimates were then used, together with the Sheth-Tormen mass function for dark matter halos to calculate the average number of halos with mass M h M b,min that each model predicts to find within the HDF survey volume. Given that at least one such source has been detected, we found that models (i) and (iii) can be ruled out at the 1σ and 1.8σ levels respectively. The confidence with which these models are ruled out increases if a less conservative estimate of the baryonic mass is taken; increasing the LVG-modeled gas masses by a factor of two to account for neglected contributions to the total baryonic mass, drives up these sigma levels to ∼ 2σ and 2.7σ respectively. Furthermore, model (iv) can now be ruled out at the 1σ level as well. We are therefore led to the conclusion that HDF 850.1 is modeled best by a collection of unvirialized molecular clouds with distinct CO and C + layers, as in PDR models. The LVG calculations for this model yield a kinetic temperature of 100 K, a velocity gradient of 1.6 km s −1 pc −1 , a molecular hydrogen density of 10 3 cm −3 , and a column density of 10 19 cm −2 . The corresponding molecular gas mass obtained using this LVG approach is MH 2 ≈ 5.16×10 10 M⊙. For this preferred model we find that the CO-to-H2 luminosity to mass ratio is α = 1.2 (K km s −1 pc 2 ) −1 , close to the value found for ULIRGs in the local universe.