CO Excitation, Molecular Gas Density, and Interstellar Radiation Field in Local and High-redshift Galaxies

Star formation in galaxies is regulated by their reservoir of molecular gas. Globally, the star formation rate (SFR) correlates with the total amount of molecular gas mass via the Kennicutt–Schmidt law (Schmidt 1959; Kennicutt 1998). Meanwhile, physical properties like density and temperature of the molecular gas also play an important role. For example, observations of different carbon monoxide (CO) rotational transition (J) lines reveal a relatively denser (${n}_{{{\rm{H}}}_{2}}\sim {10}^{3\mbox{--}4}\,{\mathrm{cm}}^{-3}$), highly excited phase of molecular gas in addition to a more diffuse (${n}_{{{\rm{H}}}_{2}}\sim {10}^{2\mbox{--}3}\,{\mathrm{cm}}^{-3}$), less excited phase (e.g., Harris et al. 1991; Wild et al. 1992; Guesten et al. 1993; Mao et al. 2000; Weiß et al. 2001, 2005; Israel & Baas 2002, 2003; Bradford et al. 2003; Bayet et al. 2004, 2006; Israel 2005, 2009a, 2009b; Israel et al. 2006, 2014, 2015; Papadopoulos et al. 2007, 2010a, 2010b, 2012; Kamenetzky et al. 2011, 2012, 2014, 2016, 2017, 2018; Zhang et al. 2014; Liu et al. 2015, hereafter L15; Daddi et al. 2015; Saito et al. 2017), while observations of rotational transition lines of high dipole moment molecules like hydrogen cyanide (HCN) reveal the densest phase of the gas (${n}_{{{\rm{H}}}_{2}}\gtrsim {10}^{3\mbox{--}4}\,{\mathrm{cm}}^{-3};$ e.g., Downes et al. 1992; Brouillet & Schilke 1993; Gao & Solomon 2004a, 2004b; Papadopoulos 2007; Shirley 2015).

In turbulent star formation theory, variations of molecular gas properties are naturally created by turbulence, which is ubiquitous in galaxies (e.g., Nordlund & Padoan 1999; Ostriker et al. 1999; Padoan & Nordlund 2002, 2011; Krumholz & McKee 2005; Krumholz & Thompson 2007; Feldmann et al. 2011; Hennebelle & Chabrier 2011; Padoan et al. 2012; Salim et al. 2015; Leroy et al. 2017; Elmegreen 2018). Turbulence generates certain gas density probability distribution functions (PDFs). At each gas density, CO molecules have different excitation conditions. By solving radiative transfer equations with the large velocity gradient (LVG) assumption (e.g., Goldreich & Kwan 1974), CO line fluxes can be calculated for each given state of gas volume density, column density, CO abundance, LVG velocity gradient, etc. The integrated CO line fluxes from all gas states give the total CO spectral line energy distribution (SLED) as observed. Therefore, CO SLED could be a powerful tracer of turbulence and of molecular gas properties.

Meanwhile, dust grains are also important ingredients of the interstellar medium (ISM), mixed with gas. They are exposed to and heated by the interstellar radiation field (ISRF), and their thermal emission dominates the (far-)infrared/(sub)millimeter part of galaxies' spectral energy distributions (SEDs). Like molecular gas, dust grains do not physically have a single state. Although observational studies sometimes approximate galaxies' dust SEDs by one or two components in modified-blackbody fitting, physical models based on assuming PDFs for the ISRF have been proposed and calculated by Dale et al. (2001), Dale & Helou (2002), Li & Draine (2002), and (Draine & Li 2007, hereafter DL07). See also subsequent applications in Draine et al. (2007, 2014), Aniano et al. (2012, 2020), Magdis et al. (2012), Daddi et al. (2015), Dale et al. (2017), and Schreiber et al. (2018).

Through the study of both CO excitation and dust SED traced mean ISRF intensity ($\left\langle U\right\rangle$) in about 20 galaxies, Daddi et al. (2015) found that the CO (5–4)/(2–1) line ratio, R₅₂, is tightly correlated with $\left\langle U\right\rangle$. This indicates that CO excitation, or its related ISM properties, is indeed sensitive to the ISRF. However, how the underlying gas density and temperature correlate with ISRF, as well as how this relates to other known correlations like the Kennicutt–Schmidt law, is still unclear.

In this work, we study the CO excitation, molecular gas density, and ISRF in a large sample of 76 (unlensed) galaxies from local to high redshift. The sample is selected from a large compilation of local and high-redshift CO observations from the literature, where we require galaxies to have both CO (2–1) and CO (5–4) detections, together with well-sampled dust SEDs. This also includes CO (5–4) observations newly presented here, from the Institute de Radioastronomie Millimétrique (IRAM) Plateau de Bure Interferometer (PdBI; now upgraded to the NOrthern Extended Millimeter Array [NOEMA]) for six starburst (SB) type galaxies at z ∼ 1.6 in the COSMOS field, which have Atacama Large Millimeter/Submillimeter Array (ALMA) CO (2–1) from Silverman et al. (2015a).

To estimate gas density and temperature from observed line ratios, we model gas density PDFs following Leroy et al. (2017) but with a new approach incorporating assumptions based on the observed correlations between the gas volume density, column density, and velocity dispersion. We propose a conversion method from the line ratio to the mean molecular hydrogen gas density $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and kinetic temperature T_kin for galaxies at global scale. ¹⁷ Our model-predicated J_u < 10 CO SLEDs also show good agreement with the current data.

The structure of this paper is as follows. Section 2 describes the sample and data. Section 3 describes the SED fitting technique for $\left\langle U\right\rangle$ and other galaxy properties. In Section 4, we present correlations between R₅₂ and various galaxy properties. Then, in Section 5, we describe details of our gas modeling and the conversion from R₅₂ to $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin, while the resulting correlations between $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$, T_kin, and $\left\langle U\right\rangle$ are presented in Section 6.2. We discuss the physical meaning of $\left\langle U\right\rangle$, the connection from the $\left\langle U\right\rangle$–and T_kin–$\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ correlations to the Kennicutt–Schmidt law, and the limitations and outlook of our study in Section 6. Finally, we summarize in Section 7.

Throughout this paper, line ratios for CO are expressed as flux-to-flux ratio, where fluxes are in units of Jy km s⁻¹. We adopt a flat ΛCDM cosmology with H₀ = 73 km s⁻¹ Mpc⁻¹, Ω_M = 0.27, and a Chabrier (2003) initial mass function (IMF).

We search the literature for CO observations of local and high-redshift galaxies and seek galaxies that have multiple CO line detections. This is not a complete search, but we have included 132 papers presenting CO observations from 1975 to 2020. ¹⁸ We require a galaxy to have one low-J CO line, CO (2–1), and one mid/high-J CO line, CO (5–4), for this work. This approach is chosen to maximize the sample size while covering most high-redshift main-sequence (MS) ¹⁹ galaxies' CO observations.

We also require multiwavelength coverage including optical, near-IR, far-infrared, and (sub)millimeter, in order to fit their panchromatic SEDs and obtain stellar and dust properties.

In this way, we build up a sample of 76 galaxies. They are divided into the following subsamples:

• 1.  
  22 "local (U)LIRG": local (ultra)luminous infrared galaxies with IR luminosity L_IR ≥ 10¹¹ L_⊙. Their high-J CO data are from the HerCULES (Rosenberg et al. 2015) and GOALS (Armus et al. 2009; Lu et al. 2014, 2015, 2017) surveys using the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) Fourier Transform Spectrometer (FTS; Naylor et al. 2010) on board the Herschel Space Observatory (Pilbratt et al. 2010) and analyzed by L15. Their low-J CO data are from ground-based observations in the literature (see references in Table 1).
• 2.  
  16 "local SFG": local star-forming galaxies, most of which have high-J CO from the KINGFISH (Kennicutt et al. 2011) and VNGS (PI: C. Wilson) surveys using Herschel SPIRE FTS, also analyzed by L15. Many of them have low-J CO mapping from the ground-based HERACLES survey (Leroy et al. 2009), while others have CO (2–1) single-pointing observations in the literature.
• 3.  
  6 "high-z SB FMOS": redshift z ∼ 1.5 SB ²⁰ galaxies from the FMOS-COSMOS survey (Silverman et al. 2015b), where CO (2–1) data are from Silverman et al. (2015a) and CO (5–4) data are newly presented in this work.
• 4.  
  4 "high-z MS BzK": z ∼ 1.5 MS galaxies from Daddi et al. (2008, 2010a, 2015), selected using BzK color criterion (Daddi et al. 2004) and representing high-redshift massive star-forming disks.
• 5.  
  4 "high-z SB SMG": z ∼ 2–6 starbursty, (sub)millimeter-selected galaxies, including GN20 (Daddi et al. 2009; Carilli et al. 2010), AzTEC-3 (Riechers et al. 2010), COSBO-11 (Aravena et al. 2008), and HFLS3 (Riechers et al. 2013).
• 6.  
  8 "high-z MS V20": high-redshift MS galaxies from Valentino et al. (2020a), with SFR within 4× the MS SFR.
• 7.  
  16 "high-z SB V20": high-redshift SB galaxies from Valentino et al. (2020a), with SFR greater than 4× the MS SFR.

Table 1. Sample of Galaxies Used in This Work with Measured and Derived Physical Properties

