Application of resolved low-$J$ multi-CO line modeling with RADEX to constrain the molecular gas properties in the starburst M82

Assuming that the gas in each sightline is uniform, the emission can be described by single values of ${n_{\mathrm{H}_{2}}},{\mathrm{T}_{\mathrm{kin}}},{N_{\mathrm{CO}}},{{X}_{12/13}},{{X}_{13/18}}$ with the assumption of an identical ${\Phi}$ for all lines. We run RADEX by iterating through the 6 parameters following the ranges and stepsizes outlined in Table 2. The radiative transfer calculations in RADEX depend solely on the ${N_{\mathrm{CO}}}/\Delta v$ ratio, resulting in a degeneracy between ${N_{\mathrm{CO}}}$ and $\Delta v$ (van der Tak et al., 2007). The RADEX calculations allow us to estimate ${N_{\mathrm{CO}}}/\Delta v$, which we then multiply by the observed CO(1–0) line width for each pixel to derive the true ${N_{\mathrm{CO}}}$ and thereby avoid this degeneracy. For the model grid, we assume a fixed $\Delta v=10$ km s^-1 for both models.²²2We note that this choice is arbitrary, since the results are identical upon normalization by the line width. These parameters generate a six-dimensional (6D) grid, amounting to $\sim$12 million points. Our modeling approach and parameter set-up follow those in Teng et al. (2022).

Theoretical and observational studies have shown that molecular clouds have a range of ${\mathrm{H}_{2}}$ densities from diffuse to dense ones, and cannot necessarily be described by a single density as we assumed above (Nishimura et al., 2019). In effort to describe the molecular gas more physically and explore the effects of a more sophisticated model, we implement a log-normal distribution for the density of ${\mathrm{H}_{2}}$. The log-normal distribution assumption follows from the literature, which analyzed high resolution observations of molecular cloud in the CO emission (Hughes et al., 2013). Other literature also employed simulations and found a similar log-normal result (see Mac Low & Klessen (2004); Elmegreen (2002); Vazquez-Semadeni (1994)). Following formulations modeling the density of a cloud, we adopt the probability density function (PDF) method used in Leroy et al. (2017a). See more details in Appendix B.

In short, the final intensities are the weighted sum over RADEX one-zone modeled intensities for different volume densities, but for a fixed temperature, abundance, and column density per unit line width. For this calculation, we sum over H₂ volume densities ranging from 10^-1 to $10^{7}$ cm^-3 in steps of 0.2 dex. Under this model, the final set of parameters governing each line of sight therefore are ${n_{0,\mathrm{H}_{2}}},{\mathrm{T}_{\mathrm{kin}}},{N_{\mathrm{CO}}},{{X}_{12/13}},{{X}_{13/18}}$ and $\sigma$, the width of the log-normal density distribution. We note, that strictly speaking, $n_{0,\rm H_{2}}$ does not correspond to the average volume density introduced for the one-zone model ($n_{\rm H_{2}}$). For the log-normal model, $n_{0,\rm H_{2}}$ more precisely corresponds to the mean volume-weighted density. For modeling purposes, we assume a beam filling factor ${\Phi}$ of 1 for all density layers to limit the degrees of freedom. We emphasize that this is a simplification, and in certain extragalactic environments, there might be indeed regions within a $4^{\prime\prime}$ beam that do not exhibit the presence of any (CO-rich) molecular gas, hence implying the need for a filling factor. However, in the case of the central region of M82, we find a beam filling factor from the one-zone model that is quite large ($\sim$0.5). The exact value will be part of a future study where we include more mm molecular line tracers.³³3A key simplification in our log-normal density model prescription lies in modeling the beam-averaged intensity as a density-weighted sum. In principle, each density layer could have a distinct beam filling factor, which likely depends on the local volume density. A more detailed modeling approach could introduce a global beam filling factor to scale the modeled intensities in equation C1. However, in this study, we refrain from adding this additional free parameter due to the limited number of available lines. While this simplification might introduce a bias, it does not impact the main conclusion: even with a global beam filling factor, the density cannot be reliably constrained above a volume density of $10^{4}\,$cm^-3. The ranges through which we iterate can be found in Table 2. More details about the implementation of this density distribution can be found in Appendix B.

We evaluate the goodness of fit of each parameter set $\theta=({n_{\mathrm{H}_{2}}},{\mathrm{T}_{\mathrm{kin}}},{N_{\mathrm{CO}}},{{X}_{12/13}},{{X}_{13/18}},{\Phi})$ or $\sigma$ (depending on model) by calculating $\chi^{2}$. To remove nonphysical parameter sets from consideration, we set a prior on the line-of-sight path length ($\ell_{\mathrm{los}}$), following Teng et al. (2022), with (where the beam filling factor term is only considered for the one-zone model)

where 500pc is based on the width of M82’s central region.

Assuming a multivariate Gaussian probability distribution, we convert the $\chi^{2}$ for each grid point to a probability. Finally, we generate marginalized 1D and 2D likelihood distributions (called PDFs for 1Ds) for every parameter and pair of parameters. These are generated by summing the probability distribution over the full range of parameter(s) except for those of interest. The resulting “1DMax solution” parameter is then selected by determining the highest 1D likelihood for each parameter. More details can be found in Appendix C.

