Molecular Gas Properties and CO-to-H₂ Conversion Factors in the Central Kiloparsec of NGC 3351

We obtained ALMA Band 3, 6, and 7 observations to capture CO, ¹³CO, and C¹⁸O lines as well as millimeter continuum emission from the central ∼2 kpc of NGC 3351 (see Table 1). These observations were designed to resolve the gas structures around the starburst ring at 1–2'' (∼50–100 pc) resolution. In this paper, we focus on the molecular line observations and their implied gas conditions. Results of the dust continuum will be included in a future work.

Table 1. ALMA Observations

Line	Project	Rest Frequency	ncrit at 20 K	Array	Beam Size	LAS	rms per Channel
	Code	(GHz)	(cm−3)		('')	('')	(${\rm{K}}\,{(2.5\,\mathrm{km}\,{{\rm{s}}}^{-1})}^{-1}$)
CO J = 1–0	(1)	115.27	2.2 × 103	12 m	2.1	11.9	0.320
CO J = 2–1	(2)	230.54	1.1 × 104	12 m+7 m+TP	1.5	⋯	0.108
13CO J = 2–1	(1)	220.40	9.4 × 103	12 m	1.8	10.4	0.036
13CO J = 3–2	(3)	330.59	3.1 × 104	12 m	1.2	7.6	0.032
C18O J = 2–1	(1)	219.56	9.3 × 103	12 m	1.8	10.4	0.036
C18O J = 3–2	(3)	329.33	3.1 × 104	12 m	1.2	7.6	0.032

Note. ALMA project codes: (1) 2013.1.00885.S, (2) 2015.1.00956.S, (3) 2016.1.00972.S. Critical densities (n_crit) are calculated from the Einstein A and collisional rate coefficients provided by the Leiden Atomic Molecular Database (Schöier et al. 2005). In optically thick cases with τ₀ = 5, n_crit of CO 1–0 and 2–1 at 20 K decrease to 6.1 × 10² and 3.1 × 10³ cm⁻³, respectively (see Shirley 2015). The rms of CO 1–0 is higher due to it being observed in a more extended array configuration, which has a lower surface brightness sensitivity.

Download table as:  ASCIITypeset image

Our Band 3 observation targets the CO 1–0 line and 100–115 GHz continuum. It uses the C43-3 configuration of the ALMA 12 m array to reach an angular resolution of 21 × 13 for the CO line, and a single 12 m pointing covered the entire central region (∼30'' in diameter). The achieved rms noise level is 0.32 K per 2.5 km s⁻¹ channel (averaged across the field of view, out to a primary beam response threshold of 0.25). We note that Band 3 has the lowest frequency, and thus it requires a more extended array configuration to achieve an angular resolution similar to the Band 6 or 7 observations. With a more extended array configuration, the surface brightness sensitivity decreases, which leads to a higher rms level for our Band 3 data compared to the other bands for a given integration time. Although a longer integration time in Band 3 observations can reduce the noise level, the rms of 0.32 K per channel is already sufficient for the secure detection of CO 1–0 over the central kiloparsec region of NGC 3351.

Our Band 6 observation targets the J = 2–1 transition of the ¹³CO and C¹⁸O lines as well as the 215–235 GHz continuum. With the C43-1 configuration of the ALMA 12 m array, it achieves an angular resolution of 18 × 12 and uses a seven pointing mosaic to cover the central 30'' × 30'' area. The achieved rms noise level is 0.036 K per 2.5 km s⁻¹ channel.

Our Band 7 observation targets the J = 3–2 transition of the ¹³CO and C¹⁸O lines and the 325–345 GHz continuum. The ALMA C40-1 configuration provides an angular resolution of 12 × 10 and a 14 pointing mosaic covers the central 30'' × 30'' area. The achieved rms noise level is 0.032 K per 2.5 km s⁻¹ channel.

To complement these observations obtained specifically for this project, we also include the PHANGS–ALMA CO 2–1 data (Leroy et al. 2021a) in our analysis to aid the imaging and provide additional constraints on the gas temperature and isotopologue ratios. This data set was observed in Cycle 3 (project code: 2015.1.00956.S) and it reaches a similar angular resolution as the other observations used in this project. It combines ALMA 12 m, 7 m, and total-power (TP) observations to ensure a complete u− v coverage. The characteristics of the PHANGS–ALMA data are described in detail in Leroy et al. (2021a).

Our Band 3, 6, and 7 observations were taken with the 12 m array alone and thus lack short-spacing information (see Table 1 for the largest angular structure (LAS) recoverable for each line). To estimate how much emission these 12 m only observations can recover, we perform a test with the PHANGS CO 2–1 data, which include coverage of all u− v scales with 12 m and 7 m arrays as well as TP observations. We imaged the 12m only PHANGS data and the 12m+7m+TP and compared the resulting maps to check flux recovery. By dividing the 12 m only image by the combined image, we find a flux recovered ratio of ∼97% over the entire observed region. The analyses presented in this paper include only the pixels within a flux recovered ratio of 1 ± 0.3.

We calibrate the raw data with the calibration scripts provided by the observatory. We use the appropriate versions of the CASA pipeline to calibrate the Cycle 2 and Cycle 5 observations as recommended by the observatory (version 4.2.2 and 5.1.1, respectively).

We then image the CO lines and continuum separately with the CASA task tclean in two steps, as inspired by the PHANGS–ALMA imaging scheme (Leroy et al. 2021b). We weight the u–v data according to the "Briggs" scheme with a robustness parameter r = 0.5, which offers a good compromise between noise and resolution. We first perform a shallow, multiscale cleaning to pick up both compact and extended emission down to a signal-to-noise ratio (S/N) of 3 throughout the entire field of view. After that, we perform a deep, single-scale cleaning to capture all remaining emission down to S/N ∼ 1. This second step is restricted to within a cleaning mask, which is constructed based on the presence of a significant CO 2–1 detection in the PHANGS–ALMA data. This way, the deep cleaning process can efficiently recover most remaining signals at their expected position−position–velocity (ppv) locations.

These calibration and imaging procedures produce high-quality data cubes for the CO lines as well as two-dimensional (2D) images for the continuum emission.

From the CO line data cubes, we derive 2D maps of integrated line intensity (moment 0), intensity-weighted velocity (moment 1), and velocity dispersion estimated using the effective line width estimator, as well as their associated uncertainties. We describe our methodology here, and note that it is very similar to that adopted in the PHANGS–ALMA pipeline (Leroy et al. 2021b, also see Sun et al. 2018).

First, we convolve all data cubes to the finest possible round beam (FWHM = 2.1''; set by the CO 1–0 resolution) to ensure that all CO line data probe the same spatial scale. Second, we estimate the local rms noise in the data cubes by iteratively rejecting CO line detection and calculating the median absolute deviation (MAD) for the signal-free ppv pixels. Third, we create a signal mask for each data cube by finding all ppv positions with an S/N > 5 detection in more than two consecutive channels, and then expand this signal mask to include all morphologically connected ppv positions with an S/N > 2 detection in more than two consecutive channels. Fourth, we combine the signal masks for all CO lines in "union" (i.e., logical OR) to generate a master signal mask. Fifth, we collapse all CO line data cubes within this master signal mask to create moment maps, and calculate the associated uncertainties based on Gaussian error propagation. Finally, we regrid all of the data products such that the new grid Nyquist samples the beam (i.e., the grid spacing equals half of the beam FWHM).