Source	Subsample	z	R52	$\mathrm{log}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$	$\left\langle U\right\rangle$	$\mathrm{log}{L}_{\mathrm{IR}}$	$\mathrm{log}{M}_{\star }$	Reference CO54	Reference CO21
Arp 193	local (U)LIRG	0.023	1.8 ± 0.3	2.4 ± 0.4	${17.0}_{+0.0}^{+2.4}$	${11.6}_{-0.0}^{+0.0}$	${10.3}_{-0.0}^{+0.0}$	L15	P14
Arp 220	local (U)LIRG	0.018	2.4 ± 0.5	2.7 ± 0.7	${20.6}_{-0.1}^{+0.4}$	${12.2}_{-0.0}^{+0.0}$	${11.0}_{-0.0}^{+0.0}$	L15	K14/G09
IRAS F17207–0014	local (U)LIRG	0.043	3.9 ± 1.0	3.5 ± 1.3	${35.0}_{-0.4}^{+0.2}$	${12.4}_{-0.0}^{+0.0}$	${11.0}_{-0.0}^{+0.1}$	L15	K14/B08/W08/P12
IRAS F18293–3413	local (U)LIRG	0.018	2.2 ± 0.1	2.6 ± 0.4	${11.1}_{-2.9}^{+0.9}$	${11.7}_{-0.1}^{+0.0}$	${9.6}_{-0.2}^{+0.2}$	L15	G93
M82	local SFG	0.001	1.7 ± 0.2	2.4 ± 0.3	${25.4}_{-0.5}^{+0.2}$	${10.6}_{-0.0}^{+0.0}$	${9.8}_{-0.0}^{+0.0}$	L15	L09
Mrk 231	local (U)LIRG	0.042	3.0 ± 0.9	2.9 ± 1.3	${50.0}_{-2.5}^{+0.0}$	${12.3}_{-0.0}^{+0.0}$	${10.9}_{-0.0}^{+0.0}$	L15	K14/P12/A07
Mrk 273	local (U)LIRG	0.038	3.4 ± 0.9	3.2 ± 1.3	${37.7}_{+0.0}^{+5.6}$	${12.0}_{-0.0}^{+0.0}$	${10.6}_{-0.2}^{+0.0}$	L15	K14/P12
NGC 0253	local SFG	0.001	2.9 ± 0.7	3.0 ± 1.0	${6.3}_{-0.5}^{+0.0}$	${10.6}_{-0.0}^{+0.0}$	${10.9}_{-0.1}^{+0.0}$	L15	K14(43.5)/H99
NGC 0828	local (U)LIRG	0.018	0.9 ± 0.4	2.1 ± 0.4	${3.5}_{-0.0}^{+0.5}$	${11.3}_{-0.0}^{+0.0}$	${11.2}_{-0.0}^{+0.0}$	L15	P12(22)
NGC 1068	local (U)LIRG	0.004	0.8 ± 0.2	2.1 ± 0.2	${5.8}_{-0.7}^{+0.2}$	${11.2}_{-0.0}^{+0.0}$	${10.5}_{-0.0}^{+0.0}$	L15	K14(43.5)/K11/B08
NGC 1266	local SFG	0.00724	2.6 ± 0.7	2.8 ± 1.0	${13.3}_{-2.0}^{+3.3}$	${10.3}_{-0.0}^{+0.1}$	${10.5}_{-0.0}^{+0.0}$	L15	K14(43.5)/A11/Y11
NGC 1365	local (U)LIRG	0.005	1.5 ± 0.4	2.3 ± 0.5	${3.5}_{-0.0}^{+0.6}$	${11.2}_{-0.0}^{+0.0}$	${10.9}_{-0.0}^{+0.0}$	L15	K14(43.5)/S95
NGC 1614	local (U)LIRG	0.016	1.4 ± 0.3	2.3 ± 0.3	${31.2}_{+0.0}^{+28.6}$	${11.6}_{-0.0}^{+0.1}$	${10.5}_{-0.2}^{+0.0}$	L15	A95(22)
NGC 2369	local (U)LIRG	0.011	2.0 ± 0.4	2.5 ± 0.5	${5.8}_{-1.9}^{+0.2}$	${11.1}_{-0.1}^{+0.0}$	${10.5}_{-0.0}^{+0.0}$	L15	A95(22)/B08
NGC 2623	local (U)LIRG	0.018	4.1 ± 0.8	3.6 ± 1.1	${19.1}_{-2.0}^{+4.3}$	${11.4}_{-0.0}^{+0.0}$	${10.3}_{-0.0}^{+0.0}$	L15	P12/W08
NGC 2798	local SFG	0.00576	1.9 ± 0.3	2.5 ± 0.4	${13.3}_{-2.6}^{+2.8}$	${10.5}_{-0.0}^{+0.0}$	${9.7}_{-0.0}^{+0.0}$	L15	L09
NGC 3256	local (U)LIRG	0.009	1.7 ± 0.2	2.4 ± 0.3	${23.9}_{-9.5}^{+0.5}$	${11.6}_{-0.1}^{+0.0}$	${10.4}_{-0.0}^{+0.0}$	L15	A95(24)/B08/G93
NGC 3351	local SFG	0.0026	0.9 ± 0.2	2.1 ± 0.2	${2.0}_{-0.3}^{+0.4}$	${9.8}_{-0.0}^{+0.0}$	${9.8}_{-0.0}^{+0.0}$	L15	L09
NGC 3627	local SFG	0.00243	1.8 ± 0.4	2.4 ± 0.6	${3.6}_{-1.0}^{+0.8}$	${10.3}_{-0.0}^{+0.0}$	${10.1}_{-0.0}^{+0.2}$	L15	L09
NGC 4321	local SFG	0.00524	0.8 ± 0.1	2.1 ± 0.2	${1.8}_{+0.0}^{+0.7}$	${10.4}_{-0.0}^{+0.0}$	${10.5}_{-0.0}^{+0.0}$	L15	L09
NGC 4536	local SFG	0.00603	1.7 ± 0.3	2.4 ± 0.4	${3.9}_{-0.1}^{+0.6}$	${10.3}_{-0.0}^{+0.0}$	${10.1}_{-0.0}^{+0.0}$	L15	L09
NGC 4569	local SFG	−0.00078	1.1 ± 0.1	2.2 ± 0.2	${1.9}_{-0.3}^{+0.2}$	${9.6}_{-0.0}^{+0.0}$	${10.5}_{-0.0}^{+0.0}$	L15	L09
NGC 4631	local SFG	0.00202	0.9 ± 0.2	2.1 ± 0.3	${2.8}_{-0.6}^{+0.4}$	${10.2}_{-0.0}^{+0.0}$	${9.4}_{-0.0}^{+0.0}$	L15	L09
NGC 4736	local SFG	0.00103	0.6 ± 0.1	2.0 ± 0.1	${4.1}_{-0.4}^{+1.5}$	${9.6}_{-0.0}^{+0.0}$	${9.8}_{-0.0}^{+0.0}$	L15	L09
NGC 4826	local SFG	0.00136	1.6 ± 0.3	2.4 ± 0.4	${3.6}_{-0.6}^{+0.8}$	${9.5}_{-0.0}^{+0.0}$	${10.4}_{-0.0}^{+0.0}$	L15	A95(28)
NGC 4945	local (U)LIRG	0.002	4.0 ± 0.8	3.6 ± 1.2	${7.0}_{-1.0}^{+0.3}$	${11.1}_{-0.0}^{+0.0}$	${9.7}_{-0.0}^{+1.1}$	L15	W04/B08(22)
NGC 5135	local (U)LIRG	0.014	1.6 ± 0.4	2.4 ± 0.4	${8.2}_{-1.2}^{+2.5}$	${11.2}_{-0.1}^{+0.0}$	${11.1}_{-0.7}^{+0.0}$	L15	P12(22)
NGC 5194	local SFG	0.002	0.7 ± 0.1	2.1 ± 0.1	${3.0}_{-0.2}^{+0.7}$	${10.2}_{-0.0}^{+0.0}$	${9.6}_{-0.0}^{+0.0}$	L15	L09
NGC 5713	local SFG	0.00633	1.0 ± 0.2	2.1 ± 0.2	${5.2}_{-1.1}^{+0.6}$	${10.4}_{-0.0}^{+0.0}$	${10.1}_{-0.0}^{+0.0}$	L15	L09
NGC 6240	local (U)LIRG	0.024	2.8 ± 0.8	2.9 ± 0.9	${20.0}_{-0.5}^{+0.2}$	${11.7}_{-0.0}^{+0.0}$	${10.8}_{-0.0}^{+0.0}$	L15	G09
NGC 6946	local SFG	0.00013	1.1 ± 0.1	2.2 ± 0.2	${4.2}_{-1.1}^{+0.4}$	${10.4}_{-0.1}^{+0.0}$	${10.3}_{-0.0}^{+0.3}$	L15	L09
NGC 7469	local (U)LIRG	0.016	1.1 ± 0.3	2.1 ± 0.2	${13.1}_{+0.0}^{+5.1}$	${11.6}_{-0.0}^{+0.0}$	${10.0}_{-0.0}^{+0.3}$	L15	P12
NGC 7552	local (U)LIRG	0.005	2.4 ± 0.5	2.7 ± 0.7	${14.0}_{-0.5}^{+0.2}$	${11.1}_{-0.0}^{+0.0}$	${10.2}_{-0.0}^{+0.0}$	L15	A95
NGC 7582	local SFG	0.005	1.5 ± 0.3	2.3 ± 0.4	${11.7}_{-0.3}^{+0.2}$	${10.9}_{-0.0}^{+0.0}$	${10.9}_{-0.0}^{+0.0}$	L15	A95
MCG +12-02-001	local (U)LIRG	0.016	1.6 ± 0.3	2.3 ± 0.3	${17.4}_{-2.1}^{+3.6}$	${11.5}_{-0.0}^{+0.0}$	${11.9}_{-1.1}^{+0.6}$	L15	K16(43.5)
Mrk 331	local (U)LIRG	0.018	2.4 ± 0.4	2.7 ± 0.7	${14.7}_{-2.1}^{+0.4}$	${11.4}_{-0.0}^{+0.0}$	${10.9}_{-1.5}^{+0.0}$	L15	K16(43.5)
NGC 7771	local (U)LIRG	0.014	1.3 ± 0.2	2.2 ± 0.3	${7.0}_{-0.5}^{+0.1}$	${11.3}_{-0.0}^{+0.0}$	${11.4}_{-0.0}^{+0.0}$	L15	K16(43.5)
IC 1623	local (U)LIRG	0.02	1.8 ± 0.3	2.4 ± 0.5	${13.2}_{-2.1}^{+0.5}$	${11.6}_{-0.0}^{+0.0}$	${9.1}_{-0.0}^{+0.0}$	L15	K16(43.5)
BzK 16000	high-z MS BzK	1.52	1.5 ± 0.3	2.3 ± 0.4	${15.2}_{-13.3}^{+33.2}$	${11.8}_{-0.0}^{+0.0}$	${11.0}_{-0.2}^{+0.0}$	D15	D15/M12
BzK 17999	high-z MS BzK	1.41	2.2 ± 0.2	2.6 ± 0.4	${14.4}_{-7.7}^{+13.7}$	${12.0}_{-0.0}^{+0.0}$	${10.7}_{-0.0}^{+0.2}$	D15	D15/M12
BzK 21000	high-z MS BzK	1.52	2.3 ± 0.2	2.6 ± 0.4	${25.2}_{-12.2}^{+5.2}$	${12.3}_{-0.0}^{+0.0}$	${11.0}_{-0.2}^{+0.1}$	D15	D15/M12
BzK 4171	high-z MS BzK	1.47	1.8 ± 0.2	2.4 ± 0.4	${16.5}_{-8.1}^{+4.5}$	${12.0}_{-0.0}^{+0.0}$	${10.7}_{-0.1}^{+0.1}$	D15	D15/M12
GN 20	high-z SB SMG	4.06	3.4 ± 0.4	3.2 ± 0.7	${35.4}_{-8.8}^{+4.9}$	${13.3}_{-0.1}^{+0.0}$	${11.2}_{-0.1}^{+0.0}$	C10	D09
AzTEC-3	high-z SB SMG	5.3	4.0 ± 0.2	3.6 ± 0.5	${120.2}_{-84.9}^{+8.9}$	${13.3}_{-0.1}^{+0.0}$	${10.8}_{-0.0}^{+0.2}$	R10	R10
COSBO-11	high-z SB SMG	1.83	3.7 ± 0.1	3.3 ± 0.5	${20.0}_{-3.0}^{+0.4}$	${12.9}_{-0.0}^{+0.0}$	${10.8}_{-0.0}^{+0.0}$	A08	A08
HFLS3	high-z SB SMG	6.34	5.9 ± 0.3	5.0 ± 0.4	${68.3}_{-0.0}^{+10.5}$	${13.7}_{-0.0}^{+0.1}$	${10.5}_{-0.4}^{+0.5}$	R13	R13
PACS-819	high-z SB FMOS	1.45	3.5 ± 0.2	3.2 ± 0.5	${27.7}_{-0.1}^{+5.6}$	${12.5}_{-0.0}^{+0.1}$	${10.7}_{-0.1}^{+0.1}$	THIS	S15
PACS-830	high-z SB FMOS	1.46	1.6 ± 0.2	2.4 ± 0.4	${24.3}_{-2.3}^{+6.7}$	${12.4}_{-0.0}^{+0.0}$	${11.0}_{-0.0}^{+0.0}$	THIS	S15
PACS-867	high-z SB FMOS	1.57	1.6 ± 0.3	2.4 ± 0.4	${2.8}_{-2.3}^{+18.6}$	${12.0}_{-0.0}^{+0.0}$	${10.8}_{-0.1}^{+0.1}$	THIS	S15
PACS-299	high-z SB FMOS	1.65	2.6 ± 0.2	2.7 ± 0.4	${28.3}_{-19.1}^{+38.5}$	${12.4}_{-0.0}^{+0.0}$	${10.1}_{-0.0}^{+0.4}$	THIS	S15
PACS-325	high-z SB FMOS	1.65	0.0 ± 3.4	⋯	${1.2}_{-0.7}^{+14.6}$	${11.8}_{-0.1}^{+0.1}$	${10.4}_{-0.1}^{+0.0}$	THIS	S15
PACS-164	high-z SB FMOS	1.65	1.9 ± 0.4	2.4 ± 0.7	${18.1}_{-15.4}^{+35.0}$	${12.5}_{-0.0}^{+0.0}$	${10.2}_{-0.2}^{+0.3}$	THIS	S15
V20-ID41458	high-z SB V20	1.29	1.8 ± 0.2	2.4 ± 0.3	${33.5}_{+0.0}^{+11.0}$	${12.5}_{-0.0}^{+0.0}$	${11.1}_{-0.0}^{+0.0}$	V20	V20
V20-ID21060	high-z SB V20	1.28	3.6 ± 0.9	3.4 ± 1.3	${51.5}_{-1.5}^{+0.5}$	${12.3}_{-0.0}^{+0.0}$	${10.0}_{-0.0}^{+0.1}$	V20	V20
V20-ID51599	high-z SB V20	1.17	2.1 ± 0.2	2.5 ± 0.4	${14.4}_{-1.9}^{+3.7}$	${12.5}_{-0.0}^{+0.0}$	${11.1}_{-0.0}^{+0.1}$	V20	V20
V20-ID30694	high-z MS V20	1.16	1.2 ± 0.2	2.2 ± 0.3	${15.0}_{-2.9}^{+5.5}$	${12.0}_{-0.0}^{+0.1}$	${10.9}_{-0.0}^{+0.2}$	V20	V20
V20-ID38053	high-z SB V20	1.15	1.3 ± 0.4	2.3 ± 0.4	${18.9}_{-0.1}^{+11.6}$	${12.0}_{-0.0}^{+0.0}$	${10.5}_{-0.0}^{+0.0}$	V20	V20
V20-ID48881	high-z SB V20	1.16	1.9 ± 0.5	2.5 ± 0.5	${42.9}_{-0.2}^{+0.5}$	${12.3}_{-0.0}^{+0.0}$	${10.6}_{-0.0}^{+0.0}$	V20	V20
V20-ID37250	high-z SB V20	1.15	1.1 ± 0.1	2.2 ± 0.2	${9.9}_{-1.1}^{+3.7}$	${12.2}_{-0.0}^{+0.0}$	${11.0}_{-0.2}^{+0.0}$	V20	V20
V20-ID44641	high-z MS V20	1.15	1.0 ± 0.3	2.2 ± 0.5	${9.4}_{-2.8}^{+2.4}$	${12.0}_{-0.1}^{+0.0}$	${11.2}_{-0.3}^{+0.0}$	V20	V20
V20-ID51936	high-z SB V20	1.4	1.9 ± 0.3	2.4 ± 0.3	${5.3}_{-0.1}^{+2.1}$	${12.0}_{-0.0}^{+0.0}$	${10.5}_{-0.0}^{+0.0}$	V20	V20
V20-ID31880	high-z SB V20	1.4	2.2 ± 0.4	2.6 ± 0.5	${20.5}_{+0.0}^{+2.8}$	${12.3}_{-0.0}^{+0.0}$	${11.0}_{-0.0}^{+0.0}$	V20	V20
V20-ID2299	high-z SB V20	1.39	3.4 ± 0.3	3.2 ± 0.5	${13.8}_{-0.4}^{+0.4}$	${12.7}_{-0.0}^{+0.0}$	${11.1}_{-0.1}^{+0.0}$	V20	V20
V20-ID21820	high-z MS V20	1.38	2.1 ± 0.4	2.5 ± 0.5	${15.7}_{-2.3}^{+7.6}$	${12.2}_{-0.0}^{+0.0}$	${11.0}_{-0.1}^{+0.0}$	V20	V20
V20-ID13205	high-z SB V20	1.27	3.0 ± 0.7	3.1 ± 1.2	${49.8}_{-10.8}^{+16.6}$	${12.3}_{-0.0}^{+0.0}$	${11.1}_{-0.2}^{+0.0}$	V20	V20
V20-ID13854	high-z MS V20	1.27	1.8 ± 0.3	2.5 ± 0.5	${20.0}_{-3.0}^{+0.4}$	${12.2}_{-0.0}^{+0.0}$	${11.1}_{-0.0}^{+0.0}$	V20	V20
V20-ID19021	high-z SB V20	1.26	1.9 ± 0.3	2.4 ± 0.4	${25.0}_{+0.0}^{+5.4}$	${12.3}_{-0.0}^{+0.0}$	${10.4}_{-0.0}^{+0.0}$	V20	V20
V20-ID35349	high-z MS V20	1.26	0.8 ± 0.2	2.1 ± 0.2	${8.2}_{-0.6}^{+4.0}$	${12.0}_{-0.0}^{+0.0}$	${11.2}_{-0.1}^{+0.0}$	V20	V20
V20-ID42925	high-z SB V20	1.6	2.1 ± 0.4	2.5 ± 0.5	${59.9}_{-18.5}^{+0.4}$	${12.7}_{-0.0}^{+0.0}$	${11.0}_{-0.0}^{+0.0}$	V20	V20
V20-ID38986	high-z MS V20	1.61	2.8 ± 0.9	2.8 ± 1.4	${19.5}_{-16.1}^{+155.1}$	${12.0}_{-0.1}^{+0.0}$	${11.1}_{-0.0}^{+0.0}$	V20	V20
V20-ID30122	high-z MS V20	1.46	2.0 ± 0.4	2.6 ± 0.6	${13.4}_{-4.2}^{+1.3}$	${12.2}_{-0.0}^{+0.0}$	${10.9}_{-0.0}^{+0.1}$	V20	V20
V20-ID41210	high-z SB V20	1.31	2.1 ± 0.2	2.5 ± 0.4	${25.0}_{-9.6}^{+0.4}$	${12.3}_{-0.1}^{+0.0}$	${10.6}_{-0.0}^{+0.0}$	V20	V20
V20-ID2993	high-z SB V20	1.19	1.4 ± 0.3	2.3 ± 0.4	${13.1}_{-3.0}^{+7.4}$	${12.2}_{-0.0}^{+0.0}$	${11.0}_{-0.2}^{+0.1}$	V20	V20
V20-ID48136	high-z MS V20	1.18	1.5 ± 0.2	2.3 ± 0.3	${14.9}_{-3.3}^{+3.0}$	${12.3}_{-0.1}^{+0.0}$	${11.1}_{-0.0}^{+0.1}$	V20	V20
V20-ID51650	high-z SB V20	1.34	2.8 ± 0.4	2.9 ± 0.6	${21.1}_{-5.0}^{+9.4}$	${12.2}_{-0.0}^{+0.1}$	${10.9}_{-0.0}^{+0.0}$	V20	V20
V20-ID15069	high-z SB V20	1.21	1.7 ± 0.6	2.4 ± 1.0	${6.1}_{-1.2}^{+2.3}$	${12.0}_{-0.0}^{+0.0}$	${10.8}_{-0.3}^{+0.1}$	V20	V20