Our two approaches to modeling emitting gas both yield 1DMax solutions for all 6 parameters over the center of the galaxy. In Figure 3, we plot the binned fits of parameters shared between the models, which include ${{X}_{12/13}},{{X}_{13/18}},{N_{\mathrm{CO}}}$, and ${\mathrm{T}_{\mathrm{kin}}}$. We bin the solutions of the 95 sightlines for which we have significant CO isotopologue emission (excluding ${\mathrm{CO}}$ (3-2)). We see similar behavior when including the ${\mathrm{CO}}$ 3-2 line, but then being limited to only 34 pixels. We make two notes: first, because we consider the effect of ${\Phi}$ in the one-zone model, the 1DMax ${N_{\mathrm{CO}}}$ value constrained is the sub-beam value, i.e the ${N_{\mathrm{CO}}}$ in the smaller beam described by ${\Phi}$. In Figure 3, we convert from the sub-beam value to the full-beam value by multiplying by ${\Phi}$. Second, the ${n_{\mathrm{H}_{2}}}$ values from the two models are not directly comparable because for the log-normal model, ${n_{\mathrm{H}_{2}}}$ is the average density rather than the singular density. It is still instructive, however, to consider their differences.

Figure 3 includes the distribution of the 1DMax solutions. The dotted vertical lines are the medians of the 1DMax solutions over all pixels, and the horizontal markers at the top of the plots indicate the weighted 16^th-to-84^th percentile range. Above the fits, we show values found in the centers of NGC3627, NGC4321, and NGC3351 which are nearby barred galaxies (Teng et al., 2022, 2023), and the M51 galaxy, which is a nearby grand design spiral galaxy (den Brok et al., 2024). The ${{X}_{13/18}}$ parameter was not constrained for in den Brok et al. (2024) and thus is left empty in the figure. We see that the two models generally agree for the medians of the ${{X}_{13/18}}$, ${N_{\mathrm{CO}}}$, and ${\mathrm{T}_{\mathrm{kin}}}$ parameters. We see that the median of the one-zone ${{X}_{12/13}}$ values is higher – but within the margin of the scatter – than that of the log-normal model. This likely reflects the fact that in the case of the one-zone model, sub-beam density variations and the resulting CO isotopologue emission are compensated by higher CO isotopologue abundance. In comparison to values found for other galaxies, we see generally more extreme conditions in the barred galaxies and less extreme conditions in M51.

We performed the Kolmogorov-Smirnov (KS) test to quantitatively evaluate the similarity between the results. The p-values from the KS test are listed in Table 3. For 3 of the 5 parameters shared by the models, the p-value was greater than 0.05, indicating that the difference between the results is not statistically significant. For the two parameters where the p-value was found to be $\ll 0.05$, namely ${\mathrm{T}_{\mathrm{kin}}}$ and ${N_{\mathrm{CO}}}$, the medians were verified to be within 0.1 dex.

For each sightline, we derive the line center optical depth of the 1-0 transition of ${\mathrm{CO}},{{}^{13}\mathrm{CO}},$ and ${\mathrm{C}^{18}\mathrm{O}}$ as follows: for each pixel, we marginalize all modeled optical depths over the whole grid to obtain a 1DMax estimate. We note that in the subsequent section, we work with the optical depths derived from the one-zone model. A similar calculation for the log-normal density model is more ambiguous, since each layer will have a different opacity.

In Figure 4, we bin the 1DMax solutions from each isotopologue. We refer to relative frequency as probability. We find (in the leftmost panel) that the ${\mathrm{CO}}$ emission is optically thick, as $\tau_{12}>1$ for all sightlines. We note that the derived values also align with those found for the centers of other nearby galaxies, including NGC 3351, 3627, and 4321, using a similar modeling approach (Teng et al., 2022, 2023). In contrast, we note that Petitpas & Wilson (2000) report lower optical depths (around 0.3–3 for ¹²CO(1-0)). However, these values correspond to measurements at 22″and also depend on the range of assumptions of kinetic temperatures and the CO isotopologue abundances. We note that their optical depths are also above unity if we select their model results which used values for kinetic temperatures (30 K) and ¹³CO/¹²CO abundance ratio (70) that are similar to what we report.

We find that the ${{}^{13}\mathrm{CO}}$ and ${\mathrm{C}^{18}\mathrm{O}}$ 1-0 emission is exclusively optically thin with $\langle\log_{10}\tau_{13}\rangle=-1.5\pm 0.1$ and $\langle\log_{10}\tau_{18}\rangle=-2.2\pm 0.1$. This aligns well with our expectations of optically thin ${{}^{13}\mathrm{CO}}$ and ${\mathrm{C}^{18}\mathrm{O}}$ emission based on the measured (2-1)/(1-0) excitation ratios. The mean integrated intensity ratio over the 95 sightlines for ¹³CO(2-1)/(1-0) is 1.62 and for C¹⁸O(2-1)/(1-0) is 2.03. Such high ratios are expected in case both the 1-0 and 2-1 transitions remain optically thin and therefore are consistent with our results from the CO line modeling.

