The CO-to-H2 Conversion Factor in the Barred Spiral Galaxy M83

We analyze the CO-to-H2 conversion factor (α CO) in the nearby barred spiral galaxy M83. We present new H i observations from the VLA and single-dish GBT in the disk of the galaxy, and combine them with maps of CO(1-0) integrated intensity and dust surface density from the literature. α CO and the gas-to-dust ratio (δ GDR) are simultaneously derived in annuli of 2 kpc width from R = 1–7 kpc. We find that α CO and δ GDR both increase radially, by a factor of ∼2–3 from the center to the outskirts of the disk. The luminosity-weighted averages over the disk are α CO = 3.14 (2.06, 4.96) M⊙pc−2[Kkms−1]−1 and δ GDR = 137 (111, 182) at the 68% (1σ) confidence level. These are consistent with the α CO and δ GDR values measured in the Milky Way. In addition to possible variations of α CO due to the radial metallicity gradient, we test the possibility of variations in α CO due to changes in the underlying cloud populations, as a function of galactic radius. Using a truncated power-law molecular cloud CO luminosity function and an empirical power-law relation for cloud mass and luminosity, we show that the changes in the underlying cloud population may account for a factor of ∼1.5–2.0 radial change in α CO.


INTRODUCTION
Molecular hydrogen (H 2 ) is the most abundant molecule in the interstellar medium (ISM) and is important for star formation and galaxy evolution.However, since the cold H 2 is not detectable in emission, molecular gas masses are often traced using emission from the second most abundant molecule CO.This assumes a constant H 2 mass to CO luminosity ratio, deemed the CO-to-H 2 conversion factor (X CO or α CO ; Bolatto et al.

2013, and references therein). 1
There are four main techniques used to determine the conversion factor: (1) the virial method, which assumes clouds are gravitationally-bound and virialized (e.g.Sanders et al. 1984;Solomon et al. 1987), (2) γ-ray emission from interactions of cosmic-rays and protons Corresponding author: Amanda M Lee amamlee@umass.edu 1 Typically, the conversion factor denoted as X CO does not include the contribution from heavier elements (e.g., helium) in mass and is expressed in units of H 2 cm −2 K • km s −1 −1 .α CO includes the heavier elements, by multiplying 1.36 to the H 2 mass, and is in units of M ⊙ pc −2 K • km s −1 −1 .
There are also new efforts to use (5) [CI] or [CII] line emission, assuming that their flux can be converted to H 2 mass (Accurso et al. 2017;Madden et al. 2020).Some other methods are less frequently used and are applied only to some specific regions (e.g., solar neighborhood) including (6) stellar extinction (Lombardi et al. 2006), and (7) X-ray shadows (Sofue & Kataoka 2016).The assumptions embedded in all these methods can be questioned, but they have given a roughly consistent α CO in the Milky Way (MW) and in nearby galaxies with solar metallicity, with a caveat that a factor of 2-3 uncertainty remains, even with all the measurements combined (Bolatto et al. 2013).
The dust-based method (3) uses measurements of the HI and dust surface densities, Σ HI and Σ d , and CO surface brightness, I CO .δ GDR is first determined in HI-dominant regions as δ GDR = Σ HI /Σ d .Assuming that this δ GDR is the same in H 2 -dominant regions, the conversion factor is calculated as α CO = (δ GDR Σ d − Σ HI )/I CO (Reach et al. 1998;Dame et al. 2001).Leroy et al. (2011) and Sandstrom et al. (2013) expanded this method and fit the HI and H 2 -dominated regions simultaneously in a large sample of nearby star-forming galaxies at kpc-scale resolutions with the fitting equation, (1) They used CO(J=2-1) emission, instead of CO(J=1-0), for I CO from single-dish observations (Leroy et al. 2009) and Σ HI from interferometric data from the Very Large Array (Walter et al. 2008).The Σ HI term provides the absolute mass scale in this equation, and thus sets the absolute scales of δ GDR and α CO .They found an average α CO marginally smaller than, but still consistent with that of the MW value, with the exception of the galaxy centers, where α CO had values of up to an order of magnitude lower.
In this paper, we present new HI 21-cm line emission observations in M83 with the Jansky Very Large Array (JVLA) and Green Bank Telescope (GBT) and produce a high-quality image of the atomic hydrogen distribution via their combination.Together with CO (J=1-0) emission data from the ALMA single-dish telescopes (Koda et al. 2020), as well as dust emission data from the Herschel Space Observatory (Bendo et al. 2012) we determine α CO and δ GDR in M83 by simultaneously fitting for them.
M83 is a close nearby barred spiral galaxy located at a distance of 4.5 Mpc (Thim et al. 2003).It has been the subject of many studies for its similarities to the Milky Way and near face-on inclination of 26 • (Koda et al. 2023).

HI DATA
New HI 21-cm data of M83 were taken with the JVLA and GBT.These data were reduced and imaged along with archival JVLA data.

New and Archival Data
We used a total of 30 measurement sets (MSs) from the JVLA, including 22 archival MSs from projects .The archival data mosaicked 10 adjacent fields around M83, and we used only the data of the central pointing on the optical disk of M83.The additional 8 new MSs from project 21B-113 (PI: Koda) observed the single central pointing with higher sensitivity at long baselines.Each MS represents an observing track starting from observations of a bandpass calibrator, and then cycling between the gain calibrator and target fields.At the frequency of the HI 21cm line (ν HI = 1420.405751MHz, λ HI = 21.1061141cm), the primary beam of the JVLA antennae has a full-width at half power (FWHP) of 29.6'.

Data Reduction
The data were reduced using version 5.6.2-2 of the Common Astronomy Software Applications package (CASA; McMullin et al. 2007;Bean et al. 2022).The MSs include 10 spectral windows (SPWs).We used only one narrow band SPW with the HI line and one broadband SPW for calibration.The archival MSs have a narrow-band SPW with 2048 channels with a channel width of 1.95 kHz (∼ 0.4 km/s) and the new MSs have 4096 channels with a 0.97 kHz (∼ 0.2 km/s) width.
The standard calibration strategy was applied as described in the CASA guide 2 , including flagging for "shadow", "clip", and "quack", and corrections for antenna positions, elevation-dependent gain variations, and opacity.The flux scale was derived using the flux of 3C286, our bandpass calibrator, measured by the observatory.Delay and bandpass solutions were applied to the entire track.We then derived time-variant gain solutions with a broad-band SPW, and applied them to both the broad and narrow-band SPWs.Some bad data with unrealistically high amplitudes or phase scatters became obvious after these calibrations.We flagged those in the raw data and repeated all the calibrations from scratch.
For imaging we used only the narrow-band SPW.We subtracted the continuum emission in the uv space using line-free channels, and regridded each MS to a common velocity grid, covering the velocity range of 200-800 km s −1 with a channel width of 5 km s −1 for HI.Imaging was performed with the CASA task "tclean".The final JVLA data cube has a beam size of 9 ′′ .89× 8 ′′ .77and a sensitivity of 0.76 mJy/beam in each channel at the field center.

