CO line emission at high-z

# CO line emission from galaxies in the Epoch of Reionization

## Abstract

We study the CO line luminosity (), the shape of the CO Spectral Line Energy Distribution (SLED), and the value of the CO-to- conversion factor in galaxies in the Epoch of Reionization (EoR). To this aim, we construct a model that simultaneously takes into account the radiative transfer and the clumpy structure of giant molecular clouds (GMCs) where the CO lines are excited. We then use it to post-process state-of-the-art zoomed, high resolution (), cosmological simulation of a main-sequence (, ) galaxy, “Althæa”, at . We find that the CO emission traces the inner molecular disk () of Althæa with the peak of the CO surface brightness co-located with that of the [C] 158 emission. Its is comparable to that observed in local galaxies with similar stellar mass. The high () gas surface density in Althæa, its large Mach number (), and the warm kinetic temperature () of GMCs yield a CO SLED peaked at the CO(7–6) transition, i.e. at relatively high-, and a CO-to- conversion factor lower than that of the Milky Way. The ALMA observing time required to detect (resolve) at 5 the CO(7–6) line from galaxies similar to Althæa is h ( h).

###### keywords:
ISM: clouds - infrared: ISM - galaxies: ISM - line: formation - galaxies: high-redshift
12

## 1 Introduction

Constraining the properties of the molecular gas in galaxies at the end () of the Epoch of Reionization (EoR) is a compelling step to understand the process of star formation in the first galaxies.

Molecular hydrogen (), the most abundant molecule in the Universe, lacks of a permanent dipole moment and its first quadrupole line has an excitation temperature ( K) significantly higher than the kinetic temperatures ( K) of giant molecular clouds (GMCs) (McKee & Ostriker, 2007). This is the reason why molecular gas in galaxies is very often traced through the detection of the rotational transitions the carbon monoxide, CO, the second most abundant molecule after . The first CO rotational transition is in fact characterised by K, with critical density (Solomon & Vanden Bout, 2005; Carilli & Walter, 2013), i.e. it is easily excited within GMCs.

Noticeably, the different excitation requirements of the various CO lines ( K, for upper state rotational quantum number ) can be exploited to constrain gas properties (e.g. density, temperature), and the gas heating mechanisms (e.g. FUV photons, X-ray photons, cosmic rays, shocks). This can be done through the analysis of the so-called CO Spectral Line Energy Distribution (CO SLED) - flux in each emission line as a function of (Kaufman et al., 1999; Meijerink et al., 2007; Obreschkow et al., 2009; Mashian et al., 2015; Rosenberg et al., 2015; Lu et al., 2017; Indriolo et al., 2017).

Searches for CO line emission at redshift have been mainly focused on the most luminous sources such as QSOs (e.g. Bertoldi et al., 2003; Maiolino et al., 2007; Wang et al., 2010; Walter et al., 2012; Combes et al., 2012; Venemans et al., 2012; Gallerani et al., 2014, 2017, for a recent review) or powerful sub-millimeter galaxies (e.g. Riechers et al., 2010a; Weiß et al., 2013; Aravena et al., 2016). On the contrary, little is known regarding the molecular gas content of high- normal star-forming galaxies – e.g. Lyman Alpha Emitters (LAEs) and/or Lyman Break Galaxies (LBGs) – which are more representative of the bulk of galaxy population at the end of EoR (e.g. Dayal et al., 2008; Dayal et al., 2009; Vallini et al., 2012). Only a handful of CO detections have been reported in (mostly lensed) LBGs at (Baker et al., 2004; Coppin et al., 2007; Riechers et al., 2010b; Livermore et al., 2012; Saintonge et al., 2013; Dessauges-Zavadsky et al., 2015, 2016; Ginolfi et al., 2017), while searches for low- CO rotational lines in LAEs at have resulted in non-detections (Wagg et al., 2009; Wagg & Kanekar, 2012).

The ALMA advent has opened new perspectives for the detection of CO rotational transitions from the EoR. It is more sensitive than the other submillimeter/millimeter facilities adopted in previous studies targeting LAEs/LGBs in the EoR and, more importantly, its receiver bands allow one to trace transitions with (high-, hereafter), from . This represents a key-point because the peak of the CO SLED in galaxies with high specific star formation rate (sSFR) is often associated with high- transitions (Mashian et al., 2015), and such galaxies are common at high- (Jiang et al., 2016). Those transitions, being more luminous than low- ones, may represent a viable option to target molecular gas at the end of EoR.

The advent of ALMA has also triggered the development of – mostly semi analytical – works that aim at modelling the CO emission signal from high-redshift () galaxies (e.g. Obreschkow et al., 2009; Lagos et al., 2012; Vallini et al., 2012; Muñoz & Furlanetto, 2013; da Cunha et al., 2013; Popping et al., 2016).

Obreschkow et al. (2009); Lagos et al. (2012) and Popping et al. (2016) combined semi analytical galaxy formation simulations with sub-grid models devised to convert the molecular mass into CO luminosities, and predict the high- evolution of the CO luminosity function (e.g. Walter et al., 2014; Decarli et al., 2016; Vallini et al., 2016). Those models account for e.g. the heating by the cosmic microwave background (CMB) (see also da Cunha et al., 2013), the Far-Ultraviolet (FUV) flux from starburst, and/or by X-rays produced by active galactic nuclei. Vallini et al. (2012) used the semi-analytical model for the molecular fraction by Krumholz et al. (2009) to infer the CO luminosity from a sample of LAEs extracted from cosmological simulations. However, the work by Obreschkow et al. (2009); Lagos et al. (2012); Vallini et al. (2012) lacked of a detailed description of the internal structure of GMCs. Muñoz & Furlanetto (2013) made a step forward in this direction. Using an analytical model, they explored the link between the GMC properties set by the turbulence, and the physics of CO emission in high- LBGs; they conclude that the CO signal could be very difficult to observe from . While Muñoz & Furlanetto (2013) catch some of the fundamental aspects of the internal structure and of the radiative transfer within the GMCs, this was at the expense of a full description of the galaxy formation process.

