The CO universe: Modelling CO emission and H 2 abundance in cosmological galaxy formation simulations
MMNRAS , 1–13 (2020) Preprint 31 August 2020 Compiled using MNRAS L A TEX style file v3.0
The CO universe: Modelling CO emission and H abundance incosmological galaxy formation simulations Shigeki Inoue , , , (cid:63) , Naoki Yoshida , , & Hidenobu Yajima Center for Computational Sciences, University of Tsukuba, Ten-nodai, 1-1-1 Tsukuba, Ibaraki 305-8577, Japan Chile Observatory, National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan Kavli Institute for the Physics and Mathematics of the Universe (WPI), UTIAS, The University of Tokyo, Chiba 277-8583, Japan Department of Physics, School of Science, The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan Research Center for the Early Universe, School of Science, The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We devise a physical model of formation and distribution of molecular gas clouds in galaxies. We use the model to predictthe intensities of rotational transition lines of carbon monoxide (CO) and the molecular hydrogen (H ) abundance. Using theoutputs of Illustris-TNG cosmological simulations, we populate molecular gas clouds of unresolved sizes in individual simulatedgalaxies, where the effect of the interstellar radiation field with dust attenuation is also taken into account. We then use thepublicly available code DESPOTIC to compute the CO line luminosities and H densities without assuming the CO-to-H conversion factor ( α CO ). Our method allows us to study the spatial and kinematic structures traced by CO(1-0) and highertransition lines. We compare the CO luminosities and H masses with recent observations of galaxies at low and high redshifts.Our model reproduces well the observed CO-luminosity function and the estimated H mass in the local Universe. About tenpercent of molecules in the Universe reside in dwarf galaxies with stellar masses lower than 10 M (cid:12) , but the galaxies aregenerally ‘CO-dark’ and have typically high α CO . Our model predicts generally lower CO line luminosities than observations atredshifts z (cid:38) Key words: methods: numerical – galaxies: evolution – ISM: molecules
Molecular gas clouds (MGCs) are the birthplaces of stars, and un-derstanding the formation and the abundance of MGCs is importantin the study of galaxy formation and evolution. Despite the ubiqui-tous existence, direct observation of H molecules in the Universe isseverely limited. Radiation from H molecules in the cold inter-stellarmedium (ISM) is not observable because of the high excitation ener-gies. Rotational transition lines of CO molecules are often used as aproxy for H in the ISM. CO is the second most abundant moleculesin the Universe, and a variety of emission lines can be observed inradio to submillimetre bands.The rotational transition line from the state of J = , and the conversion factor α CO between them is often assumed to be constant (see Section 4.1).However, it is also known that the value of α CO varies dependingon locations in a galaxy, galactic types and redshifts (e.g. Bolattoet al. 2013). Also there are many ‘CO-dark’ MGCs, where CO(1-0)emission is significantly weaker than expected for the estimated H density. The variation of α CO likely reflects differences in the localenvironments and physical properties among different MGCs, suchas their column densities, intensity of far-ultraviolet (FUV) radiation (cid:63) E-mail: [email protected] fields, dust and metal abundances (e.g. Narayanan et al. 2012; Liet al. 2018, see also Section 2).Theoretical studies have been hampered by difficulties associatedwith the thermal and chemical structure and evolution of MGCs.Since most H molecules are formed on the surface of dust grains,one needs to model the formation, growth and destruction of dustgrains that involve various physics and chemistry even at micro-scopic levels. In addition, as shown in radiative transfer simulationsof Glover et al. (2010), the H abundance can be large in centralregions of MGCs that are self-shielded against external radiationfields, which means that calculating molecular abundances needs de-tailed radiative transfer computations with resolving the small-scaleinternal structure. The conditions are simpler for CO molecules. For-mation of CO does not rely on dust grains, although it needs H (e.g.Krumholz 2017; Girichidis et al. 2020). The self-shielding effect ofCO is less important than H . In addition, metallicity gradients ingalaxies can enhance the inhomogeneity of the abundance ratios.These facts actually suggest that the abundance ratio between H andCO may not be uniform. Numerical simulations aimed at predictingaccurately the molecular abundances and line emissivities are gen-erally required to have extremely high resolutions and to implementmolecular chemistry, radiation transfer and dust models; Glover et al.(2010) and Grassi et al. (2014) study the convergence with respectto resolutions in their detailed radiative transfer simulations.Because formation of molecules does not alter the gravitational © 2020 The Authors a r X i v : . [ a s t r o - ph . GA ] A ug S. Inoue, N. Yoshida & H. Yajima assembly nor overall dynamical evolution of galaxies, often post-processing methods have been applied to outputs of simulations inorder to compute molecular abundances. A number of studies focuson a single galaxy with sufficiently high resolution to resolve giantmolecular clouds (e.g. Narayanan et al. 2012; Li et al. 2018; Valliniet al. 2018; Armillotta et al. 2020; Keating et al. 2020; Li et al.2020). These studies treat a small number of galaxies in isolated orin cosmological simulations and thus statistical quantities such as lineluminosity functions cannot be reliably determined. Also, possibleredshift evolution has not been fully addressed.Large-volume cosmological simulations are indispensable for sta-tistical studies of galaxy populations , but the typically poor mass res-olution achieved to date still hampers us from directly representingindividual MGCs. Also it is computationally expensive to performmulti-dimensional radiation transfer for a large number of galaxieseven if the mass and spatial resolutions were appropriate. Therefore,previous studies on galaxy populations utilise large-box simulationsby employing simple, empirical models or semi-analytic approaches(e.g. Obreschkow et al. 2009; Lagos et al. 2011; Popping et al. 2015;Schäbe et al. 2020; Davé et al. 2020). Popping et al. (2019) apply apost-processing method to the results of a cosmological simulationand of a semi-analytic model based on a fit to H fraction obtainedby Gnedin & Kravtsov (2011, also see Diemer et al. 2019). Theystudy in detail the evolution of H mass functions of galaxies fromredshift z = masses are in agreement withobservations of local CO-luminosity functions assuming a constant α CO . However, their results predict H masses that are significantlysmaller than high-redshift CO observations if they adopt the same α CO . As a possible reason, they argue that α CO may decrease withredshift on average in the observed galaxies.It seems that theoretical studies tend to focus on H as it is moredirectly linked to star formation (SF). Unfortunately, H is not read-ily observable, and thus modelling CO formation is necessary tomake direct comparison with observations. It is also important tostudy higher rotational transition. Although high-redshift observa-tions of CO lines are currently performed by Atacama Large Mil-limeter/submillimeter Array (ALMA), redshifted CO(1-0) lines fromdistant galaxies shift to outside the observable wavelength range ofALMA. The current high-redshift CO observations, therefore, relyon high- J transition lines, which are often converted to derive ex-pected CO(1-0) emission with assuming some certain spectral lineenergy distribution (see Section 4.2.1).In this paper, we propose a physical model of MGC formationand distribution, and compute H density and CO line luminosityin a consistent manner. We apply the post-processing method to theoutputs of a large-box cosmological simulation. Our model allowsus to calculate intensities of not only CO(1-0) but also high- J lines.We make direct comparison between the simulation and observationswithout converting CO luminosity to H density, and (dis)agreementbetween them could help us to understand differences of ISM andcloud properties among galaxies and their evolution with redshift.In the rest of the present paper, Section 2 describes the simula-tion we utilise, our modellings for gas clouds and related parameters.Section 3 gives our computations to obtain molecular abundancesand line intensities. Section 4 presents our results and comparisonwith recent survey observations for nearby and distant galaxies. Sec-tion 5 discusses agreement and disagreement between our resultsand observations. There, we address possible redshift-evolution ofmolecular clouds in the Universe. Section 6 summarises our findings. Our method is based on post-processing a large-volume cosmologicalsimulation. We aim at obtaining the large-scale distribution of galaxypopulations as well as resolving spatial and kinematic structures ofgas, stars and dark matter in individual galaxies. Unfortunately, evenstate-of-the-art cosmological simulations do not fully resolve thedistribution of stars and gas clouds at length scales of (cid:46)
100 pc.We thus resort to applying a physically motivated model to populatethe galaxies with MGCs. In Section 3, we compute the molecularabundance and atomic/molecular line emissivities using approximateradiation transfer with chemistry for individual MGCs. To this end,we also need to estimate strength of radiation affecting molecularabundances of the gas clouds, and we employ simple models forradiation fields and dust attenuation.
We utilise the data set of IllustrisTNG simulations (Nelson et al.2018a). The details of the simulations are presented on the Illus-trisTNG project web site and in related papers such as Nelson et al.(2018a), Weinberger et al. (2017) and Pillepich et al. (2018). Wespecifically use the outputs of TNG100-1 run. The simulation boxhas a comoving side length of 75 Mpc, and the mass-resolutions fordark matter and gas are 7 . . × M (cid:12) . In the simulation,dense gas cells with ρ cell > n H , SF = . − are converted to stel-lar particles according to a stochastic SF model. Therefore the stellarmass resolution is roughly the same as that of the parent gas cells.The star formation rate (SFR) is calculated as (cid:219) m star = f M m cell t SF (1)where m cell is a mass of the parent cell, the factor f M is the massfraction of cold gas (see Section 2.2), and t SF = ( G ρ cell ) − / . TheISM model of Yepes et al. (1997, see also Springel & Hernquist2003) is adopted to the star-forming gas with ρ cell > n H , SF . Type-II supernovae (SNe) are triggered immediately following the SF,and a mechanical feedback model of Springel & Hernquist (2003)is adopted to represent stellar feedback effects. Type-Ia SNe andasymptotic giant branch stars eject mass and metals into nearby gascells. The simulation also implements creation and feedback of blackholes (Springel et al. 2005) and magnetic fields (Pakmor et al. 2014).Gravitationally bound structures are identified with the friend-of-friend and SUBFIND grouping algorithms (e.g. Springel et al.2001). The total masses and line luminosities are computed for eachSUBFIND group (galaxy). We first need to determine the mass distribution of MGCs and theirphysical properties. We achieve this by using the quantities of in-dividual gas cells in the simulation as follows. In the IllustrisTNGsimulations, the ISM model of Yepes et al. (1997) is applied to gaswith ρ cell > n H , SF = . − . The dense gas is assumed to consistof cold and hot phases; the hot gas can cool and be converted tothe cold gas by thermal instability, and the cold gas triggers "unre-solved" SF and SNe. The SNe evaporate some amount of cold gasback to hot one, and the ejecta is turned into the hot phase. By con-sidering the pressure equilibrium, the two-phase ISM model yields , 1–13 (2020) he CO Universe Figure 1.
We plot the fractions of CNM as a function of the gas density.
Blue solid line: the mass fraction of CNM of a gas cell at z =
0, whichis computed by the two-phase ISM model using the same cooling functionas that of IllustrisTNG simulation.
Green dashed line: the volume fractionof CNM, which is computed by equation (2) assuming the density contrast φ = Red dotted line with the right ordinate: density enhancement factor ρ CNM / ρ cell = f M / f V . Note the logarithmic scales. a barotropic equation of state (EOS). Since this effective EOS is sig-nificantly harder than the isothermal EOS, the two-phase ISM modelin cosmological simulations tend to prevent galactic discs from frag-menting and forming giant clumps (Inoue & Yoshida 2019).The two-phase ISM model computes the mass fraction of the cold-phase gas for each gas cell, f M , as a function of density for a givencooling function. Typically f M (cid:38) . n H , SF and approaches asymptotically to f M (cid:39) φ ≡ ρ CNM / ρ WNM ∼ φ =
100 in our model. FollowingDiemer et al. (2019) and Popping et al. (2019), we assume that allthe gas in cells with ρ cell > n H , SF is neutral (atomic and molecular),whereas the gas in cells with ρ cell < n H , SF is fully or partly ionisedand contains no molecules. For the dense gas, the total cell massis given by m cell = m CNM + m WNM = ρ CNM V CNM + ρ WNM V WNM and volume V cell = m cell / ρ cell = V CNM + V WNM , where V WNM and V CNM are the volumes of WNM and CNM within the cell. We assumethat the hot and cold phases defined in the two-phase ISM modelcorrespond to the WNM and CNM in Wolfire et al., and compute thedensity and volume of the CNM. The CNM volume fraction is f V ≡ V CNM V cell = f M ( − f M ) φ + f M , (2)where f M = m CNM /( m CNM + m WNM ) . Fig. 1 shows f V as a functionof density. In spite of the dominance in mass (the blue solid line),CNM occupies a small volume at ρ cell (cid:46) cm − when φ = ρ CNM ρ cell = f M f V (3)which is shown by the red dotted line in Fig. 1. We assume thatmolecules are formed only in CNM.Next, we need to determine the physical size of a MGC. Molecules can be photo-dissociated by the inter-stellar radiation field (ISRF) butthe dissociating radiation cannot penetrate deep into the innermostregions if the gas column density is high enough for self- and/ordust-shielding against the ISRF. The degree of the shielding effectsis determined by the MGC size and density. Jeans length λ J is oftenused to approximate the characteristic size of a MGC. We assumethat a cloud in the CNM has a radius r cloud = λ J = (cid:115) γ P CNM G ρ , (4)where G is gravitational constant and γ = / P CNM = P cell . The mass and the columndensity of a single cloud are given by m cloud = ρ CNM V cloud = π ρ CNM r (5)and Σ cloud = m cloud π r = ρ CNM r cloud . (6)The above quantities of ρ CNM and Σ cloud are used as parameterswhen we compute molecular abundances and emissivities in Section3.1. Molecular line emission is powered by an internal or external ra-diation. To estimate the ISRF strength, we use a simple parametricmodel and calibrate it to match available observations. FUV radiationis most relevant for photo-dissociation of molecules, which is mainlyemitted from young massive stars. We assume that the unattenuatedISRF intensity, χ int , scales with the total SFR within a galaxy as χ int = χ (cid:12) (cid:219) M star (cid:12) yr − , (7)where χ (cid:12) is the ISRF intensity in the solar neighbourhood.To calculate the dust-attenuated ISRF, we estimate the amount ofdust within a radius characterised by the SFR distribution since FUVradiation is primarily contributed by local SF in the galaxy. The dustcolumn density is approximated to be Σ dust = f dust M metal ( < r SFR ) π r , (8)where r SFR is the three-dimensional radius within which half thetotal SFR of the galaxy is enclosed, M metal ( < r SFR ) is the total metalmass within r SFR , and we assume a constant dust-to-metal fractionto be f dust = . a = . µ m and s = . − ,the optical depth is approximately given by τ = Σ dust as . (9)Note that the dust attenuation we consider here is different fromdust-shielding within a MGC. Later in Section 3.1, we discuss thedust-shielding effect on molecular emission.In addition to the ISRF, there may exist external radiation such asthe cosmic background radiation. Its intensity χ ext is uniform but mayvary with redshift. We follow the analytic model of Puchwein et al.(2019) and set the intensity at wavelength of 1000 as log ( χ ext / χ (cid:12) ) (cid:39) MNRAS000
100 in our model. FollowingDiemer et al. (2019) and Popping et al. (2019), we assume that allthe gas in cells with ρ cell > n H , SF is neutral (atomic and molecular),whereas the gas in cells with ρ cell < n H , SF is fully or partly ionisedand contains no molecules. For the dense gas, the total cell massis given by m cell = m CNM + m WNM = ρ CNM V CNM + ρ WNM V WNM and volume V cell = m cell / ρ cell = V CNM + V WNM , where V WNM and V CNM are the volumes of WNM and CNM within the cell. We assumethat the hot and cold phases defined in the two-phase ISM modelcorrespond to the WNM and CNM in Wolfire et al., and compute thedensity and volume of the CNM. The CNM volume fraction is f V ≡ V CNM V cell = f M ( − f M ) φ + f M , (2)where f M = m CNM /( m CNM + m WNM ) . Fig. 1 shows f V as a functionof density. In spite of the dominance in mass (the blue solid line),CNM occupies a small volume at ρ cell (cid:46) cm − when φ = ρ CNM ρ cell = f M f V (3)which is shown by the red dotted line in Fig. 1. We assume thatmolecules are formed only in CNM.Next, we need to determine the physical size of a MGC. Molecules can be photo-dissociated by the inter-stellar radiation field (ISRF) butthe dissociating radiation cannot penetrate deep into the innermostregions if the gas column density is high enough for self- and/ordust-shielding against the ISRF. The degree of the shielding effectsis determined by the MGC size and density. Jeans length λ J is oftenused to approximate the characteristic size of a MGC. We assumethat a cloud in the CNM has a radius r cloud = λ J = (cid:115) γ P CNM G ρ , (4)where G is gravitational constant and γ = / P CNM = P cell . The mass and the columndensity of a single cloud are given by m cloud = ρ CNM V cloud = π ρ CNM r (5)and Σ cloud = m cloud π r = ρ CNM r cloud . (6)The above quantities of ρ CNM and Σ cloud are used as parameterswhen we compute molecular abundances and emissivities in Section3.1. Molecular line emission is powered by an internal or external ra-diation. To estimate the ISRF strength, we use a simple parametricmodel and calibrate it to match available observations. FUV radiationis most relevant for photo-dissociation of molecules, which is mainlyemitted from young massive stars. We assume that the unattenuatedISRF intensity, χ int , scales with the total SFR within a galaxy as χ int = χ (cid:12) (cid:219) M star (cid:12) yr − , (7)where χ (cid:12) is the ISRF intensity in the solar neighbourhood.To calculate the dust-attenuated ISRF, we estimate the amount ofdust within a radius characterised by the SFR distribution since FUVradiation is primarily contributed by local SF in the galaxy. The dustcolumn density is approximated to be Σ dust = f dust M metal ( < r SFR ) π r , (8)where r SFR is the three-dimensional radius within which half thetotal SFR of the galaxy is enclosed, M metal ( < r SFR ) is the total metalmass within r SFR , and we assume a constant dust-to-metal fractionto be f dust = . a = . µ m and s = . − ,the optical depth is approximately given by τ = Σ dust as . (9)Note that the dust attenuation we consider here is different fromdust-shielding within a MGC. Later in Section 3.1, we discuss thedust-shielding effect on molecular emission.In addition to the ISRF, there may exist external radiation such asthe cosmic background radiation. Its intensity χ ext is uniform but mayvary with redshift. We follow the analytic model of Puchwein et al.(2019) and set the intensity at wavelength of 1000 as log ( χ ext / χ (cid:12) ) (cid:39) MNRAS000 , 1–13 (2020)
S. Inoue, N. Yoshida & H. Yajima − . z =
0. The background intensity monotonically increases to − . z = Finally, the radiation intensity in a galaxy is modelled as χ cloud = χ int exp (− τ ) + χ ext , (10)and this is another parameter of the molecular computations in Sec-tion 3.1. Our model assumes a constant χ cloud to all the MGCs in asingle galaxy. This approximation may not represent accurately localvariations of molecular abundances in, for instance, spiral arms andinter-arm regions in a galaxy. We focus on the statistical quantitiessuch as line luminosity functions for populations of galaxies in thepresent paper. Detailed radiative transfer within clumpy galaxies willbe a subject of our future study. We use the radiation transfer code DESPOTIC (Derive the Energet-ics and SPectra of Optically Thick Interstellar Clouds, see Krumholz2014) that can compute the abundances of various chemical speciesand can predict atomic/molecular line emission from a gas cloud.Here we focus on H and CO molecules. DESPOTIC employs aspherical one-zone cloud model, and includes carbon chemistry net-work as well as various physical processes for cooling and heatingof cold ISM. For a specified parameter set that describes the physi-cal properties of a gas cloud, the code self-consistently calculates athermal and chemical equilibrium state. It then returns the gas anddust temperatures, species’ abundances and line emissivities.The basic parameters to be input to DESPOTIC are: volume andcolumn densities of a cloud and the ISRF intensity. For these, weuse the values of ρ CNM , Σ cloud and χ cloud derived in the previoussections. Non-thermal velocity dispersion σ of the cloud is calculatedwith the assumption of a marginally bound state by gravity, wherethe virial parameter α vir ≡ σ r cloud /( Gm cloud ) = σ = (cid:115) π G Σ ρ CNM . (11)Although considering non-thermal (turbulent) motions inDESPOTIC might appear inconsistent with the assumption madewhen calculating r cloud (Section 2.2), we adopt the above equationby noting that MGCs can be highly turbulent after gravitational con-traction. Federrath & Klessen (2012) demonstrate that efficient starformation is driven by compressive turbulence (see also Federrath2018). We also note that the effect of varying r cloud shall be discussedin Section 5.2. The ionization rate due to hard X-ray photons and cos-mic rays is set to ξ = − χ cloud s − per H nucleus. We assumethe dust abundance in a cloud to be proportional to its metallicity Z cloud . With the efficient mixing approximation within a gas cell,we consider Z cloud = Z cell but impose the minimum metallicity of10 − Z (cid:12) set by Population III stars (Kuhlen et al. 2012), where Z (cid:12) isthe solar metallicity. The dust abundance of the cloud is given as D cloud = D MW Z cloud Z (cid:12) , (12) The values of log ( χ ext / χ (cid:12) ) (cid:39) − . − . − . − . − . z =
1, 2, 3, 4 and 5, respectively.
Figure 2.
Normalised histograms of gas cells as functions of N cloud computedby equation (13) in the snapshots at redshifts z =
0, 1 . . where D MW is the Milky Way dust abundance. We adopt the to-tal abundance of [ C / H ] = × − Z cloud / Z (cid:12) , [ O / H ] = × − Z cloud / Z (cid:12) and [ M / H ] = × − Z cloud / Z (cid:12) for C, O andthe other refractory metals (M), respectively. These abundancesare consistent with the solar value for Z cloud = Z (cid:12) (Draine 2011).For dust grains, we set three cross sections per H nucleus: one forthermal radiation σ = × − cm − at 10 K, one for photo-electric heating σ PE = − cm − and yet another one for ISRF σ ISRF = × − cm − . The dust-gas coupling coefficient is set to α gd = . × − . The spectral index for dust thermal radiation is β = .
0. Finally, the cosmic microwave background temperature is T CMB = .
73 K at z = ρ CNM , Σ cloud , χ cloud , Z cloud and z . Since Σ cloud is a function of ρ CNM inthe two-phase ISM model used in Illustris-TNG and the cloud size r cloud = λ J / ρ CNM at a given z , we generate alook-up table of H fractions f H and line emissivities W CO as afunction of ρ CNM , χ cloud and Z cloud at each output epoch (redshift).The parameter space covers the values of the three quantities with 30,10 and 10 grids for ρ CNM , χ cloud and Z cloud in logarithmic spacing. Using the look-up table, we compute CO line intensities and H abundance for each gas cell according to the procedures described inSection 2.2. We do not consider inter-galactic components that arenot gravitationally bound to any galaxies, because such a warm/hot,diffuse gas is highly ionised and contains little molecules.The number of clouds in a gas cell is given by N cloud = V CNM V cloud = m CNM m cloud . (13)We allow this value to be less than 1, i.e. a computational cell coversonly a fraction of a MGC. Histograms in Fig. 2 illustrate the distri-bution of gas cells as functions of N cloud in the snapshots at z = . .
8. The majority of the cells have N cloud ∼ N cloud < Although these values can be directly read from the snapshot data ofIllustris-TNG, we assume the simple scaling with metallicity to reduce thenumber of input parameters. The actual abundances in the simulation do notsignificantly deviate from the scaling relations.MNRAS , 1–13 (2020) he CO Universe Figure 3.
The large-scale distribution and the structure of galaxies.
Top two panels: stellar mass (left) and CO(1-0) luminosity (right) distributions in the wholesimulation box of the TNG100-1 at z = Bottom set of panels: the five brightest galaxies in CO(1-0) luminosity in the simulation (from left to right for the firstto fifth brightest ones). From top to bottom, the panels show gas surface densities, velocity-integrated surface brightness temperatures in CO(1-0), H surfacedensities and local α CO , respectively. The galaxies are oriented in random directions. The horizontal bar on the bottom right corner in each panel indicates thephysical scale of 20 kpc. MNRAS000
Top two panels: stellar mass (left) and CO(1-0) luminosity (right) distributions in the wholesimulation box of the TNG100-1 at z = Bottom set of panels: the five brightest galaxies in CO(1-0) luminosity in the simulation (from left to right for the firstto fifth brightest ones). From top to bottom, the panels show gas surface densities, velocity-integrated surface brightness temperatures in CO(1-0), H surfacedensities and local α CO , respectively. The galaxies are oriented in random directions. The horizontal bar on the bottom right corner in each panel indicates thephysical scale of 20 kpc. MNRAS000 , 1–13 (2020) S. Inoue, N. Yoshida & H. Yajima
Figure 4.
Velocity-integrated brightness temperatures of CO(1-0) lines ofgalaxies as functions of stellar mass. The blue dots indicate our results usingTNG100-1 at z =
0. The red solid and dashed lines show the median andthe ± σ ranges of the xCOLD GASS sample with detection of the CO(1-0)lines. increases at z = .
8. By calculating a velocity-integrated line emis-sivity W in a manner described in Section 3.1, we obtain the lineintensity as I cloud = π r W . Then, the total intensity of a gas cellis simply I cell = I cloud N cloud . (14)The total H mass is m H , cell = f H m CNM , where f H is mass fractionof H with respect to all the components including hydrogen, heliumand heavy elements. We note that introducing N cloud in Equation(13) enables our model to be independent of the mass-resolution ofa simulation. We begin with testing our model prediction at z = T CO ( − ) , H surface density, and α CO = Σ H / T CO ( − ) . These CO-brightest galaxies have local values of α CO ∼ α CO in regions where CO(1-0) line is strong. Note,however, that our model does not incorporate the local variation ofISRF within a galaxy (see Section 2.3). The xCOLD GASS survey (Saintonge et al. 2017) has observed COemission of local galaxies and built a large sample that is unbiasedexcept for sampling with an equal frequency among stellar massbins between M star = and 10 . M (cid:12) . We compare our modelprediction with their observations. Fig. 4 shows galaxy-integratedbrightness temperatures of CO(1-0) lines as a function of stellar mass.The blue dots indicate results (galaxies) from our model, and the redlines delineate the median (solid) and ± σ deviations (dotted) for Figure 5.
LFs of CO(1-0) lines at z =
0. The blue solid line delineatesour result using TNG100-1. The filled circles with error bars indicate theobservations of xCOLD GASS. The yellow ones include all observed sample,but the red ones only include those with detections of CO(1-0). The yellowand red solid lines are their fittings with Schechter functions. The verticalgreen lines at log L (cid:48) CO ( − ) = . . the observed sample with CO(1-0) line detection. Clearly, our modelreproduces the correlation of L (cid:48) CO ( − ) with M star , in agreement withthe xCOLD GASS sample. The L (cid:48) CO ( − ) - M star relation is thought tocorrespond to the SF main sequence, i.e. correlation between galacticSFRs and M star . Galaxies with high SFRs are molecular-rich andtherefore bright in CO emission.Fig. 5 compares CO-luminosity functions (LFs) between ourmodel and xCOLD GASS. Again, the result is consistent with theobserved CO-LFs, with the model underpredicting slightly in therange of L (cid:48) CO ( − ) ∼ –10 K km s − pc . Below the com-pleteness limit due to the stellar mass cut of xCOLD GASS atlog L (cid:48) CO ( − ) (cid:46) .
5, the CO-LF predicted by our model decreases andsignificantly deviates from the extrapolation of the Schechter func-tions fitted to the observations. Also in the low-luminosity range, thedata points of xCOLD GASS indicate lower number density than ourmodel. This could be attributed to the stellar mass cut of xCOLDGASS, whereas our model does not impose such a lower limit ofstellar mass on our galaxy sampling.Fig. 6 shows galaxy-integrated line ratios of r ≡ L (cid:48) CO ( − ) / L (cid:48) CO ( − ) as a function of M star . All the simulated galaxiesare located in the range between log r (cid:39) − M star (cid:38) M (cid:12) . Interestingly, below the stellar masslimit of xCOLD GASS, i.e. in M star < M (cid:12) , our result predictsthat most of the low-mass galaxies have low ratios of log r (cid:39) − . α CO as a galaxy-integrated value: α CO ≡ M H / L (cid:48) CO ( − ) . We note that α CO in the bottom panels of Fig. 3 is MNRAS , 1–13 (2020) he CO Universe Figure 6.
The line luminosity ratios L (cid:48) CO ( − ) / L (cid:48) CO ( − ) against stellar mass.The colour indicates the number of galaxies in each bin normalised to thehighest value (see the colour bar on the right). The green filled circles are theobserved values of the xCOLD GASS sample. Figure 7.
Distribution of galaxy-integrated α CO and stellar mass in TNG100-1. The colour code is the same as in Fig. 6. To ensure the accuracy of α CO ,we here exclude the galaxies whose total masses of star-forming gas arelower than 10 M (cid:12) . The cyan symbol at ( log M star , log α CO ) = ( . , . ) corresponds to the values measured in the inner Galactic disc. The greenfilled circles are α CO estimated using a model of Accurso et al. (2017) for thexCOLD GASS sample. defined as the ratio of surface H mass density to CO(1-0) brightnesstemperature measured locally. Fig. 7 shows the distribution of α CO and galactic stellar masses. We find α CO is approximately constant atlog α CO (cid:39) M star ∼ –10 . M (cid:12) . However, α CO increases below M star ∼ M (cid:12) , and there is a population of massivegalaxies that have high α CO with M star (cid:38) M (cid:12) . These low- andhigh-mass galaxies with such high α CO correspond to dwarf andmassive elliptical galaxies, and their high α CO are because of theirlow metallicities and/or diffuse gas distribution (see Section 4.1.2). Figure 8.
The fractions of H (blue), H I (red) and H II (green) that reside ingalaxies above a given M star . The solid lines indicate our results. The dottedlines are the observational results of Fletcher et al. (2020) derived from thexCOLD GASS (blue dotted) and xGASS (red dotted) surveys, and the shadedregions show their observational 1 σ errors. In Fig. 7, the cyan star-shaped symbol indicates α CO = . M star = . × M (cid:12) (Licquia & Newman 2015).The green filled circles correspond to the xCOLD GASS samplewhere α CO are not observationally determined but estimated usingthe model of Accurso et al. (2017). These values of α CO appear tobe somewhat lower than the averaged values in our model althoughthese are within the range covered by our prediction.From our model, we can compute not only mass of H but alsothose of H I and H II . We assume that the star-forming gas is neutral. Insuch a gas cell, WNM is assumed to be fully atomic, and abundancesof H and H I in CNM are computed with DESPOTIC. We considerthe diffuse gas with ρ cell < n H , SF to form no molecules, and itsionised fraction is computed in the IllustrisTNG simulation. Noteagain that we do not take into account inter-galactic gas that is notbound to any galaxies. Using these quantities, we estimate the cosmicdensities of hydrogen in molecular, atomic and isonised states tobe Ω H = . × − h − , Ω H I = . × − h − and Ω H II = . × − h − , respectively. In observations, Fletcher et al. (2020)estimates Ω H = ( . ± . ) × − h − from xCOLD GASS and Ω H I = ( . + . − . )× − h − from the xGASS survey (Catinella et al.2018). Jones et al. (2018) determine Ω H I = ( . ± . )× − h − fromthe ALFALFA survey (Giovanelli et al. 2005). The values of Ω H and Ω H I in our model are thus consistent with these observationalmeasurements although Ω H is slightly above the error range of thexCOLD GASS observations. Fig. 8 shows fractions of the cosmichydrogen densities cumulated from galaxies with high M star . Thefractions of Ω H in our model is consistent with the observationsof Fletcher et al. (2020) within the error ranges although those of Ω H I appear to be somewhat higher than the observations in M star (cid:38) . M (cid:12) . As Fletcher et al. (2020) mention, our results show thatnearly ninety per cent of H in the Universe resides in galaxieswith M star > M (cid:12) , and dwarf galaxies do not host a significant The model of Accurso et al. (2017) gives α CO as a function of metallicityand offset from the SF main sequence. Although DESPOTIC also computes H II abundance in a MGC, an amountof H II is generally negligible. MNRAS000
The fractions of H (blue), H I (red) and H II (green) that reside ingalaxies above a given M star . The solid lines indicate our results. The dottedlines are the observational results of Fletcher et al. (2020) derived from thexCOLD GASS (blue dotted) and xGASS (red dotted) surveys, and the shadedregions show their observational 1 σ errors. In Fig. 7, the cyan star-shaped symbol indicates α CO = . M star = . × M (cid:12) (Licquia & Newman 2015).The green filled circles correspond to the xCOLD GASS samplewhere α CO are not observationally determined but estimated usingthe model of Accurso et al. (2017). These values of α CO appear tobe somewhat lower than the averaged values in our model althoughthese are within the range covered by our prediction.From our model, we can compute not only mass of H but alsothose of H I and H II . We assume that the star-forming gas is neutral. Insuch a gas cell, WNM is assumed to be fully atomic, and abundancesof H and H I in CNM are computed with DESPOTIC. We considerthe diffuse gas with ρ cell < n H , SF to form no molecules, and itsionised fraction is computed in the IllustrisTNG simulation. Noteagain that we do not take into account inter-galactic gas that is notbound to any galaxies. Using these quantities, we estimate the cosmicdensities of hydrogen in molecular, atomic and isonised states tobe Ω H = . × − h − , Ω H I = . × − h − and Ω H II = . × − h − , respectively. In observations, Fletcher et al. (2020)estimates Ω H = ( . ± . ) × − h − from xCOLD GASS and Ω H I = ( . + . − . )× − h − from the xGASS survey (Catinella et al.2018). Jones et al. (2018) determine Ω H I = ( . ± . )× − h − fromthe ALFALFA survey (Giovanelli et al. 2005). The values of Ω H and Ω H I in our model are thus consistent with these observationalmeasurements although Ω H is slightly above the error range of thexCOLD GASS observations. Fig. 8 shows fractions of the cosmichydrogen densities cumulated from galaxies with high M star . Thefractions of Ω H in our model is consistent with the observationsof Fletcher et al. (2020) within the error ranges although those of Ω H I appear to be somewhat higher than the observations in M star (cid:38) . M (cid:12) . As Fletcher et al. (2020) mention, our results show thatnearly ninety per cent of H in the Universe resides in galaxieswith M star > M (cid:12) , and dwarf galaxies do not host a significant The model of Accurso et al. (2017) gives α CO as a function of metallicityand offset from the SF main sequence. Although DESPOTIC also computes H II abundance in a MGC, an amountof H II is generally negligible. MNRAS000 , 1–13 (2020) S. Inoue, N. Yoshida & H. Yajima
Figure 9.
Same as Fig. 8 but for the fractions of hydrogen that resides ingalaxies below a given α CO . amount of molecules. However, it is worth reminding of the factthat the observed H mass is estimated from CO luminosity via α CO of Accurso et al. (2017); on the other hand, our model directlycomputes H mass. It is known that dwarf galaxies generally have lowSF efficiencies leading to low stellar mass to halo mass ratios (e.g.Behroozi et al. 2013a,b, 2019), and it is often attributed to intense gasoutflows by SNe due to their shallow potentials. From our result inFig. 8, we argue that the low H abundances in the low-mass galaxiescould be another cause of the low SF efficiencies of dwarfs.Fig. 9 shows the fractions of the cosmic hydrogen densitiessummed over galaxies with low α CO in our model. Nearly ninety percent of H gas in the Universe resides in galaxies with log α CO (cid:39) . .
5, and only ten per cent of H is formed in ‘CO-dark’ galaxies withlog α CO (cid:38) .
5. Although such CO-dark molecular clouds in thesegalaxies would be missed in observations, their total amount is ex-pected to be insignificant with respect to the total molecular massin the Universe. However, the CO-dark galaxies host nearly half thetotal amounts of H I and H II .As we show above, our model adopted to IllustrisTNG withour fiducial parameter settings can thus reproduce well the galaxy-integrated properties reported in the previous observational studies.We discuss parameter-dependence of our model in Section 5.1. We find that the galaxy-integrated α CO strongly correlates with theaveraged metallicities Z of star-forming gas and their column densi-ties Σ cloud , whereas the correlations with the other properties suchas ISRF χ cloud and gas fractions are less clear. This finding is con-sistent with the result of Narayanan et al. (2012), in which they haveproposed a fitting function obtained from their isolated and mergersimulations: α fitCO = . (cid:104) Z / Z (cid:12) (cid:105)(cid:104) Σ H (cid:105) . , (15)where (cid:104) Z / Z (cid:12) (cid:105) is mass-weighted mean of gaseous metallicity, and (cid:104) Σ H (cid:105) H is average of H column densities weighted by H mass over Figure 10.
Left panels: the averaged values of log α CO in our model. Rightpanels: the differences of our results from the fitting functions of Equations(15 and 16); if log ( α CO / α fitCO ) >
0, our model predicts higher α CO than the fit-ting functions. The ordinates (cid:104) Z / Z (cid:12) (cid:105) SF indicate mass-weighted metallicitiesamong star-forming gas cells within the galaxy. In the top panels, the abscissa (cid:104) Σ H (cid:105) H indicates average of H column densities weighted by molecularmass over all gas clouds in the galaxy. In the bottom panels, the abscissa (cid:104) W CO (cid:105) L CO indicates average of luminosity-weighted CO emissivities over allclouds in the galaxy. We here exclude galaxies whose masses of star-forminggas are lower than 10 M (cid:12) . all gas clouds in the units of M (cid:12) pc − . The top panels of Fig. 10 showdistribution of the ensemble averages of log α CO in our model andcomparison with the fitting function of equation (15). In the top leftpanel, the mean α CO increases with decreasing Z and Σ H . Especially,most of the metal-poor galaxies with log (cid:104) Z / Z (cid:12) (cid:105) (cid:46) − . α CO , and galaxies with (cid:104) Σ H (cid:105) H (cid:46)
10 M (cid:12) pc − have high α CO irrespective of their metallicities. In the top right panel, we showlogarithmic differences of α CO between our model and equation (15). The definition is (cid:104) Σ H (cid:105) H ≡ (cid:205) i m H , i Σ H , i / (cid:205) i m H , i , where m H , i isH mass in i -th gas cell, and Σ H , i is H column density of a gas cloud in thecell: Σ H = Σ cloud f H .MNRAS , 1–13 (2020) he CO Universe The values of α CO are consistent between our model and the fittingfunction for the galaxies with (cid:104) Σ H (cid:105) H (cid:38)
10 M (cid:12) pc − within a factorof unity. The fitting function, however, predicts significantly lower α CO than our model for the galaxies with (cid:104) Σ H (cid:105) H (cid:46)
10 M (cid:12) pc − .This may be because the fitting function is derived from a number oftheir isolated and merger simulations for relatively massive galaxieswith M baryon ∼ M (cid:12) .The fitting function of Equation (15) can be converted to an ‘ob-servable form’ (Narayanan et al. 2012), α fitCO = . (cid:104) Z / Z (cid:12) (cid:105) . (cid:104) W CO (cid:105) . CO , (16)where (cid:104) W CO (cid:105) L CO is the luminosity-weighted mean of CO(1-0) emis-sivities over all clouds in a galaxy in units of K km s − . The bottompanels of Fig. 10 show the same as the top ones but with (cid:104) W CO (cid:105) L CO forthe abscissas. The values of α CO increase with decreasing (cid:104) W CO (cid:105) L CO in our model, and equation (16) gives significantly lower α CO forgalaxies with low (cid:104) W CO (cid:105) L CO (cid:46) − . Thus, the fitting func-tions derived by Narayanan et al. (2012) appear to be accurate anduseful for CO-bright galaxies that have dense molecular clouds with (cid:104) W CO (cid:105) L CO (cid:38) − and (cid:104) Σ H (cid:105) H (cid:38)
10 M (cid:12) pc − .Our model employs a simple approximation for calculating ISRF χ cloud and dust opacity τ for gas clouds (Section 2.3). Narayananet al. (2012) use dust radiation transfer calculations (Narayanan et al.2011) for spatially resolved distribution of radiation sources suchas stars and active galactic nuclei with spectrum energy distribu-tion models. It is noteworthy that, despite the simplicity, our modelreproduces their fitting results of α fitCO for a large number of galax-ies in cosmological simulations except for the diffuse galaxies with (cid:104) Σ H (cid:105) H (cid:46)
10 M (cid:12) pc − . We post-process snapshots of TNG100-1 at high redshifts. From theagreement of our results with the xCOLD GASS survey in Section4.1.1, we can possibly apply our model to high-redshift galaxieswithout calibrating the parameters or altering our model. This may,however, be a naive expectation. For example, dust properties such astypical size a and dust-to-metal fraction f dust (Sections 2.3 and 3.1)can evolve with redshift. The inaccuracy of our modelling for ISRF χ cloud can be amplified since the ISRF can be more intense due tothe higher SFRs and lower dust abundances (lower metallicities) ofhigh-redshift galaxies. In addition, galactic morphologies are moreirregular and complex due to intense gas accretion and frequentmergers at higher redshifts. Our modelling such as assuming r cloud = λ J / z = . . We compare our results from the ASPECS surveys. Decarli et al.(2016) provide the results of their pilot observation at wavelengthsof 1 and 3 mm (band 3 and 6) for galaxies at the mean redshifts (cid:104) z (cid:105) = .
5, 1 .
0, 1 .
4, 2 . .
8. Decarli et al. (2019) present datafrom large-programme (LP) observations at 3 mm for galaxies at (cid:104) z (cid:105) = .
3, 1 .
4, 2 . .
8, which covers a wider area than the pilot survey. Redshifted CO(1-0) emission (rest-frame 2 . J lines, CO[ J -( J − J =
3, 4, 2, 3 and 4 at z = .
5, 1 .
0, 1 .
4, 2 . .
8, respectively; CO(5-4) is also observedat z = . J lines are directly calculated with DESPOTICwithout introducing any conversion factors. The effective surveyvolume of ASPECS is much smaller than the simulation box ofTNG100-1. We are thus able to estimate the cosmic variance bysampling cubic regions with side lengths of 4 .
59, 5 .
18, 8 .
30, 6 . .
48 Mpc (physical) at z = .
5, 1 .
0, 1 .
4, 2 . .
8, respec-tively. These sampling volumes are consistent with those of ASPECSpilot observation for z = .
5, 1 . z = .
4, 2 . .
8. A total of 16384 cubic regions are randomly selected in thesimulation box.Fig. 11 compares high- J CO-LFs. The blue solid lines indicateour results including all the galaxies in the simulation. Our modelappears to underpredict the CO-LFs for CO(5-4) at z = .
4, CO(3-2)at z = . z = .
8. The model CO-LFs are outsidethe ranges of uncertainties of the ASPECS pilot and LP data.For high-redshift galaxies, CO(1-0) line luminosity is estimated from the higher- J lines using the line ratios r J for an assumedline energy distribution: L (cid:48) CO [ J −( J − )] = r J L (cid:48) CO ( − ) . The ASPECSsurveys use r J estimated by Daddi et al. (2015): r J = .
76, 0 . .
31 and 0 .
23 for J =
2, 3, 4 and 5, respectively. Our model providesCO(1-0) line luminosities directly. Fig. 12 compares our model CO(1-0) LFs with those converted from the ASPECS data. We also plotthe observational results of Orellana-González et al. (2020) whichinclude brighter galaxies than ASPECS. Our model matches theobserved CO-LFs within the error ranges at z (cid:46) .
0. At z = . .
8, the modelled galaxies appear to lack galaxies brighter than L (cid:48) CO ( − ) = K km s − pc , which is clearly inconsistent with theASPECS surveys.The ASPECS observations assume α CO = . (cid:12) ( K km s − pc ) − for all observed galaxies. Accord-ingly, their galaxy-integrated H mass is estimated as M H = α CO r J L (cid:48) CO [ J −( J − )] . (17)Note the estimate is subject to uncertainties of the two conversionfactors: α CO and r J , whereas our model can calculate M H inde-pendently from CO line luminosity. Fig. 13 compares the modeland observed (estimated) H mass functions. The model agrees rea-sonably well with ASPECS and Orellana-González et al. (2020) at z = .
5, 1 . .
4. However, our model tends to overpredict M H at z >
2. This trend is actually opposite to what is found for COLFs. Overall, at high redshift, our model predicts lower L (cid:48) CO andlarger M H than ASPECS, and the discrepancy appears to be largerat higher redshift. This indicates that the average α CO in our modelincreases with redshift. In Section 5.2, we discuss possible causes ofthe variation of α CO and further investigate the inconsistency withthe observations. Orellana-González et al. (2020) estimate the CO-LFs for their compilationof various observations using correlations between the luminosities of radiocontinuum, infrared and CO emission. Orellana-González et al. (2020) also assume α CO = .000
2. This trend is actually opposite to what is found for COLFs. Overall, at high redshift, our model predicts lower L (cid:48) CO andlarger M H than ASPECS, and the discrepancy appears to be largerat higher redshift. This indicates that the average α CO in our modelincreases with redshift. In Section 5.2, we discuss possible causes ofthe variation of α CO and further investigate the inconsistency withthe observations. Orellana-González et al. (2020) estimate the CO-LFs for their compilationof various observations using correlations between the luminosities of radiocontinuum, infrared and CO emission. Orellana-González et al. (2020) also assume α CO = .000 , 1–13 (2020) S. Inoue, N. Yoshida & H. Yajima
Figure 11.
Luminosity functions of CO[ J -( J − J is consistent with that of the ASPECS surveys in each panel. The blue solid linesindicate our results including all galaxies in the whole simulation box of TNG100-1. The shaded regions with thick and thin yellow are 1 σ and 2 σ confidenceintervals of the cosmic variance computed in the volumes whose sizes are the same as the observing regions of ASPECS. The magenta dashed lines are the sameas the blue solid lines but assuming the cloud sizes to be ten times larger: r cloud = λ J (see Section 5.2). The red and green boxes indicate the observationaldeterminations of the ASPECS pilot and LP (see Decarli et al. 2016, 2019), respectively. The vertical sizes of the red boxes show the Poisson errors of the pilotobservation. The thick and thin green regions correspond to 1 σ and 2 σ confidence intervals of the LP observations. The horizontal bars with downwards arrowsindicate upper limits of the observations. Figure 12.
Same as Fig. 11 but for CO(1-0) emission. In the three panels from left, the black filled circles with error bars indicate the observational results ofOrellana-González et al. (2020) at z = .
5, 1 . .
5, respectively.
Figure 13.
Same as Figs. 11 and 12 but for H masses. The observations plotted here are the same as those in Fig. 12 but shifted horizontally by α CO = .
6. Theresults with r cloud = λ J are not shown since these are hardly different from the fiducial case. z = α CO in our model is close to those estimated for the observed galaxies(see also Fig. 7).Since our post-processing method employs several controlling pa-rameters described in Section 2, we discuss the dependence of ourresults on the parameters. In Section 2.2, we assume the density ratio MNRAS , 1–13 (2020) he CO Universe based on the two-phase ISM model with φ ≡ ρ CNM / ρ WNM = < f M / f V (cid:46)
20. Ifwe set f M = f V =
1, by effectively deviating from the ISM model, ourresults of L (cid:48) CO and M H do not systematically change. Only the scat-ter becomes large for low-mass galaxies. This is because the mostimportant parameter for molecular abundance in a cloud is Σ cloud rather than ρ CNM (e.g. Li et al. 2018). We approximate a cloud sizeto be r cloud = λ J /
2, and λ J ∝ ρ − under the pressure equilibriumbetween the CNM and WNM. Therefore, Σ cloud ∝ ρ CNM r cloud is in-dependent of the density enhancement factor, and the resulting L (cid:48) CO and M H do not change significantly. In fact, we find that increasing Σ cloud can make L (cid:48) CO considerably larger (see Section 5.2).We employ a simple model for the unattenuated ISRF χ int thatis proportional to the total SFR in a galaxy (equation 7). To testthe effect of the ISRF intensity, we alter χ int to be proportional toSFR ×( r s / r SFR ) where r s = L (cid:48) CO and M H in low-massgalaxies, and that massive galaxies are little affected. We also notethat the cosmic background radiation χ ext is generally weaker than χ int , and thus does not significantly impact our results.The dust optical depth τ is another important parameter. In equa-tion (9), the optical depth τ is proportional to an uncertain factor f dust /( as ) that involves dust-to-metal fraction, grain size and density.Decreasing f dust /( as ) can effectively lower L (cid:48) CO because of strongerISFR. If we assume a lower f dust /( as ) than the fiducial value, ourCO-LF shown in Fig. 5 becomes inconsistent with the observations.In the extreme case of τ = ∞ , our results at z = The apparent discrepancy between our model and the ASPECS ob-servations at high redshift (Section 4.2.1) may suggest evolution ofISM structure in galaxies. At z (cid:38) . L (cid:48) CO ( − ) and higher M H than the observational estimates. Com-pared with ASPECS, the luminosities of high- J lines L (cid:48) CO [ J −( J − )] are systematically lower at high redshift. The trend of producinglarge H masses appears to be opposite to the result of Popping et al.(2019), who find significantly lower M H than the ASPECS at highredshift, although their model reproduces the H mass at low-redshiftobservations. This may suggest ISM structure evolution over redshift.We identify each galaxy as a gravitational bound object, whereasradio observations capture a galaxy’s CO emission within a beamsize. This procedure may affect the total amount of L (cid:48) CO ( − ) and M H . Using the IllustrisTNG simulation, Popping et al. (2019) ex-amine the systematic ‘aperture bias’. The resulting H mass functionsare not significantly different at high redshift between the two cases ofusing the SUBFIND grouping and the 3 . τ in our model does not signif- icantly affect our CO-LFs and can only mildly reduce the difference.Even in the case of τ = ∞ , our model predicts lower L (cid:48) CO than theobservations at high redshifts. The H mass functions are hardlyaffected.It has turned out that the cloud radius is the most important quantitythat affects the CO line emission strengths. In Figs. 11 and 12, we plotCO-LFs with assuming ten times larger cloud radii as r cloud = λ J (the magenta dashed lines). These test results are in better agree-ment with the ASPECS observations. However, L (cid:48) CO ( − ) at z = . . r cloud = λ J . Overall, we find that assuming en-larged r cloud overpredicts the CO-LF at low redshift. Enlarging r cloud hardly changes the H mass functions. In determining λ J (equation4), our model does not take into account contributions by turbulentor magnetic pressure. Although the turbulent velocities in localgalaxies are typically lower than the ISM sound velocity, it is knownthat star-forming galaxies at high redshift are highly turbulent withvelocity dispersions of ∼
100 km s − (e.g. Förster Schreiber et al.2009). Zhou et al. (2017) show illustrative comparison of velocitydispersions between star-forming galaxies observed at low and highredshifts. Ignoring the turbulent pressure can thus underestimate λ J for high-redshift galaxies. Intriguingly, large turbulent velocities arethought to produce large, massive gas clouds (e.g. Dekel et al. 2009;Dessauges-Zavadsky et al. 2019). A possible prescription to matchour results to the observations may be to introduce a fudge factor toincrease r cloud depending on redshift.It may be possible that state-of-the-art cosmological hydrodynam-ics simulations are still unable to reproduce gas properties at highredshifts. The IllustrisTNG simulation reproduce a variety of statis-tical properties of observed galaxies such as stellar mass function,mass-metallicity relation and redshift-evolution of the cosmic SFRdensity (e.g. Vogelsberger et al. 2013; Genel et al. 2014; Nelson et al.2018b; Naiman et al. 2018; Pillepich et al. 2018). However, it hasnot been fully examined whether the properties of the local ISM arecompatible with real galaxies. Inoue & Yoshida (2019) have demon-strated that clumpiness of high-redshift disc galaxies can stronglydepend on the ISM model (EOS of star-forming gas) assumed insimulations while global properties such as stellar and gas massesand SFRs are almost unchanged. Davé et al. (2020) analyse by post-processing various cosmological simulations with the same methodto calculate L (cid:48) CO ( − ) and M H . They find that the results do not con-verge between different sets of simulations. Hayward et al. (2020)post-process the IllustrisTNG and original Illustris simulations toreproduce number counts of submillimetre galaxies (SMGs) at theredshift z =
2. They find that utilising IllustrisTNG significantly un-derpredicts the number of bright SMGs although their result usingthe original Illustris simulation is consistent with observations. Theyargue that this is because galaxies of IllustrisTNG have lower dustmasses (metallicities) and SFRs than those of the original Illustris.Bright SMGs generally have quite high SFRs in observations, andthis fact implies that such SMGs are expected molecular-rich. Ourresults predicting the low CO luminosities may stem from the samereason. If the computation of λ J includes the turbulent and magnetic pressure, γ P CNM in equation (4) is replaced with γ P CNM + ρ CNM ( v + v ) , where v turb and v A are turbulent and Alfvén velocities of the CNM (e.g. Federrath& Klessen 2012). MNRAS , 1–13 (2020) S. Inoue, N. Yoshida & H. Yajima
We utilise the IllustrisTNG simulation and populate unresolved gasclouds whose sizes are approximated as ‘thermal’ Jeans lengths ofCNM. Adopting DESPOTIC, our method can compute not onlyCO(1-0) but also higher- J lines, calculate H mass.For galaxies at the redshift z =
0, we can reproduce the LF of CO(1-0) obtained by xCOLD GASS. Our values of α CO are consistent withthose estimated for the observed galaxies. Although we assume thesimple model to approximate the radiation fields, the distribution ofmodel α CO is in agreement with the results of the more detailedmodel of Narayanan et al. (2012). We find that about ten per cent ofH in the Universe resides in galaxies with M star (cid:46) M (cid:12) . Thesedwarfs have significantly low molecular abundances, which may bea reason for their low SF efficiencies.For high-redshift galaxies, our method underpredicts CO lumi-nosities. We find that we can mitigate the discrepancy of the CO-LFsif we enlarge the cloud sizes by a factor of ten, which corresponds toassuming a larger λ J than the ‘thermal’ Jeans lengths in high-redshiftgalaxies. Highly turbulent states of galaxies observed in high-redshiftUniverse is expected to suppress gravitational collapse in small scalesand lead to the formation of large and massive clouds. Thus, our re-sults imply the redshift-evolution of ISM properties in molecular-richand star-forming galaxies.Our method enables direct comparison between simulations andobservations. It forms the basis for a wealth of future studies, includ-ing mock observations using simulations and evaluation for potentialbiases in kinematic analyses using gas densities and velocity disper-sion estimated from measurements of CO intensities and line widths. ACKNOWLEDGEMENTS
This study was supported by World Premier International ResearchCenter Initiative (WPI), NAOJ ALMA Scientific Research GrantNumber 2019-11A, MEXT, Japan and by SPPEXA through JSTCREST JPMHCR1414. SI receives the funding from KAKENHIGrant-in-Aid for Young Scientists (B), No. 17K17677, and HY re-ceives the funding from Grant-in-Aid for Scientific Research (No.17H04827 and 20H04724) from the Japan Society for the Promotionof Science (JSPS). The numerical computations presented in this pa-per were carried out on the analysis servers, the general-purpose PCcluster and Cray XC50 at Center for Computational Astrophysics,National Astronomical Observatory of Japan.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Accurso G., et al., 2017, MNRAS, 470, 4750Armillotta L., Krumholz M. R., Di Teodoro E. M., 2020, MNRAS, 493, 5273Behroozi P. S., Wechsler R. H., Conroy C., 2013a, ApJ, 762, L31Behroozi P. S., Wechsler R. H., Conroy C., 2013b, ApJ, 770, 57Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488,3143Bolatto A. D., Wolfire M., Leroy A. K., 2013, ARA&A, 51, 207Catinella B., et al., 2018, MNRAS, 476, 875Daddi E., et al., 2015, A&A, 577, A46 Davé R., Crain R. A., Stevens A. R. H., Narayanan D., Saintonge A., CatinellaB., Cortese L., 2020, arXiv e-prints, p. arXiv:2002.07226Decarli R., et al., 2016, ApJ, 833, 69Decarli R., et al., 2019, ApJ, 882, 138Dekel A., Sari R., Ceverino D., 2009, ApJ, 703, 785Dessauges-Zavadsky M., et al., 2019, Nature Astronomy, 3, 1115Diemer B., et al., 2019, MNRAS, 487, 1529Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium.Princeton University PressFederrath C., 2018, Physics Today, 71, 38Federrath C., Klessen R. S., 2012, ApJ, 761, 156Fletcher T. J., Saintonge A., Soares P. S., Pontzen A., 2020, preprint (astro-ph/2002.04959), p. arXiv:2002.04959Förster Schreiber N. M., et al., 2009, ApJ, 706, 1364Genel S., et al., 2014, MNRAS, 445, 175Giovanelli R., et al., 2005, AJ, 130, 2598Girichidis P., et al., 2020, Space Sci. Rev., 216, 68Glover S. C. O., Federrath C., Mac Low M. M., Klessen R. S., 2010, MNRAS,404, 2Gnedin N. Y., Kravtsov A. V., 2011, ApJ, 728, 88Grassi T., Bovino S., Schleicher D. R. G., Prieto J., Seifried D., SimonciniE., Gianturco F. A., 2014, MNRAS, 439, 2386Hayward C. C., et al., 2020, arXiv e-prints, p. arXiv:2007.01885Heyer M., Dame T. M., 2015, ARA&A, 53, 583Inoue S., Yoshida N., 2019, MNRAS, 488, 4400Jones M. G., Haynes M. P., Giovanelli R., Moorman C., 2018, MNRAS, 477,2Keating L. C., et al., 2020, arXiv e-prints, p. arXiv:2001.08197Krumholz M. R., 2014, MNRAS, 437, 1662Krumholz M. R., 2017, Star Formation. World Scientific Publishing Co. Pte.Ltd., doi:10.1142/10091Kuhlen M., Krumholz M. R., Madau P., Smith B. D., Wise J., 2012, ApJ,749, 36Lagos C. D. P., Baugh C. M., Lacey C. G., Benson A. J., Kim H.-S., PowerC., 2011, MNRAS, 418, 1649Li Q., Narayanan D., Davè R., Krumholz M. R., 2018, ApJ, 869, 73Li Y., Gu M. F., Yajima H., Zhu Q., Maji M., 2020, MNRAS, 494, 1919Licquia T. C., Newman J. A., 2015, ApJ, 806, 96Naiman J. P., et al., 2018, MNRAS, 477, 1206Narayanan D., Krumholz M., Ostriker E. C., Hernquist L., 2011, MNRAS,418, 664Narayanan D., Krumholz M. R., Ostriker E. C., Hernquist L., 2012, MNRAS,421, 3127Nelson D., et al., 2018a, arXiv e-prints, p. arXiv:1812.05609Nelson D., et al., 2018b, MNRAS, 475, 624Obreschkow D., Croton D., De Lucia G., Khochfar S., Rawlings S., 2009,ApJ, 698, 1467Orellana-González G., et al., 2020, MNRAS, 495, 1760Pakmor R., Marinacci F., Springel V., 2014, ApJ, 783, L20Pillepich A., et al., 2018, MNRAS, 473, 4077Popping G., Behroozi P. S., Peeples M. S., 2015, MNRAS, 449, 477Popping G., et al., 2019, ApJ, 882, 137Puchwein E., Haardt F., Haehnelt M. G., Madau P., 2019, MNRAS, 485, 47Saintonge A., et al., 2017, ApJS, 233, 22Schäbe A., Romano-Díaz E., Porciani C., Ludlow A. D., Tomassetti M., 2020,arXiv e-prints, p. arXiv:2003.04329Springel V., Hernquist L., 2003, MNRAS, 339, 289Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328,726Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776Vallini L., Pallottini A., Ferrara A., Gallerani S., Sobacchi E., Behrens C.,2018, MNRAS, 473, 271Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L.,2013, MNRAS, 436, 3031Weinberger R., et al., 2017, MNRAS, 465, 3291Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 1995, ApJ,453, 673Yepes G., Kates R., Khokhlov A., Klypin A., 1997, MNRAS, 284, 235MNRAS , 1–13 (2020) he CO Universe Zhou L., et al., 2017, MNRAS, 470, 4573This paper has been typeset from a TEX/L A TEX file prepared by the author. MNRAS000