Observations
GBT observations were also made (project 20A-432).We used the single-element L-band dual-polarization receiver and the VErsatile GBT Astronomical Spectrom-Figure 1.An optical B-V-Hα composite image of M83 taken at the Cerro Tololo Inter-American Observatory and downloaded from the NASA/IPAC Extragalactic Database (Meurer et al. 2006;Cook et al. 2014).The white box indicates the 11.7 ′ × 11.7 ′ area observed in CO(1-0) (Koda et al. 2020).The coverage in HI is significantly larger than this area, with the VLA primary beam size of 29.6 ′ (FWHM) and GBT coverage even larger.Radial contours (cyan) show the radii from R=1-7 kpc at a 1 kpc interval with the position angle and inclination of 225 deg and 26 deg, respectively (Koda et al. 2023).
eter (VEGAS).We employed one of the eight banks in the narrowest bandwidth, 11.72 MHz with 32768 channels.The corresponding velocity coverage is about 2,400 km s −1 at a 0.08 km s −1 channel width.The beam size at the full width half maximum (FWHM) is ∼ 8.7 ′ at the HI frequency/wavelength.
A total of 5.25 hours was allocated over two sessions in May 2020.Each session started from receiver tuning, and telescope pointing and focus corrections using 1311-2216.A region of about 1.2 • × 1.5 • around M83 was mapped with the On-The-Fly (OTF) mapping technique in the RA and DEC directions, with a step size of 3.6 ′ (= 1 2 λ HI /D, between OTF scans, where D = 100 m is the dish diameter).Each OTF mapping in RA and DEC required 26 and 22 scans to cover the whole region, and each scan took 72 and 90 seconds, respectively.After every 4-5 scans, we integrated an OFF position of (RA,DEC)=(13:30:00.0,-29:51:19.0)for 30 seconds.We obtained a total of 6 OTF maps (3 in RA and 3 in DEC).We observed 3C286 for flux calibration.

Data Reduction
The data were reduced in a standard way with the GBTIDL software and pipeline.Time and elevationdependent amplitude variations were calibrated using the noise diode in the telescope, ON/OFF sky measurements, and T sys measurements.The calibrated spectra were gridded spatially to produce a data cube using the GaussBessel kernel function.We adopted spatial and velocity pixel scales of 1.47 ′ and 0.75 km/s.
In order to correct for relative flux variations, we made 12 data cubes for the two polarizations of the 6 OTF maps.We measured their flux ratios and inverted the ratio matrix to derive relative flux scaling factors.Their range was small and was 0.98-1.02(peak-to-peak) after setting the median to be unity.We applied these factors to the data and made the integrated data cube.
The absolute flux scale was calibrated using 3C286.This calibrator had been monitored and confirmed to be stable (Perley & Butler 2013) with the flux density of S ν = 14.97 Jy at 1.420 GHz (Perley & Butler 2017).

Combination of JVLA and GBT Data
The JVLA and GBT cubes were combined with the "feather" task in CASA.We regridded the GBT cube, applied the JVLA primary beam to it, combined it with the JVLA cube, and applied the JVLA primary beam correction to the combined cube.The final JVLA+GBT data cube has a beam size of 9 ′′ .89× 8 ′′ .77and a sensitivity of about 1.1 mJy/beam in each 5 km s −1 channel at the field center.This HI data covers a large area, and in this paper, we focus only on the disk of M83 (Figure 1).
We applied a mask to the final cube for generating an integrated intensity map for a 11'×11' box around the galaxy center (the region analyzed in Section 2.4; see Figure 2).We made the mask in two steps.First, we utilized the clean component cube from the JVLA imaging (Section 2.1.2),smoothed it with a Gaussian kernel with a FWHM of 60 ′′ (about the width of the HI spiral arms), and took regions with fluxes > 5 mJy/beam in the smoothed cube to make the first mask with pixel values of either 0 or 1.The first mask encloses significant emission in the JVLA imaging.Second, we expanded this mask by smoothing it with a Gaussian kernel with a FWHM of 90 ′′ and by taking the regions > 0.1.We start from the clean component cube because, by definition, it includes only the flux components identified as emission.This operation is to include any extended emission from the GBT data.The parameters in these operations were chosen manually by looking at the emission distributions in the cubes, and conservatively, so that the final mask may include blank sky but does not miss any emission.We used this mask also to calculate a noise map using the standard noise propagation based on the number of velocity pixels in each spatial pixel.Within the 11'×11' region, we calculated the total flux with and without applying the mask in the cube, and found that the ratio is ∼ 1.09.Thus, the mask likely encloses all of the emission.

Comparisons with Previous HI Data
To check for consistency with previous HI data, we acquired archival HI data cubes and integrated intensity (moment 0) maps of M83 from the THINGS survey (The HI Nearby Galaxy Survey Walter et al. 2008), which achieved a high resolution of 15 ′′ .16× 11 ′′ .44 with natural weighting without single-dish data, and the LVHIS survey (The Local Volume HI Survey; Koribalski et al. 2018), which achieved a lower spatial resolution of ∼ 85 ′′ .8× 60 ′′ .8but with single-dish data for flux recovery.While the THINGS survey also distributes data at a resolution of 10 ′′ .40× 5 ′′ .60 with uniform weighting without single-dish data, we used the naturally weighted data for comparison, since it typically reproduces the emission distribution more accurately.We make two comparisons: (1) between their archival data cubes and our cube to evaluate the quality of the data and (2) between their distributed moment 0 maps and our map to understand the difference in derived α CO (section 6.1).For (1) and (2) we smoothed and regridded our data to the resolution and pixel scale of the archival data.For the flux comparisons, we consider the ratio of total fluxes in the 11'×11' region shown in Figure 2.
For the first comparison, the ratio between the total flux of the archival data cubes and ours, without applying the mask, is THINGS/GBTVLA ∼ 0.39 and LVHIS/GBTVLA ∼ 0.92, where GBTVLA denotes the data from this work.After applying the mask, the ratios are THINGS/GBTVLA ∼ 0.65 and LVHIS/GBTVLA ∼ 0.84.Because of negatives in the THINGS data, the THINGS/GBTVLA ratio with and without the mask vary, and depend on the mask.
For the second comparison, we compare the archival moment 0 maps to our map because in Section 6.1, we will compare our results with those from previous studies, which used the distributed moment 0 map (Foyle et al. 2012;Jiao et al. 2021;Yasuda et al. 2023).These moment 0 maps were made with some masks applied, although we do not have access to those masks.The ratio between the total fluxes within the 11'×11' region is THINGS/GBTVLA ∼ 0.39 and LVHIS/GBTVLA ∼ 0.82. Figure 2 shows the ratio maps.The THINGS data show the galactic structures in HI (e.g., spiral arms) clearly at the high resolution, but suffer from missing flux, due to the lack of single-dish data.The fraction of recovered flux varies systematically as a function of the galactic structures.For example, it is as high as ∼80% along the spiral arms, but decreases to ∼20% in the interarm regions.This systematic variation is expected since interferometric data without single-dish data tend to miss extended emission.On the other hand, the LVHIS data is relatively more consistent with our new map, though it has a much lower resolution.The LVHIS HI fluxes also appear consistent with recent HI observations in M83 (Eibensteiner et al. 2023, see their Figure 3).
We see that the THINGS/GBTVLA ratio varies significantly depending on what mask was applied/not applied.Overall, the THINGS flux appears lower than that of our new data and analysis.

CO(1-0) AND INFRARED DATA
We compare the distributions of atomic hydrogen (in HI, discussed above), molecular hydrogen (in CO), and dust (in far-infrared).
We adopt the CO(J=1-0) integrated intensity map (and noise map) from Koda et al. (2020).This was obtained via observations with ALMA's Total Power (TP) array, which consists of four single-dish telescopes.The map has a low spatial resolution of 56.6 ′′ (FWHM; ∼ 1.2 kpc).We use the integrated intensity map in the main beam temperature in units of K • km s −1 .As a consistency check, we compare this map to the CO(1-0) map obtained from the Nobeyama 45-m telescope (Kuno et al. 2007) in Figure 2. The Nobeyama 45-m data is higher in resolution (15 ′′ ) but is limited in spatial coverage (radial extent R ∼ 4 kpc).The intensity weighted average ratio between the Nobeyama CO(1-0) and ALMA TP maps within R = 0-4 kpc is ∼ 0.8.
As a tracer of the dust distribution, we use the Herschel SPIRE 500µm image (Bendo et al. 2012;Foyle et al. 2013).The final calibration uncertainties in the dust map include a systematic uncertainty of ∼5% and random uncertainty of ∼2% (Bendo et al. 2012).
The image is obtained from the NASA/IPAC data archive for the Very Nearby Galaxy Survey (VNGS) (VNGS Team 2020)3 .