This product creation scheme yields a coherent set of moment maps and uncertainty maps for all of the CO lines, which have matched resolution and consistently cover the same ppv footprint.

We use the non-LTE radiative transfer code RADEX (van der Tak et al. 2007) to model the observed line intensities under various combinations of H₂ density (${n}_{{{\rm{H}}}_{2}}$), kinetic temperature (T_k), CO column density per line width (N_CO/Δv), CO/¹³CO (X_12/13) and ¹³CO/C¹⁸O (X_13/18) abundance ratios, and the beam-filling factor (Φ_bf). The beam-filling factor is a fractional area of the beam covered by the emitting gas, and thus it should be ≤1. Tests of the recoverability of model parameters with RADEX fitting have shown that modeling the isotopologue abundance ratio as a free parameter leads to more accurate results than assuming a fixed ratio (Tunnard & Greve 2016). With six measured lines, we construct two different models with a one-component or two-component assumption on gas phases. The setups for these models are described separately in Sections 4.2 and 4.3.

The input to RADEX includes a molecular data file specifying the energy levels, statistical weights, Einstein A coefficients, and collisional rate coefficients of each specific molecule. The data files we use for CO, ¹³CO, and C¹⁸O are from the Leiden Atomic and Molecular Database (Schöier et al. 2005) with collisional rate coefficients taken from Yang et al. (2010). To run RADEX, we directly input the T_k, ${n}_{{{\rm{H}}}_{2}}$, molecular species column density (N), and line width (Δv). The input T_k and ${n}_{{{\rm{H}}}_{2}}$ are used to set collisional excitation, while N and Δv are the radiative transfer parameters that determine the optical depth and escape probability for the modeled molecular emission line.

RADEX assumes a homogeneous medium and uses radiative transfer equations based on the escape probability formalism (Sobolev 1960) to find a converged solution for the excitation and radiation field. Initially, RADEX guesses the population of each energy level under an LTE assumption and then computes for each transition the resultant optical depth at the line center (τ₀) by

$$\begin{eqnarray}&&{\tau }_{0}=\displaystyle \frac{{A}_{{ij}}}{8\pi {\tilde{\nu }}^{3}}\ \displaystyle \frac{N}{{\rm{\Delta }}v}\left({x}_{j}\displaystyle \frac{{g}_{i}}{{g}_{j}}-{x}_{i}\right)\,,\end{eqnarray} \tag{ 3 }$$

where A_ij is the Einstein A coefficient, $\tilde{\nu }$ is the wavenumber, and x_i and g_i are the fractional population and statistical weight, respectively, in level i. The optical depth then determines the escape probability (β) for a uniform sphere:

$$\begin{eqnarray}&&\beta =\displaystyle \frac{1.5}{\tau }\left[1-\displaystyle \frac{2}{{\tau }^{2}}+\left(\displaystyle \frac{2}{\tau }+\displaystyle \frac{2}{{\tau }^{2}}\right){e}^{-\tau }\right]\,,\end{eqnarray} \tag{ 4 }$$

which estimates the fraction of photons that can escape the cloud (Osterbrock & Ferland 2006). This directly constrains the radiation field, and thus a new estimation of level population and excitation temperature, leading to new optical depth values. This procedure is done iteratively until convergence is reached.

We first run RADEX with CO to build up a three-dimensional (3D) grid with varying N_CO/Δv, ${n}_{{{\rm{H}}}_{2}}$, and T_k. To construct a model grid with X_12/13 and X_13/18 axes, we calculate the column densities of ¹³CO and C¹⁸O that correspond to each varying N_CO/Δv and abundance ratios, and then re-run RADEX with the calculated ${N}_{{13}_{\mathrm{CO}}}/{\rm{\Delta }}v$ and ${N}_{{{\rm{C}}}^{18}{\rm{O}}}/{\rm{\Delta }}v$ values. The RADEX output of the three molecules predicts a Rayleigh–Jeans equivalent temperature (T_RJ) for each line under the varied parameters. Assuming a Gaussian profile, the predicted line fluxes can then be estimated using a Gaussian integral:

$$\begin{eqnarray}&&\int {T}_{\mathrm{RJ}}\,{\rm{d}}v=\displaystyle \frac{\sqrt{\pi }}{2\sqrt{\mathrm{ln}2}}\,{T}_{\mathrm{RJ}}\ {\rm{\Delta }}v,\end{eqnarray} \tag{ 5 }$$

where Δv is the FWHM line width. Finally, we expand the grid with an additional axis of beam-filling factor (Φ_bf) by multiplying the predicted line integrated intensities with the varying Φ_bf. This assumes the same beam-filling factor among all six emission lines. Note that since we are essentially fitting N_CO/Δv with RADEX (see Equation (3)), the estimation of N_CO would vary with line widths in different regions. Thus, when N_CO values are needed in further analyses, we multiply N_CO/Δv with the observed CO 1–0 line width for each pixel. This assumes all six lines to have the same line width, which is consistent with our assumption of single-component gas in Section 4.2. We have checked that adopting the CO 2–1 line width also results in a similar N_CO, as the line widths of 1–0 and 2–1 agree to within 50%. The high-S/N emission in ¹³CO also yields similar line widths to those from the low-J CO emission—e.g., the mean CO 1–0/¹³CO 2–1 line width ratio is 1.07 ± 0.21. We do not consider the low-S/N pixels because the line widths are very uncertain at low S/N and therefore do not provide a strong constraint.

To evaluate the goodness of fit for each parameter set $\theta =({n}_{{{\rm{H}}}_{2}},{T}_{{\rm{k}}},{N}_{\mathrm{CO}}/{\rm{\Delta }}v,{X}_{12/13},{X}_{13/18},{{\rm{\Phi }}}_{\mathrm{bf}})$, we compute the χ² values at each point in the model grid:

$$\begin{eqnarray}&&{\chi }^{2}(\theta )=\sum _{i=1}^{n}\displaystyle \frac{{\left[{S}_{i}^{\mathrm{mod}}(\theta )-{S}_{i}^{\mathrm{obs}}\right]}^{2}}{{\sigma }_{i}^{2}}\,,\end{eqnarray} \tag{ 6 }$$