The goal of this work is to assess the feasibility of detecting CO lines from typical LBGs at high- by simultaneously capturing the full cosmological context of high-redshift galaxy formation, and the radiative transfer from the outer photodissociation regions (PDRs) (Hollenbach & Tielens, 1999) up to the fully molecular inner part of GMCs. To this aim, we construct a physically motivated model that simultaneously takes into account the radiative transfer and the clumpy structure of GMCs. We then use it to post-process state-of-the-art zoomed cosmological simulations (Pallottini et al., 2017b, a). This type of sub-grid approach has been already shown to be an optimal strategy to obtain predictions and insights on the luminosity of, e.g., [C] 157.7 line emission tracing the neutral diffuse gas and dense PDRs (Vallini et al., 2013, 2015; Pallottini et al., 2015) in galaxies at the end or EoR. It allows one, on the one hand, to properly treat the small physical scales ( pc) of clumps in GMCs, and the complex network of chemical and physical reactions in the PDR layer and in the fully molecular parts of GMCs. On the other hand, it benefits from the high-resolution hydrodynamical simulations by obtaining a proper description, down to scales of pc, of the ISM density, turbulence level, metal enrichment, and radiation field into which GMCs are embedded.

The paper is structured as follows: in Sec. 2 we describe how we implement the CO line emission calculation taking into account the internal density structure of GMCs. In Sec. 3 we validate the model against local observations. Finally, in Sec. 4, we apply the model to high-resolution cosmological simulations to compute the CO emission from normal star-forming galaxies at the end of EoR. We draw our conclusions in Sec. 5.

## 2 Model outline

The model strategy adopted in this work is summarised in Fig. 1, where we outline its modular structure. The first part (see Sec. 2.1) deals with the analytical description of the internal density structure of GMCs and its time evolution. The second part concerns the radiative transfer performed to compute the CO line emission once the density field is established (see Sec. 2.2). The sub-grid model is optimised for implementation in cosmological simulations able to approximately resolve GMC scales ( pc, see McKee & Ostriker, 2007).

### 2.1 The internal structure of GMCs

Fragmentation of gas clouds has been intensively studied in the past, and both numerical and analytical models (e.g Vazquez-Semadeni, 1994; Krumholz & McKee, 2005; Padoan & Nordlund, 2011; Hennebelle & Chabrier, 2011, 2013; Kim et al., 2003; Wada, 2008; Tasker & Tan, 2009; Federrath & Klessen, 2013) have shown that the density field of an isothermal, non-gravitating, turbulent gas of mean density is well described by a log-normal probability distribution function (PDF). The volume-weighted PDF () can be written as:

 Missing or unrecognized delimiter for \left (1)

where is the mean cloud density, and the volume-averaged value of the logarithm of the density is related to by (e.g Ostriker et al., 2001; Federrath & Klessen, 2013). The latter depends on the sonic Mach number () through the following relation:

 σ2=ln(1+b2M2), (2)

where parametrises the kinetic energy injection mechanism (often referred to as forcing) driving the turbulence (, see Molina et al., 2012). Throughout this paper we assume .

When self-gravity becomes important, the probability of finding dense regions increases, and a power-law tail develops on the high-density side of the PDF. The occurrence of the power-law tail is confirmed both theoretically (e.g. Krumholz & McKee, 2005; Hennebelle & Chabrier, 2011; Padoan & Nordlund, 2011; Federrath & Klessen, 2013), and observationally via molecular line detections (e.g. Goldsmith et al., 2008; Goodman et al., 2009; Schneider et al., 2016) or dust extinction measurements carried out in nearby GMCs (e.g. Kainulainen et al., 2009; Lombardi et al., 2015; Stutz & Kainulainen, 2015; Schneider et al., 2016). Treating properly the high-density tail of the density PDF is pivotal when computing the emission of high- CO rotational lines as they have high critical densities ( for ), and trace the densest regions of GMCs (e.g. Carilli & Walter, 2013, and references therein).

In our model we describe the time evolution of the density PDF of self-gravitating GMCs via the formalism developed by Girichidis et al. (2014). Assuming a pressure-free collapse, these authors provide a set of analytical equations to calculate the functional form of the PDF at any given time, , given the initial . According to their model, the high-density tail of the PDF quickly asymptotes to a power-law, consistently with observations. Details on the calculation of from the initial are given in Appendix A.

As the initial PDF, in this paper we take , where is defined in eq. 1. The density at which starts to deviate from the lognormal, , moves with time to lower densities. This is due to the fact that collapse at a given density can only manifest itself after approximately a free-fall time, (see eq. 17), implying that the collapse of lower density parcels takes a longer time. Girichidis et al. (2014) find that the tail is well defined above a density that evolves with time according to the following relation:

 ρtail(t)ρ0≈0.2(ttff(ρ0))−2.0. (3)

This the reason why (see the schematic in Fig. 1), the density PDF of a GMC is ultimately function of two parameters: (i) the mach number () – which affects the lognormal part of the distribution –, and (ii) the ratio () – that determines the point at which the PDF significantly deviates from the lognormal distribution. Note that the additional time-dependence of the density PDF due to FUV-photoevaporation (Gorti & Hollenbach, 2002; Vallini et al., 2017; Decataldo et al., 2017) is neglected in this work.

Let now consider a GMC of radius , volume , and mean density3 (), characterised by a fixed Mach number   and an evolutionary time . The normalisation of the volume-weighted density PDF must satisfy the following relation:

 Vtot=∫PV(ρ|M,t/tff)dρ, (4)