METHODS
We derive the CO-to-H 2 conversion factor (α CO ) and gas-to-dust mass ratio (δ GDR ) from the observations of Σ d , Σ HI , and I CO using eq.( 1).Σ HI and the molecular gas surface density Σ H2 , include the masses of heavier elements (e.g., helium).Figure 3 shows the HI, CO, and far-infrared intensity images that we use for the analysis.
Here we describe how to evaluate Σ HI , Σ H2 , and Σ d from the observed intensities, and how to use them for fitting to derive α CO and δ GDR .

HI and H 2 Surface Densities
(2) The coefficient, α CO , is the conversion factor and will be determined by fitting.Following the definition by Sandstrom et al. (2013), α CO includes the mass of helium (a factor of 1.36) and is in units of where ing the factor of 1.36 to account for helium) is derived theoretically in an optically-thin condition under local thermodynamic equilibrium (LTE).Middle: LVHIS/GBTVLA.Right: NRO/TP.GBTVLA and TP denote the HI and CO(1-0) data used in this study, while the archival HI and CO(1-0) data are denoted by THINGS and LVHIS, and NRO respectively.The map resolutions are at the lower resolution between the two used in the ratio (see text).Radial contours at R = 4 and R = 7 kpc are overlaid in magenta.In the 11'×11' region for the left two panels (HI data), the ratio of the total fluxes are THINGS/GBTVLA ∼ 0.39 and LVHIS/GBTVLA ∼ 0.82 (see text).For the rightmost panel (CO data), the intensity-weighted average ratio is ∼ 0.83 within R ≤ 4 kpc.

Dust Surface Density
galaxies over a wide wavelength range of dust emission from 3.6 to 500µm.They derived Σ d maps of the galaxies, as well as maps of other parameters, such as the dust fraction in polycyclic aromatic hydrocarbons (PAHs), background stellar light intensity, and infrared luminosity from the dust.So far, their study provides the most accurate Σ d maps.Sandstrom et al. (2013) adopted their maps for the analysis of α CO .
M83 is not among the sample of Aniano et al. (2020).However, Aniano et al. (2020) compared their Σ d maps with 500µm luminosity surface density maps (Σ νLν ) and found the relation, (4) They concluded that Σ νLν reproduces their Σ d maps to an accuracy of ∼ 0.07 dex over almost three decades of Σ νLν and we accordingly adopt an uncertainty of 20% for Σ d .Since this is larger than the uncertainties in the 500µm data mentioned in Bendo et al. 2012, we assume that 20% is the overall uncertainty.We also assume this includes the uncertainty from sky fluctuation including background sources, since even in the outermost region considered in our analysis, most pixels (98%) have intensities at least a factor of 5 higher than the typical sky fluctuation around M83.
They found that their Σ d can be estimated more reliably with this equation, than with the more widelyused SED-fitting with an optically thin modified blackbody approximation.A similar relation is obtained with galaxy-integrated values by Groves et al. (2015).Hence, we adopt this equation and convert the Herschel 500µm image into a Σ d map.
The Σ d map inherits the resolution and pixel scale of the Herschel 500µm image, which are 39.0 ′′ and 12.0 ′′ .We smooth this image with a Gaussian to the CO(1-0) resolution of 56.6 ′′ .The gas-to-dust ratio (GDR) in Aniano et al. (2020) is ∼ 140 in the conditions of the local interstellar medium (ISM).Thus, we measure the GDR relative to this value in the following sections.The Aniano et al. (2020) GDR is from the mass of atoms depleted from the gas phase in the local ISM: the hydrogen-to-dust mass ratio (including PAHs) is ∼ 103 (=1/0.0097;Draine & Hensley 2021, see their Table 2), and by taking into account the mass of helium (a factor 1.36), their GDR is ∼ 140 (Bruce Draine, in private communication).

Finding Solutions
We search for the solutions of eq. ( 1) within annuli of 2 kpc width centered on the galactic center from R = 1 to 7 kpc in 1 kpc increments (Section 5.1), as well as for the whole disk area of R =1-7 kpc (Section 5.2).The central R < 1 kpc is excluded in the analysis.The annuli are defined on the disk plane with the position angle of 225 deg and inclination of 26 deg (Koda et al. 2023).There is some redundancy between adjacent annuli.We adopt two methods to find the solutions.
The first method (the χ 2 method) is to minimize the χ 2 of the equation (referred to as "plane fits" in the Σ HI -I CO -Σ d space by Leroy et al. 2011).The second method (the ∆ log δ GDR method) is to minimize ∆ log δ GDR , the scatter in log δ GDR .Leroy et al. (2011) and Sandstrom et al. (2013) tested both methods and found no significant differences in results.They favored the ∆ log δ GDR method as they found it more stable in fitting.Here, we use both methods with a slight preference to the χ 2 method, as it gives a convenient estimation of errors, assuming that all errors in the analysis follow the normal distribution.
We use only the pixels with a signal-to-noise (S/N) ratio above 3 in all of the maps.95% of the pixels in the radial range of R = 1-7 kpc are used.The top panels of Figure 4 show the S/N maps for I CO and Σ HI .The pixels with S/N≤ 3 are blanked (white in these maps), and all of them are in the outermost annulus.The bottom panel shows a correlation between the I CO and Σ HI S/N ratios.The S/N ratio of Σ HI levels around ∼20-30, while that in I CO reaches higher values, because the observed Σ HI saturates around 10 M ⊙ pc −2 , above which the gas phase is generally molecular (Wong & Blitz 2002;Bigiel et al. 2008).In the lower S/N range (≲ 20), the S/N ratios in I CO and Σ HI are more or less comparable (i.e., one does not significantly dominate the other in a systematic way).Sandstrom et al. (2013) thoroughly discussed the capabilities and limitations of these fitting methods, including their dependence on the S/N of the tracers.The uncertainties here are small enough that the fitting methods, e.g. the ∆logδ GDR scatter minimization approach, should give unbiased results.

The χ 2 Minimization
Using eqs.( 1) and ( 2), a solution of (α CO , δ GDR ) is derived with the observed variables of (Σ HI , I CO , Σ d ).We rewrite eq. ( 1) as, (5) and define χ 2 as σ ΣHI , σ ICO , and σ Σ d are the errors in Σ HI , I CO , and Σ d respectively.The summations are over all image pixels, i, within a region of interest in the galaxy.The χ 2 values are calculated on a 2-dimensional grid of (α CO , δ GDR ).The minimum χ 2 is searched in this space.Figure 5a shows an example for an annulus of R = 3-5 kpc.If the errors follow the normal distribution, the χ 2 values can be converted to confidence intervals to assess the uncertainties (errors).The figure shows three contours at ∆χ 2 ≡ χ 2 − χ 2 min = 1, 4, and 9, which correspond to the confidence levels (CLs) of 68, 95, and 99.7% when they are projected to the α CO or δ GDR axis.The results of the fitting are discussed in Section 5.1.