where S^mod and S^obs are the modeled and observed line integrated intensities with uncertainty σ_i , and n = 6 since we consider six lines. To include both the measurement and calibration uncertainties, σ_i can be written as ${\sigma }_{i}^{2}={\sigma }_{\mathrm{noise},i}^{2}+{\sigma }_{\mathrm{cal}}^{2}$, where σ_noise,i is simply the measurement uncertainty of the ith observed line intensity and σ_cal represents the calibration uncertainty. For ALMA Band 3, 6, and 7 observations, a flux calibration uncertainty of ≳5% is suggested for point sources (Bonato et al. 2018, 2019), and an additional few percent should be added for extended sources (e.g., Tunnard & Greve 2016; Leroy et al. 2021b). Thus, we adopt a 10% calibration uncertainty (i.e., σ_cal = 0.1 × flux) for each line and assume independent calibrations. Although there should exist correlations between the calibration uncertainties for lines observed in the same spectral setup with the same calibrator, assuming independent calibrations is a less stringent constraint and therefore more conservative for the modeling.

The best-fit parameter set can be determined by selecting the combination with the lowest χ² value. By assuming a multivariate Gaussian probability distribution, the χ² value for each parameter set can be converted to the probability as

$$\begin{eqnarray}&&P({S}^{\mathrm{obs}}| \theta )=\displaystyle \frac{1}{Q}\,\exp \left[-\displaystyle \frac{1}{2}{\chi }^{2}(\theta )\right]\,,\end{eqnarray} \tag{ 7 }$$

where Q² = Π_i (2π σ_i ). We calculate this probability for each grid point in the model parameter space, and we make our grid space wide enough to cover all grid points with reasonably high likelihood. Next, we generate the marginalized 1D or 2D likelihood distributions for any given parameter(s) by summing the joint probability distribution over the full range of all other parameters except the one(s) of interest. With such marginalized 1D likelihood distributions, we find the "1DMax" solutions that give the highest 1D likelihood in each parameter. We also determine the 50th percentile values (i.e., medians) of the cumulative 1D likelihoods and compare them with the 1DMax values to better understand the probability distributions. This procedure of computing the χ² values and the likelihood distributions is applied to every pixel of the galaxy image in order to generate maps of the galaxy showing the derived physical parameters. In Sections 4.2 and 4.3, we present these results based on a one-component model with six free parameters and a two-component model with eight free parameters, respectively. For reproducibility, we release all of the source code and parameters in a GitHub repository. ²⁰

Assuming that the gas along each line of sight is uniform and that the emission can be described by a single physical condition (${n}_{{{\rm{H}}}_{2}}$, T_k, N_CO/Δv, X_12/13, X_13/18) with an identical filling factor Φ_bf across all lines, we create intensity model grids for each observed transition. As listed in Table 3, the model can be represented by a six-dimensional (6D) grid with $\mathrm{log}({n}_{{{\rm{H}}}_{2}}\,[{\mathrm{cm}}^{-3}])$ varied from 2 to 5 in steps of 0.2 dex, T_k from 10 to 200 K in steps of 0.1 dex, N_CO/Δv from 10¹⁶/15 to 10²¹/15 ${\mathrm{cm}}^{-2}\ {(\mathrm{km}\,{{\rm{s}}}^{-1})}^{-1}$ in steps of 0.2 dex, X_12/13 from 10 to 200 in steps of 10, X_13/18 from 2 to 20 in steps of 1.5, and Φ_bf from 0.05 to 1 in steps of 0.05. The range of N_CO/Δv covers a reasonable ${N}_{{{\rm{H}}}_{2}}$ range of ∼10²⁰–10²⁵ cm⁻² with typical CO/H₂ abundance ratios of ∼10⁻⁴, and all of the parameter ranges were optimized from representative pixels such that the shapes of the 1D likelihood distributions are well covered for typical bright and faint regions in the galaxy.

Table 3. Model Grid Parameters

Parameter	Range	Step Size
	One-component Models
$\mathrm{log}({n}_{{{\rm{H}}}_{2}}\,[{\mathrm{cm}}^{-3}])$	2.0–5.0	0.2 dex
$\mathrm{log}({T}_{{\rm{k}}}\,[{\rm{K}}])$	1.0–2.3	0.1 dex
$\mathrm{log}({N}_{\mathrm{CO}}\,[{\mathrm{cm}}^{-2}])$	16.0–21.0	0.2 dex
X12/13	10–200	10
X13/18	2–20	1.5
Φbf	0.05–1.0	0.05
Δv [km s−1]	15.0	⋯
	Two-component Models
$\mathrm{log}({n}_{{{\rm{H}}}_{2}}\,[{\mathrm{cm}}^{-3}])$	2.0–5.0	0.25 dex
$\mathrm{log}({T}_{{\rm{k}}}\,[{\rm{K}}])$	1.0–2.0	0.1 dex
$\mathrm{log}({N}_{\mathrm{CO}}\,[{\mathrm{cm}}^{-2}])$	16.0–19.0	0.25 dex
$\mathrm{log}({{\rm{\Phi }}}_{\mathrm{bf}})$	−1.3 to −0.1	0.1
X12/13	25	⋯
X13/18	8	⋯
Δv [km s−1]	15.0	⋯

Download table as:  ASCIITypeset image

We apply a prior on the path length of the molecular gas along the line of sight. The idea of applying priors on line-of-sight lengths to avoid unrealistic solutions was also adopted in Kamenetzky et al. (2014) and is supported by Tunnard & Greve (2016). The line-of-sight path length is calculated by ${{\ell }}_{\mathrm{los}}={N}_{\mathrm{CO}}{(\sqrt{{{\rm{\Phi }}}_{\mathrm{bf}}}\ {n}_{\mathrm{CO}})}^{-1}$, where ${n}_{\mathrm{CO}}={n}_{{{\rm{H}}}_{2}}\,{x}_{\mathrm{CO}}$ and x_CO is the CO/H₂ abundance ratio for which we assume a typical value of 3 × 10⁻⁴ for starburst regions (Lacy et al. 1994; Ward et al. 2003; Sliwa et al. 2014). The $\sqrt{{{\rm{\Phi }}}_{\mathrm{bf}}}$ factor can be interpreted as the one-dimensional (1D) filling factor along the line of sight to match the area filling factor Φ_bf, and the same equation was also adopted by Ward et al. (2003) and Kamenetzky et al. (2012, 2014). With an angular resolution of ∼100 pc, we require all parameter sets in our model grids to have a line-of-sight path length ℓ_los ≤ 100 pc. The 100 pc path length is also consistent with the vertical height (thickness) of CO observed in the central few kiloparsec of the Milky Way and other nearby edge-on galaxies (Yim et al. 2014; Heyer & Dame 2015). We have checked that relaxing the constraint to 200 pc results in similar best-fit and 1DMax solutions, and thus the result is not sensitive to factor of ∼2 changes in the line-of-sight constraint. The change in the gas thickness along the line of sight due to galaxy inclination would not alter the results. We note that all of the parameters for calculating ℓ_los are the varied inputs of our model except for the constant x_CO and the observed line width to obtain N_CO from N_CO/Δv. In other words, we adopt ℓ_los < 100 pc as a prior by setting a zero probability to all of the parameter sets that violate this constraint, and then investigate the resultant likelihood distributions. This prior generally rules out solutions with N_CO ≳ 10¹⁹ cm⁻² or ${n}_{{{\rm{H}}}_{2}}\lesssim 3\times {10}^{2}$ cm⁻³, as either a high column density or a low volume density would lead to large ℓ_los. We also find that the primary effect of allowing larger path lengths would be to shift the solutions to higher N_CO while ${n}_{{{\rm{H}}}_{2}}$ and T_k stay roughly constant, and we do not think that such large N_CO values and path lengths are reasonable solutions.