which allows us to associate to each density a typical length scale,

 ri=(Vtot∫ρi+δiρi−δiPV(ρ|M,t/tff)dρ)1/3, (5)

and column density . The CO emission per unit volume from each density element can be computed as follows:

 lCO,J(ρi)=1VtotεCO,J(ni,Ni,Z,G0)4πr2i, (6)

where is the CO emissivity of the transition as a function of the gas element (i) density, , (ii) column density, , (iii) metallicity, , and (iv) Far-Ultraviolet (FUV) flux, , in the Habing band () normalised to that in the solar neighbourhood ( Habing, 1968). The total CO emission4 from the GMC is then:

 LtotCO,J=∫lCO,J(ρ)PV(ρ|M,t/tff)dρ. (7)

It is useful to express eq. 7 in terms of :

 LtotCO,J=∫LCO,J(s)ds, (8)

where .

We calculate the value of for the first 9 rotational transitions of the CO molecule with the cloudy  version c13.03 (Ferland et al., 2013).

cloudy includes reactions involving molecules containing H, He, C, N, O, Si, S, and Cl atoms (see Appendix A of Abel et al., 2005, for details on the molecular network). Details on the CO network, including CH, CH, OH, OH, , , , and are presented in Ferland et al. (1994). The majority of reaction rates come from the UMIST 2000 database (Le Teuff et al., 2000). The treatment of the formation and dissociation of the molecule is outlined in Shaw et al. (2005) and it accounts for formation via gas-phase reactions, and on the dust grain surface (Cazaux & Tielens, 2004). The local grain properties (temperature, charge) at each point in the cloud are computed self-consistently. cloudy also treats the primary and secondary cosmic-ray (CR) ionisation processes. We adopt the default CLOUDY prescriptions for the CR ionization rate background (; CRIR, hereafter) (Indriolo et al., 2007). The secondary ionisation rate is (Glassgold & Langer, 1974). As the CRIR is a fundamental parameter that may have strong effects on the gas temperature and chemical composition at high densities, we discuss the impact of our assumption on the CO luminosity inferred with our model in the Appendix B.

We adopt a 1D geometry, assuming a gas slab of density , and fixed metallicity . The spectral energy distribution (SED) of the radiation field impinging on the slab surface is calculated using the stellar population synthesis code starburst99 (Leitherer et al., 1999), assuming a Salpeter Initial Mass Function in the range , and considering a continuous star formation mode with an age of the stellar population of 10 Myr. We adopt the Geneva standard evolutionary tracks (Schaller et al., 1992) with metallicity , and Lejeune-Schmutz stellar atmospheres which incorporate plane-parallel atmospheres and stars with strong winds (Lejeune et al., 1997; Schmutz et al., 1992). The SED is normalised so that the Habing flux varies in the range .

We adopt the gas-phase abundances (, , , )5 provided in cloudy  for the ISM of the Milky Way. The abundances are an average of those measured in the cold and warm diffuse phases of the galactic ISM by Cowie & Songaila (1986) and Savage & Sembach (1996).

We implement the Weingartner & Draine (2001) grain size distribution (; GSD) which has been shown to reproduce (among others) the SMC extinction curve. This study provides a set of analytical equations (see eqs. 4–6, and the best fit parameters in their Tab. 3) that (i) imposes a smooth cutoff for sizes greater than a threshold , (ii) controls the steepness of this cutoff, and (iii) allows for a change in the slope () below .

As we are interested in obtaining predictions on the CO luminosity in high- galaxies, we set the floor temperature of our fiducial simulations to that of CMB at , i.e. K. The line emission from molecular gas is affected in two ways by the CMB: (i) the higher leads to an increase of the line excitation, and thus of the line luminosities; (ii) the background against which the line is measured also increases (e.g Obreschkow et al., 2009; da Cunha et al., 2013). To test the behaviour of our model on local observations, and to study the impact of the CMB temperature on our predictions, we run also a set of cases using the present-day CMB temperature, (see Sec. 3.2 for a discussion on this point).

We run a total of cloudy  simulations varying (in 0.5 dex steps) in the range , in , and in , scaling the gas-phase abundances and the dust-to-gas ratio with the metallicity of each specific model.

Such parameter space brackets the plausible range of GMC parameters relevant to high- in galaxies. The code computes the radiative transfer through the slab up to a hydrogen column density . This stopping criterion is chosen to fully sample the molecular part of the illuminated slab, typically located at (McKee & Ostriker, 2007). The output of each run is the CO line emissivity of each transition:

 εCO,J=εCO,J(ni,Ni,Z,G0) (9)

which enters in eq. 6 and hence in eq. 7.

## 3 Model validation

By adopting the procedure described in the previous sections we compute the luminosity of the first 9 CO rotational transitions. The link between the time evolution of the density PDF, and the resulting CO luminosity is illustrated in Figure 2. In the upper panel, we plot volume-weighted density PDF at (lognormal, ) and the resulting after (solid line). We choose to maximise the mass ( of the total mass of the GMC, Girichidis et al., 2014) in the tail, hence highlighting more clearly the effect of the tail appearance on the CO emission. The fiducial cloud (see Tab. 1) is characterised by: (), , , and illuminated by a FUV field with . In the lower panel we plot the corresponding (see eq. 8) at (solid lines) and at (dashed lines). The CO(1–0) emission is boosted at high densities, once the tail has developed at .

### 3.1 The Mvir−L′CO relation

Studies of resolved molecular clouds find that GMCs are in approximate virial equilibrium, (e.g. Larson, 1981; Solomon et al., 1987; Bolatto et al., 2008) and, because of that, they obey scaling relations, ofter referred to as Larson laws. Those relations ultimately link the size (), the velocity dispersion (), and the CO luminosity ()6of GMCs (Solomon et al., 1987):

 σ≈0.7R0.5kms−1 L′CO≈130σ5 Kkms−1pc2 (10) L′CO≈25R2.5 Kkms−1pc2.