The log δGDR Scatter Minimization
A solution can also be derived by minimizing the scatter, ∆ log δ GDR , in log δ GDR over a region of interest.For a fixed α CO , we calculate a map of and evaluate the scatter ∆ log δ GDR .This is repeated for a range of log α CO ∈ (−1.0, +2.0), hence for a wide range of α CO ∈ (0.1, 100.0), to find the α CO that minimizes the scatter.Following Sandstrom et al. (2013), the scatter and mean of log δ GDR are calculated as the scale and location of the robust statistics, respectively.Sandstrom et al. (2013) showed, in their Figure 1c, that the minimization of ∆ log δ GDR is equivalent to adopting a constant δ GDR (i.e., the inverse of their dust-to-gas ratio -DGR) independent of gas phase (i.e., I CO /Σ HI , the ratio between I CO and Σ HI ). Figure 5b shows an example for the same region (R = 3-5 kpc) as in Figure 5a.The color in the background shows ∆ log δ GDR as a function of α CO : each α CO corresponds to one value of ∆ log δ GDR , and hence, ∆ log δ GDR changes only in the horizontal direction.In addition, an average log δ GDR at each α CO can be evaluated with eq. ( 7).We adopt the biweight location to approximate the average following Sandstrom et al. (2013), which is shown as the white line in Figure 5b.Along this line, the point at the minimum of ∆ log δ GDR , in the color background, is the solution.The result of the fitting is discussed in Section 5.1.

Radial Gradients
We search for solutions in each annulus with a width of 2 kpc and increment of 1 kpc.The width is selected so that it is greater than the 56.6 ′′ resolution, which is 1.23 kpc at the distance of M83.Adjacent annuli overlap and are not independent from each other, but are useful in finding a radial trend.The central 1 kpc in radius is excluded as Σ HI cannot be measured accurately, due to the bright continuum emission at the galactic center, and HI absorption.
Table 1 summarizes the derived (α CO , δ GDR ) solutions for each fitting method and annulus.Overall, the (α CO , δ GDR ) values obtained between the fitting methods for the same annulus are similar.The values within the parentheses show the minimum and maximum α CO and δ GDR of the ∆χ 2 = 1 region, which would correspond to the 68% CL (i.e., 1σ) if all errors in this analysis follow the Gaussian statistics entirely, and gives some gauge of the error in this analysis.Under this treatment, at least for the χ 2 method, the overall errors in α CO in each annulus is about a factor of 2.
Figure 6 shows the radial variations of α CO and δ GDR .All the solutions suggest that both α CO and δ GDR increase as a function of radius, independent of the fitting method.If we take the results from the χ 2 method (column 4), α CO monotonically increases, by a factor of 2-3, from α CO ∼1.82 (1.04-3.53 at 68% CL) in R = 1-3 kpc to ∼5. 19 (3.83-6.90) δ GDR also increases by a factor of 2 from ∼ 95 (70.2, 149) in R = 1-3 kpc to ∼ 217 (191,250) in R = 5-7 kpc (column 5).Again, both methods show the same radial trend in δ GDR .Eq. ( 1) indicates that the absolute value of δ GDR is sensitive to the normalization of the scale of Σ d , which depends on the dust model.For example, if Σ d is overestimated by, say, a factor of 2, δ GDR would be underestimated by the same factor.Hence, the absolute value of δ GDR carries such systematic uncertainties.
Here, we assumed that α CO and δ GDR are constant in each annulus.There is discussion on the metallicitydependence of α CO and δ GDR (Arimoto et al. 1996;Israel 1997;Leroy et al. 2011;Schruba et al. 2012;Rémy-Ruyer et al. 2014;Hunt et al. 2015).The metallicity in M83 has been measured in several studies (e.g., Bresolin et al. 2009;Pilyugin et al. 2014;Grasha et al. 2022), and their absolute values are not necessarily consistent with each other, often due to different calibrations.However, the azimuthal variation, which is seen as the scatter at each radius, is consistently small in all the studies.Therefore, it seems safe to assume that the metallicity is approximately constant in each annulus, and thus that α CO and δ GDR do not suffer much from metallicity-dependent variations within each annulus.We find a change of ∆(12 + log(O/H)) = -0.06within one annulus using the radial gradients discussed in Section 6.2.1.

Galaxy-wide Averages
The average α CO and δ GDR give the representative values of the whole galaxy.Since both of these parameters vary systematically along the galactic radius, their (3) Total CO luminosity in corresponding annulus with units [K • km s −1 ] • pc 2 .(4)( 6) CO-to-H2 conversion factor including the contribution of helium (a factor of 1.36 The two values within the parenthesis are the minimum and maximum values within χ 2 = χ 2 min + 1. (6)(7) Gas-to-dust mass ratios.
averages depend on how the annuli are weighted in averaging.Here we calculate the galaxy-average values by applying different weightings on the annulus-based solutions (Table 1) or by applying the fitting methods in Section 4.3 to the whole galaxy.The best averaging scheme may depend on the purpose of usage.For example, to derive the total molecular gas mass over a galaxy, the luminosity-weighted average α CO may be the best.
We first calculate the averages from the annulus-based analysis (Table 1).With the χ 2 method, the averages weighting by luminosity and area are 3.14 (2.06 -4.96) and 3.96 (2.73 -5.81), respectively, with the values within the parenthesis indicating the 68% confidence points.With the ∆logδ GDR method, the average α CO are 3.13 and 3.93 for the luminosity and area weightings, respectively.Since the luminosity is concentrated more towards the inner galaxy, the luminosity-weighted galaxy average tends to weight the inner part of the disk more, and results in a smaller α CO .On the other hand, there is more area in the outer galaxy, and the areaweighted average skews more towards the values of the outer disk, leading to the larger α CO .Nevertheless, all these averages are consistent within an expected factor of 2 uncertainty.
The luminosity and area weighted δ GDR using the χ 2 method are 137 (111, 182) and 167 (140, 209), while the ∆logδ GDR method gives 136 and 164.The δ GDR increases with radius, and the luminosity weighting weighs the inner disk more, while the area weighting weighs the outer disk more.
We also apply the χ 2 and ∆ log δ GDR methods to the whole area (R = 1-7 kpc; excluding the R < 1 kpc area).The results are in the last row of Table 1.The χ 2 fitting gives α CO = 6.03 with ∆χ = 1 range of 5.19-6.90.The minimization of log δ GDR also gives a very similar value of α CO = 5.95.These values are almost a factor of 2 larger than the annulus-based averages, and are likely due to the non-linearity between I CO and the molecular gas surface density Σ H2 , calculated as the difference between the dust-based total gas and atomic gas surface densities (i.e., Σ H2 = δ GDR Σ d − Σ HI ; see Section 6.2).In other words, the assumed linear fit equation (5) may break down on the scales of the entire galactic disk in M83, as most environmental parameters (such as metallicity, cloud populations, dominant gas phasemolecular vs. atomic) can vary across the disk.This is a potential limitation of this method, but in our annulus-based analysis in Section 5.1, this limitation is mitigated.Bolatto et al. (2013) summarized the various measurements of α CO in the MW and recommended α MW = α CO = 4.35M ⊙ pc −2 [K km s −1 ] −1 with uncertainties of ±30% within the inner MW disk (R=1-9 kpc) and a factor of 2 for "normal" galactic disks with solar metallicity.Taking into account our own uncertainties, the derived α CO averages are consistent within the α CO in min for the χ 2 minimization method.The contours are at ∆χ 2 = 1.0, 4.0, and 9.0, which correspond to the confidence levels of 1, 2, and 3σ when they are projected to the αCO or δGDR axis, if the noises in the images and analysis follow a Gaussian distribution.Panel (b) shows the scatter in log δGDR for the ∆ log δGDR minimization method.The white solid line traces the average value of log δGDR derived when the scatter ∆ log δGDR is minimized at a fixed αCO.

Previous Measurements in M83
In this study, we derived α CO and δ GDR , simultaneously in M83.Previous studies of this galaxy derived α CO by fixing δ GDR (Jiao et al. 2021), derived δ GDR by fixing α CO (Foyle et al. 2012), or derived α CO and δ GDR together (Yasuda et al. 2023).1) plotted at the midpoint of each annulus.2021) estimated the CO-to-H 2 and CIto-H 2 conversion factors in M83, also using the dustbased method as we do in this study (eq.1), but with THINGS HI and Nobeyama CO(1-0) data (Walter et al. 2008;Kuno et al. 2007) and the Draine & Li (2007) dust model.They fixed δ GDR , assuming three forms of metallicity-dependent δ GDR , and derived α CO finding an increase from ∼1.0 to ∼3.2 from the central region to their CO detection radius limit of R ∼ 4 kpc (their Figure 2).These appear consistent with the results derived in Section 5, but with caution.As demonstrated in Figure 2, the ratios of the archival HI and CO(1-0) fluxes used in their study over our fluxes are, on average, only 40 and 80%, respectively.These ratios also vary radially and azimuthally, e.g. between the arm and interarm regions.If we simplistically take the average differences and compare them with eq. ( 1), their I CO /Σ HI ratio would be a factor of two larger than ours, which should result in twice smaller α CO than ours.The apparent consistency may indicate that the dust-based method, which both their and our studies employed, has an intrinsic uncertainty of a factor of two.Yasuda et al. (2023) analyzed the same HI and CO(1-0) data as Jiao et al. (2021) and adopted the SED fitting by Casey (2012) for Σ d .They derived both α CO and δ GDR together.They obtained an average of α CO = 0.84 for R ≲ 4.7 kpc (0.60 -1.23, within the 68% confidence level in the α CO -δ GDR plane), and a radial increase in α CO by about a factor of two from 0.34 (0.18 -0.80) to 0.78 (0.43 -1.42) between the inner (R < 1.7 kpc) and outer (> 1.7 kpc, up to R ∼ 4.7 kpc) disk.The radial increase is consistent with the factor of two increase in our results, from 1.82 to 3.91 (Table 1).The values of their α CO are about twice smaller than our results, which can be explained by the flux errors in their HI and CO data (as discussed above).Yasuda et al. (2023) also found δ GDR =73 (59-94) on average, with an increasing trend from 33 (21-66) in R < 1.7 kpc to 72 (55-102) in > 1.7 kpc.These are about three times smaller than results (Table 1).This discrepancy could also be explained by the factor of 2.5 smaller HI flux in their data, which reduces δ GDR by the same factor (see eq. 1).