Figures 5 and 6 are two corner plots (Foreman-Mackey 2016) showing the marginalized likelihood distributions of a bright pixel at $({\rm{R}}.{\rm{A}}.,{\rm{decl}}.)=({10}^{{\rm{h}}}{43}^{{\rm{m}}}57\mathop{.}\limits^{s}9,11^\circ 42{\rm{^{\prime} }}19\buildrel{\prime\prime}\over{.} 5)$ in the northern contact point and a faint pixel at (${10}^{{\rm{h}}}{43}^{{\rm{m}}}57\buildrel{\rm{s}}\over{.} 8,11^\circ 42^{\prime} 15\buildrel{\prime\prime}\over{.} 5$) in the gap region between the ring and the nucleus. The 1D likelihoods of N_CO, T_k, ${n}_{{{\rm{H}}}_{2}}$, and X_13/18 are well constrained in both pixels. By contrast, X_12/13 is loosely constrained compared to other parameters, and the peaks in both 1D likelihoods (i.e., 1DMax, indicated by the dashed lines) are consistently at a lower X_12/13 of ∼20–30. This leads to a deviation between the 1DMax solutions and the medians (median X_12/13 ∼ 90–100, indicated by the dotted lines). Even though the 1D likelihoods of X_12/13 along each sight line may be less constrained than other parameters, we will show that the 1DMax solutions of X_12/13 are consistent over the central ∼1 kpc region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that some parameters show correlations in their 2D likelihoods, including ${N}_{\mathrm{CO}}-{n}_{{{\rm{H}}}_{2}}$, ${T}_{{\rm{k}}}-{n}_{{{\rm{H}}}_{2}}$, and N_CO−Φ_bf. Below the critical densities of the optically thin lines (see Table 1), the emissivity depends on the density. Therefore, the anticorrelation between N_CO and ${n}_{{{\rm{H}}}_{2}}$ is expected, as a larger column density is needed to produce the same emission if we decrease the volume density with all other parameters fixed. Similarly, as ${n}_{{{\rm{H}}}_{2}}$ decreases, T_k has to increase to produce the same intensity if all other parameters remain unchanged. The inversely correlated N_CO−Φ_bf may imply a well-constrained total mass of CO, as M_mol ∝ N_CO Φ_bf (see Equation (9)). In Figure 7, we show how well the best-fit (i.e., highest likelihood in the grid) solutions for the two corresponding pixels in Figure 5 and 6 match the constraints given by the observed line fluxes. The contour values are set as the ranges of observed line intensities ± 1σ uncertainties, which include the measurement uncertainty from the data and a 10% flux calibration uncertainty. Note that the best-fit solutions of X_12/13 also deviate from the 1DMax solutions at ∼20–30.

Zoom In Zoom Out Reset image size

We present the pixel-by-pixel maps and regional statistics of all six 1DMax parameters in Figure 8 and Table 4, respectively. It is clear that the CO column densities are highest in the center and the ring, while the average H₂ densities are similar in all regions. Although the average temperatures from Table 4 are also similar in each region, Figure 8(b) reveals a higher kinetic temperature of ∼30–60 K near the nucleus and both contact points. High-resolution studies at pc scales toward the inner ∼500 pc (Central Molecular Zone, CMZ) of the Milky Way suggested a dominant gas component with a temperature of ∼25–50 K (Longmore et al. 2013; Ginsburg et al. 2016; Krieger et al. 2017), which is consistent with the temperatures we find near the nucleus and contact points. In addition, we find that these peaks in the T_k map cover several circumnuclear star-forming regions identified by previous UV, Hα, near-infrared, or radio observations (Colina et al. 1997; Planesas et al. 1997; Hägele et al. 2007; Linden et al. 2020; Calzetti et al. 2021), and thus higher temperatures could be related to the heating from young stars. We also find consistent 1DMax values of X_12/13 ∼ 20–30 and X_13/18 ∼ 6–10 inside the central ∼1 kpc region. Unlike for X_12/13, the 1DMax solutions for X_13/18 are consistent with the medians and are very well constrained over the central region. Both the 1DMax X_12/13 and X_13/18 abundance ratios are consistent with those found in the central kiloparsec of the Milky Way (Langer & Penzias 1990; Wilson & Rood 1994; Milam et al. 2005; Wouterloot et al. 2008; Areal et al. 2018). We note that only the "center" pixels defined in Figure 2(c) have S/N > 3 in the C¹⁸O images, and thus a pixel-by-pixel estimation beyond this region could be subject to larger uncertainties. Therefore, to examine our pixel-based results in the arm regions, we conduct similar analysis using stacked spectra in Section 5.2.

Zoom In Zoom Out Reset image size

Table 4. The 1DMax Solutions from One-component Modeling

Region		$\mathrm{log}\left({N}_{\mathrm{CO}}\,\tfrac{15\,\mathrm{km}\,{{\rm{s}}}^{-1}}{{\rm{\Delta }}v}\right)$	$\mathrm{log}{T}_{{\rm{k}}}$	$\mathrm{log}{n}_{{{\rm{H}}}_{2}}$	X12/13	X13/18	Φbf
		(cm−2)	(K)	(cm−3)
Whole	Median	17.8	1.2	3.2	30	8.0	0.20
	Mean	17.46	1.35	3.46	51.46	10.43	0.23
	Std. Dev.	1.13	0.41	0.76	58.31	6.48	0.19
Center	Median	18.4	1.3	3.2	20	6.5	0.25
	Mean	18.30	1.34	3.33	22.42	7.38	0.26
	Std. Dev.	0.59	0.22	0.40	14.64	1.66	0.16
Ring	Median	18.4	1.3	3.2	20	6.5	0.20
	Mean	18.21	1.33	3.38	22.94	7.12	0.26
	Std. Dev.	0.65	0.22	0.43	16.66	1.74	0.18
Arms	Median	16.2	1.1	3.2	40	15.5	0.15
	Mean	16.83	1.35	3.48	78.56	12.38	0.21
	Std. Dev.	0.99	0.51	0.93	67.51	7.86	0.20

Note. The statistics are determined for the ensemble of 1DMax solutions across the pixels. The standard deviation does not reflect the uncertainties in the 1D likelihoods of each pixel.

Download table as:  ASCIITypeset image