These relations are valid under the assumption that molecular gas is dominating the mass enclosed in the cloud radius and therefore the virial mass () is a good measure of the traced by CO.7

The virial equilibrium implies that , thus eqs. 3.1 translate into the following:

 Mvir≈39L′0.81COM⊙, (11)

where we note that the slope of the mass-luminosity relation, , is the one found by Solomon et al. (1987). Subsequent studies (e.g. Bolatto et al., 2008; Bolatto et al., 2013, and references therein) have confirmed that .

In Figure 3 we plot as resulting from our model. More precisely, we calculate the CO(1–0) luminosity fixing the mass (), radius (), considering the at , i.e. when of the GMC mass is in the PDF tail, and the rest in the lognormal distribution. The Mach number is selected so that it satisfies the following condition:

 Mcs=√GMGMCRGMC. (12)

We set the sound speed of the GMC to . For the RT we adopt the cloudy  runs at (i.e. those with ) for different metallicities . Our results are in nice agreement with observations, and a linear fit between and returns a slope for , respectively. The CO luminosity decreases with decreasing at fixed , implying that the CO-to- conversion factor (see Sec. 4.4) increases for lower , as noticed also by e.g. Wolfire et al. (2010); Glover & Mac Low (2011); Narayanan et al. (2012); Bolatto et al. (2013). The model results at enclose the observations by Pineda et al. (2009); Wong et al. (2011) of GMCs in the LMC ( Rolleston et al., 2002; Chevance et al., 2016; Lee et al., 2016).

### 3.2 The impact of the high-density tail

In this Section we discuss how the high-density power-law tail affects the luminosity of the various CO transitions. To do that we compare the pure log-normal density field for the fiducial GMC (()), and a distribution which includes the power-law tail after (()).

In Figure 4 we plot the ratio as a function of the upper rotational quantum number . We separately address the case and , and we assume two different values for the Mach number (). In all cases, and for all the rotational transitions, . More precisely, for the pure lognormal density distribution can account only for () of the CO emission model including the tail. If , (, ) and (, ). Both at and , we note a clear decreasing trend of with , highlighting the strong contribution of the dense tail gas to the high- lines that have increasingly high excitation temperatures and critical densities (Carilli & Walter, 2013). Hence, for a given density distribution, the emission from high- CO lines is boosted for warmer kinetic temperatures.

This is exactly what we obtain at fixed Mach numbers, where the drop of at high- transitions is steeper at than at . This is because the mean kinetic temperature of the molecular gas (i.e. for ) is lower at (), than at (). The increase of the Mach number from to boosts , i.e. the emission of the lognormal-distributed gas. The reason is that for large Mach numbers the lognormal distribution becomes wide (eq. 2) and hence a non-negligible fraction of the gas is compressed in high density, albeit not in gravitationally bound, structures. Such turbulent density enhancements are transient, and might not survive long enough time to allow formation of CO and molecules (see later in Sec. 4.3).

### 3.3 The CO SLED

Observations of the CO SLEDs can be used as a tool to constrain the properties of molecular gas and to link them to the star formation process, both from single GMCs (e.g. Pon et al., 2016; Lee et al., 2016; Indriolo et al., 2017) and on galactic scales (e.g. Rosenberg et al., 2015; Mashian et al., 2015; Lu et al., 2017; Pozzi et al., 2017). Any model that aims at predicting and interpreting the CO emission must be able to reproduce the relative strength of the various lines under different ISM conditions. In Fig. 5 we investigate how the resulting CO SLED from the fiducial cloud (Tab. 1) is influenced by variations of (i) the mean density, , (top left) (ii) FUV Habing flux, , (top right) and the density PDF shape parametrised by (iii) (bottom left) and  (bottom right). Given that we will compare our results with local observations of starburst galaxies in the local Universe (Mashian et al., 2015), in this Section we adopt cloudy runs with present-day CMB temperature. For the reasons explained later, we finally analyse the effects of metallicity variations separately.

#### Mean density effect

As already mentioned in Sec. 1 different CO rotational transitions trace gas with different properties. While low- rotational transitions () arise from diffuse (), cold () molecular ISM, higher- transitions are excited in denser () and warmer () gas (Kaufman et al., 1999). In Fig. 5 we plot the resulting CO SLEDs, normalised to the CO(1–0) transition, for different mean density ; the values of the other parameters are kept fixed to the fiducial ones. As expected, the peak of the CO SLED rise and shifts from with increasing .

This does not come as a surprise, as the population of high- CO levels – set by the competition of collisional excitation and radiative de-excitation – increases with increasing mean gas density. This eventually boosts the emission of high- CO lines, and shifts the peak of the CO SLED towards larger (see e.g. Weiß et al., 2007; da Cunha et al., 2013; Narayanan & Krumholz, 2014).