Jiao et al. (
Foyle et al. ( 2012) fixed α CO and derived δ GDR on spatial scales of ∼ 0.8kpc.They used a modified blackbody fit for Σ d and applied both a constant (MW value) and metallicity-dependent α CO .They used THINGS HI and CO(3-2) emission from the JCMT for the atomic and molecular gas, respectively.From their Figure 10, the radial δ GDR profile looks roughly constant, with δ GDR values of ∼ 80 and ∼ 60 when using a constant and metallicity-dependent α CO (their equation 6) respectively.These values could be corrected upwards, possibly by a factor of two, if the factor of 2 flux error in the THINGS HI data is taken into account (see eq. 5).To that end, the δ GDR derived in this study are consistent with their values.

Possible Origin of Radial Gradients
We find that α CO and δ GDR systematically increase with galactic radius (Section 5.1).Such variations are often attributed to the dependence on metallicity (Wilson 1995;Arimoto et al. 1996;Israel 1997;Leroy et al. 2011), and we discuss their possible metallicitydependence in Section 6.2.1.At the same time, many galactic parameters vary radially (e.g., stellar density and radiation field, gravitational potential, strength of galactic structures, cosmic ray intensity, gas amount, and phase, as well as metallicity), and the α CO and δ GDR could correlate with any of those radially-variable parameters.It is valuable to discuss other possible origins of the α CO and δ GDR variations.In Section 6.2.2, we discuss another possibility: a dependence on the underlying molecular cloud populations.Identifying the exact origin is beyond the scope of this work, and here we aim to investigate a possible origin.
Figure 7 presents a clue to understanding the radial gradient of α CO .This plot shows molecular and atomic gas surface densities vs.I CO .The top panel is for the disk (R =1-7 kpc), and the subsequent panels show the same data, but for each 2 kpc-width annulus separately.If α CO is constant, then we expect a straight line for the molecular gas, with α CO as the slope.The blue points "δ GDR Σ d − Σ HI " represent the molecular gas surface density (equivalent to Σ H2 when δ GDR is calibrated; see eq. 1).Again, α CO = Σ H2 /I CO is the ratio between the y-axis and the x-axis in this plot and is the CO-to-H 2 conversion factor, when the y-intercept is zero.For comparison, the figure also includes the total HI+H 2 gas surface density "δ GDR Σ d " (yellow) and Σ HI (magenta).
The blue points show non-linear (curved) distributions.The black solid lines are 2-d polynomial fits to the blue data points for R =1-7 kpc (top panel) and are the same in all panels for reference.The curved distributions indicate that the slope α CO decreases with increasing I CO .At the same time, these panels also show that the range of I CO that each annulus occupies declines radially, and hence, α CO increases with galactic radius.
This suggests that if the non-linear dependence of α CO on I CO reflects intrinsic gas properties, the range of I CO that each radius covers determines the average α CO in that radius.The non-linearity thus could explain the radial increase of α CO . 5he above discussion relies on an assumption that the other parameters, especially Σ d , are determined accurately.It is also possible that the determination of Σ d could suffer from some systematic errors even after significant efforts, due to the inherent difficulty in constraining dust properties.The relation between I CO and δ GDR Σ d (note that δ GDR is constant in the plots) is already non-linear in a wide range of I CO , which is inherited by the molecular gas surface density "δ GDR Σ d − Σ HI " (blue).

Metallicity Dependence
Metallicity, as well as many other parameters, is likely to change radially, and α CO and δ GDR could depend on any of those parameters (Arimoto et al. 1996 The three include the HI gas surface density, "ΣHI" (magenta), the total gas surface density traced by dust, "δGDRΣ d " (yellow), and the H2 surface density calculated as a difference of the dust-based total gas and HI, "δGDRΣ d − ΣHI" (blue).All the surface densities include the contribution of He mass (a factor of 1.36).The panels correspond to, from top to bottom, the radial ranges of R =1-7 (whole area), 1-3, 2-4, 3-5, 4-6, and 5-7 kpc.For this example, we set δGDR = 137, the luminosity weighted average for δGDR.The plots for δGDR=167 and 221, which correspond to the values for the area weighted average and over the entire disk (R = 1-7 kpc), show similar curved shapes.The black solid line is a 2-d polynomial fit to the blue points in the top panel, and between the x and y axes, it is y = −0.055x 2 + 4.43x − 1.54.Note that the slope of "δGDRΣ d − ΣHI" (blue points) vs. ICO is αCO.
1997; Leroy et al. 2011;Schruba et al. 2012;Rémy-Ruyer et al. 2014;Hunt et al. 2015;Accurso et al. 2017;Madden et al. 2020;Chiang et al. 2021) (see Appendix B).There are also theoretical explanations on the metallicity dependence of α CO (Wolfire et al. 2010;Gong et al. 2020;Hu et al. 2022).Recent works have emphasized the importance of metallicity for molecular gas mass estimations even within a single galaxy.Evans et al. (2022) demonstrated that the α CO variations with metallicity have large effects on star formation efficiency estimates in the MW.Metallicity determinations with optical line emission have systematic uncertainties.Previous studies have measured the metallicity in M83 with roughly consistent results (e.g., Bresolin et al. 2009;Pilyugin et al. 2014;Grasha et al. 2022), but with some differences due to their different calibrations.Hence, we focus primarily on the radial changes (differentials) in α CO due to radial metallicity changes, and secondarily on their absolute values.
The absolute value of α CO suffers from uncertainties in the calibration of metallicity measurements.From R = 2 to 6 kpc, the B07 metallicity calibration gives an increase of α CO = 5.66 to 6.80 for the W95 α COmetallicity relation, and α CO = 8.45 to 11.10 for the A96 relation.The KD02 calibration gives α CO = 2.70 to 3.24 for the W95 relation, and α CO = 2.80 to 3.68 for the A96 relation.Table 3 summarizes these α CO in columns 2-5.The α CO values from the KD02 calibration are consistent with the results of this study, while those from the B07 calibration are about a factor of 2 larger (see Table 1).Figure 8 plots the metallicity-dependent α CO variations with respect to our results.The bottom panel is the same as the top, but shows only the α CO from the KD02 calibration.

Cloud Population Dependence
As discussed in Section 6.2 (see also Figure 7), α CO appears to depend on the range of I CO that each annulus traces, and hence depends on the range of Σ H2 .It therefore seems natural to consider variations of the physical conditions of the gas, i.e., the average properties of the underlying molecular cloud populations, as another potential origin of the α CO variation.Dickman et al. (1986) calculated α CO as an average over a cloud population (see Appendix C).We extend this idea and examine an effect of varying cloud populations on α CO .While our 56.6 ′′ (1.23 kpc) resolution does not resolve individual molecular clouds, it has been known that the cloud population in M83 varies radially: massive (luminous) clouds are populated more in the inner disk than in the outer disk (Freeman et al. 2017, Hirota et al. in prep.).
The dynamical state of molecular clouds is often characterized with the virial parameter α vir , which is the ratio between the kinetic energy K and gravitational potential energy U of cloud (Bertoldi & McKee 1992), A cloud is gravitationally-bound when α vir ≤ 2 (i.e., K ≤ |U |), and is in virial equilibrium when α vir = 1 (2K = |U |, the virial theorem).Hence, gravitationallybound clouds should show α vir =1-2.For a spherical cloud with the mass M , line-of-sight velocity dispersion σ v , and density profile ρ(r) ∝ r −ξ truncated at r = R, we have U = −(3 − ξ)/(5 − 2ξ)(GM 2 /R) and K = 3/2(M σ 2 v ).Thus, with a definition of the virial mass Molecular clouds in the Milky Way (MW) show central concentrations (ξ > 0) but do not have as steep a profile as the isothermal sphere (ξ < 2).If we adopt ξ = 1 as done by Solomon et al. (1987), M vir = 9/2(Rσ 2 v /G). 6The true cloud mass M is often represented by the luminous mass M lum (= α CO L CO ) using the cloud CO luminosity L CO , and hence, 6 We take this definition of M vir so that the virialized state is represented by α vir = 1, since this seems to be a natural definition of the virial parameter.Some studies define M ′ vir ≡ 5Rσ 2 v /G independent of ξ and in turn their virial parameter is α ′ vir ≡ M ′ vir /M = (5/3)((3−ξ)/(5−2ξ))α vir .In their case, the virialized state is α ′ vir =1.00, 1.11, and 1.67 for ξ= 0, 1, and 2, respectively.Solomon et al. (1987) analyzed a population of molecular clouds with masses > 10 4 M ⊙ in the MW and found a power-law relation of M vir = 39L 0.81 CO with M vir and L CO in units of M ⊙ and K km s −1 pc 2 , respectively.Using the same data but conducting a separate analysis, Scoville et al. (1987) obtained a linear relation, M vir = 7.9L CO , for clouds with > 10 5 M ⊙ .In a generalized form, this empirical relation, along with eq. ( 10), is We adopt a CO reference luminosity of L ref = 10 5 K • km s −1 pc 2 .The relations found by Solomon et al. (1987) and Scoville et al. (1987) translate to β = −0.19 and β = 0, respectively.In the case of β ̸ = 0, eq. ( 11) means that α CO and/or α vir change with the cloud luminosity.It has been discussed in two ways, by assuming either α CO or α vir to be constant.Chevance et al. ( 2022) assumed a constant α CO and discussed that α vir , the dynamical state of clouds, changes with cloud luminosity, and hence with mass.Bolatto et al. (2013) adopted α vir = 1 (all clouds are in virial equilibrium) and explained that α CO decreases with cloud luminosity and mass.If α vir is constant, the average α CO in a region should depend on the underlying cloud population.
In order to evaluate the average α CO in each annulus, we adopt the simplest form of commonly used luminosity functions, a truncated power-law cloud luminosity function (Sanders et al. 1985;Solomon et al. 1987;Rosolowsky 2005;Roman-Duval et al. 2010;Colombo et al. 2014) defined between a maximum and minimum cloud luminosity, L max and L min where Φ(L CO ) dL CO is the number of clouds with luminosity L CO in a volume and Φ max is the number density of clouds with L max .We adopt a luminosity function instead of a mass function because a mass function would require assuming a fixed conversion factor.The total luminosity and mass in an annulus are L tot = L CO Φ dL CO / Φ dL CO and M tot = α CO L CO Φ dL CO / Φ dL CO , respectively.With eq. ( 11) and constant α vir , we derive the average ⟨α CO ⟩ = M tot /L tot for the cloud population, Note that α ref is in the units of α CO , and α vir is a dimensionless parameter.L min /L max is small (∼ 10 −2 , Solomon et al. 1987;Scoville et al. 1987), and the last term within the bracket is around 1. Thus, in the limit where L min /L max → 0, ⟨α CO ⟩ depends on the luminosity of the brightest and hence most massive cloud.Hirota et al. (in prep.)generated a catalog of molecular clouds in M83.In their analysis, ∼91% of the total CO luminosity within the mapped area was assigned to molecular clouds.They fitted a truncated powerlaw cloud luminosity function and obtained (γ, L max ) in each annulus (Table 2).Using their cloud catalog, we determine (α ref , β) of eq. ( 11) in Appendix A. These parameters are listed in Table 4 in the appendix.In this catalog, small clouds are not entirely spatially resolved.Hirota et al. (in prep.)uses their peak brightness temperature T peak as a rough measure to separate likelyresolved clouds (T peak > 2 K) from likely-unresolved ones (< 2 K) (Appendix A).Therefore, we obtained (α ref , β) for all clouds (denoted "M83 all ") and likelyresolved clouds with T peak > 2 K ("M83 >2K ").As an example of unambiguously-resolved clouds, we also use the clouds in the MW ("MW"; Solomon et al. 1987).For consistency, we made our own fit to the MW clouds (Appendix A) and use the MW relation for the M83 clouds as well.
With these parameters, we derive the populationaveraged ⟨α CO ⟩ for clouds in M83.The results in all cases are listed in Table 3 (columns 6-11).The ⟨α CO ⟩ obtained in each annulus reproduces the radially increasing trend (Table 3).If we adopt the case of gravitationally-bound clouds (α vir = 2) as an example, the ⟨α CO ⟩ increases radially from α CO = 3.84 to 7.58 using (α ref , β) obtained for all M83 clouds (M83 all ), α CO = 3.12 to 5.61 for those of the likely-resolved clouds in M83 (M83 >2K ), and α CO = 1.74 to 2.62 for those of the MW clouds (MW), from the center (R = 1-3 kpc) to the outer parts (R = 5-7 kpc).These increases are similar to the increase from α CO = 1.0-3.5 to 3.8-6.9derived in this study (Table 1).Figure 8 compares the α CO predicted by this cloud population-based model with the results in this study (Section 5).
We should note that because ⟨α CO ⟩ ∝ 1/α vir (eq.13) the absolute values increase by a factor of two if we adopt α vir = 1 (all clouds are in virial equilibrium; see Table 3 columns 6-8), instead of α vir = 2 (the clouds are merely gravitationally-bound; columns 9-11).This uncertainty stems from our lack of knowledge on the exact dynamical state of molecular clouds.
Here, we tested the idea that α CO may depend on the underlying cloud population and its variations, and showed that the cursory estimation can explain the radial increase of α CO .There are of course some caveats.For example, we did not consider a potential dependence of the cloud luminosity function (eq.12) on metallicity.It is possible that previous studies might have already unintentionally treated the possible dependence on cloud population as a dependence on other parameters, since the variation of cloud population itself could depend on other parameters, e.g., dynamical environment, metallicity, etc.These dependencies make it difficult to isolate the origin of the α CO variations.Put simply, if the properties of molecular clouds vary, e.g., as in eq. ( 11), for any reason, the α CO averaged over a cloud population should also change accordingly.

SUMMARY AND CONCLUSIONS
We present new HI observations in the disk of M83 combining JVLA and single-dish GBT data.In the region we analyzed in this study, the new observations recover a factor of 2 more flux than the distributed THINGS moment 0 map used in previous studies.We compare these HI observations to CO(J=1-0) and a dust surface density map derived with the conversion equation from a 500 µm image (Aniano et al. 2020).We derive (α CO , δ GDR ) simultaneously with the dust-based approach, i.e., using dust masses as a tracer of total gas masses via two methods: (1) χ 2 minimization (planefitting) and (2) log δ GDR minimization.These gave similar results.The dust-based approach requires accurate Σ HI measurements as Σ HI determines the absolute scale of δ GDR , which calibrates the scale for α CO .
We find that α CO and δ GDR both increase as a function of radius, by a factor of ∼ 2-3 from R = 1-3 to 5-7 kpc.These radial changes are consistent with the ones found in previous studies, when the errors in flux calibration in the previous studies are taken into account.As representative values for the galactic disk, we also average the annulus values and derive, e.g., luminosityweighted averages of α CO = 3.14 (2.06, 4.96) and δ GDR = 137 (111,182).These are similar to the values measured in the Milky Way.
the observed nonlinearity between Σ d and Σ H2 (assuming it primarily arises from gas properties, and not potential systematic uncertainties in dust), and the tracing of different ranges of I CO in each annulus.We test this possibility by assuming constant α vir for clouds, and using a truncated power-law cloud luminosity function and an empirical dependence of α CO on cloud CO luminosity.We show that the cloud averaged ⟨α CO ⟩ can explain a factor of 1.5 − 2.0 change in α CO , compared to the factor of 2 − 3 change from the annulus-based fits.Beyond the relative radial changes, the absolute values of α CO predicted by the metallicity and cloud populationdependences are consistent with the results we obtained, if we adopt a particular calibration for the metallicity measurements and α vir value (Figure 8).This is suggestive, but must be revisited in the future, as there is no strong evidence to support the particular calibration nor value.
The exact origin of α CO variations in nearby galaxies remains difficult to disentangle as most environmental parameters vary with galactic radius and correlate with each other.While the dependence on cloud population , with the results in this study from the χ 2 method (black diamonds).The results from the log ∆δGDR minimization method are similar to those derived from the χ 2 method, and are omitted for clarity (see Figure 6).The bottom panels are the same as the top, but only for a smaller αCO range around our results, including the αCO predicted (bottom left) using the metallicity relations under the KD02 calibration, and (bottom right) if the clouds are marginally bound with αvir = 2.The αCO predicted based on this particular calibration for the metallicity gradient or αvir = 2 for clouds are consistent with the αCO derived in this study, within uncertainty.
may in fact embed dependencies on dynamical environment, disk potential, metallicity, etc., it is important to recognize that α CO could depend on cloud population, especially as we move towards more spatially resolved analyses of molecular clouds.Similarly, in terms of cloud population and their properties, we discussed only about the dependence of α CO on cloud luminosity.α CO is also known to vary with cloud temperature T and density ρ (i.e., α CO ∝ ρ 1/2 /T , ee Appendix C) (Dickman et al. 1986;Solomon et al. 1987;Heyer et al. 2001).These parameters could also depend on each other.

A. THE CONVERSION FACTOR AND VIRIAL PARAMETER OF MOLECULAR CLOUDS IN M83
For the discussion in Section 6.2.2, we obtain the relationship among the conversion factor (α CO ), virial parameter (α vir ), and CO luminosity (L CO ) of molecular clouds in M83, using a new cloud catalog (Hirota et al. in prep.).Detailed discussions on this catalog and its limitations (e.g., cloud identification, measurements, and possible biases) are given in Hirota et al. in prep..In particular, the clouds at low L CO are unlikely to be resolved in their observations at a 40 pc resolution, since their counterparts in the MW are typically smaller.
Figure 9 shows a correlation between the α vir α CO and L CO of individual molecular clouds in M83.We adopt α vir α CO = M vir /L CO = 9/2(Rσ 2 v /G)/L CO , where R and σ v are the radius and velocity dispersion of the cloud and G is the gravitational constant.For comparison, we also plot the clouds in the MW (Solomon et al. 1987).
The lines are the fits to all M83 clouds (dashed line), M83 clouds with a peak temperature of T peak > 2 K (solid), and MW clouds (dot-dashed line).Hirota et al. (in prep.)uses the peak brightness temperature T peak as a rough measure to separate likely-resolved clouds (T peak > 2 K) from likely-unresolved ones (< 2 K), although the separation is hardly perfect.They adopted T peak for this separation since it is among the direct observables in their data cube.Cloud radius may be an alternative parameter for this separation; however, while it can be measured, it is uncertain for clouds with sizes close to or smaller than the beam size.Galactic clouds appear to have T peak > 2 K when their surface brightness is averaged over entire cloud areas (e.g., the Taurus molecular cloud; see Goldsmith et al. 2008).Hence, the clouds with T peak < 2 K at the 40 pc resolution are likely suffering from a beam dilution effect (i.e., likely-unresolved).In addition, the majority of the clouds in Hirota et al. in prep.with T peak > 2 K have masses of ≳ 2 × 10 5 M ⊙ when they are calculated from L CO , assuming the Galactic α CO .Using the scaling relations by Solomon et al. (1987), these masses correspond to cloud diameters of ∼40 pc, which are about their resolution.The clouds with T peak < 2 K are typically less massive and are likely smaller.More detailed discussions are in Hirota et al. (in prep.).
We adopt orthogonal distance regression (ODR) for the fitting.The slopes and y-intercepts of the lines are listed in Table 4.We note that Hirota et al. in prep.will derive and discuss these relations more carefully by taking into account possible biases in the catalog., with a range of power-law indices (η = -0.7 to -3.4; Schruba et al. 2012;Hunt et al. 2015;Accurso et al. 2017;Madden et al. 2020) .While there are not many alternative ways to constrain α CO ) and in the MW (Solomon et al. 1987).The linear fits are made with the orthogonal distance regression, for all M83 clouds (dashed line), M83 clouds with Tp > 2 K (solid), and MW clouds (dot-dashed line).The peak brightness temperature T peak is a crude measure to separate resolved clouds (T peak > 2 K) from unresolved ones (< 2 K), but this separation is enough for the purpose of the study presented in the main text.These fits may be subject to further refinement, but are sufficient for the discussion in Section 6.2.2.Hirota et al. (in prep) will present more accurate results with possible biases taken into account.
especially in low-metallicity environments, many of these studies adopt assumptions that could be questioned (e.g., the universality of the star formation process/Kennicutt-Schmidt relation independent of environments, the accuracy of chemical models that assume geometry, density, and temperature of the gas).There are also theoretical works that analyzed simulations to understand the metallicity-dependence of α CO (Gong et al. 2020;Hu et al. 2022), which reproduced the observed α CO -metallicity dependences with power-law indices of -0.80 and -0.70.Hence, instead of adopting a single relation directly, we explore a range of indices encompassing suggested α CO − Z relations in the literature.For this, we adopt the power-law model α CO = α MW Z Z⊙ η with η = 0.0, −1.0, −2.0 and −3.0.We use the metallicity gradient for M83 and calibrations described in section 6.2.1.In Table 5, we show the results from this model as a function of radius.In Figure 10, we compare them with the χ 2 results obtained in section 5.The absolute value of the predicted α CO depends on the metallicity calibration adopted, but η = −2.0 or −3.0 (the steepest slopes among the suggested) with the B07 calibration, may potentially explain the factor of ∼2-3 radial variation in α CO derived in section 5.For η = −2.0 or −3.0 with the B07 calibration to match in absolute value with the χ 2 results, α MW would have to be ∼2 times smaller.
i.e., a dependence of α CO on cloud size (radius) under constant T .With this α CO -R relation, they averaged α CO over the distribution of R in the MW cloud population, and obtained α CO ∼ 4. We extended this idea, applied it to varying cloud populations, and explained the variations of α CO with galactic radius in M83.We used [1] an empirical relation of eq. ( 11) under the assumptions equivalent to [2] and [3] -in our case, we set α vir =1 or 2 and used L CO measured by integrating T over a cloud area (equivalent to the R 2 term) and velocity (σ v ).By starting from eq. ( 11) instead of eq.(C1), we did not need to adopt the assumption [4] of constant T .As a further note, from assumptions [1]-[3], the volume density of cloud is ρ ∝ M vir /R 3 ∝ R 2(γ−1) .With eq.(C1), we can derive independent of γ.
Aniano et al. (2020) used the Draine & Li (2007) dust model, and fitted it to SEDs of a large set of nearby