While single-component modeling has been widely used to determine physical properties of molecular gas in various galaxies (e.g., Topal et al. 2014; Sliwa et al. 2017; Imanishi et al. 2018; Teng & Hirano 2020), it is likely that the molecular line emission comes from multiple gas components with different physical conditions (e.g., Ginsburg et al. 2016; Barnes et al. 2017; Krieger et al. 2017). Some studies have modeled the gas with smoothly varying distributions of temperatures and densities (e.g., Leroy et al. 2017; Bisbas et al. 2019; Liu et al. 2021), while others modeled the emission with a warm and cold or dense and diffuse component (e.g., Kamenetzky et al. 2014; Schirm et al. 2014; Liu et al. 2015). In addition, the theoretical expectation for isolated and virialized molecular clouds and suggestively some observational studies of galaxy centers imply a two-phase molecular gas with a high-α_CO dense component and a low-α_CO diffuse component (e.g., Papadopoulos et al. 2012a; Leroy et al. 2015). Israel (2020) also conducted a two-phase analysis toward the centers of ∼100 galaxies with single dish measurements, although NGC 3351 was not included. To test how multiple gas components could change our results, we construct a two-component model of N_CO/Δv, T_k, ${n}_{{{\rm{H}}}_{2}}$, and Φ_bf. This is primarily to see if we get dramatically different results if we change our assumptions about the gas components.

Based on the consistent 1DMax values of X_12/13 and X_13/18 from the one-component modeling and their consistency with the values of the Milky Way center, we assume fixed isotopologue abundances of X_12/13 = 25 and X_13/18 = 8 for both gas phases. The two-component model can thus be described by eight parameters, namely, four parameters for each gas phase. We remind the readers that even if the 1DMax values of X_12/13 are consistently at ∼25 in the one-component modeling, the 1D likelihoods of X_12/13 are loosely constrained; the X_13/18 likelihoods, by contrast, are well constrained at 6–10 (see Figures 5 and 6). Jiménez-Donaire et al. (2017) have observed the central kiloparsec of NGC 3351 at an 8'' resolution with ¹³CO and C¹⁸O 1–0 and suggested an X_13/18 value that is also consistent with our one-component modeling as well as the Galactic Center value. Nevertheless, we note that applying fixed, Galactic Center like abundance ratios to the two-component modeling may be a rather strict assumption. We emphasize that the primary goal of our two-component model is to test if the solutions will differ substantially after altering the assumption of the numbers of gas components, and this goal can be achieved by comparing the solutions of N_CO/Δv, T_k, or ${n}_{{{\rm{H}}}_{2}}$ while fixing the isotopologue abundances to remain consistent with the one-component modeling.

Due to large grid size of the eight-dimensional (8D) model, we slightly adjust the parameter ranges and step sizes (see Table 3): N_CO/Δv ranges from 10¹⁶/15 to 10¹⁹/15 ${\mathrm{cm}}^{-2}\,{(\mathrm{km}\,{{\rm{s}}}^{-1})}^{-1}$ in steps of 0.25 dex, T_k ranges from 10 to 100 K in steps of 0.1 dex, $\mathrm{log}({n}_{{{\rm{H}}}_{2}}\,[{\mathrm{cm}}^{-3}])$ ranges from 2 to 5 in steps of 0.25 dex, and $\mathrm{log}({{\rm{\Phi }}}_{\mathrm{bf}})$ ranges from −1.3 to −0.1 in steps of 0.1 dex. We lower the upper limit of T_k to 100 K in the two-component model grids because the sensitivity to distinguish different temperatures above 100 K is very weak with J = 3–2 as the highest transition. In addition, the one-component modeling results show that T_k is well below 100 K across the whole region. We also lower the upper limit of N_CO/Δv to 10¹⁹/15 because a higher value is normally ruled out by the line-of-sight constraint of ℓ_los < 100 pc as mentioned in Section 4.2.

The 8D model grid consists of all combinations of either two components from the four-dimensional (4D) parameter space. By summing up the predicted intensities of both components, each grid point has a prediction of the total integrated intensity. Since the two components are interchangeable under summation, nearly half of the grid points are redundant, and we thus require the first component to have higher T_k than the second component. This not only allows us to distinguish between warm and cold or dense and diffuse components, but also avoids getting the exact same marginalized likelihoods for both components due to symmetric model grids. Similar to the 6D one-component modeling, we include a 10% flux calibration uncertainty and the measurement uncertainties for fitting.

When deriving the values of N_CO, we assume both components to have the same line widths as the observed CO 1–0 line width. This assumption has a negligible effect on our α_CO results in Section 5, as N_CO and Φ_bf are being fit simultaneously to determine α_CO values, which means that the two parameters can adjust themselves to fit best with the line intensities. In addition, we mostly observed a single-peaked spectrum that can be well represented by a single line width (see Figure 4), so assuming the same line widths as the measured line width is likely the most feasible thing we could do in the two-component analysis. If both components have Gaussian spectra, this assumption would imply that both components are at the same velocity.

Figure 9 shows the 1DMax solutions from the two-component modeling, where the component with a higher ${n}_{{{\rm{H}}}_{2}}$ is shown in the top row. Regional statistics are also summarized in Table 5 with the denser component being presented as the first component. We note that the separation by density is made for visualization purposes and has no role in the subsequent calculation of α_CO. In general, the dense components correspond to lower T_k and N_CO values than the diffuse components, which is similar to the correlations observed in the one-component result within the probability distribution of a single pixel. From the regional averages of both components (see Table 5), we find that the differences in N_CO, ${n}_{{{\rm{H}}}_{2}}$, and Φ_bf between the two components are all less than ∼0.5 dex. Moreover, the difference between their masses, which are proportional to their N_CO and Φ_bf (see Equation (9)), is even found to be less than ∼0.3 dex. This means that the two components have similar physical properties. We also find these properties to be similar to those suggested by the one-component model. For example, the 1DMax ${n}_{{{\rm{H}}}_{2}}$ distribution from the two-component model suggests a density of ∼2–3 × 10³ cm⁻³ in the center, which is consistent with ∼2 × 10³ cm⁻³ suggested by the one-component model. These are also similar to the average gas density of ∼5 × 10³ cm⁻³ found in the Milky Way's CMZ (e.g., Longmore et al. 2013). The only substantial difference between our two-component solutions is the kinetic temperature. While the denser component has a similar range of temperature and region-by-region variations as observed in the one-component model (see Figure 8(b)), the other component shows an evidently higher kinetic temperature of ∼100 K over the whole region. These temperatures are quantitatively similar to previous two-component studies toward the GMCs in our Galactic Center, which suggest a dominant component by mass with T_k ∼ 25–50 K and a less-dominant warm component with T_k ≳ 100 K (e.g., Huettemeister et al. 1993; Krieger et al. 2017). We thus conclude that the dominant component in our two-component model is well represented by the one-component model alone, but there might be a secondary, warmer component unaccounted for in the one-component model, as previously seen in the Galactic CMZ.

Zoom In Zoom Out Reset image size

Table 5. The 1DMax Solutions from Two-component Modeling

Region		$\mathrm{log}\left({N}_{1}\,\tfrac{15\,\mathrm{km}\,{{\rm{s}}}^{-1}}{{\rm{\Delta }}v}\right)$	log T1	log n1	$\mathrm{log}{{\rm{\Phi }}}_{1}$	$\mathrm{log}\left({N}_{2}\,\tfrac{15\,\mathrm{km}\,{{\rm{s}}}^{-1}}{{\rm{\Delta }}v}\right)$	log T2	log n2	$\mathrm{log}{{\rm{\Phi }}}_{2}$
		(cm−2)	(K)	(cm−3)		(cm−2)	(K)	(cm−3)
