Spritz is sparkling: simulated CO and [CII] luminosities

We present a new prediction of the luminosity functions of the [CII] line at 158 $\mu$m, of the CO lines from J=0 to J=24, and of the molecular gas mass density up to z=10, using the Spectro-Photometric Realisations of Infrared-selected Targets at all-z (SPRITZ) simulation (Bisigello et al. 2021). We update the state-of-the-art phenomenological simulation SPRITZ to include both the CO ($J\leq24$) and [CII] line luminosities. This has been performed using different empirical and theoretical relations to convert the total infrared luminosity (or star formation rate) to [CII] or CO luminosity. The resulting line luminosity functions have been compared for validation with a large set of observations available in the literature. We then used the derived CO and [CII] line luminosities to estimate the molecular gas mass density and compare it with available observations. The CO and [CII] luminosity functions presented here are well in agreement with all the available observations. In particular, the best results for [CII] are obtained deriving the [CII] luminosity directly from the star formation rate, but considering a dependence of this relation on the gas metallicity. For all the CO luminosity functions, the estimates favoured by the data are derived considering different relations, depending on the ionisation mechanism dominating each galaxy, i.e. star formation or active galactic nuclei, and, moreover, deriving the $J\geq4$ CO lines directly from the [CII] luminosity. However, further data are necessary to fully discriminate between models. Finally, the best agreement with observations of the molecular gas mass density are derived by converting the [CII] luminosity to H2 mass, using a [CII]-to-H2 conversion ~130 $\rm M_{\odot}/{\rm L}_{\odot}$. All the line luminosity functions, useful for planning and interpreting future observations, are made publicly available.