For high- CO lines, large gas temperatures are required to populate the corresponding rotational levels. An increase of the FUV fluxes at the GMC surface produces warmer gas and can ultimately boost the peak of the normalised CO SLED. Recall that CO(1–0) traces gas that is colder than that in which high- CO lines are excited. Note however that, in spite of a large (four dex) variation in , the CO SLED peak increases only by times (upper right panel of Fig. 5. As emphasised also by Kaufman et al. (1999), the temperature in the CII/CI/CO transition layer – directly influencing the CO emission – has a weak dependence on the Habing flux. A larger also forces such boundary to move deeper in the cloud.

#### Density PDF shape effects

In the lower left panel of Fig. 5 we plot the CO SLED dependence on the cloud evolutionary time . As discussed in Fig. 4, the emergence of the high-density tail in the more evolved stages boosts the emission of all CO lines, and more noticeably of the high- ones. This explain the shift-and-increase effect of the CO SLED peak at larger . A similar effect is also produced by large Mach numbers (lower right panel). As already explained a larger   causes a increase in the standard deviation of the of lognormal density distribution, thus allowing the gas to achieve larger densities with a non-negligible probability. This ultimately enhances the emission of high- CO lines.

For comparison, in Fig. 5 we plot in grey the range covered by the observed CO SLEDs (up to ) of five nearby starburst galaxies extracted from the sample of Mashian et al. (2015). We do not expect that any of our single cloud model can, alone, reproduce the observed CO SLED on global galactic scales. In fact, the galaxy-integrated CO SLEDs results from the overlap of many different GMCs which have different illumination, density, and temperature conditions.

Nevertheless the comparison of our simulated CO SLEDs with those observed allows us to highlight three points:

• Even though it is out of our scope to reproduce in detail the CO SLED shapes of the galaxies in the sample, the predictions for reproduce well the observed trends.

• The model under-predicts the luminosity of lines. This does not come as a surprise, as we have not included shocks in our treatment. Shocks in GMCs are known (Pon et al., 2012, 2015; Pon et al., 2016) to dissipate their energy primarily through CO rotational transitions. In particular, lines come often from shocked gas, and are typically brighter than those predicted by PDR models.

• Our model assumes a fixed CRIR and this might affect the shape of the CO SLED. For example, the CO SLED of NGC 253, in which the CRs flux is times the MW one (Bradford et al., 2003), is not reproduced by any combination of parameters shown in Fig 5. In the Appendix B we address in detail the effect of the CRIR variation on the CO line luminosity.

#### Metallicity effect

The analysis of the CO SLED presented in Fig. 5 has been performed at solar metallicity, as we were comparing our results with observations of local starburst galaxies whose . This is likely not the case of high- sources ( Pallottini et al., 2014), which are the primary target of this work. In what follows we will discuss the impact of variations on the CO SLEDs and we will use the observations of multiple CO emission lines in the N159W region of the Large Magellanic Cloud (LMC) (Lee et al., 2016) as benchmark.

N159W is one of the three prominent GMCs ( Fukui et al., 2015) in the N159W complex in the LMC ( kpc, Schaefer, 2008) whose physical conditions have been extensively studied at multiple wavelengths (Lee et al., 2016, and references therein). Recently, Lee et al. (2016) presented a coherent analysis of the CO and fine structure ([CII] 158, [OI] 145, [CI] 370) line emission to assess the properties of the molecular gas in N159W. N159W is a perfect target to test our model at low metallicities for two reasons: (i) the size8 of the observed region ( corresponding to ) is comparable to the diameter of our fiducial GMC (), and (ii) Lee et al. (2016) performed a PDR and shock analysis of N159W CO SLED up to the CO(12–13) transition with which we will compare our results.

In Fig. 6 we plot the CO SLED as a function of assuming two different FUV fluxes at the GMC surface (). We note that the luminosity of all CO lines decreases with . The effect of the variation (solid and dotted lines in Fig. 6) is negligible at , but it increases at lower metallicities. This is because, as pointed out by e.g. Chevance et al. (2016), a low metal (and dust) abundance results in less shielding. The FUV photons penetrate deeper into the cloud producing thicker PDRs and smaller CO cores. In Fig. 6 we show that there is a positive correlation between ,   and/or . This is in line with what discussed in the previous sections for the runs. In their PDR analysis of metal FIR lines, Lee et al. (2016) found that the best-fit is obtained for . Moreover, from the study of the CO-traced molecular gas they conclude that and . The PDR models used in their analysis, however, fail to explain the CO observations. Their conclusion is that the CO-emitting gas is excited by something other than UV photons, possibly shocks. Our model, including the internal density structure of the GMC is instead successful in matching the data. In fact, we find that models with and can reproduce the CO SLED up to if and .

## 4 Galaxy simulations

We apply the above CO-emission model to post-process a recently produced zoom-in simulation described in Pallottini et al. (2017a). Below we briefly summarise its main features.

Starting from cosmological initial conditions, we have used a modified version of the Adaptive Mesh Refinement code ramses (Teyssier, 2002) to carry out a zoom-in simulation of a dark matter (DM) halo of mass . In the zoomed-in region the gas has mass resolution of , and dynamics is followed down to spatial scales of . Stars are formed from molecular hydrogen, whose abundance is computed on the fly using the non-equilibrium chemistry code krome (Grassi et al., 2014) which is conveniently coupled to our customised version of ramses. The thermal and turbulent energy content of the gas is modelled following to e.g. Agertz & Kravtsov (2015). As detailed in Pallottini et al. (2017b), stellar feedback includes supernovae, winds from massive stars and radiation pressure. Stellar energy inputs and chemical yields depend both on time and stellar populations; the feedback prescription accounts for energy losses inside the GMC (albeit the density structure modelling introduced here has not yet been included).

The selected DM halo hosts “Althæa”, a galaxy characterised by a stellar mass , a mean gas surface density , and a at . Althæa features a SFR-stellar mass relation compatible with observations (e.g. Jiang et al., 2016), is in agreement with the Schmidt-Kennicutt relation (Krumholz et al., 2012), and has a [C] emission , slightly lower than the one expected from the local [C]-SFR relation (De Looze et al., 2014) and compatible with some high-z galaxy upper-limits (Schaerer et al., 2015).

### 4.1 Luminosities of individual simulated clouds

As outlined in Fig. 1, and discussed in Sec. 2, our CO model needs , , , , , and as inputs. In Fig. 7 we plot the mass weighted PDF of - relation in Althæa, considering a snapshot of the simulation. In Althæa the density PDF peaks at . In addition there is a small fraction of very dense gas () – the low   diagonal stripe. This part of the PDF describes the virtually metal-free gas in which production proceeds via gas-phase reactions rather than on dust grain surfaces and can survive only if self-shielded by a high density. The Mach number has a relatively wide distribution (see the inset of the Figure) with a pronounced peak at . This high level of turbulence is mostly supported by momentum injection associated with radiation pressure onto dust around massive stellar clusters and supernova explosions. It is a consequence of the very high star formation rate per unit area which in Althæa is about 1000 times higher than in the Milky Way.

To guide the following interpretation, in Fig. 8 we plot as a function of  and   the luminosity of the brightest transition (indicated by the white contours) from a single GMC characterised by and . Considering that and , we expect a typical from such a cloud, which has a mass of . The CO SLED peaks at CO(7–6) transition for all if , otherwise, if the peak shifts towards higher with decreasing metallicity.

### 4.2 CO emission from Althæa

Figure 9 shows the morphology of the CO(1–0) and CO(7–6) emission in Althæa. We select these two lines among all the CO rotational transitions because: (i) the CO(1–0) enters in the calculation of the CO-to- conversion factor, , which we will discuss in detail in Sec. 4.4) and, (ii) the CO(7–6) is the most luminous transition of the Althæa CO SLED (see Sec. 4.3) and it is observable with ALMA from .