Figure 2 .
Figure2.Ratio maps between archival HI and CO(1-0) data, and the data used in this study.Left: THINGS/GBTVLA.Middle: LVHIS/GBTVLA.Right: NRO/TP.GBTVLA and TP denote the HI and CO(1-0) data used in this study, while the archival HI and CO(1-0) data are denoted by THINGS and LVHIS, and NRO respectively.The map resolutions are at the lower resolution between the two used in the ratio (see text).Radial contours at R = 4 and R = 7 kpc are overlaid in magenta.In the 11'×11' region for the left two panels (HI data), the ratio of the total fluxes are THINGS/GBTVLA ∼ 0.39 and LVHIS/GBTVLA ∼ 0.82 (see text).For the rightmost panel (CO data), the intensity-weighted average ratio is ∼ 0.83 within R ≤ 4 kpc.

Figure 3 .
Figure 3.The IHI, ICO, and Σ d maps from left to right, shown at their original resolutions (top) and smoothed to the resolution of the CO(1-0) map (56.6 ′′ ) (bottom).Radial contours from R=1-7 kpc are overlaid in white.We exclude the central R < 1 kpc in the analysis, since it suffers from HI absorption.

Figure 4 .Figure 5 .
Figure 4. Comparison of the CO and HI signal-to-noise (S/N) ratios.Top: The S/N maps of ICO (left) and ΣHI (right), with radial contours R = 1 − 7 kpc overplotted in magenta.Only the pixels used in the analyses (i.e., S/N> 3) are shown.We do not show a map in Σ d since we set it constant (S/N=5, see Section 4.2) everywhere.Note the ICO S/N color map does not show its full range of S/N.Bottom: Pixel-to-pixel correlation between the ΣHI and ICO S/N.