Whole	Median	17.00	1.1	3.75	−0.6	17.00	2.0	3.25	−0.8
	Mean	17.01	1.27	3.55	−0.68	17.17	1.74	3.00	−0.79
	Std. Dev.	0.85	0.31	0.74	0.42	1.00	0.36	0.64	0.40
Center	Median	17.75	1.3	3.50	−0.5	18.25	2.0	3.25	−0.7
	Mean	17.80	1.38	3.61	−0.55	18.11	1.90	3.15	−0.71
	Std. Dev.	0.53	0.25	0.48	0.30	0.45	0.24	0.39	0.26
Ring	Median	17.75	1.4	3.75	−0.5	18.00	2.0	3.25	−0.7
	Mean	17.83	1.42	3.65	−0.57	18.11	1.89	3.16	−0.74
	Std. Dev.	0.58	0.26	0.94	0.34	0.47	0.27	0.39	0.29
Arms	Median	16.25	1.0	3.50	−0.6	16.25	1.9	3.25	−0.8
	Mean	16.38	1.13	3.36	−0.71	16.39	1.67	2.88	−0.77
	Std. Dev.	0.45	0.25	0.83	0.48	0.55	0.37	0.70	0.47

Note. The denser component is represented as the first component, i.e., N₁, T₁, n₁, and Φ₁.

Download table as:  ASCIITypeset image

In this subsection, we describe the unified procedure of how we generate and deal with the marginalized likelihoods and then derive plausible solutions for parameters like α_CO. This is fundamental to the results presented in Section 5.

In Sections 4.2 and 4.3, we derived the marginalized 1D likelihoods for each input parameter. In both cases, we have a full 6D or 8D grid, with each grid point having an associated probability. Since our input parameters are sampled uniformly, their marginalized 1D likelihood distribution is simply a probability-weighted histogram, where the probability value of each grid point is reflected proportionally in the number of counts toward each bin. However, to derive the marginalized likelihoods of parameters that are functions of the intrinsic grid parameters (e.g., the modeled line intensities or α_CO), an additional normalization on the histograms would be needed due to the possibility of irregular sampling on the derived parameter grid space. For instance, given the grids of modeled line integrated intensities, we could determine 1D likelihoods for each line intensity with some chosen bin spacing and range. However, since multiple combinations of parameters in the 6D or 8D grid space could produce the same integrated intensity, and some line intensities may only be produced by a small subset of the parameter combinations, the resulting 1D likelihoods will have an additional nonuniform weighting due to the starting grid. To properly deal with the grid irregularity, we can generate a uniformly weighted histogram for line intensities under the same parameter range and bin size, treating that as a normalization. To obtain the marginalized 1D likelihood of any derived parameter from the original grid, we divide the probability-weighted histogram by the uniformly weighted histogram to normalize out any grid irregularity. The marginalized likelihoods of α_CO, which is a function of N_CO, Φ_bf, and the CO 1–0 line intensity, can also be derived in this way (see Section 5).

Once the 1D likelihood distributions of a parameter are obtained, we can determine either the 1DMax solutions by finding the values that correspond to the peaks, or the medians, and ±1σ values by finding the 16th and 84th percentiles of the cumulative 1D likelihood distribution. Then, maps of the 1DMax solutions or medians ± 1σ can be generated by iterating over each pixel. While we could look individually into the likelihood distributions at each pixel, comparing between these maps is a more simple and feasible way to get a sense of the probability distribution in different regions of the galaxy. With this procedure, it is important to note that our solutions for α_CO in Section 5 do not rely on the individually derived N_CO and Φ_bf distributions presented in Sections 4.2 and 4.3, since those parameters are being fit simultaneously within the full grid before we obtain the marginalized likelihoods of α_CO.

Using the method described in Section 4.4, we derive the marginalized 1D likelihoods of $\mathrm{log}({\alpha }_{\mathrm{CO}})$ for each pixel with $\mathrm{log}({\alpha }_{\mathrm{CO}})$ ranging −2.5 to 2.5 in steps of 0.125. Figures 10 and 11 present the α_CO variations derived from the one-component and two-component modeling, respectively. The marginalized 1DMax α_CO shown in Figure 10(a) reveals an average of ∼2 M_⊙ (K km s⁻¹ pc²)⁻¹ in the inner 1 kpc region, while the α_CO along the arms is a factor of ∼10 lower than the center. Figure 10(b) shows how the median and 1DMax α_CO vary with the galactocentric radius of NGC 3351. Both the median and 1DMax α_CO values exhibit a slow decrease from the nucleus to the ring and a significant drop at a radius of ∼10'' (∼0.5 kpc), which is where the ring and the arms intersect. Interestingly, while α_CO remains roughly constant at a radius of ∼10–20'', it starts to increase beyond a radius of ∼20''. We find that the ≳20'' region corresponds to the curved-up feature in the southern arm (see Figure 2(c)). This region is called the "transverse dust lane" in Leaman et al. (2019), where they suggested the feature is evidence for an outflow pushing gas/dust away from the nuclear ring via stellar feedback processes.

Zoom In Zoom Out Reset image size

Overall, all of the pixels in our observed region show at least factor of ∼2–3 (0.3–0.5 dex) lower α_CO values than the Milky Way average at the assumed x_CO of 3 × 10⁻⁴. If x_CO were set to 10⁻⁴, α_CO would increase by a factor of three. In such a case, the central 1 kpc region would have a nearly Galactic α_CO, whereas the α_CO is still substantially lower in the inflow arms since a global change in x_CO would not change the relative drop in α_CO. Note that if there were spatial variations of x_CO across this region, it could alter our α_CO variation results. Leaman et al. (2019) estimated the circular rotational velocity to be 150–200 km s⁻¹ on the circumnuclear ring of NGC 3351, implying a gas rotation period of 15–20 Myr. Since this orbital timescale is much shorter than the stellar evolution timescale that can cause a difference in chemical abundances, we expect the abundances of carbon and oxygen to be well mixed in the central kpc of NGC 3351, and thus sub-kpc scale variations of x_CO should be small if x_CO is essentially determined by the availability of carbon and oxygen. Despite such expectation, we emphasize that disentangling x_CO variations and α_CO is infeasible under current modeling and measurements. Therefore, our estimated α_CO distribution can be affected if there is any x_CO variation on sub-kpc scales.