As discussed in Section 4.1, we find that the constraints between the one-zone and log-normal model results do not significantly differ for the volume density, ${n_{\mathrm{H}_{2}}}$, of H₂ and the CO isotopologue abundances.

In principle, we expect the log-normal distribution to more closely trace the true physical conditions of the molecular gas since it adopts a range of densities rather than a single value, as in the one-zone model. Inspecting the individual corner plots, however, we see that there is a sharp cut off for the ${n_{\mathrm{H}_{2}}}$ PDF roughly near $10^{4}$ cm^-3, with the values above being weakly or not constrained. Similarly, for the log-normal distribution, while the volume density seems to be constrained around $10^{4}$ cm^-3, we only get a loose constraint on a wide $\sigma$ PDFs. This implies that low $-J$ CO lines are unable to effectively trace higher gas densities that are most likely present in starburst galaxies. This is consistent with expectations based on the critical density and optical depth effects (Shirley, 2015). Indeed, the CO isotopologues used in the study have critical densities on the order of $10^{3}-10^{4}$ cm^-3, which decrease further by roughly another order of magnitude in the optically thick gas.

Finally, we also note that the two model results show similar distributions of derived molecular gas properties (see Figure 3), despite the simplified assumptions employed for both prescriptions. This suggests that the major limitation is inherent to the general low-$J$ CO line ratio modeling approach and only to a lesser degree to simplified assumptions, such as a uniform gas layer (for the one-zone model), or no global beam filling factor for the log-normal density model.

There are several indications that support the presence of high density ($n_{\rm H_{2}}{>}10^{4}$ cm^-3) and warmer ($T_{\rm kin}{>}50\,$K) molecular gas in the central region of M82:(i) the molecular gas in starbursts is observed to show densities that are much higher than those found in typical disk galaxies (Jackson et al., 1995; Iono et al., 2007). (ii) observations also suggest that the molecular gas in starbursts is warmer (Bradford et al., 2003; Rangwala et al., 2011). Such high temperatures and densities are consistent with the high (2-1)/(1-0) line ratios observed for the CO isotopologues in M82 (see Section 4.2).

Studies using different molecular gas tracers than the low-$J$ CO isotopologues indeed report higher density and temperature gas in the central region of M82. For instance, Panuzzo et al. (2010) rely on higher-$J$ CO lines observed with Herschel, including the CO(13-12) and report temperatures in the range of 400–800 K, an order of magnitude higher than the average temperatures we find using the CO line modeling approach. Using the formaldehyde lines as a thermometer, Mühle et al. (2007) also report temperatures of around 200 K for the central region. Previous studies also found densities above 10⁴ cm^-3. For instance, using higher critical density tracers, including HCN, CS, and HCO⁺, Naylor et al. (2010) report average molecular gas densities not in excess of 10^4.4 cm^-3.

Therefore, our results suggest that even with sophisticated modeling approaches incorporating numerous low-$J$ CO isotopologue lines, it remains crucial to include higher excitation energy lines (e.g., higher-$J$ CO lines beyond the 2–1 transition) and high critical density transitions (e.g., multiple HCN transitions) to achieve robust constraints in extreme environments, such as the central region of M82. Indeed, previous Galactic and extragalactic studies focusing on higher-$J$ lines or lines of higher critical density do recover higher volume densities. For instance, studies using Herschel SPIRE/FTS and PACS spectra covering higher-$J$ lines of up to $J_{\rm up}{=}30$, do recover higher temperatures (${>}100\,$K) and volume densities (${\geq}10^{5}$ cm^-3) in local active galaxies (Pereira-Santaella et al., 2013, 2014; Mashian et al., 2015) and nearby starbursts (NGC 6240 and Arp 193; Papadopoulos et al., 2014). Similarly, studies including dense gas tracers, such as the rotational HCN and HCO⁺ lines, also recover and constrain the volume densities and temperatures of a hot (${>}100\,$K) and dense ($>10^{5}$ cm^-3) molecular gas component in extreme conditions found in ULIRGs (e.g., Mashian et al., 2015; Imanishi et al., 2023a, b) or other local starbursts (Krips et al., 2008; Greve et al., 2009). These prior findings reinforce our conclusion that relying on low-$J$ CO isotopologue transitions alone can lead to an incomplete characterization of molecular gas conditions in extreme environments, such as the central region of M82, and underscore the importance of expanding the range of molecular tracers to capture the varied physical conditions in starbursts.

The intense star formation and high molecular gas densities in M82 serve as a compelling nearby analog to the molecular gas conditions in high-redshift galaxies ($z\gtrapprox 1$; e.g., Weiß et al., 2005; Casey et al., 2014). Our findings underscore the need for caution when interpreting low-$J$ CO lines in studies of high-$z$ systems. Analyses relying solely on low-$J$ transitions risk mis-characterizing the density and excitation properties of molecular gas in these environments. To improve the accuracy and diagnostic power of molecular line studies, future work should prioritize the inclusion of multi-line, multi-tracer modeling and higher-density tracers, particularly in starburst galaxies, both locally and at high redshift.