Figure 6 .
Figure 6.Radial variations of αCO and δGDR (Table1) plotted at the midpoint of each annulus.

Figure 7 .
Figure7.Plots of three surface densities as a function of ICO.The three include the HI gas surface density, "ΣHI" (magenta), the total gas surface density traced by dust, "δGDRΣ d " (yellow), and the H2 surface density calculated as a difference of the dust-based total gas and HI, "δGDRΣ d − ΣHI" (blue).All the surface densities include the contribution of He mass (a factor of 1.36).The panels correspond to, from top to bottom, the radial ranges of R =1-7 (whole area), 1-3, 2-4, 3-5, 4-6, and 5-7 kpc.For this example, we set δGDR = 137, the luminosity weighted average for δGDR.The plots for δGDR=167 and 221, which correspond to the values for the area weighted average and over the entire disk (R = 1-7 kpc), show similar curved shapes.The black solid line is a 2-d polynomial fit to the blue points in the top panel, and between the x and y axes, it is y = −0.055x 2 + 4.43x − 1.54.Note that the slope of "δGDRΣ d − ΣHI" (blue points) vs. ICO is αCO.

Figure 8 .
Figure 8.Comparison of the αCO predicted by the metallicity dependent relations (left column) and by the cloud populationbased model (right column) as labeled (see text), with the results in this study from the χ 2 method (black diamonds).The results from the log ∆δGDR minimization method are similar to those derived from the χ 2 method, and are omitted for clarity (see Figure6).The bottom panels are the same as the top, but only for a smaller αCO range around our results, including the αCO predicted (bottom left) using the metallicity relations under the KD02 calibration, and (bottom right) if the clouds are marginally bound with αvir = 2.The αCO predicted based on this particular calibration for the metallicity gradient or αvir = 2 for clouds are consistent with the αCO derived in this study, within uncertainty.