The 1DMax α_CO distributions derived from the two-component modeling results, as shown in Figure 11(a), suggest a similar spatial variation to that determined by the one-component results, while the α_CO values in the center are generally lower than those in Figure 10(a). Figure 11(b) shows the median and 1DMax solutions determined from the marginalized α_CO likelihoods and their variation with galactocentric radius. The α_CO solutions inside the central 1 kpc region are shifted ∼0.3–0.5 dex downward compared with Figure 10(b). By contrast, the α_CO solution in the inflow arms is consistent with Figure 10(b), and we can see a clear distinction between α_CO values in the center and inflow regions, which are separated at the 10'' or 0.5 kpc radius. We note that the possible cause for a ∼0.3–0.5 dex shift of α_CO solutions in the center may originate from the assumption of a fixed X_12/13 at 25. The 2D correlations in Figure 5 and 6 show that X_12/13 = 25 corresponds to a lower N_CO/Δv and a higher Φ_bf compared with the 1DMax solutions, which may altogether result in a 0.3–0.5 dex offset in α_CO. In general, both our one- and two-component models result in lower than Galactic α_CO values for the entire observed region and significantly (≳10×) lower α_CO values in the inflow arms. Furthermore, by comparing with the derived T_k distributions shown in Figure 8(b), there is no sign that lower α_CO values are associated with higher temperatures (see discussion in Section 6). In Section 5.2, we will discuss possible scenarios that cause the α_CO variation and such low α_CO values in the arms.

Zoom In Zoom Out Reset image size

To compare the derived α_CO with those observed previously at kiloparsec scales, we compute the intensity-weighted mean α_CO over our entire observed region of ∼2 kpc. The intensity-weighted mean α_CO is calculated as the ratio of total molecular mass to the total luminosity over the whole region. Based on Equation (10), this can be done with the full grids of the modeled N_CO, Φ_bf, and intensity of CO 1–0 (S^mod) at every pixel. For each pixel, we first perform a likelihood-weighted random draw from the model grids to obtain N_CO, Φ_bf, and the corresponding modeled I_CO(1−0) value. Next, we sum up the values over the pixels to derive the total molecular mass and CO intensity, and then calculate the resulting intensity-weighted mean α_CO. We then perform another random draw and repeat 2000 times. Finally, the mean and standard deviation over the 2000 measurements of α_CO are obtained. This approach is a combination of Monte Carlo sampling with likelihood weighting of the model grids, and is similar to the "realize" method used by Gordon et al. (2014) for dust spectral energy distribution (SED) fitting.

The intensity-weighted mean α_CO of the entire region is 1.79 ± 0.10 and 1.11 ± 0.09 M_⊙ (K km s⁻¹ pc²)⁻¹ based on the one- and two-component models, respectively. Since the sum of our I_CO(1−0) and I_CO(2−1) over the entire observed region are similar, weighting either by the CO 1–0 or 2–1 intensities gives approximately the same intensity-weighted mean α_CO. With an angular resolution of ∼40'', Sandstrom et al. (2013) found an α_CO(2−1) of ${1.0}_{-0.3}^{+0.4}$ M_⊙ (K km s⁻¹ pc²)⁻¹ in the central 1.7 kpc of NGC 3351 using dust modeling with CO 2–1 intensities. This is similar to the intensity-weighted mean α_CO values derived from both our one- and two-component models. However, since they assumed a constant CO 2–1/1–0 ratio (R₂₁) of 0.7, they determined a lower mean α_CO(1−0) value of ${0.7}_{-0.20}^{+0.27}$ M_⊙ (K km s⁻¹ pc²)⁻¹ is appropriate for CO 1–0 intensities. While R₂₁ = 0.7 may be a good approximation across galaxy disks (e.g., den Brok et al. 2021), studies of nearby galaxy centers have suggested a higher R₂₁ (≳0.9) in the central kiloparsec region of disk galaxies (Israel 2020; Yajima et al. 2021). As presented in Table 2, the mean R₂₁ over our entire observed region is 1.0, and both our one- and two-component models predict a mean R₂₁ higher than 0.9. Thus, if a higher R₂₁ of >0.9 had been adopted by Sandstrom et al. (2013), their α_CO estimates with CO 1–0 would be in good agreement with our modeling results.

As pixel-based analyses in the arm regions may have higher uncertainties due to low S/N, we stack the spectra from inflow regions to recover low brightness emission applying the approach by Schruba et al. (2011). First, we set the moment 1 velocities of CO 2–1 as the reference velocity (i.e., v = 0) for each pixel and shift the spectra from all other lines to v = 0. Next, we extract the spectra in the v = ±100 km s⁻¹ range and sum them up to obtain a stacked spectrum for each line. Figure 12 presents the shifted and stacked spectra over the arms (defined in Figure 2(c)), overlaid with the best-fit Gaussian function. The spectra of the ¹³CO lines show much improvement in S/N after stacking. However, the spectrum of C¹⁸O 3–2 is still noisy, as most pixels in the arms are below 1σ in C¹⁸O 3–2 (see Figure 1) and therefore do not have reliable measurements even after stacking. Thus, we should be aware that C¹⁸O 3–2 is not a significant detection in the stacked spectrum of inflow arms, and it is included in the fitting as an upper limit.

Zoom In Zoom Out Reset image size

We measure the observed integrated flux of each line and use it as input to the multiline modeling technique described in Section 4.2. Figures 13 and 14 show the likelihood distributions and the ${T}_{{\rm{k}}}-{n}_{{{\rm{H}}}_{2}}$ slice where the best-fit parameter set lies. The 1DMax and best-fit physical conditions to the stacked spectra of inflow arms are found at (N_CO/Δv, T_k, ${n}_{{{\rm{H}}}_{2}}$, X_12/13, X_13/18, Φ_bf) = ($\tfrac{{10}^{16}}{15}$ ${\mathrm{cm}}^{-2}\ {(\mathrm{km}\,{{\rm{s}}}^{-1})}^{-1}$, 10 K, 10^3.6 cm⁻³, 30, 14, 0.05) and ($\tfrac{{10}^{16.6}}{15}$ ${\mathrm{cm}}^{-2}\ {(\mathrm{km}\,{{\rm{s}}}^{-1})}^{-1}$, 10 K, 10^4.4 cm⁻³, 30, 14, 0.3), respectively. Both solutions imply a lower N_CO, lower T_k, and higher ${n}_{{{\rm{H}}}_{2}}$ in the arms compared with those in the central ∼1 kpc. As shown in Figure 14, the constraints from five of the six lines have similar shapes, while CO 1–0 is the main constraint that pushes the solution to the lower-right corner. Therefore, even if we exclude the constraints from the C¹⁸O lines that have low S/N in the arms, the best-fit solution would remain, which means that C¹⁸O lines are not the key data in the arm regions. The line width determined by the best-fit Gaussian to the stacked CO 1–0 spectrum is Δv_FWHM = 25.18 ± 0.58 km s⁻¹. By substituting the best-fit or 1DMax solutions of N_CO/Δv and Φ_bf, and the fitted CO 1–0 line width into Equation (10), the α_CO of inflow arms is estimated to be 0.01–0.1 M_⊙ (K km s⁻¹ pc²)⁻¹. This matches our pixel-based estimation in the arms despite with the lower S/N, and thus supports the idea that the inflow regions have substantially lower α_CO values than the center/ring.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figures 3(d) and (e), we find that the CO/¹³CO and CO/C¹⁸O line ratios in the inflow arms are ∼2× higher than in the center. This means that we observe more CO emission in the arms compared to the total molecular gas mass, which by definition indicates a lower α_CO value. One possible scenario that could cause such high line ratios are high CO/¹³CO and CO/C¹⁸O abundance ratios. It is likely that the circumnuclear star formation activities have enriched ¹³C or ¹⁸O in the center but not in the inflows. This could lower the X_12/13 or X_13/18 abundances in the center and cause lower CO/¹³CO and CO/C¹⁸O line ratios than those in the arms. However, both the best-fit and 1DMax solutions based on the stacked spectra of inflow arms suggest X_12/13 ∼ 30 and X_13/18 ∼ 10, which are similar to those found in the center region as shown in Figures 8(d) and (e). Note that the 1D likelihood of X_12/13 is very well constrained for the arm regions with stacking (see Figure 13), which is different from the loosely constrained X_12/13 in the individual pixels of the center region. Even if the center pixels have a higher X_12/13 (e.g., median X_12/13 ∼ 90 in Figure 5 and 6), this would contradict the expectation of star formation enrichment and would instead enhance CO/¹³CO and CO/C¹⁸O line ratios in the center. Therefore, changes in CO isotopologue abundances may not be the main reason that causes higher line ratios in the inflow regions.