Note. Only a few selected key columns are shown here. The full sample table has more columns including galaxy properties of DL07 warm- and cold-dust luminosities, AGN luminosities, and offset from the MS, which are used in Figure 2. The full machine-readable table is available at 10.5281/zenodo.3958271.

References THIS = this work (see Appendix A); L15 = Liu et al. (2015); P14 = Papadopoulos et al. (2014); K14 = Kamenetzky et al. (2014); G09 = Greve et al. (2009); E90 = Eckart et al. (1990); I14 = Israel et al. (2014); B08 = Baan et al. (2008); W08 = Wilson et al. (2008); P12 = Papadopoulos et al. (2012); G93 = Garay et al. (1993); L09 = Leroy et al. (2009); B06 = Bayet et al. (2006); A07 = Albrecht et al. (2007); H99 = Harrison et al. (1999); B92 = Braine & Combes (1992); K11 = Kamenetzky et al. (2011); A11 = Alatalo et al. (2011); Y11 = Young et al. (2011); S95 = Sandqvist et al. (1995); A95 = Aalto et al. (1995); W04 = Wang et al. (2004); K16 = Kamenetzky et al. (2016); D15 = Daddi et al. (2015); M12 = Magnelli et al. (2012); D09 = Daddi et al. (2009); W12 = Walter et al. (2012); C10 = Carilli et al. (2010);R10 = Riechers et al. (2010); A08 = Aravena et al. (2008); R13 = Riechers et al. (2013); S15 = Silverman et al. (2015a); V20 = Valentino et al. (2020a).

Download table as:  ASCIITypeset images: 1 2

Our sample is shown in Table 1, where references for the CO (2–1) and CO (5–4) observations are provided. We note that there are also additional interesting galaxies observed in these CO lines: for example, strongly lensed galaxies (e.g., Yang et al. 2017; Harrington et al. 2018), or galaxies that have observations of different CO lines (e.g., Boogaard et al. 2019). As the sample we compiled in this work already covers a large variety of galaxy types (e.g., MS/SB, local/high-redshift), we chose not to further include these data for simplicity and consistency. Applying our method to an extended sample of galaxies could be the subject of a future study.

In the following, we present more details about the CO and multiwavelength photometry data for subsamples.

2.1. Local (U)LIRGs and SFGs

2.2. High-z SB FMOS Galaxies with New PdBI Observations

2.3. High-z MS BzK Galaxies

2.4. High-z SB SMGs

2.5. High-z MS and SB Galaxies from V20

The well-sampled SEDs from optical to far-IR/millimeter allow us to obtain accurate dust properties by fitting them with SED templates. Particularly, since dust grains do not have a single temperature in a galaxy, the mean ISRF intensity, $\left\langle U\right\rangle$, has been considered to be a more physical proxy of dust emission properties (DL07). It represents the 0–13.6 eV intensity of interstellar UV radiation in units of the Mathis et al. (1983) ISRF intensity (see Draine et al. 2007).

The $\left\langle U\right\rangle$ parameter has advantages in describing mixture states of ISRF over using a single or several dust temperature values to describe galaxy dust SEDs. In DL07 dust models, the majority of dust grains are exposed to a minimum ambient ISRF with intensity ${U}_{\min }$, while the rest are exposed to the photon-dominated region (PDR) ISRF, with intensities ranging from ${U}_{\min }$ to ${U}_{\max }$ in a power-law PDF (in mass). The mass fraction of the latter dust grain population ("warm dust" or "PDR dust") is expressed as f_PDR in this work and is a free parameter in the fit. ${U}_{\min }$ is another free parameter, while ${U}_{\max }$ is empirically fixed, as well as the power-law index (see more detailed introduction in Draine et al. 2007, 2014; Aniano et al. 2012, 2020). As pointed out by Dale & Helou (2002), such a physically driven dust model actually fits the mass distribution of molecular clouds (Stutzki 2001; Shirley et al. 2002; Elmegreen 2002). Based on this model, DL07 generated SED templates that can then be used for fitting by other works using their own SED fitting code.

In this work, we use our own-developed SED fitting code, MiChi2, ²⁶ providing us the flexibility in combining multiple SED components and choosing SED templates for each component. Comparing with popular panchromatic (UV-to-millimeter/radio) SED fitting codes, e.g., MAGPHYS (da Cunha et al. 2008, 2015), LePhare (Arnouts et al. 1999; Ilbert et al. 2006), and CIGALE (Noll et al. 2009; Ciesla et al. 2015; Boquien et al. 2019), our code fits SEDs well and produces similar best-fitting results (see Appendix C). Our code also performs χ²-based posterior probability distribution analysis and estimates reasonable (asymmetric) uncertainties for each free or derived parameter (e.g., Figure 1).

Zoom In Zoom Out Reset image size

Our code can also handle an arbitrary number of SED libraries as the components of the whole SED. For example, we use five SED libraries/components representing stellar, AGN, DL07 warm dust, DL07 cold dust, and radio emissions (see below). Our code samples their combinations in the five-dimensional space, then generates a composite SED (after multiplying the model with the filter curves), and then fits to the observed photometric data and obtains χ² statistics. The post-processing of the χ² distribution provides the best-fit and probability range of each physical parameter in the SED libraries (following Press et al. 1992, chapter 15.6).

Details of the five SED libraries/components are as follows:

• 1.  
  Stellar component: for high-redshift (z > 1) SFGs, we use the Bruzual & Charlot (2003) code to generate solar metallicity, constant star formation history (SFH), and Chabrier (2003) IMF SED templates and then apply the Calzetti et al. (2000) attenuation law with a range of E(B − V) = 0.0−1.0 to construct our SED library. For local galaxies, we use the FSPS (Conroy et al. 2009; Conroy & Gunn 2010a, 2010b) code to generate solar -metallicity, τ-declining SFH, Chabrier (2003) IMF SED templates (also with the Calzetti et al. 2000 attenuation law), as this generates a larger variety of SED templates that fit local galaxies better.
• 2.  
  Mid-IR AGN component: we use the observationally calibrated AGN torus SED templates from Mullaney et al. (2011). They cover 6–100 μm in wavelengths and can fit type 1, type 2, and intermediate-type AGNs as demonstrated by Mullaney et al. (2011).
• 3.  
  DL07 warm dust component for dust grains exposed to the PDR ISRF with intensity ranging from ${U}_{\min }$ to ${U}_{\max }={10}^{7}$ in a power-law PDF with an index of −2 (updated version; see Draine et al. 2014; Aniano et al. 2020). The fraction of dust mass in polycyclic aromatic hydrocarbons (PAHs) is described by q_PAH. The contribution of such warm dust to total ISM dust in mass is described by f_PDR in this work (i.e., the γ in DL07). Free parameters are ${U}_{\min }$, q_PAH, and f_PDR.
• 4.  
  DL07 cold-dust component for dust grains exposed to the ambient ISRF with intensity of ${U}_{\min }$. The ${U}_{\min }$ and q_PAH of the cold dust are fixed to be the same as the warm dust in our fitting.