Figure 9 .
Figure 9.Comparison of αvirαCO and LCO of individual molecular clouds in M83 (Hirota et al. in prep.) and in the MW(Solomon et al. 1987).The linear fits are made with the orthogonal distance regression, for all M83 clouds (dashed line), M83 clouds with Tp > 2 K (solid), and MW clouds (dot-dashed line).The peak brightness temperature T peak is a crude measure to separate resolved clouds (T peak > 2 K) from unresolved ones (< 2 K), but this separation is enough for the purpose of the study presented in the main text.These fits may be subject to further refinement, but are sufficient for the discussion in Section 6.2.2.Hirota et al. (in prep) will present more accurate results with possible biases taken into account.

Figure 10 .
Figure10.αCO in M83 as predicted by adopting a power-law metallicity relation for a range of indices η (orange, green, pink, and blue) suggested by observational studies in the literature, and the results in this study from the χ 2 method (black diamonds).The left and right panels use the metallicity calibrations from(Bresolin 2007, B07)  and(Kewley & Dopita 2002,  KD02).

Table 1 .
Results using ΣHI, ICO, and Σ d at 56.6 ′′ resolution with the χ 2 and ∆ log δGDR minimization methods.Radial range of aperture for fit in kpc.(2) Number of independent beams in aperture.

Table 2 .
Parameters of the cloud luminosity function.