The high line ratios observed in CO/¹³CO and CO/C¹⁸O 2–1 could also be explained if the optical depths of CO changes significantly. Figure 2(b) reveals a broader line width in the arms, which may indicate a more turbulent environment or the existence of shear in these bar-driven inflows. This could lead to lower optical depths in the gas inflows. As shown in Figure 15, the optical depths derived from the 1DMax physical conditions are found to be lower in the arms compared to the center, which is likely a combined effect of low N_CO and high Δv (see Equation (3)). Since higher velocity dispersions and consequently lower optical depths allow a larger fraction of CO emission to escape, these may explain why α_CO is low in the inflow arms.

Zoom In Zoom Out Reset image size

In optically thin regions, the average α_CO is ${0.04}_{-0.03}^{+0.13}$ M_⊙ (K km s⁻¹ pc²)⁻¹ based on the marginalized 1DMax α_CO map shown in Figure 10(a). The optically thin regions were determined by selecting pixels having τ < 1 in Figure 15. Bolatto et al. (2013) derived an LTE equation to determine α_CO as a function of the CO/H₂ abundance (x_CO) and excitation temperature (T_ex):

$$\begin{eqnarray}&&{\alpha }_{\mathrm{CO}}\approx \displaystyle \frac{1.6\times {10}^{19}}{4.5\times {10}^{19}}\,\displaystyle \frac{{10}^{-4}}{{x}_{\mathrm{CO}}}\,\displaystyle \frac{{T}_{\mathrm{ex}}}{30\,{\rm{K}}}\,\exp \left(\displaystyle \frac{5.53\,{\rm{K}}}{{T}_{\mathrm{ex}}}-0.184\right),\end{eqnarray} \tag{ 11 }$$

where 4.5 × 10¹⁹ is the factor to convert ${X}_{\mathrm{CO}}={N}_{{{\rm{H}}}_{2}}/{I}_{\mathrm{CO}(1-0)}$ ${\mathrm{cm}}^{-2}\ {({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1})}^{-1}$ to α_CO (see Equation (10)). By plugging our modeled T_k (Figure 8(b)) and the assumed x_CO of 3 × 10⁻⁴ into Equation (11), we get an average α_CO of ${0.12}_{-0.07}^{+0.16}$ M_⊙ (K km s⁻¹ pc²)⁻¹ in the optically thin regions. This overlaps well with the average derived from our modeled α_CO distribution, showing that the solutions from multiline modeling are physically reasonable. We have also compared the excitation temperatures between CO and ¹³CO predicted by the best-fit physical conditions of the averaged spectra in each region. In the arm regions, the predicted excitation temperatures of CO and ¹³CO 2–1 lines are at 9.0 and 8.8 K, implying that the conditions in the arms are very close to LTE. In addition, the excitation temperatures of ¹³CO in the center/ring are lower than those of CO by a factor of ∼1.5–2, indicating a clear departure from LTE that may originate from radiative trapping affecting the CO excitation. Radiative trapping lowers the effective critical density for the CO lines, leading to differences in the excitation for CO and ¹³CO. We also note that temperature inhomogeneities in the gas will tend to bias CO lines to higher temperatures since warmer gas has a higher luminosity and, unlike ¹³CO or optically thin tracers, not all CO emission from the full line of sight contributes to the measured integrated intensity due to optical depth effects. It is also important to note that our derived values of α_CO vary inversely with x_CO. Thus, any increase/decrease to the assumed x_CO of 3 × 10⁻⁴ would lower/raise the numerical values of α_CO.

To compare our modeling results with other observations, we also acquired ALMA data of ¹³CO and C¹⁸O 1–0 with an angular resolution of 8'' toward the center of NGC 3351. These data were analyzed in Jiménez-Donaire et al. (2017) and Gallagher et al. (2018). Since the resolution is lower than the matched 21 resolution in our analysis, we only extract regional averages in the center, ring, and arms for comparison. First, we input the integrated intensities measured from the stacked spectra of each region (e.g., Figure 4) into the six-line modeling described in Section 4.2. This gives us a full probability grid for each region. Then, we generate the modeled intensity grids of ¹³CO and C¹⁸O 1–0 using RADEX, such that the modeled intensity at each grid point has an associated probability. Finally, we obtain the predicted 1D likelihoods for both lines using the method described in Section 4.4.

Figure 16(a) shows the predicted 1D likelihoods for the center region, covering an intensity range of 0 to 15 K km s⁻¹ in steps of 1 K km s⁻¹ for ¹³CO 1–0 and 0 to 3 K km s⁻¹ in steps of 0.25 K km s⁻¹ for C¹⁸O 1–0. The dotted lines mark measured intensities and the shaded regions show the ±1σ uncertainties that include the measurement uncertainty and a 10% calibration uncertainty. Following Schruba et al. (2011), we estimate the measurement uncertainties using the signal-free part of the spectra and the fitted line widths. In the center region, both the observed intensities are within ∼50% of the predicted 1DMax intensities. The 1D likelihoods and observed intensities for the ring region are similar to those in the center. However, Figure 16(b) shows that both the predicted and observed emission in the arms are fainter than the center/ring. For the arms, the 1D likelihoods range from 0 to 2 K km s⁻¹ in steps of 0.1 K km s⁻¹ for ¹³CO and 0 to 0.2 K km s⁻¹ in steps of 0.01 K km s⁻¹ for C¹⁸O. As there is no detection of C¹⁸O 1–0 in the arms, the dotted line in the right panel of Figure 16(b) marks the 1σ upper limit, which is also within 50% of the 1DMax intensity. Notably, the observed ¹³CO 1–0 intensity is a factor of ∼2 higher than the prediction. This may indicate a more optically thick CO or a lower X_12/13 abundance ratio in the arms compared with our modeling results. It could also imply that the excitation temperature of ¹³CO is underestimated in our models.

Zoom In Zoom Out Reset image size