As extensively discussed in Pallottini et al. (2017a), Althæa features a clearly defined, even though rather perturbed, spiral disk of radius kpc, embedded in a lower density () medium. The CO emission traces the disk where indeed most of the mass resides. Both the CO(1–0) and the CO(7–6) maps show an enhanced emission clump along the spiral arms. The peak of the CO(1–0) and CO(7–6) surface brightnesses are coincident (, and , respectively). Not surprisingly, they are co-located with a high-density () clump where also the 17.03 and [C] lines reach their maximum surface brightness (Pallottini et al., 2017a).

To understand how the total CO(1–0) luminosity of Althæa compares with local observations, in the upper panel of Fig. 10 we plot the SFR relation for galaxies in the COLD GASS (Saintonge et al., 2011), and ALLSMOG (Cicone et al., 2017) samples, as well as the location of Althæa in the same plane. The points are colour-coded in . We also highlight the evolution of the star-forming “main-sequence” (MS) from (Renzini & Peng, 2015) to (Speagle et al., 2014; Jiang et al., 2016).

Althæa nicely falls on the MS at (Speagle et al., 2014). Its SFR relation is in agreement with that recently found by Jiang et al. (2016) in LAEs and LBGs characterised by stellar ages . As shown in the lower panel of Fig. 10, the CO(1–0) luminosity of Althæa is (), and it is comparable to that of galaxies with the same stellar mass at .

The specific star formation rate () of Althæa is higher than that of MS-galaxies at , i.e. its SFR is larger than that of galaxies with comparable (see the color code of the Althæa symbol in the lower panel of Fig. 10). Therefore the CO(1–0) luminosity per unit SFR in Althæa is lower than that of MS galaxies.

### 4.3 The CO SLED in Althæa

The CO SLED is a unique tool to infer the properties of molecular gas, even if the simultaneous effect of the various parameters (e.g. the gas density, the FUV field, the gas metallicity) often makes its interpretation challenging. This is the reason why, in Sec. 3.3, we separately discussed the effect of each of the relevant parameters entering in our modelling. Here, we will refer to the results of that analysis to interpret the CO SLED of Althæa plotted in Fig. 11.

The CO SLED of Althæa includes the effects of CMB background radiation at . To isolate the effects of the CMB, we have also produced a case in which the CMB temperature is fixed at its present-day () value. The latter peaks at and is consistent with the range spanned by the observed CO SLEDs in a sample of local starburst galaxies (Mashian et al., 2015). The higher average GMC temperature due to the warmer CMB at , shifts the SLED peak at ; in addition, the galaxy is about 2 times more luminous in the brightest line available if located at compared to the case.

According to da Cunha et al. (2013) the observed flux of the -th line () against the CMB over the intrinsic one can be expressed as:

 SobsJSintrinsicJ=1−Bν[TCMB(z)]Bν[Tex], (13)

where is the rest frame frequency of the transition, and is the black body spectrum at temperature . Using this equation, the ratio of the observed versus intrinsic CO(7–6) flux from (substituting K and ) is .

Then, we can estimate the ratio of the observed CO(7–6) from over that at . This can be done by computing the intrinsic CO(7-6) luminosity in LTE (Obreschkow et al., 2009, see eq. 4,5), using () for (). This yields .

Qualitatively the CO(7–6) suppression due to the CMB (eq. 13) is compensated by the increased temperature of GMCs at high-. However the ratio as estimated in LTE () is slightly lower than what we find () from our model (Fig. 11). The origin of such discrepancy is due to a combination of the following factors: (i) cloudy calculations show that the excitation temperature inside clouds is not uniform but has a varying spatial profile; (ii) the predicted CO luminosity in Althæa arises from a collection of emitting clouds with different densities; (iii) the presence of an external stellar radiation field, , induces deviations from the LTE regime implicitly assumed by eq. 13.

In Fig. 11 the Althæa SLEDs are represented with shaded areas which highlight the variation of the CO line luminosity as a function of . We let  vary in the range , that causes an increase of by a factor . The impact of the tail is relatively small because the average Mach number of molecular cells in the simulation is (Fig. 7), which implies (see eq. 2). In this case the median density of the mass-weighted PDF is (see e.g. Tab. 1 in Girichidis et al., 2014). Given that the typical number density of molecular cells is , this means that 50% of the mass actually resides in structures with densities above even without including the the contribution of the tail. Given that for Althæa the  parameter has a minor impact on the CO luminosity, in what follows we fix .