• 5.  
  Radio component: a simple power-law with index −0.8 is assumed. Our code has the option to fix the normalization of the radio component at rest-frame 1.4 GHz to the total IR luminosity L_{IR(8–1000 μm)} (integrating warm- and cold-dust components only) via assumptions about the IR−radio correlation (e.g., Condon et al. 1991; Yun et al. 2001; Ivison et al. 2010; Magnelli et al. 2015) when galaxies lack sufficient IR photometric data and display no obvious radio excess due to AGNs (e.g., Liu et al. 2018). As radio is not the focus of this work, we only use the simple power-law assumption for illustration purposes.

Note that we do not balance the dust-attenuated stellar light with the total dust emission. This has the advantage of allowing for optically thick dust emission that is only seen in the infrared. Our fitting then outputs χ² distributions for the following parameters of interest (see bottom panels in Figure 1):

• 1.  
  Stellar properties, including stellar mass M_⋆, dust attenuation E(B − V), and light-weighted stellar age.
• 2.  
  AGN luminosity L_AGN, integrated over the AGN SED component.
• 3.  
  IR luminosities for cold dust (L_{IR,cold dust}), warm dust (L_{IR,PDR dust}) and their sum (L_{IR,total dust}).
• 4.  
  Mean ISRF intensity $\left\langle U\right\rangle$, minimum ISRF intensity ${U}_{\min }$, and the mass fraction of warm/PDR-like dust in the DL07 model f_PDR.

In Figure 1 we show two examples of our SED fitting. Best-fit parameters and their errors are also listed in our full sample table (Table 1 online version).

To verify our SED fitting, we also fit our high-z galaxies' SEDs with MAGPHYS and CIGALE (see more details in Appendix C). We find that for most high-z galaxies the stellar masses and IR luminosities agree within ∼0.2–0.3 dex. The IR luminosities are more consistent than stellar masses among the results of three fitting codes, with a scatter of ∼0.2 dex. In several outlier cases, our code produces more reasonable fitting to the data (e.g., AzTEC-3, Arp220, NGC 0253), which is likely because we do not have an energy balance constraint in the code. Our code has no systematic bias against CIGALE, but there is a noticeable trend that MAGPHYS fits slightly larger stellar masses than the other two. A possible reason is the use of the Charlot & Fall (2000) double attenuation law in MAGPHYS (see Lo Faro et al. 2017) rather than the Calzetti et al. (2000) attenuation law in our MiChi2 and CIGALE fitting.

Given the general agreement between our code and CIGALE/MAGPHYS, and to be consistent within this paper, we fit all SEDs with our MiChi2 SED fitting code with the five SED libraries as mentioned above.

We use our SED fitting results and the compiled CO data to study the empirical correlation between the CO (5–4)/CO (2–1) line ratio R₅₂ and the mean ISRF intensity $\left\langle U\right\rangle$. This correlation physically links molecular gas and dust properties together, supporting the idea that gas and dust are generally mixed together at large scales and exposed to the same local ISRF.

In Figure 2 we correlate R₅₂ with various galaxy properties derived from our SED fitting. Panel (a) shows a tight correlation between R₅₂ and the ambient ISRF intensity ${U}_{\min }$, and panel (b) confirms the tight correlation between R₅₂ and $\left\langle U\right\rangle$ that was first reported by Daddi et al. (2015). Panels (c) and (d) show that CO excitation is also well correlated with galaxies' dust luminosities, but not with their stellar masses. In panels (e)−(g), we show that R₅₂ exhibits no correlation with f_PDR and mid-IR AGN fraction, while a very weak correlation seems to exist between R₅₂ and the offset to the MS SFR, SFR/SFR_MS. In each panel, the Pearson correlation coefficient P is computed and shown in the lower right corner. These correlations, or lack thereof, demonstrate that R₅₂ or mid-J CO excitation is indeed mostly driven by dust-related quantities, i.e., L_IR, $\left\langle U\right\rangle$, and ${U}_{\min }$.

Zoom In Zoom Out Reset image size

Our best-fitting R₅₂–$\left\langle U\right\rangle$ correlation is close to the one found by Daddi et al. (2015), yet somewhat shallower than that. Valentino et al. (2020a) also reported a shallower slope of the R₅₂–$\left\langle U\right\rangle$ correlation, given that the high-z V20 sample is used in both their and this work. Indeed, subsamples behave slightly differently in Figure 2. While local SFGs and local (U)LIRGs are scattered well around the average R₅₂–$\left\langle U\right\rangle$ correlation line, high-z MS and SB galaxies from the FMOS and V20 subsamples tend to lie below it. Given the varied S/N of IR data as reflected by the $\left\langle U\right\rangle$ error bars, the majority of those high-z galaxies do not have a high-quality constraint on $\left\langle U\right\rangle$. High-z sample selections for CO observations are usually also biased to high-z IR-bright galaxies. Therefore, it is difficult to draw a conclusion about any redshift evolution of the R₅₂–$\left\langle U\right\rangle$ correlation with the current data set.

From panel (f) of Figure 2, we can see that there are several galaxies showing a high AGN-to-ISM dust luminosity ratio (note that the AGN luminosity is integrated over all wavelengths, while the IR luminosity is only DL07 warm+cold dust integrated over 8–1000 μm). The three galaxies with L_AGN,allλ /L_{IR,8–1000μm} ≳ 0.9 are V20-ID38986, V20-ID51936, and V20-ID19021, from high to low, respectively. They all clearly show power-law shape SEDs from the near-IR IRAC bands to mid-IR MIPS 24 μm and PACS 100 μm. ²⁷ However, their R₅₂ do not tend to be higher. This likely supports that these mid-J (J_u ∼ 5) CO lines are not overwhelmingly affected by AGNs.

We note that the correlations in Figure 2 are not the only ones worth exploring. R₅₂ also correlates with dust mass in a way similar to $\left\langle U\right\rangle$ but with larger scatter, and $\left\langle U\right\rangle$ can be considered as the ratio of L_IR/M_dust; therefore, here we omit the correlation with M_dust. Daddi et al. (2015) also investigated how SFR surface density (Σ_SFR), star formation efficiency (SFR/M_gas), gas-to-dust ratio (δ_GDR), and massive star-forming clumps affect $\left\langle U\right\rangle$ and R₅₂. Their results support the idea that a larger fraction of massive star-forming clumps with denser molecular gas compared to the diffuse, low-density molecular gas is the key for a high CO excitation (as proposed by the simulation work of Bournaud et al. 2015). Therefore, to understand the key physical drivers of CO excitation, information on molecular gas density distributions is likely the most urgently required.

CO line emission in galaxies arises mainly from the cold molecular gas, and CO line ratios/SLEDs are sensitive to local molecular gas physical conditions, i.e., volume density ${n}_{{{\rm{H}}}_{2}}$, column density ${N}_{{{\rm{H}}}_{2}}$, and kinetic temperature T_kin. These properties typically vary by one to three orders of magnitude within a galaxy, e.g., as seen in observations as reviewed by Young & Scoville (1991), Solomon & Vanden Bout (2005), Carilli & Walter (2013), and Combes (2018) and references therein, and also in modeling and simulations, e.g., by Krumholz & Thompson (2007), Glover & Clark (2012), Smith et al. (2014b), Narayanan & Krumholz (2014), Bournaud et al. (2015), Glover et al. (2015), Glover & Smith (2016), Popping et al. (2016, 2019), Renaud et al. (2019b, 2019a), and Tress et al. (2020).

In practice, studies of the CO SLED at a global galaxy scale or at subkiloparsec scales usually require the presence of a relatively dense gas component (${n}_{{{\rm{H}}}_{2}}\sim {10}^{3\mbox{--}5}\,{\mathrm{cm}}^{-3};$ T_kin ≳ 50–100 K) in addition to a relatively diffuse gas (${n}_{{{\rm{H}}}_{2}}\sim {10}^{2\mbox{--}3}\,{\mathrm{cm}}^{-3};$ T_kin ∼ 20–100 K), via non-local thermodynamic equilibrium (non-LTE) LVG radiative transfer modeling (e.g., Israel et al. 1995; Mao et al. 2000; Israel & Baas 2001, 2002, 2003; Weiß et al. 2001, 2005; Bradford et al. 2003; Zhu et al. 2003; Bayet et al. 2004, 2006, 2009; Papadopoulos et al. 2007, 2008, 2010a, 2010b, 2012; Hailey-Dunsheath et al. 2008, 2012; Panuzzo et al. 2010; Rangwala et al. 2011; Papadopoulos et al. 2012; Spinoglio et al. 2012; Meijerink et al. 2013; Pereira-Santaella et al. 2013; Rigopoulou et al. 2013; Greve et al. 2014; Kamenetzky et al. 2014, 2016, 2017; Lu et al. 2014, 2017; Rosenberg et al. 2014a, 2014b; Schirm et al. 2014, 2017; Zhang et al. 2014; Israel et al. 2015; Liu et al. 2015; Mashian et al. 2015; Wu et al. 2015; Yang et al. 2017; Valentino et al. 2020a). A third state that is mostly responsible for J_u ≳ 10 CO lines is also found in the case of AGNs (e.g., van der Werf et al. 2010; Rangwala et al. 2011; Spinoglio et al. 2012) or mechanical heating (e.g., Rosenberg et al. 2014a). Therefore, a mid-J-to-low-J CO line ratio like R₅₂ reflects not only the excitation condition of a single gas state but also the relative amount of the denser and warmer to the more diffuse gas component.

Leroy et al. (2017) have conducted pioneer modeling of the sub-beam gas density PDF to understand line ratios of CO isotopologue and dense gas tracers. The method includes constructing a series of one-zone clouds, performing non-LTE LVG calculation, and compositing line fluxes by the gas density PDF. They demonstrated that such modeling can successfully reproduce observed isotopologue or dense gas tracers to CO line ratios. Inspired by this work, we present in this section similar sub-beam density–PDF gas modeling to study the CO excitation and propose a useful conversion from R₅₂ observations to $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin for galaxies at global scales.

5.1. Observational Evidences of Gas Density PDF

5.2. Sub-beam Gas Density PDF Modeling

We thus assume that the line-of-sight volume density of molecular gas in a galaxy follows a lognormal PDF, with a small portion of lines of sight following a power-law PDF at the high-density tail. Representative PDFs are shown in Figure 3. Each PDF samples the ${n}_{{{\rm{H}}}_{2}}$ from 1 to 10⁷ cm⁻³ in 100 bins in logarithmic space. For each ${n}_{{{\rm{H}}}_{2}}$ bin, the height of the PDF is thus proportional to the number of sight lines with a density of ${n}_{{{\rm{H}}}_{2}}$. We assume that the CO line emission surface brightness from each line of sight can be computed from an equivalent "one-zone" cloud with a single ${n}_{{{\rm{H}}}_{2}}$, ${N}_{{{\rm{H}}}_{2}}$, T_kin, velocity gradient, and CO abundance. Thus, the total CO line emission surface brightness is the sum of all sight lines in the PDF.

Zoom In Zoom Out Reset image size

The shape of the gas density PDF is described by the following parameters: the mean gas density of the lognormal PDF $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$, the threshold density of the power-law tail ${n}_{{{\rm{H}}}_{2},\mathrm{thresh}}$, the width of the lognormal PDF, and the slope of the power-law tail. We model a series of PDFs by varying the $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ from 10^2.0 to 10^5.0 cm⁻³ in steps of 0.25 dex and ${n}_{{{\rm{H}}}_{2},\mathrm{thresh}}$ from 10^4.0 to 10^5.25 cm⁻³ in steps of 0.25 dex, to build our model grid, which can cover most situations observed in galaxies. The slope of the power-law tail is fixed to −1.5, which is an intermediate value as indicated by simulations (Federrath & Klessen 2013), also previously adopted by Leroy et al. (2017). The width of the lognormal PDF is physically characterized by the Mach number of the supersonic turbulent ISM (see Padoan & Nordlund 2002, 2011 and Equation (5) of Leroy et al. 2017): $\sigma \approx 0.43\sqrt{\mathrm{ln}(1+0.25{{ \mathcal M }}^{2})}$, which ranges typically from 4 to 20 in star-forming regions as shown by simulations (e.g., Krumholz & McKee 2005; Padoan & Nordlund 2011). Here we adopt a fiducial Mach number of 10, as done previously by Leroy et al. (2017). Note that a high Mach number of ∼80 is also found in merger systems and SB galaxies (e.g., Leroy et al. 2016). It corresponds to a lognormal PDF width 1.56× our fiducial value and marginally affects the CO excitation in a similar way to a higher $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$. Thus, for simplicity in this work we fix the Mach number and allow $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ to vary.