Introduction
The molecular phase of the interstellar medium (ISM) is the birthplace of stars, and therefore it plays a central role in galaxy evolution (see e.g. the review by Tacconi et al. 2020).The direct detection of molecular hydrogen (H 2 ) in galaxies is hampered by the fact that this molecule, lacking a permanent dipole moment, possesses no corresponding dipolar rotational transition.The lowest energy transitions of H 2 are the purely rotational quadrupole lines that require high temperatures (T > 500−1000 K) in order to be excited (Bolatto et al. 2013).For this reason carbon monoxide (CO), which is the most abundant molecule after H 2 and is easily excited even in cold molecular clouds, is usually used to estimate the molecular gas mass (M H 2 ) via the CO(1−0) emission, assuming a CO-to-H 2 conversion factor, α CO (e.g.Bolatto et al. 2013;Decarli et al. 2019;Riechers et al. 2019).
CO detections at high-z are almost exclusively reported in rare, highly star-forming, sub-millimetre galaxies (e.g.Jarugula et al. 2021;Dye et al. 2022) and quasars (e.g.Carniani et al. 2019;Pensabene et al. 2021), albeit with some exceptions (D'Odorico et al. 2018;Pavesi et al. 2019).With the Atacama Large Millimetre/submillimetre Array (ALMA), the detection of CO from intermediate redshifts (z ≈ 1−2) has become feasible for normal star-forming galaxies as well (Valentino et al. 2020), but at z > 5 it remains extremely timeconsuming (e.g.Vallini et al. 2018).This is primarily due to: (i) the overall lower metallicity and dust abundance of early galaxies, resulting in CO being easily dissociated (e.g.Madden et al. 2020); and (ii) the effect of the increased cosmic microwave background (CMB) temperature that represents a stronger background against which CO lines are observed (da Cunha et al. 2013).
Given the difficulties of observing the CO emission in faint galaxies beyond the Local Universe, several works (e.g.Keating et al. 2020) have focused on the exploitation of the global CO line emission signal from unresolved sources using the so-called line intensity mapping (LIM) technique.Models of the expected LIM signal require the derivation of the line luminosity-halo mass relation, which has often been obtained through hydrodynamical simulations and semianalytical models (SAMs; e.g.Lidz et al. 2011;Gong et al. 2011;Mashian et al. 2015;Li et al. 2016;Sun et al. 2019;Yang et al. 2022).However, this relation can also be inferred from the observed CO luminosity function (LF, Padmanabhan 2018) using an abundance matching technique analogous to that assumed for the stellar mass-halo mass relation (e.g.Behroozi et al. 2010).Nevertheless, an extensive use of this approach has been hampered, up to now, due to the sparse CO observations available and the resulting huge uncertainties regarding the evolution of the CO LFs.
CO fails, sometimes, in tracing the whole H 2 mass, particularly in low-metallicity galaxies where the reduced dust content results in a deeper penetration of far-ultraviolet photons, which are able to dissociate the CO but not the self-shielded H 2 .The H 2 thus survives outside the CO region (e.g.Gnedin & Draine 2014) in the so-called CO-dark clouds (e.g.Wolfire et al. 2010) and can instead be efficiently traced by another (much brighter) proxy of cold gas, namely the [C ii] line at 158 µm.[C ii] is the major coolant of the cold diffuse medium (Wolfire et al. 2003), and dense photodissociation regions (PDRs, Hollenbach & Tielens 1999;Wolfire et al. 2022, for a recent review) associated with molecular clouds.Most importantly, it is now routinely detected in large samples of galaxies at z > 4−5, such as those targeted by the ALMA Large Program to INvestigate CII at Early Times (ALPINE; Le Fèvre et al. 2020;Gruppioni et al. 2020;Yan et al. 2020;Loiacono et al. 2021) and the Reionization Era Bright Emission Line Survey (REBELS; Bouwens et al. 2022).
The [C ii] line is a reliable tracer of the total molecular gas mass (Zanella et al. 2018;Madden et al. 2020;Vizgan et al. 2022), and thus is a fundamental tool for following the cosmic evolution of the fuel in star formation, and helps to better understand how the gas supply in galaxies has moderated the star formation rate (SFR) across the history of the Universe.
The study of the molecular gas mass density is fundamental for understanding the physical processes that are driving the change in the star formation rate density (SFRD) occurring at cosmic noon (e.g.Madau & Dickinson 2014).Indeed, it is still a matter of debate whether this is due to a lack of cold gas supply, or a lower efficiency in converting the gas into stars, or to the presence of strong outflows, preventing the infall of new cold material.The simplest scenario would expect the SFRD to mirror the cold gas evolution, as gas is being consumed by star formation (e.g.Driver et al. 2018).To further improve our understanding on this topic, it is desirable to complement targeted studies with blind measurements to derive the CO (or [C ii]) LF at different cosmic epochs.The recent ALMA Spectroscopic Survey in the Hubble Ultra-Deep Field (ASPECS; Walter et al. 2016;Decarli et al. 2019;Boogaard et al. 2020) was designed exactly for this purpose.
At the same time, SAMs and empirical models of [C ii] and CO LFs have started providing predictions for the LF evolution, and, most importantly, a framework within which the upcoming data can be interpreted (Obreschkow et al. 2009;Lagos et al. 2012;Vallini et al. 2016;Popping et al. 2019a).However, the majority of these models have difficulty in reproducing the bright end of the observed CO LFs at z > 1 (Decarli et al. 2019;Riechers et al. 2019), similarly to what has been observed for other related quantities, such as the total infrared (IR) LF (e.g.Gruppioni et al. 2015;Alcalde Pampliega et al. 2019) or the dust mass (e.g.Pozzetti & Mannucci 2000;Calura et al. 2017;Magnelli et al. 2020;Pozzi et al. 2021).
An alternative approach is based on the exploitation of empirical relations to associate the nebular line emission with dark-matter halos, as was recently done by Chung et al. (2020) and Bethermin et al. (2022).The work by Chung et al. (2020) is based on the halo-galaxy connection presented by Behroozi et al. (2019), which includes the observed UV LFs as constraints, while Bethermin et al. (2022) adopt the stellar mass functions (see also Béthermin et al. 2017).In this paper, we consider a similar empirical approach by extending the work presented in Vallini et al. (2016) to derive the evolution of the [C ii] and CO LFs, together with the molecular gas mass density.
In particular, our work uses different constraints with respect to Chung et al. (2020) and Bethermin et al. (2022), as it is based on the state-of-the-art Spectro-Photometric Realisations of Infrared-selected Targets at all-z (Spritz; Bisigello et al. 2021, hereafter B21) simulation, which uses as input the observed IR LF (Gruppioni et al. 2013) and is not linked to any dark-matter only simulation.
The paper is organised as follows.The Spritz simulation is described in detail in Sect.2, while in Sect. 3 we list all the relations considered to include CO and [C ii] in Spritz.We compare the CO and [C ii] LFs with observations available in the literature in Sect. 4. In Sect. 5 we focus on the molecular gas mass, describing its derivation in Spritz and comparing it with observations, and we finally report our conclusions in Sect.6.We consider a ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.27 and Ω Λ = 0.73, and a Chabrier initial mass function (Chabrier 2003).

The Spritz simulation
The CO and [C ii] line luminosities reported in this paper were obtained from the Spritz simulation, described in detail by B21, which includes elliptical galaxies, dwarf irregulars, star-forming galaxies, and active galactic nuclei (AGN).Spritz is derived starting from a set of galaxy stellar mass functions (GSMFs) and LFs, mainly in the IR.Then, a spectral energy distribution (SED) template is assigned to each simulated galaxy, allowing us to make predictions for several past, current, and future facilities covering different ranges of the electromagnetic spectrum, from the X-ray to the IR.By using the SED templates and a set of empirical and theoretical relations, we derived all the main physical properties for each simulated galaxy, including stellar masses and line luminosities.We now focus on the part of the simulation relevant for this work.
As mentioned before, all the simulated galaxies were extracted from a set of GSMF or LFs.First, we included the IR LFs, as derived from Herschel observations by Gruppioni et al. (2013).These LFs were estimated for different galaxy populations, including normal star-forming galaxies (hereafter spirals), starburst galaxies (SBs), and two composite systems (SF-AGN and SB-AGN).The latter describe two populations with AGN components that are, however, not the dominant source of power, except in the mid-IR and, partially, in the X-ray.In particular, SF-AGN contain intrinsically faint AGN hosted by star-forming galaxies, while SB-AGN host bright, but heavily obscured AGN hosted by SBs.All of these LFs are extrapolated at z > 3, where Herschel observations are not available, by assuming a constant characteristic luminosity (L * ) and by decreasing the number density at the knee (Φ * ) as ∝(1 + z) k Φ .For the power-law exponent k Φ , we considered a range of values from −4 to −1 to span different possible scenarios.
Second, part of the simulated galaxies were extracted starting from the IR LF of AGN-dominated systems (un-obscured AGN 1 and obscured AGN 2) derived by B21 starting from Herschel observations, complemented by far-ultraviolet observations up to z = 5 (e.g.Croom et al. 2009;McGreer et al. 2013;Ross et al. 2013;Akiyama et al. 2018;Schindler et al. 2019).The LF is described by a modified Schechter function and its evolution was extrapolated at z > 5 following the observations at lower z (i.e. A193, page 2 of 18 Third, we also extracted galaxies from the K-band LF of elliptical galaxies (hereafter, Ell), estimated by averaging the LFs by Arnouts et al. (2007), Cirasuolo et al. (2007), and Beare et al. (2019).At z > 2 the LF was extrapolated by keeping the characteristic luminosity constant and decreasing the number density at the knee as ∝(1 + z) −1 .This assumption has little impact, as the number density of elliptical galaxies at z = 2 is already quite low, and it will be tested with future observations.Finally, we included the GSMF of dwarf irregular galaxies (hereafter, Irr) derived by Huertas-Company et al. (2016) and complemented by the local GSMF of irregular galaxies observed in the Galaxy And Mass Assembly (GAMA) survey (Moffett et al. 2016).More details are given in Appendix A.
As mentioned before, for each simulated galaxy we derived the main physical properties (e.g.stellar mass, star formation rate, stellar metallicity, and luminosities).In particular, for the majority of them, the total IR luminosity L IR was directly taken from the best template used to derive the observed IR LF.For Irr and Ell, L IR was instead obtained starting from either the galaxy stellar mass or the K-band luminosity, assuming an SED template (Polletta et al. 2007;Bianchi et al. 2018, for dwarf elliptical and irregulars galaxies, respectively).
The IR component of the SFR was derived from the L IR assuming the Kennicutt (1998a) conversion, while the UV component was derived from the luminosity at 1600 Å not corrected by dust absorption, assuming the Kennicutt (1998b) relation.The 1600 Å luminosity was derived directly from the SED template associated with each simulated galaxy, which depended on the galaxy population to which it belongs, and was taken from a set of 35 empirical templates (Polletta et al. 2007;Rieke et al. 2009;Gruppioni et al. 2010;Bianchi et al. 2018).These empirical templates are of low-z galaxies, but they represent a good description of galaxies observed by Herschel up to z = 3.5 (Gruppioni et al. 2013).The same procedure applied to derive the UV luminosity was performed to obtain the galaxy stellar mass, as each template was normalised to 1 M .Finally, the stellar metallicity was derived from the galaxy stellar mass assuming the mass-metallicity relation by Wuyts et al. (2014): where the asymptotic metallicity is Z 0 = 8.69 and the power-law slope at low metallicity is γ = 0.40.The Spritz simulation is in agreement with a large set of observations, ranging from number counts at different wavelengths to the total GSMF (see Sect. 4 in B21).Of particular interest for this work is the agreement with the observed IR LF at z ∼ 5 (Gruppioni et al. 2020), which was not included as input in the simulation.In particular, among the different high-z extrapolations tested, the best agreement is present when assuming that the number density at the knee of the IR LF evolves as ∝(1 + z) −1 for spirals, SBs, SF-AGN, and SB-AGN.This agreement supports the validity of the extrapolation performed at z > 3.
On the contrary, some tensions are present between Spritz and the observed UV LF (see Fig. 20 in B21).In particular, the bright end (M 1600 Å < −22.5 at z = 0.5 and −21.0 at z = 1.5) of the galaxy UV LF is overestimated at z ≤ 1.5, while the faint end is underestimated, particularly at z > 2.
When looking at the SFR-M * plane in Spritz, star-forming galaxies correctly populate the galaxy main sequence, whose normalisation increases with increasing redshift, as also visible in observations (e.g.Noeske et al. 2007;Speagle et al. 2014;Bisigello et al. 2018).SBs are correctly placed above the main sequence, while elliptical galaxies, by construction, are placed below.However, at z > 4, the simulation does not have galaxies with a specific SFR high enough to account for the observed SBs (Caputi et al. 2017).The mentioned discrepancies in the faint end of the UV LF and SFR-M * plane can be due either to the absence of a particularly dust-poor galaxy population, which has not been previously observed by Herschel or included in the dwarf irregular galaxy population, or to a limitation on the set of templates.We analyse the impact of these discrepancies on our results in Appendix B.
In Spritz, the spatial distribution of galaxies, which is fundamental for line intensity mapping, is included starting from the observed two-point correlation function using the algorithm by Soneira & Peebles (1978).All of the details of this procedure are reported in B21.Briefly, galaxies are distributed following an angular correlation function w(θ) = A w θ 1−γ , with a power-law slope δ = γ − 1 = 0.7 as suggested by observations (e.g.Wang et al. 2013).At the same time, the spatial correlation length r 0 , of which A w represents the angular projection, has a dependence on stellar mass, as derived from observations (Wake et al. 2011;Hatfield et al. 2016): with the stellar mass break log 10 (M * break /M ) = 10.6, a low-mass slope of k M,1 = 0.0959 ± 0.0003 and a high-mass slope of k M,2 = 0.181 ± 0.006.The full procedure is repeated, splitting the mock catalogue on different redshift slices.
In the Spritz workflow, the spatial distribution derived as a function of stellar mass, as just mentioned, is obtained after assigning physical properties to each simulated galaxy.Therefore, a reader, if preferred, can ignore the included method and match the Spritz catalogue to a dark-matter simulation, using a stellar-to-halo mass relation (e.g.Girelli et al. 2020).
More information on Spritz and its comparison with observations is available in B21.

[C II]
In this and the following section, we summarise the relations considered to include the [C ii] and CO emission lines in Spritz.
In particular, for the first one, we assumed three different methods to derive the expected [C ii] emission (Fig. 1).
First, we included the empirical relation derived by Gruppioni et al. (2016) starting from a local sample of Seyfert galaxies observed with Herschel: IR is the component of the IR luminosity due to starformation.In their work, they verify that this relation does not change with the AGN fraction, and thus it can be used either in sources with no or little AGN contribution, or in AGNdominated objects.
Second, we considered the empirical relation between SFR and L [C ii] proposed by De Looze et al. ( 2014) using a broad sample of galaxies, including dwarfs.In particular, they derived different relations for the overall galaxy sample and for different sub-samples of galaxies, namely metal-poor dwarf galaxies, starforming galaxies or SBs, composite or AGN sources, and galaxies at high-z (z > 0.5).The relations are in the form of: A193, page 3 of 18 we converted L SF IR to SFR using the relation by Kennicutt (1998a).In the bottom right, we show the uncertainties associated with each relation.
whose α and β depend on redshift and galaxy type (see Table 3 in De Looze et al. 2014).
Finally, we included the relation by Vallini et al. (2015, hereafter V15) derived by post-processing a radiative-transfer cosmological simulation at z ∼ 6.6.The post-processing was constructed to obtain, on a cell-by-cell basis, the [C ii] emission from dense PDRs and from the diffuse neutral gas (Wolfire et al. 2003).In particular, the [C ii] emission from V15 is well described by the following analytical relation, depending on metallicity (Z) and SFR: (5) For this relation we considered a scatter similar to that derived by De Looze et al. ( 2014) for star-forming galaxies, that is 0.27 dex, and a larger one of 0.37 1 dex for galaxies at z > 5 (Carniani et al. 2018).

CO lines
The CO LF evolution at different J is a great tool for shedding light on the underlying physical properties of different galaxy populations across cosmic time (see e.g.Obreschkow et al. 2009;Lagos et al. 2012;Vallini et al. 2016;Popping et al. 2019a).The relative luminosity of different CO lines, also referred to as CO spectral line energy distribution (CO SLED), in fact gives unique insights into the gas density and temperature, and into the heating mechanisms acting in the ISM of galaxies (e.g.Meijerink et al. 2007;Rosenberg et al. 2015;Pozzi et al. 2017;Mingozzi et al. 2018;Talia et al. 2018; 1 Derived from the observed scatter of 0.48 dex and considering an error on the SFR of 0.3 dex. Table 1.Relations used to include CO lines in Spritz for different J and galaxy populations.As for the [C ii] line, we considered different recipes to include CO J → (J − 1) rotational transitions in Spritz.We assumed different CO excitation depending on the galaxy population, as summarised in Table 1 and described in detail in the following sections.In Appendix C we show the CO SLEDs associated with each galaxy type using the different relations.
The reader can refer to Fig. C.1 to understand which relation is preferable in order to minimize discontinuities in the CO SLED.We assumed that elliptical galaxies do not emit in CO, as few of them have been observed in CO, even in the Local Universe (∼22%; Young et al. 2011) and, even through stacking, they generally show a low gas fraction (<1% in the Local Universe and <8% at z ∼ 1.5; Magdis et al. 2021).

J < 14
In order to derive the CO J < 14 luminosities2 , we considered different relations from the literature, either based on the ratio of the CO J → (J − 1) transitions with respect to the CO(1−0), or with respect to the [C ii] line.Moreover, we also tested relations between the IR and the CO luminosities.In particular, to derive the CO(1−0) of star-forming galaxies with or without a non-dominant AGN component (i.e.spirals, SBs, SF-AGN, Irr), we considered the relation by Sargent et al. (2014): This was obtained in a sample of 131 galaxies at z < 3 with M * > 10 10 M , and has a scatter of 0.21 dex.
For the same galaxy populations, but from J = 2 to J = 8, we considered the recent L CO(J→(J−1)) /L CO(1−0) ratios by Boogaard et al. (2020), which were derived from observations of 22 star-forming galaxies up to z = 3.6, as part of the ASPECS survey (Walter et al. 2016).The ratios for J from 2 to 6 are reported in Table 2 for galaxies above and below z = 2, except for J = 7 and 8, for which a single value is present, corresponding to observations at z > 2. To avoid discontinuities between z < 2 and z > 2, we interpolated between the two different ratios at z = 1.5−2.5.To estimate the CO ratios, we used the CO(1−0) luminosity derived using the previously mentioned relation by Sargent et al. (2014).
Some studies have shown that the CO(1−0)-SFR relation is different for low-metallicity galaxies (e.g.Cormier et al. 2014).Therefore, we applied a correction to the CO(1−0) luminosities for galaxies with sub-solar metallicity (i.e. 12 + log 10 (O/H) < 8.7), following the results derived by Hunt et al. (2015): For the galaxy populations with a dominant AGN component (i.e.SB-AGN, AGN1, and AGN2) we considered the relation between CO(1−0) and the total IR luminosity by Greve et al. (2014): with a scatter of 0.27 dex.The relation was derived from a sample of 62 local ultra-luminous infrared galaxies (ULIRGs), but consistent results were obtained including AGN-dominated systems (Greve et al. 2014).We also included the relations presented in the same paper to convert the total IR luminosity to the luminosity of the J = 2−13 CO transitions.We also explored another possibility for deriving the 1 < J < 14 transitions for the same galaxy populations, namely we considered the CO(J → (J − 1))/CO(1−0) ratios from a sample of 35 local AGNs (L 2−10 keV ≥ 10 42 erg s −1 ; Esposito et al. 2022).In particular, we derived the CO(J → (J − 1))/CO(1−0) ratios after cross-matching the sample from Esposito et al. (2022) with the one from Gruppioni et al. (2016) to identify galaxies with an AGN fraction at 5−40 µm above and below 40% (eight and seven objects, respectively).This threshold on the AGN fraction separates, in Spritz, the SF-AGN ( f AGN < 40%) from the other AGN populations (i.e.SB-AGN, AGN1, and AGN2).We then derived the median CO(J → (J − 1))/CO(1−0) ratios for both AGN populations (Table 3) and we normalised such ratios to the CO(1−0) derived in Eq. ( 8).
In Liu et al. (2015), the FIR (40−400 µm) luminosity of 167 local galaxies with Herschel spectroscopic observations is related to the CO luminosity from J = 4 to J = 12, as: with the values of N and A reported in Table 4.In Spritz we applied these relations to each simulated galaxy, without separating for galaxy population and redshift.We note that beyond J = 10, no observed LFs are available for comparison.
For J ≥ 4, we included the CO luminosities derived considering the CO(J → (J − 1))/[C ii] ratio from Rosenberg et al.
(2015) and considering the three classes of objects presented on the same paper.In particular, the first class (c1) does not require any mechanism in addition to the UV-heating from starformation to reproduce the observed CO ladder, while the third class (c3) includes galaxies with an AGN-component and probably requires mechanical heating in addition to UV-heating to describe its excited CO ladder.The second class (c2) simply indicates an intermediate case, where it is not possible to discriminate which heating mechanism dominates the CO ladder.
The CO(J → (J − 1))/[C ii] ratio should be more stable than Notes. (a) We calculated this value averaging the two ratios at J − 1 and J + 1.
The ratios for the different classes are reported in Table 5.In this work we first considered the three extreme cases where all galaxies behave as a single class, and then we examined the case where spirals, SBs and Irr are in c1, SF-AGN are in c2, and SB-AGN, AGN1, and AGN2 are in c3.For the [C ii] luminosities, we fixed the ground luminosities to be those derived considering the relation by Vallini et al. (2015, see Sect.(10) We also included this relation in Spritz with a scatter of 0.18 dex, as reported in the reference paper, without distinguishing between different galaxy populations or redshifts.

J ≥ 14
At present, no relations are available in the literature to derive the CO transitions with J ≥ 14 starting from the SFR or the IR luminosity.Therefore, we decided to adopt observed ratios as reference for our simulated galaxies.In particular, we considered the CO(J → (J − 1))/CO(1−0) (Table 6) estimated from observations by Mashian et al. (2015), using NGC 6240 as a reference for SB-AGN, Mrk 231 for AGN1 and AGN2, and M 82 for SBs.These galaxies are also among the templates included in Spritz to derive photometry and physical properties for the same galaxy populations.
The CO SLED of NGC 6240 is detected up to J = 24, while Mrk 231 and M 82 are detected up to J = 20 and J = 18, respectively.Beyond these transitions, only upper limits are available, and we therefore considered no CO emission.The CO SLED of Mrk 231 has been previously studied (van der Werf et al. 2010;Vallini et al. 2019), showing that the excitation of the CO J > 8 lines cannot be completely reproduced considering the PDR emission only, but also requires an XDR component created by the X-rays from the accretion onto the central black hole.Moreover, the CO emissions for J ≥ 13 are completely dominated by the emission coming from the XDR.Given the absence of such a source of high X-ray excitation in the spiral, SF-AGN and dwarf populations, we assumed that their CO J ≥ 14 transitions are negligible.

Comparison with CO and [C ii] observations
In this section we compare some observed [C ii] and CO LFs with the LFs in Spritz derived considering all the relations previously discussed.All LFs are made publicly available3 .We also report the 1σ confidence intervals associated with each LF.These intervals were obtained taking into account the uncertainties associated with the observed LFs or GSMFs used as inputs in the simulation (see Sect. 2) and the errors associated with the relations used to derive the [C ii] or CO line luminosities.
The LFs were derived considering the full Spritz catalogue (i.e.L IR > 10 6 L ) and the entire sky, as the observed LFs considered in the following sections for the comparison were corrected for incompleteness and include cosmic variance in their errors.

[C II]
In Fig. 2 2020) are derived from a sample of UV-selected galaxies; for this reason, they may be affected by observational biases.
Given the absence of spectroscopic instruments covering the wavelength range between 160 and 600 µm, no observed  At z > 4 the observed values show a significant dispersion and all the relations, which mainly differ at L [C ii] > 10 9.5 L , are broadly consistent with the observations.The LFs reported in the figure correspond to the flattest high-z extrapolation (i.e.Φ * ∝ (1 + z) −1 ) included in Spritz, but we also report, as an example, the LF derived considering the relation by Gruppioni et al. (2016) and a number density at the knee (Φ * ) decreasing as ∝(1 + z) −4 .The latter LF is well below the observed values, showing that the data are consistent with the first extrapolation (i.e.Φ * ∝ (1 + z) −1 ), as also observed for the total IR LFs (B21).
In Fig. 3, we split the [C ii] LF, which were derived considering the relation by V15, into the different contributions of the single galaxy populations.In this way, we can appreciate that, in the Local Universe, the [C ii] LF is dominated at all luminosities by spiral galaxies, while at z > 4 it is dominated by dwarf irregular galaxies at L [C ii] < 10 8.5 L and by SBs and SB-AGN at brighter luminosities.These two populations have specific SFRs (sSFRs) ranging from log(sSFR/yr −1 ) = −8.8 to −8.1, and therefore they are considerably above the main sequence at low-z, but they are in the main sequence at z = 5−7.Indeed, following the parametrisation by Speagle et al. (2014), the main sequence at z = 5 corresponds to a log(sSFR/yr −1 ) = −8.5 to −7.9, depending on the stellar mass, with an observed scatter of 0.3 dex.This is also consistent with the results by Faisst et al. (2020), who found that ALPINE sources at z = 4−6 are star-forming galaxies on the main sequence.We remind the reader that, as explained in Sect.2, the templates associated with each galaxy population do not evolve with redshift.

Going into further detail, in Spritz the faint end of the [C ii]
LF moves from being dominated by spirals to being dominated by dwarf irregulars around z ∼ 1 (Fig. 4).On the other hand, the contribution of the SBs and the SB-AGN populations to the bright end of the [C ii] LF becomes dominant at z ∼ 1 and z ∼ 2.8, respectively.However, it is evident that more observations are needed between z = 0.3−4 to verify these predictions.
Given LF in different luminosity regimes.

CO lines
Here we compare the observations by Saintonge et al. (2017), Riechers et al. (2019), and Decarli et al. (2019, 2020), all corrected for incompleteness in the respective works, with the different relations (see   from a mass-selected sample of a galaxy (M * > 10 9 M ) at z < 0.05, while the other two works are based on blind line observations.In this section we analyse only the CO transi-tions with J ≤ 10, as no observed LFs are available for higher J values.To facilitate the comparison between the different models, we report in Table 7 the total χ 2 derived by comparing each model with the available observations4 , taking into account the observational errors.

Low J
In Figs. 5 and 6, we report the results for J = 1 to 3. The Spritz LFs are in agreement with the observations for these J values at the available redshifts, with no strong difference between the results derived with the different relations (∆χ 2 ≤ 1).The main discrepancy between models and observations is for J = 2 at z ∼ 1.4,where the Spritz LFs are lower than the observed LFs by ∼0.5 dex, but they are still consistent within the errors.This difference is possibly linked to the light offset present at z = 1.0−1.7 between the parametric description of the redshift evolution of the IR LF included as input in Spritz and the observed Herschel IR LF (see Appendix E).As for the [C ii] LF, the slope of the density extrapolation Φ * ∝ (1 + z) −4 in Spritz is discarded by the data, as it leads to a significant underestimation (i.e. up to 0.8 dex) of the observed J = 2 LF at z 6 (see dotted black line in Fig. 6).
In the same Figs.5 and 6, we also report the models by Popping et al. (2016Popping et al. ( , 2019a)), which are close to our predictions for the CO(1−0) at z < 0.3 and CO(3−2) at z ∼ 0.5.However, at z > 2 the agreement between the data and our predictions is much better than with those of the Popping et al. (2016Popping et al. ( , 2019a) ) Table 7. χ 2 derived comparing the different models with CO observations.models, particularly at the bright end.To investigate this discrepancy further, it would be interesting to verify if those models are also missing, or at least under-predict, the most massive and dusty galaxies, and therefore underestimate the bright end of the IR LF.

Mid-and high-J
In Fig. 7 we report the comparison between the results from Spritz and the observed CO LFs from J = 4 to J = 10.We tested different relations in Spritz to estimate the mid-J CO LFs, as described in Sect.3.2 and summarised in Table 1.
The faint-end slopes of the CO(4−3) LFs are similar for all the considered relations and they are within 0.5 dex, even for J > 4. The main differences between the considered predictions are, therefore, on the knee position and the bright-end slope of the LFs.The CO(4−3) observational data cover a narrow luminosity range around the knee at z ∼ 1, with all the considered relations slightly below the observations, but still consistent within the error bars.This may be linked to the light offset between the input Spritz IR LF and the observed one at z = 1.0−1.7 (see Appendix E).At z = 3.5−3.9 the Spritz and the observed LFs have different shapes, with the knee of the latter being at lower luminosities and at higher densities, and the bright-end slope being steeper than the bright-end slopes predicted by Spritz with any relation.This discrepancy has no obvious explanation, and more observations over a larger luminosity range (ASPECS observations correspond to a single independent luminosity bin) and at additional redshfits are needed to investigate it further.
For CO(5−4) and CO(6−5), the observed LFs are generally higher than the different relations included in Spritz.One exception is the LF derived using the relation by Liu et al. (2021), which includes a dependence of the CO(5−4)/CO(2−1) ratio on the IR luminosity.This indicates that the different CO transitions may have a different dependence on the IR luminosity.However, the observed underestimation of the CO(5−4) LF may be linked to the observed underestimation of the CO(1−0) LF at the same redshift (i.e.z = 1.2−1.6,see previous section and Appendix E).No CO(1−0) observations are available at z = 1.7−2.1 to further investigate the offset observed for CO(6−5), but we highlight that the observed IR LF at these redshifts are perfectly reproduced in Spritz, and therefore we would expect to equally reproduce the CO(1−0) LF (assuming the CO(1−0)-L IR does not strongly evolve with redshift).Finally, for transitions with J > 6, all of the predictions are consistent with the observations within the error bars.
For comparison, the model by Popping et al. (2016) generally follows the observed LFs for J = 4 to 6, with some A193, page 9 of 18 discrepancies at z > 2, which is present only for J = 4.As mentioned in the previous section, the same model also strongly differs from the observed LFs for J < 3 at the same redshifts.
Overall, looking at all the high-J LFs and the available observations, the best agreement for Spritz is obtained by considering the CO(J → (J − 1))/[C ii] ratios for the different classes presented in Rosenberg et al. (2015), as they generally result in smaller χ 2 values than the other relations (see Table 7).For the classes by Rosenberg et al. (2015), with increasing redshift or J values, the inclusion of only high-excitation classes (c2 and c3) seems to be preferred by the observations.However, more observations, for instance those of J > 7 at z < 1, are necessary to distinguish whether the better agreement with the highexcitation classes is driven by the redshift evolution (i.e. the excitation increases with redshift) or by the J value (i.e. the excitation increases with J value), or a combination of both.At the same time, more observations of the bright end (L CO(4−3) > 10 10 K km s −1 pc 2 ) are necessary to further discriminate between the different relations considered.

Molecular gas mass
As previously mentioned, CO can be efficiently used to trace the molecular content of a galaxy.However, sometimes H 2 may survive outside the CO regions, in the so-called CO-dark clouds (e.g.Wolfire et al. 2010), and instead can be efficiently traced by [CII] emission (e.g.Zanella et al. 2018;Madden et al. 2020;Wolfire et al. 2022).For this reason, to calculate the molecular gas mass of each simulated galaxy in Spritz, we decided to use both proxies.
In particular, on the one hand, we derived the molecular gas mass directly from the CO(1−0) by considering the Milky Way value α CO = 4.3 M (K km s −1 pc 2 ) −1 for normal star-forming galaxies (Spiral and SF-AGN) and the value derived for ULIRGs (α CO = 0.86 M (K km s −1 pc 2 ) −1 ) for the most active galaxies (intense star-formation or nuclear activity, i.e.SBs, SB-AGN, AGN1, and AGN2).For dwarf galaxies we considered a metallicity-dependent CO-to-H 2 conversion factor, as derived by Madden et al. (2020), namely α CO = 10 0.58 × (Z/Z ) −3.39 , by taking into account the contribution of the CO-dark clouds.As reference, we considered the CO(1−0) luminosity estimated from the IR emission using the relation by Sargent et al. (2014) for star-forming galaxies, with the correction by Hunt et al. (2015) for galaxies with sub-solar metallicity, and by Greve et al. (2014) for AGN-dominated systems (see Sect. 3.2).
On the other hand, we calculated the molecular gas mass from the [C ii] line luminosity as M H 2 = 10 2.12 (L [C ii] /L ) 0.97 (Madden et al. 2020), without any variation among the different galaxy populations.In this case we considered the [C ii] luminosity estimated using the relation by Vallini et al. (2019), which was the one with the best agreement with the observations (see Sect. 4.1).
To validate the molecular gas mass included in Spritz, derived either from CO or from [C ii] for each simulated galaxy, we estimated the cosmic evolution of the molecular gas mass and compared it with available observations (Fig. 8).
On the one hand, when comparing the molecular gas mass density derived from the predicted CO(1−0) luminosity, it is evident that using the Milky Way α CO value for all the galaxies leads to under-predict the molecular gas mass at z > 0.5 by ∼0.5 dex.This underestimation increases up to 1 dex if we assign a lower α CO value to the most active galaxies (intense star-formation or nuclear activity, i.e.SBs, SB-AGN, AGN, and AGN2), but it is balanced if we include an α CO that varies with metallicity in dwarf galaxies.Therefore, the best option to convert the CO into molecular gas mass seems to be represented by the use of a different α CO for normal star-forming galaxies, active systems, and dwarf galaxies.
On the other hand, a single α C ii value seems to be enough to reproduce the observed shape of the molecular gas mass density, showing a peak around z = 2, then decreasing at lower and higher redshifts.A value of α C ii ∼ 130, as proposed by Madden et al. (2020), is necessary to reproduce the observed normalisation with the Spritz simulation.
For a proper comparison with observations, it is necessary to apply the observational limits of the ASPECS survey5 , as their molecular mass density is derived without extrapolating their CO LF.For example, once the ASPECS observational limits have been applied, the model by Popping et al. (2019a) is a factor of two to three lower than the observations (Popping et al. 2019b).Similar, or even larger, discrepancies are present in Spritz when deriving the molecular gas from the CO, as seen in Fig. 9.These discrepancies are present because the majority of dwarf galaxies are below the ASPECS observational limits.Taking into account the fact that the ASPECS H 2 values were derived from the CO(1−0) to CO(4−3) line luminosities, the discrepancies present in the H 2 may arise from the differences observed for the J = 2 and J = 4 transitions (see Sect. 4.2).
The molecular gas mass density derived from [C ii] is, instead, consistent with the observations, even after applying the observational limits of the ASPECS survey.Observational data show a light decrease in the molecular gas mass at z > 2, while in Spritz the contribution of SBs and SB-AGN (see Fig. 4) keep the molecular gas density almost constant.However, the big uncertainties do not allow this to be investigated further.Therefore, we can conclude that with Spritz it is possible to estimate a reliable molecular gas mass density starting from [C ii], at least up to z = 4.

Summary and conclusions
In this work we used the state-of-the-art Spritz simulation to predict the [C ii] and CO LFs at different redshifts, as well as the molecular gas mass density.In particular, we considered differ- In addition, we included predictions for the LFs of CO with J = 14−24, using the CO SLED of galaxies observed by Mashian et al. (2015) as templates.However, no observed LFs are currently available to test the predictions for CO transitions with J ≥ 14.
For [C ii], at z < 0.3 the best result is obtained by considering the relation by V15, which takes into account not only a dependence of [C ii] on SFR, but also on metallicity.The relation is also consistent with observations at high-z, where, however, more observations at L [C ii] > 10 10 L are necessary to unambiguously discriminate between the considered relations.Future IR spectroscopic observations (from space, given the Earth's atmospheric transmission), covering wavelengths shorter than those sampled by ALMA, will be essential in order to explore intermediate redshifts and provide valuable constraints for the different relations.
A193, page 11 of 18 (2020, dash-dotted orange line).We also report the H 2 mass density derived from the CO, using the CO(1−0) luminosity derived from Greve et al. (2014) and Sargent et al. (2014), and considering different α CO values for star-forming galaxies, ULIRG-like galaxies, and dwarfs (solid blue line).We show the molecular gas mass density obtained by considering a single α CO value for all the galaxies (dotted blue line), or two different α CO values for star-forming and ULIRG-like galaxies (dashed blue line).We report observational results by Decarli et al. (2016Decarli et al. ( , 2019Decarli et al. ( , 2020))    For mid-and high-J CO transitions (i.e.J > 3), the best results are obtained by considering the CO/[C ii] ratios derived for the different classes by Rosenberg et al. (2015), while for the CO(5−4) LF, the relation by Liu et al. (2021), including a further dependence of the CO(5−4)/CO(2−1) ratio on IR luminosity, provides one of the best results, when compared with the available observations.However, all relations are generally consistent with each other in the faint end, and more observations at luminosities brighter than ∼10 10.5 K km s −1 pc 2 are necessary to unambiguously discriminate between the predictions of the different models.
Finally, we integrated the CO and [C ii] LFs, after converting them to molecular gas masses through different recipes, to obtain an estimate of the molecular gas mass density at different redshifts.The evolution of the molecular gas mass density is correctly reproduced by Spritz over the whole redshift range where observations are available (i.e.0 < z < 4), in particular by deriving the H 2 mass directly from the [C ii] LF.
We conclude that the Spritz simulation can be used to predict the evolution of both the [C ii] and CO luminosity, as well as that of the molecular gas mass.This work constitutes a useful reference for any future sub-millimeter and millimeter observations, and strongly outlines the need for a future far-IR spectroscopic instrument that covers the huge gap between the past Herschel observations (λ ≤ 210 µm) and the current ones with ALMA (λ ≥ 300 µm).This would be fundamental to obtain A193, page 12 of 18 statistical samples of galaxies over a continuous redshift range and derive a better understanding of the quantities discussed in this paper.
The CO SLED of dwarf irregular galaxies is generally smooth when considering the relations by Greve et al. (2014), Liu et al. (2015), and Boogaard et al. (2020), while the CO SLED by Rosenberg et al. (2015) is five to ten times higher than the others.This is at least partially due to the different metallicity dependence considered when deriving the [C ii] and CO luminosity.In general, further observations of the CO SLED of dwarf galaxies are necessary to disentangle the different models.
Moving to galaxies with an AGN component, the CO SLED of SF-AGN follows two separated trends as the CO SLEDs by Boogaard et al. (2020), Greve et al. (2014), andLiu et al. (2015) are between 5 and 20 times lower than the CO SLEDs by Esposito et al. (2022) and Rosenberg et al. (2015).We remind the reader that the relations by Esposito et al. (2022) and Rosenberg et al. (2015) are specific for objects hosting a lowluminosity AGN, and should therefore be, therefore better suited to describe SF-AGN, while the other relations are broadly valid for star-forming galaxies.Moreover, for galaxies at high redshifts (i.e.z > 1.5), which should be more contaminated by AGN, the L CO(J−(J−1)) /L CO(1−0) ratio by Boogaard et al. (2020) increases, becoming closer to the ratios by Esposito et al. (2022).
Finally, AGN1, AGN2, and SB-AGN show similar CO SLEDs with a large scatter (i.e. up to a factor of 13) between the different relations.For J > 14, the CO SLEDs by Mashian et al. (2015) are completely in agreement with the CO SLEDs by Esposito et al. (2022) for AGN1 and AGN2, while they are a factor of five lower than the same CO SLED for SB-AGN.The CO SLED derived using results from Greve et al. (2014) is always the steeper one, showing the smallest L CO(J−(J−1)) /L CO(1−0) values.
Overall, switching from one relation to another to estimate the CO luminosity at different J values may produce, in some cases, large discontinuities in the estimated CO SLED.In the future, observations of statistical samples of galaxies over wide J ranges may be used to improve the CO SLED included in Spritz.
, we report the Spritz [C ii] LF compared with observational results by Capak et al. (2015), Yamaguchi et al. (2017), Hemmati et al. (2017), Decarli et al. (2020), Yan et al. (2020), and Loiacono et al. (2021).Results by Capak et al. (2015), Hemmati et al. (2017), and Yan et al. (2020) are corrected for incompleteness, while the correction in the work by Loiacono et al. (2021) is limited, given that it is based on a single detection.The works by Hemmati et al. (2017) and Loiacono et al. (2021) are based on direct [C ii] observations, based on a blind Herschel survey and on an ALMA serendipitous detection, respectively.On the other hand, the results by Capak et al. (2015) and Yan et al. (

[
C ii] LFs are currently available to fill the gap between the Local Universe and z ∼ 4.5.In the same figure, for a direct comparison with our predictions, we also report the [C ii] LFs derived from the SAMs by Popping et al. (2016) and Lagache et al. (2018), taking into account that the latter is valid at 4.7 ≤ z ≤ 8.In the Local Universe, Spritz can reproduce the observed values only when considering the relation by V15, which includes a dependence on both the SFR and metallicity.It is worth noticing that the relation by V15 is in agreement with the relation by De Looze et al. (2014) for solar metallicity, once the respective uncertainties have been considered.

Fig. 2 .
Fig. 2. [C ii] LF derived with Spritz assuming different relations (see Sect. 3.1).We compare the results with [C ii] observations in the Local Universe by Hemmati et al. (2017), at z ∼ 4.45 by Yan et al. (2020) and Loiacono et al. (2021), at z ∼ 5.5 by Capak et al. (2015), at z ∼ 6.3 by Yamaguchi et al. (2017), and at z ∼ 6.9 by Decarli et al. (2020).We also report the model predictions of Popping et al. (2016, dashed magenta lines) and Lagache et al. (2018, dotted red lines) at z ≥ 4.7.The shaded areas (same colours as the solid lines) show the uncertainties of the considered [C ii] relations and the 1σ errors on the Spritz input LFs and GSMFs.
the urgency for future far-IR probes covering the gap between Herschel (z = 0) and ALMA (z > 4) [C ii] observations, in Appendix D we report, for reference, the [C ii] LF derived by Spritz, considering the relation by V15, up to z = 10.Similarly, we also report the area and depth necessary to sample the [C ii]

Fig. 4 .
Fig. 4. Contribution to the [C ii] luminosity density of the different galaxy populations (see legend) included in Spritz.We report two different luminosity ranges: 10 7 L < L [C ii] ≤ 10 8.5 L (top) and L [C ii] > 10 8.5 L (bottom).The grey shaded areas show the redshift range where observations are available.We highlight that SBs and SB-AGN have sSFRs ranging from log(sSFR/yr −1 ) = −8.8 to −8.1, which is above the main sequence at low-z, but on it at z = 5−7.

Fig. 8 .
Fig. 8. Molecular gas mass density derived by Spritz from [C ii], using the [C ii] luminosity derived from V15 and the relation by Madden et al.
,Saintonge et al. (2017),Riechers et al. (2019), and Magnelli et al. (2020, from  the evolution of the dust mass density).The shaded areas include the 1σ uncertainties on the input LF, the high-z extrapolations, and the H 2 conversions.For the sake of clarity, we report uncertainties only for the results obtained from [C ii] and from CO, assuming three different α CO values for star-forming galaxies, ULIRG-like galaxies, and dwarfs.The uncertainties associated with other results derived from CO, under different α CO assumptions, are comparable.

Fig. 9 .
Fig. 9. Molecular gas mass density derived by Spritz, by considering the ASPECS observational limits.The results are derived from [C ii]
For low-J values (J ≤ 3), the CO LFs in Spritz, derived using the relations bySargent et al. (2014), Greve et al. (2014), Boogaard et al. (2020), or Esposito et al. (2022) for different galaxy populations, are in good agreement with the observations.The only discrepancy occurs for the CO(2−1) transition at z = 1.2−1.6,where the Spritz LF is slightly below the observed data, but still consistent within the errors.

Fig
Fig. C.1.Median CO SLED of the different galaxy types included in Spritz: spirals, SBs, Irr, SF-AGN, SB-AGN, and AGN1.The CO SLED of AGN2 is identical to that of AGN1 and is not shown.Different symbols indicate the relations considered to derive the CO luminosity for different J.The dash-dotted cyan lines show the CO SLEDs by Mashian et al. (2015) for J < 14, which are not included in this work, but are shown for consistency.Shaded areas show the one σ variation of the CO SLED for galaxies with 10 10 ≤ L IR /L ≤ 10 12 .

Table 2 .
Boogaard et al. (2020)0) ratios byBoogaard et al. (2020)used for J = 2−8 for spirals, SF-AGN, SBs, and Irr.Notes.The first column indicates the J value, and second and third columns indicate the L CO(J→(J−1)) /L CO(1−0) ratios at z < 2 and z > 2, respectively.For observational limitations, a single ratio is given for J = 7 and 8, and it is used at all redshifts.
Notes.The last column shows the 1σ scatter around the relations.

Table 1
Saintonge et al. (2017)to estimate the CO luminosities.The work bySaintonge et al. (2017)is based on a representative sample of galaxies extracted starting A193, page 7 of 18