However, we re-emphasise that, as discussed in Sec. 3.1, this statement holds true only in chemical equilibrium, i.e. if the timescale for (and CO) formation (, with Jura, 1975) is shorter than the typical lifetime of turbulent density enhancements. A simple estimate can be obtained as follows. Let us consider a clump with , corresponding to a typical scale (see eq. 5). At such density, the H formation timescale is (see Liszt, 2007; Glover et al., 2010, for a detailed calculation). Let us assume that the clump lifetime is the eddy turnover time (). Using Larson’s law, we know that the turbulent velocity scales as hence, . Substituting , and we obtain . This simple estimate shows that it can be difficult for (and CO) to form in a purely turbulent environment and gravitational collapse might be needed to keep the overdense regions bound. Finite lifetime of clumps and non-equilibrium chemistry is typically not considered in CO emission calculations that explore a wide range of physical conditions (e.g. Kazandjian et al., 2016), while it is accounted for in single cloud simulations (e.g. Glover et al., 2010; Shetty et al., 2011). This is an interesting aspect that is worth to be investigated in future work.

### 4.4 CO-to-H2 conversion factor

The so-called CO-to- conversion factor is defined as the ratio of the molecular gas mass to the CO(1–0) line luminosity:

 αCO≡(MH2L′CO)M⊙(Kkms−1pc2)−1. (14)

Observationally, is determined by combining independent measurements of gas mass with the detection of the CO(1–0) line. There are three methods to infer the molecular mass and/or column density: (1) assume that GMCs are in virial equilibrium, and derive the mass from the CO line width (e.g. Larson, 1981; Solomon et al., 1987), (2) assume a constant dust-to-gas ratio and use the dust continuum emission, possibly combined to HI measurements, to infer , e.g. Pineda et al. (2008, in resolved Galactic GMCs) and Leroy et al. (2011, for extragalactic studies); Magdis et al. (2011, for extragalactic studies); Sandstrom et al. (2013, for extragalactic studies); (3) through -ray emission induced by cosmic ray interactions with molecules (e.g. Padovani et al., 2009). In the MW disk the conversion factor is fairly constant, (Bolatto et al., 2013) on large (kpc) scales.

In recent years, a number of observational studies have provided evidence for at least two physical regimes where departs from the MW value. The first deviation is observed in sources with high-surface density and/or warm molecular gas such as mergers, and starburst galaxies where the conversion factor is lower than MW one (e.g. Yao et al., 2003; Tacconi et al., 2008; Papadopoulos et al., 2012). To first order high gas temperatures, and gas surface densities (), yield brighter CO emission at fixed molecular mass, thus decreasing the conversion factor (Narayanan et al., 2011, 2012).

A deviation in the opposite direction is instead observed in low-metallicity galaxies. In these sources is larger than in the MW (e.g. Bolatto et al., 2008; Leroy et al., 2011). This is due to the low C and O abundances and, most importantly, to the low dust-to-gas ratios which prevent an efficient shielding against FUV dissociating photons9. The reduced CO abundance then produces a fainter luminosity. The actual value of in high-redshift () galaxies is far from being firmly constrained. On the one hand, high- galaxies are more compact, dense, and star-bursting than the MW; all these facts lead to lower values. On the other hand, if their metallicity is sub-solar, one would expect higher ratios.

In Fig. 12 we plot the map of Althæa. The map is obtained by dividing the CO(1–0) surface brightness by the corresponding mass surface density. The mass of molecular hydrogen in each cell is:

 MH2,cell=fH2,GMCMgas,cell, (15)

where is the molecular fraction obtained with the sub-grid model (see Appendix C for details). The mean value (standard deviation) of the CO-to- in Althæa is . Even though  and have similar relative variations in the disk of Althæa, i.e. , most of the dispersion of is due to the density fluctuation. The is rather constant trough the disk and this is in agreement with what found by Sandstrom et al. (2013) that pointed out that the radial profile of in spiral galaxies is mostly flat. However, Sandstrom et al. (2013) found also a slight decrease of the conversion factor towards the center of the observed galaxies. In Althæa we do not see such a trend.

As pointed out previously, can be influenced by the gas surface density , by the strength of the Habing field, by the gas metallicity , and also by the CRIR.

Recently, Clark & Glover (2015); Glover & Clark (2016) performed detailed numerical simulations of turbulent molecular clouds that cover a wide range of metallicities, strength of the interstellar radiation field (ISRF), and cosmic ray ionisation rate, to investigate the impact of these parameters on the conversion factor. They find that increases with decreasing metallicities, in agreement with observations in the nearby Universe (e.g. Bolatto et al., 2008; Leroy et al., 2011). Increasing the ISRF and the CRIR produces different effects on the well-shielded clumps, where CO survives effectively, and in the diffuse interclump medium, where CO is instead dissociated effectively. Hence, even for high values of the ISRF, the integrated intensity from dense clumps increases owing to the heating of the gas. This is especially relevant in the case of model clouds characterised by high-density and high turbulent velocity dispersion, for which the conversion factor drops close to the MW value, even when the local ISRF is 100 times the fiducial (i.e. ) value.

The warm temperature of the molecular gas, sustained by the CMB at , and the high level of turbulence and gas surface densities of Althæa, are likely the two main reasons pushing the Althæa CO-to- conversion factor below that of the MW, despite the sub-solar metallicity, , and featured by Althæa. As a caveat, we also note that in our model we totally neglect any type of stellar feedback on GMCs (e.g. Gorti & Hollenbach, 2002; Vallini et al., 2017; Decataldo et al., 2017), which may affect the density field of GMCs, and ultimately increase the conversion factor.