5.3. One-zone Gas Cloud Calculation

For a given gas density PDF, each ${n}_{{{\rm{H}}}_{2}}$ bin is composed of the same "one-zone" gas clouds for which we will compute the line surface brightness. A one-zone cloud has a single volume density ${n}_{{{\rm{H}}}_{2}}$, column density ${N}_{{{\rm{H}}}_{2}}$, gas kinetic temperature (T_kin), CO abundance [CO/H₂], and velocity gradient dv/dr. Note that although an equivalent cloud size r is implied from the ratio of ${N}_{{{\rm{H}}}_{2}}$ and ${n}_{{{\rm{H}}}_{2}}$, given that the calculation is in 1D, r should not be taken as a physical cloud size. Also note that in our study we do not model the 3D distribution of one-zone models; therefore, any radiative coupling between one-zone models along the same line of sight cannot be accounted for. This is likely a minor issue for star-forming disk galaxies given their thin disks (a few hundred parsecs; Wilson et al. 2019) and systematic rotation that separates molecular clouds in the velocity space for inclined disks, but the actual effects need to be studied by detailed numerical simulations (e.g., Smith et al. 2020; Tress et al. 2020).

Here we use RADEX (van der Tak et al. 2007) to compute the 1D non-LTE radiative transfer. For a given ${n}_{{{\rm{H}}}_{2}}$, we loop ${N}_{{{\rm{H}}}_{2}}$ from 10²¹ to 10²⁴ cm⁻², and r is then determined by

$$\begin{eqnarray}&&{N}_{{{\rm{H}}}_{2}}=2\times r\times {n}_{{{\rm{H}}}_{2}}=6\times {10}^{18}\times \displaystyle \frac{r}{\mathrm{pc}}\times \displaystyle \frac{{n}_{{{\rm{H}}}_{2}}}{{\mathrm{cm}}^{-3}}\ \ ({\mathrm{cm}}^{-2}).\end{eqnarray} \tag{ 1 }$$

We also loop over T_kin values of 25, 50, and 100 K, while we fix [CO/H₂] = 5 × 10⁻⁵, a reasonable guess for star-forming clouds (e.g., Leung & Liszt 1976; van Dishoeck & Black 1987, 1988), although it varies from cloud to cloud and depends on chemistry (e.g., Sheffer et al. 2008). Note that there is one additional free parameter to set, i.e., the LVG velocity gradient dv/dr, or the line width FWHM ΔV, or the velocity dispersion σ_V . They are related to each other by

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}V & = & 2\times r\times {dv}/{dr}\ \ (\mathrm{km}\,{{\rm{s}}}^{-1})\\ {\sigma }_{V} & = & {\rm{\Delta }}V/(2\sqrt{2\mathrm{ln}2})\ \ (\mathrm{km}\,{{\rm{s}}}^{-1}).\end{array}\end{eqnarray} \tag{ 2 }$$

To determine these quantities and effectively reduce the number of free parameters while being consistent with observations, we use an empirical correlation between ${N}_{{{\rm{H}}}_{2}}$, r, σ_V , and the virial parameter α_vir. α_vir describes the ratio of a cloud's kinetic energy and gravitational potential energy (e.g., Bertoldi & McKee 1992) and can be written as $\tfrac{5{\sigma }_{V}^{2}r}{{fGM}}$, where σ_V and r are introduced above, G is the gravitational constant, M is the cloud mass, and f is a factor to account for the lack of balance between kinetic and gravitational potential (see Equation (6) of Sun et al. 2018). Observations show that clouds are not always virialized, i.e., α_vir is not always unity. Based on ∼60 pc CO mapping of 11 galaxies in the PHANGS-ALMA sample, Sun et al. (2018) reported the following correlation in their Equation (13) (helium and other heavy elements are included; see also Equation (2) in the review by Heyer & Dame 2015):

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{\mathrm{vir}} & = & 5.77\times {\left(\displaystyle \frac{{\sigma }_{V}}{\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{2}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}}{{{\rm{M}}}_{\odot }{\mathrm{pc}}^{-2}}\right)}^{-1}{\left(\displaystyle \frac{r}{40\mathrm{pc}}\right)}^{-1}\\ & = & \ 5.77\times {\left(\displaystyle \frac{{\sigma }_{V}}{\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{2}{\left(\displaystyle \frac{{N}_{{{\rm{H}}}_{2}}}{1.55\times {10}^{20}{\mathrm{cm}}^{-2}}\right)}^{-1}{\left(\displaystyle \frac{r}{\mathrm{pc}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 3 }$$

They find α_vir ≈ 1.5–3.0 with a 1σ width of 0.4–0.65 dex. For simplicity and also with the idea of focusing primarily on the effect of gas density, we adopt a constant α_vir of 2.3. As shown in later sections, this is already sufficient to explain the observed CO line ratios/SLEDs by our modeling. But note that more comprehensive descriptions of α_vir can be achieved in simulations and can be compared with the results from this work to better understand how a changing α_vir could affect CO line ratio prediction.

Figure 4 presents how R₅₂ changes with the gas densities of one-zone cloud models for four representative redshifts where the cosmic microwave background (CMB) temperatures are different. We repeat our calculations for three representative T_kin as labeled in each panel. The comparison shows that T_kin significantly affects the R₅₂ line ratio, especially at low densities and at low redshifts. Note that due to the constant α_vir assumption, for a given ${n}_{{{\rm{H}}}_{2}}$, Equation (3) implies that σ_V ∝ r and that dv/dr is not varying with ${N}_{{{\rm{H}}}_{2}}$. Thus, the actual choices of ${N}_{{{\rm{H}}}_{2}}$ (or r) for each single one-zone model will not affect the modeling of R₅₂ (and of the optical depth τ).

Zoom In Zoom Out Reset image size

In addition, our modeling is also able to produce reasonable line optical depths (τ) and [C i]/CO line ratios, as presented in Appendix D.

5.4. Converting R₅₂ to $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin with the Model Grid

We compute the global line surface brightness by summing one-zone line surface brightnesses at each ${n}_{{{\rm{H}}}_{2}}$ bin according to the gas density PDF. With our assumptions, there are only four free parameters: the mean gas density of the lognormal PDF $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$, the threshold density of the power-law tail ${n}_{{{\rm{H}}}_{2},\mathrm{thresh}}$, gas kinetic temperature T_kin, and redshift. Their grids are described in Section 5.2.

In Figure 5, we present the predicted R₅₂ as a function of the four free parameters. R₅₂ increases smoothly with $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin, while ${n}_{{{\rm{H}}}_{2},\mathrm{thresh}}$ does not substantially alter the R₅₂ ratio, as indicated by the color-coding. The minimum R₅₂ at the lowest density (${\mathrm{log}}_{10}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle /{\mathrm{cm}}^{-3}\sim 2$) is nearly doubled from redshift 0 to 6 owing to the increasing CMB temperature, but such a redshift effect is less prominent (<×1.5) both at higher density (${\mathrm{log}}_{10}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle /{\mathrm{cm}}^{-3}\gt 3$) and for higher T_kin.

Zoom In Zoom Out Reset image size

In Figure 6, we further show the full CO SLEDs at J_u = 1–9 from our model grid and compare them with a subsample of galaxies with multiple CO transitions at various redshifts. These galaxies are displayed in panels where the $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ is closest to their R₅₂-derived $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ (see below). Our modeling can generally match these CO SLEDs given certain choices of $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin. Yet we caution that this is not a thorough comparison, and our model grid might not fit entirely well the CO SLED shape owing to our simplifying assumptions of fixed Mach number and power-law tail slope or α_vir. While this work only focuses on R₅₂ with the simplest assumptions, the model predictions seem overall already quite promising for the whole CO SLEDs and can be further improved in future works.

Zoom In Zoom Out Reset image size

Based on the model grid, we describe below a method to determine the most probable $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$, T_kin, and ${n}_{{{\rm{H}}}_{2},\mathrm{thresh}}$ ranges for a given R₅₂ and its error in galaxies with known redshift. This is done with a Monte Carlo approach. We first interpolate our 4D model grid to the exact redshift of each galaxy using Python scipy.interpolate.LinearNDInterpolator and then resample the 3D model grid to a finer grid, perturb the R₅₂ given its error over a normal distribution for 300 realizations, and find the minimum χ² best fits for each realization. Finally, we combine best fits to obtain posterior distributions of $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$, T_kin, and ${n}_{{{\rm{H}}}_{2},\mathrm{thresh}}$ and determine their median, 16th percentile (L68), and 84th percentile (H68). This fitting method is coded in our Python package co-excitation-gas-modeling that is made publicly available.

We note that although there is a single input observation (R₅₂) whereas there are three parameters to be determined ($\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$, T_kin, and ${n}_{{{\rm{H}}}_{2},\mathrm{thresh}}$), our method still produces reasonable results. In fact, our method is able to take into account the internal degeneracy between $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin inside model grids, thus obtaining reasonable probability ranges. Figure 7 shows the fitted $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin for our galaxy sample, resulting in a nonlinear trend between $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin. The galaxy-wide mean pressures of gas can also be calculated as ${T}_{\mathrm{kin}}\times \left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and are found to agree with estimates in local galaxies (Kamenetzky et al. 2014).

Zoom In Zoom Out Reset image size

In Figures 8 and 9, we present correlations between the R₅₂-fitted $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin, respectively, and various galaxy properties, similarly to what is presented in Figure 2 for R₅₂. We discuss them in detail in the next sections (Section 6.2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1. The Underlying Meaning of $\left\langle U\right\rangle$: A Mass-to-light Ratio for Dust

By definition, $\left\langle U\right\rangle$ is the mass-weighted ISRF intensity created by UV photons from stars in a galaxy. As indicated by the DL07 model and many of its applications, e.g., Draine et al. (2007, 2014), Aniano et al. (2012, 2020), Dale et al. (2012, 2017), Magdis et al. (2012, 2017), Ciesla et al. (2014) and Schreiber et al. (2018), $\left\langle U\right\rangle$ is actually a mass-weighted, average mass-to-light ratio for the mixture of dust grains in a galaxy. It is driven by the young stars emitting most of the UV photons, but it also reflects the mean distance between young stars and interstellar dust and the efficiency of UV photons heating the dust. For a given DL07 ISRF distribution power-law index (=−2) and ${U}_{\max }$ (=10⁷; Draine et al. 2014), $\left\langle U\right\rangle$ is proportional to the ratio between L_IR and M_dust, with a coefficient P₀ ≈ 138 from this work, where P₀ represents the power absorbed per unit dust mass in a radiation field U = 1:

$$\begin{eqnarray}&&\begin{array}{l}{L}_{\mathrm{IR},8\mbox{--}1000\mu {\rm{m}}}={P}_{0}\cdot \left\langle U\right\rangle \cdot {M}_{\mathrm{dust}}\\ {\rm{where}}\ {P}_{0}\approx 120\mbox{--}150\ ({\rm{mean}}=138).\end{array}\end{eqnarray} \tag{ 4 }$$

Note that the P₀ factor is calibrated to be equal to 125 in Magdis et al. (2012) owing to a slightly different ${U}_{\max }={10}^{6}$, a small 10% systematic difference.

$\left\langle U\right\rangle$ is also positively linked to dust temperature, but it depends on how dust temperature is defined. For example, Draine et al. (2007) find that T ≈ 17 · U^1/6 (K) for dust grains with sizes greater than 0.03 μm whose blackbody radiation peaks around 160 μm. Schreiber et al. (2018) calibrate the light-weighted dust temperature ${T}_{\mathrm{dust}}^{\mathrm{light}}=20.0\cdot {U}^{1/5.57}\ ({\rm{K}})$ (and mass-weighted ${T}_{\mathrm{dust}}^{\mathrm{mass}}=0.91\cdot {T}_{\mathrm{dust}}^{\mathrm{light}}$) by fitting Wien's law to each elementary Galliano et al. (2011) template.

Studies of T_dust and $\left\langle U\right\rangle$ have shown that dust (ISRF) is warmer (stronger) for increasing IR luminosity from local SFGs to (U)LIRGs (e.g., Hwang et al. 2010; Symeonidis et al. 2013; Herrero-Illana et al. 2019) and increases with redshift for the majority of MS galaxies (e.g., Magdis et al. 2012; Magnelli et al. 2014; Béthermin et al. 2015; Schreiber et al. 2018). Some observations show colder dust temperatures in a few among the most extreme SB systems (e.g., Lisenfeld et al. 2000; Jin et al. 2019; Cortzen et al. 2020). These are likely due to the presence of high dust opacity at shorter wavelengths, which makes the dust SED apparently colder. Observations of SMGs also show colder dust temperatures in some of the less luminous ones. This phenomenon is likely driven by the fact that (sub)millimeter selection favors cold-dust galaxies whose SEDs peak closer to (sub)millimeter wavelengths (e.g., Chapman et al. 2005; Kovács et al. 2006; Symeonidis et al. 2009, 2011; Hwang et al. 2010; Magdis et al. 2010; Magnelli et al. 2010).

There is also an interesting finding that for extreme SB galaxies with SFR/SFR_MS > 4 their $\left\langle U\right\rangle$ seem to not evolve with redshift (e.g., Béthermin et al. 2015), while $\left\langle U\right\rangle$ in MS galaxies does evolve with redshift, and extrapolation. This suggests that $\left\langle U\right\rangle$ in MS galaxies might become stronger than those in extreme SB galaxies at z > 2.5, which seems at odds with the expectation. Yet this finding might also be limited by sample size and selection method, templates used for SED fitting, and the dust optically thin assumption in DL07 templates (e.g., Jin et al. 2019; Cortzen et al. 2020).

Combining with the results from this work, the T_dust or $\left\langle U\right\rangle$ trends are easier to understand when correlating them with molecular gas mean density and temperature. We propose a picture in which the general increase of dust temperature and ISRF is mainly due to the increase in cold molecular gas temperature, due to either higher CMB temperature at higher redshifts or more intense star formation and feedback. While the mean molecular gas density has a weaker, nonlinear trend driving $\left\langle U\right\rangle$ in most galaxies, merger-driven compaction could strongly increase gas density and hence boost $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$, T_kin, and $\left\langle U\right\rangle$ in a small number of SB galaxies. Such an increase in gas density creates more contrast at lower redshifts owing to the general decrease of the cosmic molecular gas density and CMB temperature. This could explain why $\left\langle U\right\rangle$ is more different between MS and SB galaxies at lower redshifts.

6.2. Density- or Temperature-regulated Star Formation? The $\left\langle U\right\rangle$–$\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and $\left\langle U\right\rangle$–T_kin Correlations

6.3. Implication for Star Formation Law Slopes

The star formation law is known as the correlation between gas mass (or surface density) and SFR (or surface density) and can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{SFR} & = & A\cdot {M}_{{{\rm{H}}}_{2}}^{N}\quad \mathrm{or}\\ {{\rm{\Sigma }}}_{\mathrm{SFR}} & = & A\cdot {{\rm{\Sigma }}}_{\mathrm{gas}}^{N},\end{array}\end{eqnarray} \tag{ 5 }$$

where A is the normalization and N is the slope. After the initial idea presented by Schmidt (1959), Kennicutt (1998) first systematically measured the star formation law to be ${{\rm{\Sigma }}}_{\mathrm{SFR}}\propto {{\rm{\Sigma }}}_{\mathrm{gas}}^{1.4\pm 0.15}$ based on observations of nearby spiral and SB galaxies, where Σ_gas is the mass surface density of atomic plus molecular gas, and Σ_SFR is the SFR surface density traced by Hα and/or L_IR. This Kennicutt–Schmidt law with N ≈ 1.4 has been extensively studied in galaxies with Σ_gas ∼ 1–10⁵ M_⊙ pc⁻² and is widely used in numerical simulations (see reviews by Kennicutt & Evans 2012; Carilli & Walter 2013).

However, the actual slope N of the star formation law has been long debated. High-resolution (subkiloparsec-scale) observations in nearby spiral galaxies revealed that atomic gas does not correlate with SFR, whereas only molecular gas traces SFR, and N is close to unity in these galaxies (e.g., Wong & Blitz 2002; Bigiel et al. 2008; Leroy et al. 2008, 2013; Schruba et al. 2011). Meanwhile, from local SFGs to (U)LIRGs, observations suggest that N is superlinear, ranging from ∼1 to ∼2 (e.g., Kennicutt 1998; Yao et al. 2003; Gao & Solomon 2004a; Shetty et al. 2013, 2014b, 2014a; de los Reyes & Kennicutt 2019; Wilson et al. 2019). Furthermore, Daddi et al. (2010b) and Genzel et al. (2010) found that high-redshift MS and SB galaxies follow two parallel sequences in the star formation law (${M}_{{{\rm{H}}}_{2}}$–SFR) diagram, each with substantial breadth, and both with N ∼ 1.1–1.2 but with a 0.6 dex mean offset in normalization. Thus, why local SFG regions show a linear star formation law, while high-z SB galaxies have a much higher $\mathrm{SFE}\equiv \mathrm{SFR}/{M}_{{{\rm{H}}}_{2}}$, is still to be understood.

Here we decompose the star formation law into $\left\langle U\right\rangle$ and $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ to gain some insights. First, it is known that the dust-obscured SFR can be traced by the IR luminosity (e.g., Kennicutt 1998; Kennicutt & Evans 2012) as SFR = L_IR/C_IR, where ${C}_{\mathrm{IR}}\sim {10}^{10}\,[{L}_{\odot }\,{({M}_{\odot }{\mathrm{yr}}^{-1})}^{-1}]$ assuming a Chabrier (2003) IMF. Second, as mentioned in the previous section, $\left\langle \langle U\right\rangle \,={P}_{0}^{-1}\cdot {L}_{\mathrm{IR}}/{M}_{\mathrm{dust}}$ . Third, we use the gas-to-dust ratio δ_GDR ≡ M_gas/M_dust to link gas to dust mass. This ratio varies with metallicity (e.g., Israel 1997; Leroy et al. 2007, 2011; Bolatto et al. 2013; Sandstrom et al. 2013; Rémy-Ruyer et al. 2014, 2015), and also note that the definition of gas in δ_GDR is atomic plus molecular gas. We include an additional molecular hydrogen fraction ${f}_{{{\rm{H}}}_{2}}\equiv {M}_{{{\rm{H}}}_{2}}/{M}_{\mathrm{gas}}$ to the gas-to-dust ratio, having finally ${M}_{\mathrm{dust}}={M}_{{{\rm{H}}}_{2}}\cdot {({f}_{{{\rm{H}}}_{2}}{\delta }_{\mathrm{GDR}})}^{-1}$. Fourth, we consider MS galaxies to be disks with radius r and height h and assume that $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ is the global mean gas density; thus, the molecular gas mass can be expressed as the product of the volume and the mean molecular gas density: ${M}_{{{\rm{H}}}_{2}}=\pi \,\cdot \left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle \cdot {r}^{2}\cdot h$ . And fifth, we ignore atomic gas and only consider the molecular gas star formation law.

Then, we rewrite the star formation law equation as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{SFR} & = & A\cdot {M}_{{{\rm{H}}}_{2}}^{N}\\ & \Rightarrow & \ \left\langle U\right\rangle \cdot {M}_{\mathrm{dust}}=A\cdot {P}_{0}^{-1}\cdot {C}_{\mathrm{IR}}\cdot {M}_{{{\rm{H}}}_{2}}^{N}\\ & \Rightarrow & \ \left\langle U\right\rangle \cdot {\left({f}_{{{\rm{H}}}_{2}}{\delta }_{\mathrm{GDR}}\right)}^{-1}=A\cdot {P}_{0}^{-1}\cdot {C}_{\mathrm{IR}}\cdot {M}_{{{\rm{H}}}_{2}}^{N-1}\\ & \Rightarrow & \ \left\langle U\right\rangle \cdot {\left({f}_{{{\rm{H}}}_{2}}{\delta }_{\mathrm{GDR}}\right)}^{-1}\\ & = & \ A\cdot {P}_{0}^{-1}\cdot {C}_{\mathrm{IR}}\cdot {\left(\pi \cdot \left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle \cdot {r}^{2}\cdot h\right)}^{N-1}.\end{array}\end{eqnarray} \tag{ 6 }$$

Taking the logarithm of both sides, and assuming that $\mathrm{log}\left\langle U\right\rangle$, $\mathrm{log}({f}_{{{\rm{H}}}_{2}}{\delta }_{\mathrm{GDR}})$, and $\mathrm{log}({r}^{2}\,h)$ are functions of $\mathrm{log}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$, we have

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}\left\langle U\right\rangle -\mathrm{log}({f}_{{{\rm{H}}}_{2}}{\delta }_{\mathrm{GDR}})\\ \ =\ \mathrm{log}(A\,{P}_{0}^{-1}\,{C}_{\mathrm{IR}})\,+\,(N-1)\left[\mathrm{log}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle +\mathrm{log}(\pi {r}^{2}h)\right]\\ \ \Rightarrow \ N=\displaystyle \frac{\tfrac{d\mathrm{log}\left\langle U\right\rangle }{d\mathrm{log}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle }-\tfrac{d\mathrm{log}({f}_{{{\rm{H}}}_{2}}{\delta }_{\mathrm{GDR}})}{d\mathrm{log}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle }}{1+\tfrac{d\mathrm{log}({r}^{2}h)}{d\mathrm{log}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle }}+1.\end{array}\end{eqnarray} \tag{ 7 }$$

Therefore, the star formation law slope N depends on how $\left\langle U\right\rangle$, ${f}_{{{\rm{H}}}_{2}}{\delta }_{\mathrm{GDR}}$ (metallicity), and r² h (galaxy size) change with $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$, which can further be described by the differentials $\tfrac{d\mathrm{log}\left\langle U\right\rangle }{d\mathrm{log}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle }$, $\tfrac{d\mathrm{log}({f}_{{{\rm{H}}}_{2}}{\delta }_{\mathrm{GDR}})}{d\mathrm{log}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle }$, and $\tfrac{d\mathrm{log}({r}^{2}h)}{d\mathrm{log}\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle }$, respectively. These differentials strongly depend on galaxy samples. When studying subkiloparsec regions in local SFGs, if the ISRF, metallicity, and galaxy size are similar among these SFGs, N is close to 1. When studying a sample including both SFGs and (U)LIRGs, $\left\langle U\right\rangle$ increases by a factor of a few tens with $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ changing from 10² to 10⁴ cm⁻³, and r decreases by a factor of a few with $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$, as (U)LIRGs are usually smaller and more compact (while the scale height h seems constant; e.g., Wilson et al. 2019). As for the ${f}_{{{\rm{H}}}_{2}}{\delta }_{\mathrm{GDR}}$ term, because ${f}_{{{\rm{H}}}_{2}}$ increases with metallicity while δ_GDR decreases with it, their product ${f}_{{{\rm{H}}}_{2}}{\delta }_{\mathrm{GDR}}$ likely does not change much. Therefore, N can be much higher than 1. The overall effect is that the star formation law does not have a single slope, yet the overall N is about 1–2.

6.4. Limitations and Outlook

In this work, we compiled a comprehensive sample of galaxies from local to high redshift with CO (2–1) and CO (5–4) detections and well-sampled IR SEDs. This includes our new IRAM PdBI CO (5–4) observations of six z ∼ 1.5 COSMOS SB galaxies. With this large sample, we measure their mean ISRF intensity $\left\langle U\right\rangle$ from dust SED fitting (Section 3) and their mean molecular gas density $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ converted from R₅₂ = S_CO(5-4)/S_CO(2-1) line ratios based on our density–PDF gas modeling (Section 5). Our results can be summarized as follows.

• 1.  
  We confirm the tight $\left\langle U\right\rangle$–R₅₂ correlation first reported by Daddi et al. (2015) and find that $\left\langle U\right\rangle$, ${U}_{\min }$, and L_IR all strongly correlate with R₅₂, while stellar mass, AGN fraction, and the SFR offset to the MS all show weaker or no correlation with R₅₂ (Figure 2).
• 2.  
  We conduct density–PDF gas modeling to connect the mean molecular gas density $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and kinetic temperature T_kin to the observable CO line ratio R₅₂. Based on this, we provide a Monte Carlo method (and a Python package co-excitation-gas-modeling) to compute $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin's probability ranges using our model grid for any given J_u = 1–10 CO line ratio (and for CO SLED as the next step; see, e.g., Figure 6).
• 3.  
  We find that both $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin increase with $\left\langle U\right\rangle$, with T_kin having a tighter correlation with $\left\langle U\right\rangle$.
• 4.  
  Based on these correlations, we propose a scenario in which the ISRF in the majority of galaxies is more directly regulated by the gas temperature and nonlinearly by the gas density. A fraction of SB galaxies have gas densities larger by more than one order of magnitude with respect to MS galaxies and are possibly in a merger-driven compaction stage (Sections 6.2 and 6.1).
• 5.  
  We link the $\left\langle U\right\rangle$–$\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ correlation to the Kennicutt–Schmidt star formation law and discuss how the star formation law slope N can be inferred from the $\left\langle U\right\rangle$–$\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ correlation slope and other galaxy properties versus $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ correlations. We find that N ∼ 1–2 can be inferred from the trends of how $\left\langle U\right\rangle$ and galaxy size change with $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ in different galaxy samples (Section 6.3).

Our study demonstrates that ISRF and molecular gas are tightly linked to each other, and density–PDF gas modeling is a promising tool for probing detailed ISM physical quantities, i.e., molecular gas density and temperature, from observables like CO line ratios/SLEDs.

We thank the anonymous referee for helpful comments. D.L., E.S., and T.S. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 694343). F.V. acknowledges support from the Carlsberg Foundation Research Grant CF18-0388 "Galaxies: Rise and Death". G.E.M. and F.V. acknowledge the Villum Fonden research grant 13160 "Gas to stars, stars to dust: tracing star formation across cosmic time" and the Cosmic Dawn Center of Excellence funded by the Danish National Research Foundation under then grant No. 140. Y.G.'s research is supported by the National Key Basic Research and Development Program of China (grant No. 2017YFA0402700), National Natural Science Foundation of China (grant Nos. 11861131007, 11420101002), and Chinese Academy of Sciences Key Research Program of Frontier Sciences (grant No. QYZDJSSW-SLH008). S.J. acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) under grant AYA2017-84061-P, co-financed by FEDER (European Regional Development Funds). A.P. gratefully acknowledges financial support from STFC through grants ST/T000244/1 and ST/P000541/1. We thank A. Weiss and C. Wilson for helpful discussions. This work used observations carried out under project No. W14DS with the IRAM Plateau de Bure Interferometer (PdBI). IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). This work used observations carried out under project 17A-233 with the National Radio Astronomy Observatory's Karl G. Jansky Very Large Array (VLA). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Our MiChi2 SED fitting code is publicly available at https://ascl.net/code/v/2533. Our Python package co-excitation-gas-modeling for computing $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ and T_kin from CO line ratios is publicly available at https://pypi.org/project/co-excitation-gas-modeling. And our SED fitting figures as shown in Figure 1 and the full Table 1 are publicly available at 10.5281/zenodo.3958271.

We present the sample table and CO (5–4) imaging of our PdBI observations in Table 2 and Figure 10. The observations are described in Section 2.2.

Zoom In Zoom Out Reset image size

CenA. We excluded this galaxy because its CO (2–1) and CO (5–4) data only cover the center of the galaxy, and significant correction is needed for recovering the entire galaxy. For example, Kamenetzky et al. (2014) applied a correction factor of 1/0.48, where 0.48 is the beam-aperture-to-entire-galaxy fraction denoted as "BeamFrac" and reported in the full Table 1 (available online), to convert the CO (2–1) observed at the galaxy center with a beam of 22'' (Eckart et al. 1990) to a beam of 43'' for their study. They derived this factor based on SPIRE 250 μm image aperture photometry. Based on PACS 70–160 μm (as presented in L15), we obtain correction factors of 1/0.088 and 1/0.185 from a 22'' and 43'' beam to the entire galaxy, respectively. Thus, the 22''-to-43'' correction factors in two works fully agree (0.088/0.185 ≈ 0.48). Despite the good agreement in beam-related correction among these works, we caution that using far-infrared data to correct CO (2–1) is very uncertain, as low-J CO lines do not linearly correlate with far-infrared emission (L15).

M83. We excluded this galaxy in this work as well. Israel & Baas (2001) reported a CO (2–1) line flux of 261 ± 15 K km s⁻¹ (5501 ± 316 Jy km s⁻¹) within a 22'' beam (with SEST 15 m) at the M83 galaxy center, Lundgren et al. (2004) reported 98.1 ± 0.8 K km s⁻¹ (2068 ± 17 Jy km s⁻¹) within a 22'' aperture (with JCMT 15 m) at the same center, and Bayet et al. (2006) reported 67.4 ± 2.2 K km s⁻¹ (2721 ± 88 Jy km s⁻¹) within a 305 beam (with CSO 10.4 m) also at the center position. Kamenetzky et al. (2014) adopted the Bayet et al. (2006) line flux and applied a factor of 1/0.76 correction to obtain the line flux within a 43'' beam. This correction factor agrees with L15. However, if we want to obtain the entire flux for M83, we will need to correct the 43'' flux by a factor of 1/0.232 based on the Herschel PACS aperture photometry in L15. We caution that the uncertainty in such a correction is large, and the CO (2–1) line fluxes at the galaxy center in the literature are already inconsistent by a factor of two.

NGC 0253. At the galaxy center position, the reported CO (2–1) line fluxes are 6637 ± 996 Jy km s⁻¹ within a 12'' beam (Bradford et al. 2003), 10,684 ± 1602 Jy km s⁻¹ within a 15'' beam (Bradford et al. 2003), 17,757 ± 3551 Jy km s⁻¹ within a 21'' beam (Bayet et al. 2004), 24,428 ± 2686 Jy km s⁻¹ within a 23'' beam (Bayet et al. 2004), 33,800 ± 3200 Jy km s⁻¹ within a 435 beam (Kamenetzky et al. 2014), and 34,300 ± 3600 Jy km s⁻¹ within a 435 beam (Kamenetzky et al. 2014; corrected from the original beam in Harrison et al. 1999). The beam-to-entire-galaxy fraction, "BeamFrac," is 0.518 from 435 to the entire galaxy based on L15. These fluxes are roughly consistent, and the "BeamFrac"-based correction factor is only a factor of two; thus, we take the last 435 beam flux and obtain 66,216 ± 15010 Jy km s⁻¹ as the CO (2–1) flux for the entire NGC 0253, where we added a 0.2 dex uncertainty to the 435 beam flux error. We caution that, even with the additional uncertainty, the flux error might still underestimate the true uncertainty, which includes original flux calibration and measurement error in Harrison et al. (1999), correction from original beam to 435 by Kamenetzky et al. (2014), and correction from 435 beam to entire galaxy by L15.

NGC 0891. We excluded this galaxy in this work. Braine & Combes (1992) reported a CO (2–1) line flux of 86 ± 6 K km s⁻¹ (1974 ± 138 Jy km s⁻¹) within a convolved 23'' beam at the galaxy center position. Baan et al. (2008) converted the same line brightness temperature from Braine & Combes (1992) to a flux of 381 ± 26 Jy km s⁻¹, which, however, is lower than our converted value in parentheses and is possibly mistaking the original 12'' beam for calculation, while the brightness temperature that Braine & Combes (1992) reported has been convolved to 23'' beam as mentioned in their Table 1 caption. We note that the correction factor from 12'' or 23'' to the entire NGC 0891 is as large as ∼10, e.g., the "BeamFrac" from 169 to entire galaxy is 0.112 as measured by L15. Thus, it is too uncertain to consider this galaxy in this work.

As for CO (1–0), Braine & Combes (1992) reported a line flux of 96 ± 5 K km s⁻¹ (551 ± 28 Jy km s⁻¹) within a 23'' beam at the galaxy center position. This can be corrected to the entire galaxy scale as 3908 ± 204 Jy km s⁻¹ based on L15. Gao & Solomon (2004a, 2004b) reported a line flux of 35.5 ± 5 K km s⁻¹ (963 ± 136 Jy km s⁻¹) within a 50'' beam (with FCRAO 14 m) and a global scale integrated flux of 3733.7 Jy km s⁻¹. They are consistent within errors.

NGC 1068. Braine & Combes (1992) reported a CO (2–1) line flux of 240 ± 10 K km s⁻¹ (5488 ± 229 Jy km s⁻¹) within a convolved 23'' beam at the galaxy center position. Baan et al. (2008) converted the same line brightness temperature from Braine & Combes (1992) to a flux of 1967.2 ± 80 Jy km s⁻¹, which is also inconsistent with our converted value (in parentheses) and possibly due to the mistaking of the original 12'' beam in their calculation. Papadopoulos et al. (2012) reported a CO (2–1) line flux of 11300 ± 2200 Jy km s⁻¹ within the inner 40'' of NGC 1068 (originally from Papadopoulos & Seaquist 1999). Kamenetzky et al. (2011) reported a CO (2–1) line flux of 8366 ± 19 Jy km s⁻¹ within a beam of 30'' (with CSO 10.4 m), which is then corrected to 43''-beam flux of 11,700 ± 1100 Jy km s⁻¹ by Kamenetzky et al. (2014). Kamenetzky et al. (2014) also cited the Baan et al. (2008) flux and reported a 43''-beam flux of 12,600 ± 2500 Jy km s⁻¹ converted from a 12'' beam. But note that Baan et al. (2008) might have mistaken a 12'' beam for the calculation. If we directly correct the Braine & Combes (1992) 23''-beam flux to a 43'' beam, it is 8669 Jy km s⁻¹, which, however, is 30% smaller than that in Kamenetzky et al. (2014). Meanwhile, if we correct the Kamenetzky et al. (2011) flux from 30'' beam to 43'' beam, it is 10,542 Jy km s⁻¹, consistent with both Kamenetzky et al. (2014) and Papadopoulos et al. (2012). Given that the difference is only about 30%, in this work we adopt the average of these fluxes, i.e., 10170 Jy km s⁻¹ for a 43'' beam, or 15,551 Jy km s⁻¹ corrected to the entire galaxy scale (based on L15 BeamFrac).

For CO (1–0), we perform our own photometry using the Nobeyama 45 m COAtlas Survey data (Kuno et al. 2007) and obtain a flux of 5228 Jy km s⁻¹. This is 40% higher than the global scale line flux of 3651.1 Jy km s⁻¹ measured by Gao & Solomon (2004b) using FCRAO mapping observations, but closer to the line flux of 4240 Jy km s⁻¹ within a 43'' beam reported by Kamenetzky et al. (2014), which is citing Baan et al. (2008) and originally also from Gao & Solomon (2004b).

NGC 1365. NGC 1365 were observed at two positions by Herschel SPIRE FTS, one at northeast (NGC 1365-NE) and one at southwest (NGC 1365-SW). They have similar CO (5–4) within 10%, but the IR luminosities within each aperture differ by 25%. This means that our aperture-based beam-to-entire-galaxy correction has at least 25% uncertainty (same in the independent analysis of the similar method by Kamenetzky et al. 2014). For CO (2–1) we use the same Sandqvist et al. (1995) SEST 15 m (24'' beam) data as in Kamenetzky et al. (2014) and correct it to the entire galaxy scale to match our corrected CO (5–4).

NGC 1614. CO (2–1) is from Aalto et al. (1995), observed with SEST 15 m (22'' beam; η_mb = 0.5, ∫T_mb dv = 56 ± 2 K km s⁻¹ or line flux 1180 ± 42 Jy km s⁻¹). We correct from 22'' beam to the entire galaxy with a BeamFrac of 0.792 (L15). Meanwhile, note that Wilson et al. (2008) reported an interferometric integrated CO (2–1) flux of 670 ± 7 Jy km s⁻¹ (synthesized beam 37 × 33). The discrepancy of about 50% is likely due to the missing flux of the interferometry (see Wilson et al. 2008).

NGC 2369. CO (2–1) is from Aalto et al. (1995), observed with SEST 15 m (22'' beam; η_mb = 0.5, ∫T_mb dv = 74 ± 2.4 K km s⁻¹, or line flux 1560 ± 51 Jy km s⁻¹). Meanwhile, note that Baan et al. (2008) reported 959.4 ± 14.3 Jy km s⁻¹, which is originally from Garay et al. (1993) also with SEST 15 m (∫T_mb dv = 46.8 ± 0.7 K km s⁻¹, with η_mb = 0.54). The reason for this factor of two discrepancy is unclear. Here we take their average (1259.7 Jy km s⁻¹) and correct from the 22'' beam to the entire galaxy with a BeamFrac of 0.808 (L15).

NGC 2623. Wilson et al. (2008) reported an interferometric integrated CO (2–1) flux of 267 ± 8 Jy km s⁻¹ observed with SMA. Papadopoulos et al. (2012) cited this flux in their study and discussed that this flux is unlikely affected by missing flux.

NGC 3256. Aalto et al. (1995) reported a CO (2–1) flux of ∫T_mb dv = 314 ± 8 K km s⁻¹ (6619 ± 169 Jy km s⁻¹) observed with SEST 15 m (22'' beam; η_mb = 0.5). Meanwhile, Baan et al. (2008) reported 2980.7 ± 14.3 Jy km s⁻¹, which is originally from Garay et al. (1993), also observed with SEST 15 m (∫T_mb dv = 145.5 ± 0.7 K km s⁻¹, with η_mb = 0.7). Similar to NGC 2369, the reason for the factor of two to three discrepancy is unclear. We take their average (4799.85 Jy km s⁻¹) and correct from the 22'' beam to the entire galaxy with a BeamFrac of 0.744 (L15).

NGC 3351. We obtain CO (2–1) and CO (1–0) line fluxes for the entire galaxy with our own photometry as 2681 and 1138 Jy km s⁻¹, respectively, to the HERACLES data and the Nobeyama 45 m COAtlas Survey (Kuno et al. 2007) data. Uncertainties contributed by the noise in the moment-0 maps are about 6% of the measured fluxes. Note that Braine & Combes (1992) observed a CO (2–1) and CO (1–0) flux of about 642 and 97 Jy km s⁻¹, respectively, convolved to a 23'' beam. Leroy et al. (2009) reported a CO (2–1) luminosity of $0.78\times {10}^{5}\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{arcsec}}^{2}$, or a line flux of 2808 Jy km s⁻¹, for the entire galaxy, consistent with ours. Usero et al. (2015) reported a CO (1–0) flux of about 210 Jy km s⁻¹ within a 213 beam at the central position.

NGC 3627. Similar to NGC 3351, we obtain the CO (2–1) and CO (1–0) line fluxes for the whole galaxy via our photometry using the HERACLES and the NRO45 m COAtlas data, respectively. We measured 9219 and 7366 Jy km s⁻¹, respectively. Note that Gao & Solomon (2004b) reported a global CO (1–0) flux of 4477 Jy km s⁻¹, which is about 40% lower than ours.

NGC 4321. Similar to NGC 3351 and NGC 3627, the global CO (2–1) and CO (1–0) line fluxes are obtained as 9088 and 2251 Jy km s⁻¹, from the HERACLES and the NRO45 m COAtlas data, respectively. Note that Braine & Combes (1992) observed a CO (1–0) flux of 445 Jy km s⁻¹ within a 23'' beam, which can be corrected to a consistent entire galaxy flux of 2280 Jy km s⁻¹ by a BeamFrac of 0.195 (L15), while Komugi et al. (2008) observed a CO (1–0) flux of 174 Jy km s⁻¹ within a 16'' beam, which is somehow lower than others.

NGC 4945. Wang et al. (2004) observed the central position of NGC 4945 with SEST 15 m and obtained a CO (2–1) flux of ∫T_mb dv = 920.9 ± 0.6 K km s⁻¹ (19,412 ± 12.6 Jy km s⁻¹, for point-source response in a 22'' beam). Baan et al. (2008) cited the same Wang et al. (2004) CO (2–1) flux as 18,878.5 ± 12.3 Jy km s⁻¹, which is consistent with our conversion. Curran et al. (2001) also observed the central position of NGC 4945 with SEST 15 m. They reported a CO (2–1) flux of ∫T_mb dv = 740 ± 40 K km s⁻¹, about 20% lower than that of Wang et al. (2004). As discussed in Wang et al. (2004), the reason for the discrepancy is unclear, but this shows that the uncertainty in the CO (2–1) flux at the galaxy center is at least 20%. We take the average (17,505 Jy km s⁻¹) in this work and estimate the entire galaxy CO (2–1) flux to be 31,770 Jy km s⁻¹ based on a BeamFrac of 0.551 (L15) from the 22'' beam.

NGC 6946. The global CO (2–1) and CO (1–0) line fluxes are obtained as 36,296 and 11,454 Jy km s⁻¹, from the HERACLES and the NRO45 m COAtlas data, respectively. This is in good agreement with the global scale CO (1–0) flux of 11,400.5 Jy km s⁻¹ reported by Gao & Solomon (2004b) using NRAO 12 m mapping data.

We performed additional MAGPHYS (da Cunha et al. 2008, 2015; Battisti et al. 2019) and CIGALE (Burgarella et al. 2005; Noll et al. 2009; Ciesla et al. 2014, 2015; Boquien et al. 2019; Yang et al. 2020) SED fitting to verify our MiChi2 SED fitting results. We use the updated MAGPHYS version with high-z extension (http://www.iap.fr/magphys/download.html), and CIGALE version 2020.0 (2020 June 29) (https://cigale.lam.fr/download/). We modified the MAGPHYS FORTRAN source code to allow for longer photometry filter names and larger filter number. MAGPHYS and CIGALE require a list of preset filters, for which we choose the following list: GALEX FUV and NUV, KPNO MOSAIC1 u, CFHT MegaCam u band, SDSS ugriz, Subaru SuprimeCam BVriz, GTC griz, VISTA VIRCAM Y, J, H, K_s , HST ACS F435W/F606W/F755W/F814W and WFC3 F125W/F140W/F160W, Spitzer IRAC ch1/2/3/4, IRS PUI 16 μm and MIPS 24 μm, Herschel PACS 70/100/160 and SPIRE 250/350/500 μm, SCUBA2 450/850 μm, VLA 3/1.4 GHz, and pseudo-880/1100/1200/2000 μm filters. Other photometry data like submillimeter interferometry data (e.g., from ALMA) and some optical data are ignored. Note that in our MiChi2 fitting these bands without a known filter curve are automatically used with a pseudo-delta-function filter curve.

The current MAGPHYS code does not include the fitting of a mid-IR AGN SED component, although such an extension has been used nonpublicly in some studies (Chang et al. 2015). MAGPHYS has preset stellar libraries, dust attenuation laws, and dust libraries; therefore, there is no need to adjust any parameters, except that we run MAGPHYS only for z > 0.03 galaxies, as MAGPHYS computes the luminosity and mass properties with the luminosity distance, which does not match the physical distance at a very low z.

CIGALE has the capability of including a mid-IR AGN component, as does our MiChi2 code. The current version of CIGALE uses AGN emission models computed from physical modeling of AGN torus by Fritz et al. (2006). It has much more freedom than the observationally derived AGN templates by Mullaney et al. (2011) used by MiChi2. However, this can also easily overfit the data when there are only a few broadband photometry data points at mid-IR ∼8–100 μm. For our fitting with CIGALE, we fix several AGN parameters based on the fitting results of SB galaxies in Fritz et al. (2006), r_ratio = 60, beta = −1.0, gamma = 6.0, opening_angle = 140.0, and let the following parameters vary: tau = 1.0, 3.0, 6.0, psy = 0.001, 10.100, 20.100, 30.100, and fracAGN = 0.0, 0.2, 0.4, 0.6.

For the stellar component in our CIGALE fitting of high-redshift galaxies, we use a constant SFH as in our MiChi2 fitting. This is achieved by adding sfhperiodic into the CIGALE sed_modules and setting type_bursts = 2, delta_bursts = 200, tau_bursts = 200. To allow the fitting of a range of stellar ages, we set age = 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000 for the bc03 SED module. And we adopt the Calzetti et al. (2000) dust attenuation law as in our MiChi2 fitting by adding dustatt_modified_starburst to the SED modules and setting E_BV_lines to 0.0 to 2.0 in steps of 0.2. Meanwhile, for local galaxies (z < 0.03) in our sample whose stellar ages are generally older, a constant SFH stellar component cannot fit the stellar SED well. Thus, we adopt the exponentially declining SFH sfh2exp in CIGALE and set tau_main = 200, 500, 1000, 2000, 4000 and age = 200, 500, 1000, 2000, 4000, 6000, 8000, 10000. We turn off the burst model by setting f_burst = 0.0.

The dust template used in CIGALE fitting is also the same as used in our MiChi2 fitting, i.e., the Draine et al. (2014) updated DL07 templates. ${U}_{\min }$ (umin) is set to vary from 1.0 to 50, and f_PDR (gamma) from 0.0 to 1.0. We also set lim_flag = True to allow CIGALE to analyze the photometry upper limits (to achieve this, we need to flip the sign of the flux errors for the photometry with an S/N < 3). Then, each galaxy has about 696,960 models fitted. In comparison, in MAGPHYS in general 13,933 optical models and 24,999 IR dust models are fitted for each galaxy.

In Figure 11, we compare the fitted dust 8–1000 μm luminosities and stellar masses from the three SED fitting codes. In the left panel of Figure 12 we compare the fitted $\left\langle U\right\rangle$ from MiChi2 and CIGALE, as MAGPHYS does not have the same Draine & Li (2007) library. Note that not all fittings show a reasonable χ², as can be seen in the right panel of Figure 12, where the histograms of reduced-χ² are shown for the three fitting codes. CIGALE fittings in general have a higher reduced-χ², which means poorer fitting than MiChi2, whereas MAGPHYS produces slightly better fittings than MiChi2. However, both MAGPHYS and CIGALE have a number of very poor/failed fitting cases that have reduced-χ² ≳ 10 and are the outlier data points in Figures 11. The threshold reduced-χ² ∼ 8–10 is empirically estimated after visually examining the SED fitting results. ²⁹ There are only two sources that exhibit reduced-χ² > 10 in MiChi2 fitting, and their IR-to-millimeter are actually well fitted, leaving the stellar part poorly constrained (CenA and NGC 0253). In comparison, there are 12 poor/failed cases in CIGALE fitting, and 4 in MAGPHYS fitting. The main reason for these poor/failed cases is likely the energy balance forced in MAGPHYS and CIGALE. In these cases, the stellar part of the SED fitting gives a dust attenuation that cannot fully balance the far-IR/millimeter emission, and this is also mentioned in other studies of extremely dust-obscured high-redshift galaxies (e.g., Casey et al. 2017; Miettinen et al. 2017b, 2017a; Simpson et al. 2017). Except for these poor/failed fittings, the fitted IR luminosities and stellar masses reasonably agree within about 0.3 dex.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the comparison of $\left\langle U\right\rangle$ in Figure 12, we excluded the 12 sources with reduced-χ² > 10 in CIGALE fitting. Although most sources have consistent $\left\langle U\right\rangle$, a small number of sources do not have consistent $\left\langle U\right\rangle$, and they mostly come from the V20 subsample. This is mainly because they have very poor IR photometry except for one or two submillimeter interferometry photometries. But we chose to skip these interferometry photometries in our CIGALE fitting tests owing to the filter setting. Adding a fake filter in CIGALE and rerunning the fitting for each of these sources are required in order to fit the submillimeter interferometry photometry, and then more consistent results are expected. Therefore, these comparisons show that for sources with good reduced-χ² and photometry data MiChi2 and CIGALE have similar constraints on $\left\langle U\right\rangle$.

We present the predicted [C i] (³ P₁ − ³ P₀) (hereafter [C i] (1–0)) and CO (1–0) line optical depths and the [C i] (1–0)/CO (1–0) line ratio in surface brightness unit (${R}_{\mathrm{CICO}}^{{\prime} }$) in Figures 13–15. The optical depths shown in Figure 13 agree with normal conditions where CO (1–0) is optically thick, while [C i] (1–0) is roughly optically thin or has τ ∼ 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 14 and 15 show ${R}_{\mathrm{CICO}}^{{\prime} }$ as a function of one-zone cloud molecular gas density and mean molecular gas density of the composite PDF, respectively. Similarly to Figures 4 and 5, we show our prediction at four redshifts, z = 0, 1.5, 4, and 6, and with three representative T_kin = 25, 50, and 100 K. ${R}_{\mathrm{CICO}}^{{\prime} }$ increases with ${n}_{{{\rm{H}}}_{2}}$ or $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$ but strongly decreases with T_kin at intermediate ${n}_{{{\rm{H}}}_{2}}\sim {10}^{3-4}\,{\mathrm{cm}}^{-3}$. Future study of this line ratio with our gas modeling will shed light on how to better constrain T_kin and $\left\langle {n}_{{{\rm{H}}}_{2}}\right\rangle$.