Finally, we compare our inferred conversion factor with the one resulting from the best fit formula, , obtained by Narayanan et al. (2012). They simulate the hydrodynamic evolution of both isolated and merging disk galaxies finding that . If we substitute and we obtain in agreement within a factor with our results.

### 4.5 ALMA observability

In Sec. 4.3 we have shown that the peak of the Althæa CO SLED coincides with the CO(7–6) line. As CO transitions with fall in the ALMA bands from , in what follows, we use the ALMA Sensitivity Calculator to compute the observing time required to detect the CO(7–6) line from (a galaxy similar to) Althæa. The total CO(7–6) luminosity of Althæa is , i.e. of its [C] luminosity. If we assume the line width to be equal to that of [C] observations in LBGs at ( Maiolino et al., 2015; Pentericci et al., 2016; Knudsen et al., 2016; Bradač et al., 2017), this yields a peak flux density 10. The ALMA full-array observing time required to detect (resolve over ) the CO(7–6) line with a signal-to-noise ratio S/N=5 is h ( h), i.e. it would be challenging but doable within the maximum observing time () allowed for ALMA regular programs. This is in line with what found by Muñoz & Furlanetto (2013), which pointed out that ALMA observations of high- CO can be performed with reasonable integration times only in those galaxies which, like Althæa, are chemically evolved () and UV-bright ().

## 5 Conclusions

In this paper we have studied the CO emission properties of galaxies at the end of the Reionization Epoch. First, we have developed a semi-analytical model that, given the internal density field of a GMC and the impinging FUV flux () at the cloud surface, computes the CO line emission. The density PDF of the GMC is set both by turbulence (parametrised by the Mach number ) and self-gravity (parametrised by ). The radiative transfer is performed with cloudy and includes CMB radiation. The model takes the mean gas density , , ,   as inputs, and it returns the luminosity of the first 9 CO rotational transitions. We validated the model with local observations, demonstrating its capability to reproduce:

• the - relation

• the observed CO excitation in local SB galaxies

• the CO SLED of the low-metallicity () molecular cloud N159W in the LMC.

We used our validated CO emission model to post-process a cosmological zoom-in simulation of a prototypical galaxy (, ), Althæa. We have made the following additional assumptions: (i) we adopt a constant CRIR (cfr. Appendix B for a discussion on the effects of the CRIR variation); (ii) the dust-to-gas ratio scales linearly with metallicity; (ii) all clumps within a GMCs are exposed to the same (external) FUV flux predicted by the simulation; (iv) photoevaporation of clumps affecting the density field and the lifetime of GMCs has been neglected.

The key results can be summarised as follows:

• The CO emission traces the innermost disk of Althæa where most of the mass resides.

• The CO(1–0) luminosity of Althæa is comparable to that of main sequence galaxies with similar stellar mass at . As the MS evolves with redshift, and the sSFR increases, this means that the CO(1–0) per unit SFR in Althæa is lower than measured in local galaxies.

• The CO-to- conversion factor is , with little spatial variation throughout the disk. The dispersion is primarily introduced by density variations in the Althæa disk.

• The maximum of the CO(1–0) and CO(7–6) surface brightnesses are colocated in the disk (, and , respectively).

• The suppression of the observed CO luminosity due to the CMB at high- is compensated by the increased temperature of GMCs. The net result is both an increase of the CO SLED excitation, and a shift of the peak at higher .

• The peak of the Althæa CO SLED coincides with the CO(7–6) transition, and i.e. of the [C] luminosity. This is due to the relatively high surface density of Althæa and to the warm temperature () of the GMCs. To resolve the CO(7–6) line with a S/N=5 an ALMA observing time of is required.

## Acknowledgments

We thank G. Ucci, D. Cormier, R. Maiolino, M. Bothwell for useful discussions. We thank the anonymous referee for the valuable feedback that increased the clarity of the paper. AF acknowledges support from the ERC Advanced Grant INTERSTELLAR H2020/740120. ES acknowledges support from the Israeli Science Foundation under Grant No. 719/14.

## Appendix A Time evolution of the density PDF of a GMC

Following Girichidis et al. (2014), we first consider a homogeneous sphere, with initial density , collapsing under its own gravity. The density at any later time can be well approximated by

 ρ=ρ0[1−(ttff)2]−2, (16)

where

 tff≡√3π32Gρ0 (17)

is the free-fall time of the sphere.

Due to the conservation of the probability density, one can calculate the evolved density PDF as

 PV(ρ,t)=PV(ρ0,0)dρ0dρ, (18)

where needs to be expressed as a function of and inverting eq. (16).

## Appendix B Effect of varying the cosmic ray ionisation rate

In our simulations we assume (Indriolo et al., 2007). This value is x higher than the CRIR () generally adopted as default in many PDR calculations (Glover & Clark, 2012; Bisbas et al., 2015, e.g.). The actual CRIR value in Althæa and in general in high- galaxies is highly uncertain, owing to e.g. the unknown magnetic field in the ISM of high- galaxies.

A simple estimate of the CRIR can be obtained assuming a linear scaling with the SFR, as the main source of CRs is the Fermi acceleration in supernova (SNe) remnants, and the rate of SNe is related to the rate at which stars form. In Althæa the SFR, thus, a linear CRIR-SFR scaling would imply i.e. at most a factor of 5 greater than the value adopted in our cloudy simulations.

To test the impact of the CRIR variation on the CO (and CI) abundance, gas temperature, and CO line emissivity, we run a set of cloudy models with CRIR in the range , keeping and fixed to the fiducial values of Althaea. In Fig. 13, we show the variation of the and ratios (upper left panel), gas temperature (upper right) at , CO(1–0) emissivity (bottom left), CO(7–6) emissivity (bottom right), that we obtain in two different set of CLOUDY runs at fixed gas number density and . The and ratios are obtained by averaging the and returned by cloudy as a function of the depth (