Issue 
A&A
Volume 600, April 2017



Article Number  A27  
Number of page(s)  14  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201629300  
Published online  23 March 2017 
The molecular gas mass of M 33
^{1} Laboratoire d’Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, Allée Geoffroy SaintHilaire, 33615 Pessac, France
email: pierre.gratier@ubordeaux.fr
^{2} Institut de Radioastronomie Millimétrique (IRAM), 300 rue de la Piscine, 38406 SaintMartind’ Hères, France
^{3} Department of Physics, 4181 CCIS, University of Alberta, Edmonton, AB T6G 2E1, Canada
^{4} Unidad de Astronomía, Fac. Cs. Básicas, Universidad de Antofagasta, Avda. U. de Antofagasta 02800, Antofagasta, Chile
^{5} Department of Astronomy, University of Massachusetts – Amherst, Amherst, MA 01003, USA
^{6} Observatoire de Paris, LERMA (CNRS: UMR 8112), 61 Av. de l’Observatoire, 75014 Paris, France
^{7} Instituto de Radioastronoma Milimtrica (IRAM), Av. Divina Pastora 7, Nucleo Central, 18012 Granada, Spain
^{8} MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
^{9} Astron. Dept., King Abdulaziz University, PO Box 80203, 21589 Jeddah, Saudi Arabia
^{10} Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, The Netherlands
^{11} CSIRO Astronomy and Space Science, Australia Telescope National Facility, PO Box 76, Epping, NSW 1710, Australia
^{12} Tata Institute of Fundamental Research, Homi Bhabha Road, 400005 Mumbai, India
^{13} Instituto de Astrofísica de Canarias, via Láctea S/N, 38205 La Laguna, Spain
^{14} Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Spain
^{15} KOSMA, I. Physikalisches Institut, Universität zu Köln, Zülpicher Strasse 77, 50937 Köln, Germany
^{16} SRON Netherlands Institute for Space Research, Landleven 12, 9747 AD, Groningen, The Netherlands
^{17} Kapteyn Astronomical Institute, University of Groningen, The Netherlands
Received: 13 July 2016
Accepted: 25 August 2016
Do some environments favor efficient conversion of molecular gas into stars? To answer this, we need to be able to estimate the H_{2} mass. Traditionally, this is done using CO observations and a few assumptions but the Herschel observations which cover the farIR dust spectrum make it possible to estimate the molecular gas mass independently of CO and thus to investigate whether and how the CO traces H_{2}. Previous attempts to derive gas masses from dust emission suffered from biases. Generally, dust surface densities, H i column densities, and CO intensities are used to derive a gastodust ratio (GDR) and the local CO intensity to H_{2} column density ratio (X_{CO}), sometimes allowing for an additional COdark gas component (K_{dark}). We tested earlier methods, revealing degeneracies among the parameters, and then used a sophisticated Bayesian formalism to derive the most likely values for each of the parameters mentioned above as a function of position in the nearby prototypical low metallicity (12 + log (O/H) ~ 8.4) spiral galaxy M 33. The data are from the IRAM Large Program mapping in the CO(2–1) line along with highresolution H i and Herschel dust continuum observations. Solving for GDR, X_{CO}, and K_{dark} in macropixels 500 pc in size, each containing many individual measurements of the CO, H i, and dust emission, we find that (i) allowing for CO dark gas (K_{dark}) significantly improves fits; (ii) K_{dark} decreases with galactocentric distance; (iii) GDR is slightly higher than initially expected and increases with galactocentric distance; (iv) the total amount of dark gas closely follows the radially decreasing CO emission, as might be expected if the dark gas is H_{2} where CO is photodissociated. The total amount of H_{2}, including dark gas, yields an average X_{CO} of twice the galactic value of 2 × 10^{20} cm^{2}/ K km s^{1}, with about 55% of this traced directly through CO. The rather constant fraction of dark gas suggests that there is no large population of diffuse H_{2} clouds (unrelated to GMCs) without CO emission. Unlike in large spirals, we detect no systematic radial trend in X_{CO}, possibly linked to the absence of a radial decrease in CO line ratios.
Key words: ISM: general / galaxies: individual: M 33 / submillimeter: ISM / radio lines: ISM / Local Group / ISM: structure
© ESO, 2017
1. Introduction
Recent work has shown that largescale star formation in galaxies is strongly linked to the molecular gas reservoir, in particular the dense molecular gas, and less so to the total amount of neutral gas (H_{2} + H i) (Kennicutt & Evans 2012; Lada et al. 2012). If we are to understand what affects the relationship between molecular gas and star formation, we need to be able to measure the amount of molecular gas at all positions within the disk of galaxies, ideally down to the scale of individual starforming regions. In lowmetallicity objects, we are very far from such an understanding. The cosmic starformation rate density rises rapidly with redshift (Madau & Dickinson 2014), suggesting that either or both the molecular gas content and the starformation efficiency (mass of stars formed per unit time and unit H_{2} mass) also increase while the fraction of metals decreases with redshift (Combes 2013). This is such that what we learn about local star formation at subsolar metallicities may be useful to better interpret observations of the young universe. The small Local Group spiral galaxy M 33 has a halfsolar metallicity and is near enough (840 kpc, Galleti et al. 2004) to resolve Giant Molecular Clouds (GMCs) and has an inclination (i = 56°) that makes the position of the clouds in the disk well defined (in contrast to e.g. M 31).
Fig. 1
Dust surface density [g/cm^{2}] maps of M 33 at 25′′ resolution: (left) for a constant β = 2, (right) radially variable β (2−1.3) as derived in Tabatabaei et al. (2014). The ellipsed correspond to a galactocentric radius of 7 kpc. 
Description variables.
The whole bright stellar disk of M 33 (up to a radius of ~7 kpc) was recently observed in the CO(2–1) line down to a very low noise level (Druard et al. 2014; Gratier et al. 2010a) using the IRAM 30 m telescope on Pico Veleta. The singledish CO(2–1) data do not suffer from missing flux problems which is an essential asset to the understanding of the entire molecular phase in the galactic disk. M 33 is a chemically young galaxy with a high gas mass fraction and as such represents a different environment in which to study cloud and star formation with respect to the Milky Way. As the average metallicity is subsolar by only a factor of two and the morphology remains that of a rotating disk, M 33 represents a stepping stone towards lower metallicity and less regular objects. Measuring the link between CO and H_{2} is particularly important given the evidence that the conversion of H_{2} into stars becomes more efficient at lower metallicities (Gardan et al. 2007; Gratier et al. 2010a; Druard et al. 2014; Hunt et al. 2015).
With the advent of high resolution dust maps in the Herschel SPIRE and PACS, and Spitzer MIPS and IRAC bands it is possible to determine reliable dust column densities with spatial resolution close to the size of individual GMCs in M 33 (see Kramer et al. 2010; Braine et al. 2010; Xilouris et al. 2012). Under the assumption of local independence of the gastodust ratio (GDR) with respect to the H_{2}/H i fraction, it is possible to determine the local CO intensity to H_{2} column density ratio (X_{CO}).
A simplified global version of such an approach has been applied in Braine et al. (2010, Fig. 4). A more sophisticated method based on maximizing correlation between dust column density structure and that of the gas as derived from H i and CO through an optimal X_{CO} factor has recently been proposed and successfully demonstrated by Leroy et al. (2011) and Sandstrom et al. (2013).
However, these methods have biases and/or degeneracies which will be studied in Sects. 3 and 4, in particular they often do not consider a possible contribution from CO dark molecular gas. In this work, the dust, CO, and H i data covering the disk of M 33 are analyzed using existing these methods along with simulations to quantify bias and degeneracy. A new Bayesian approach is then used and tested in order to calculate the GDR and X_{CO} for any position but also the amount of potential CO dark gas, unseen in H i or CO. All the methods take as a basic assumption that any gas not traced by CO, or potentially optically thick H i, contains dust with similar properties as in the gas traced by CO and H i. This is common to all other studies using dust emission.
2. Data
The CO data are from the recently completed CO(2–1) survey of M 33, which now covers the bright optical disk at high sensitivity (Druard et al. 2014; Gratier et al. 2010b; Gardan et al. 2007). The H i data are from Gratier et al. (2010b). In both lines, we use the datasets produced at 25′′ resolution. The dust surface density is estimated from the Herschel observations (Kramer et al. 2010; Boquien et al. 2011; Xilouris et al. 2012), using the 100, 160, 250, and 350 micron flux densities convolved when necessary to a resolution of 25′′ (see Fig. 1). Thus, the linear spatial resolution at which this study is carried out is 100 pc.
In Fig. 1 (left panel), we show the dust surface density estimated from the SPIRE 250 and 350μm fluxes, using the ratio of these two bands to define the temperature, and assuming a dust opacity of κ = 0.4(ν/ 250 GHz)^{2} cm^{2} per gram of dust (Kruegel & Siebenmorgen 1994), or κ_{350} = 4.7 cm^{2} g^{1} at 350μm.
It is now clear that the dust emissivity index, traditionally designated β, is not necessarily 2 as has generally been assumed. In particular Tabatabaei et al. (2014) have shown that β is variable and lower in M 33 (β = 2−1.3 from the center to the outer disk). However, without being able to calibrate the value κ at the wavelength of interest, it is difficult to be sure of the constant (0.4 above for the dust opacity) as extrapolations have generally assumed β = 2. If the intrinsic β of the dust grains is less than 2, then using β = 2 will result in an underestimate of the temperature and thus an overestimate of the dust mass (compare the two panels of Fig. 1). In this context, a more accurate but more complex means of deriving the dust surface density has been tested. Tabatabaei et al. (2014) find a link between the galactocentric distance and β in M 33 (their Fig. 3). This β(r)is used to derive dust temperatures over the disk of M 33.
In a similar way as in Braine et al. (2010), we then take pixels with H i column density measurements and dust temperatures but no CO emission and compute the median dust crosssection (σ_{dust}) per Hatom: σ_{dust} = S_{ν}/ (B_{ν,T}N_{H}), where S_{ν} is the dust emission and B_{ν,T} the Planck black body emissivity for a frequency ν and a temperature T. At submillimeter wavelengths the dust emission is optically thin. This yields a crosssection per Hatom which naturally varies with radius, much like the metallicity (Magrini et al. 2009). Using σ_{dust}(r), we calculate the total H (i.e., cold, neutral hydrogen gas: H i + H_{2}) column density. The dust opacity is N_{H}σ_{dust} = κΣ_{dust} = S_{ν}/B_{ν,T(β)}, and the dust surface density Σ_{dust} = S_{ν}/ (B_{ν,T(β)}κ). For κ_{350} as above, the dust surface density can be computed for all points in M 33, as shown in Fig. 1 (right panel), such that the difference with respect to Fig. 1 (left panel) is that the temperature is computed with a radially varying β. The values of β are below 2 in M 33 (Tabatabaei et al. 2014) so the temperatures are higher. Since the Planck function B_{ν,T(β)} increases with T, the dust surface density in Fig. 1 (right panel) is lower, particularly in the outer disk where β is lower.
In this work, we only discuss hydrogen content and do not include helium. As helium is present in both the atomic and molecular phases in equal proportion, this does not affect the calculations. As in many other works, we use the term GDR to refer to the hydrogen to dust mass ratio.
3. Dustderived H_{2} versus CO intensity
A simple approach is to take the preexisting map of the H_{2} column density based on Herschel and H i data from Braine et al. (2010) where N(H_{2}) is estimated from the dust and H i emission as N(H_{2}) = (N(H)−N(H i))/2, as in their Fig. 4.
In this case, the variables are X_{CO} and, potentially, a COdark gas column density designated . Figure 2 shows the scatter plots for a sample of three radial bins – 0 kpc < r < 1 kpc, 1 kpc < r < 2 kpc, and 4 kpc <r < 5 kpc. These radii show progressively the transition from an H_{2} dominated ISM, to approximate H i–H_{2} equality between radii 1 and 2 kpc, to the H i dominated outer regions.
Thick red lines show the binning of the scatterplot in 0.5 K wide intervals. The cloud of points are fit by two lines, one assuming N(H_{2}) = X_{CO} × I_{CO} (light red line) and in green. As described by Dickman et al. (1986) a X_{CO} ratio is an average over many different clouds so it cannot be expected to characterize all clouds, or all of our data points.
Fig. 2
Fit of dustderived N(H_{2}) as a function of I_{CO} for data in radial intervals between 0 and 1 kpc (top), 1−2 kpc (middle), and 4−5 kpc (bottom). No cut in intensity has been applied. The color scale indicates the density of points and the thick red histogram shows the N(H_{2}) data averaged in bins of 0.5 K km s^{1}. The thin green line shows an affine fit between N(H_{2}) and I_{CO}; the corresponding fit results are printed in green. The thin red line is a linear fit without an offset; the corresponding fit results are printed in red. Blue cross: average value of the plotted data. 
Figure 2 shows the relationship between the dustderived H_{2} column density and I_{CO} for three radial intervals chosen to represent the inner and outer regions, respectively H_{2} dominated, slightly H i dominated (1−2 kpc), and strongly H i dominated with weak CO emission. From the inner to outer regions, the X_{CO} factor increases, as could be expected given that there is a metallicity gradient and a decline in CO emission (Gratier et al. 2010b) and cloud temperature (Gratier et al. 2012).
The lines without a systematically overestimate the H_{2} mass at moderate and high I_{CO} and both fits overestimate N(H_{2}) at high I_{CO}. There is no physical reason to expect a constant offset () but it appears that there is gas whose dust emission is detected but is not seen in CO – this could be optically thick H i, molecular gas where CO has not formed or is photodissociated, low density H_{2} clouds, or unexpectedly large quantities of ionized gas.
4. LeroySandstrom method
4.1. Prior discussion on the gastodust ratio (GDR)
The GDR is likely wellconstrained by the metallicity, at least for metallicities reasonably close to solar. The solar metallicity is about Z = 0.0142 by mass (Asplund et al. 2009, Sect. 3.1.2). Assuming the standard hydrogentodust mass ratio of 100 (Draine & Li 2007, Table 3), the total gas/dust mass ratio is M(H + He + gas  phasemetals) /M(dust), assuming H and He to be negligible contributors to the dust mass. From Asplund, M(H) = 0.7154 and M(He) = 0.2703, and denoting the gasphasemetalfraction as Z_{gas}, we define the hydrogen gastodust mass ratio as GDR = (0.7154 + 0.0142 Z_{gas})/(0.0142(1−Z_{gas})). Helium adds just under 40% to this number. For GDR = 100, the typical Galactic value, the gasphasemetalfraction Z_{gas} = 0.49 and 51% of the metals are in the dust phase. This value is reasonably robust; for a solar composition, if GDR = 100 ± 20 then 50 ± 10% of the metals are in the gas phase.
What about lower metallicity environments? Since dust condenses from the gas in AGB stellar winds (Gielen et al. 2010) and super nova remnants (Matsuura et al. 2011), one expects that when there is less dust and less metals, the gasphase metal fraction will tend to be higher. At very low metallicities, except for very dense environments, the GDR should be higher than the relation given above due to the difficulty in forming dust grains and mantles sufficiently quickly such that evaporation or destruction processes do not reduce the dust mass (RémyRuyer et al. 2014).
4.2. Method and application to M 33
Developed in Leroy et al. (2011) and later extended and applied to the HERACLES/KINGFISH data in Sandstrom et al. (2013), the idea is that the dust emission can be expressed as the sum of the emission from the atomic and molecular components, implicitly assuming that the contribution from the ionized gas is negligible. The latter assumption is likely appropriate and is also common to other studies (1)where α_{CO} is a surface density conversion factor from I_{CO} to Σ_{H2}. Equating the righthand terms gives us the relation equivalent to Sandstrom et al. (2013, Eq. (3)). In order to allow for some form of CO dark gas, we allow for an additional term, such that the basic equation becomes (2)The procedure is fairly simple: the α_{CO} – K_{dark} space is explored on a regularly spaced grid and, for each couple (α_{CO}, K_{dark}), the dispersion in log (GDR) over the ensemble of pixels is computed. The best fit parameters (α_{CO}, K_{dark}) are chosen as the ones that minimize the log (GDR) dispersion, similar to what was done in Leroy et al. (2011). Sandstrom et al. (2013, their appendix) later studied the influence of different methods to identify the best solution finding robust results over the different methods and settling to using a minimization of the (robust) standard deviation of the logarithm of the GDR. Our maps of M 33 cover an area of several thousand beams. This enables us to look for variations, in particular radial variations, of GDR, α_{CO}, and K_{dark}. Figure 3 shows this space for three radial intervals in M 33, with a minimum computed assuming that a single value for each of the three parameters GDR, α_{CO}, and K_{dark} is appropriate. The best fits are shown as a function of radius in Fig. 4 where the same procedure is applied to concentric elliptical rings sampling 1 kpc in radius.
Fig. 3
Scatter in log (GDR) as a function of X_{CO} and K_{dark}. The color scale and solid white contours indicate the amplitude of the scatter in log (GDR) as measured by the standard deviation for varying X_{CO} and K_{dark} offsets. Radii between 0 and 1 kpc (top), 1−2 kpc (middle), and 4−5 kpc (bottom). The white cross corresponds to the minimum scatter (i.e., best fit). The contours correspond to constant scatter values and give an indication of the uncertainties and degeneracies. The white lines correspond to constant GDR values of 100 (solid), 150 (dashed), 200 (dotted), 250 (dashdotted). 
Fig. 4
Average values for K_{dark}, X_{CO}, and GDR derived for 1 kpc radial bins using the LeroySandstrom method. The black histogram shows results derived with the variable beta (Tabatabaei et al. 2014) prescription and the red used the β = 2 to determine dust temperatures. 
Figure 3 shows that a very broad region of α_{CO} – K_{dark} space yields similar quality fits but that a prior on GDR would help break this degeneracy. The radial behavior shown in Fig. 4 appears somewhat unphysical as the metallicity gradient necessarily yields an increasing GDR and would be expected to also yield X_{CO} increasing with radius.
If we assume that K_{dark} = 0, then we see from Fig. 3 (horizontal line where K_{dark} = 0) that the fit is clearly poorer than the best fit. The same is true for the individual radial bins. The physical interpretation of K_{dark} is far from straightforward. The same procedure has been applied but with a filter only accepting pixels with I_{CO} > 2σ. The result is essentially the same: the slope of the ellipses decreases steadily with radius, showing how difficult it is to measure X_{CO} in the outer regions. The radial variation of the parameters with radius is shown in Fig. 4.
The somewhat more complicated nature of the L–S method (3 parameters: α_{CO}, GDR and K_{dark}) and the broad degeneracies prompted us to explore the effect of noise on typical values (Sect. 4.4) and the recoverability of input parameters using realistic simulated data (Sect. 4.3).
4.3. Recoverability
In order to check the recoverability of the parameters, we have created simulated dust observations with known parameters α_{CO}, GDR and K_{dark}. The I_{CO} and I_{H i} used are the observed values for M 33 to maintain the right correlation between these two quantities. Simulated observations are created following Eq. (6). Noise is then added to each observable quantity I_{CO}, I_{H i} and Σ_{dust}.
We then create the same figures as in Sect. 4.2. The figures are not shown because they are indistinguishable in shape from those in Sect. 4.2 (Figs. 3 and 4). This is not surprising as the data are the same. However, we can add many mock runs of the noise and examine how the biases are affected by differing noise levels and intensity cuts.
Figure 5 shows the result of 200 sets of trial data based on the inner kpc. Input parameters are X_{CO} = 4 × 10^{20} cm^{2}/ (K km s^{1}), GDR = 150 and K_{dark} = 5 M_{⊙}/ pc^{2}, indicated as red lines.
It is immediately clear that the optimization (i.e., the lowest log (GDR) dispersion in Fig. 3) favors lowvalued solutions, with “optimal” values clearly below the input. Even in this high signaltonoise (S/N) region, X_{CO} is underestimated by 25% as is K_{dark} and the GDR by half as much. The GDR is less affected because the H i column density is not modified by K_{dark} or X_{CO} but contributes close to half of the GDR.
Two variants were tested as well. Although a K_{dark} was present in the input parameters, we test the values obtained if it is assumed that K_{dark} = 0, as in Eq. (3) of Sandstrom et al. (2013). In this case, the GDR is underestimated, presumably because more dust is present (as a K_{dark} was injected) than what is seen in H i or CO. Near the center, (Fig. 5) X_{CO} is underestimated (see middle row) but at larger radii the situation is different (cf. next paragraph). If metallicity measurements are reliable, then the GDR is quite constrained (Sect. 4.1). The top row shows the values for X_{CO} and K_{dark} if the true GDR is injected. If a prior on GDR is injected, then we approximately recover X_{CO} and K_{dark}. The dispersion in the histograms is rather small, showing that the results do not depend on the number of realizations.
In the H i dominated outer regions, Fig. 6 shows the same biases as before except that X_{CO} is overestimated when K_{dark} is forced to zero. The prior on GDR again helps recover the input values with reasonable fidelity. There is only weak CO emission at these radii so the constraint on X_{CO} is weak. We therefore made a test excluding values where I_{co}< 2σ. The differences with respect to the input parameters are somewhat less severe (compare Figs. 6 and 7). For the inner kpc, excluding values below 2σ makes no difference because virtually all of the values exceed the threshold.
Fig. 5
Histogram of recovered values for the generative model including noise in all three observables X_{CO}, GDR, and K_{dark}. Bottom row: recovering the 3 parameters. Middle row: recovering only α_{CO}and GDR even though K_{dark} is present in the data. Top row: same as bottom row but imposing the correct value of GDR. This figure is for the central kpc of M 33. Input values are in red. 
Fig. 7
Same as Fig. 6 but only considering pixels where I_{CO} > 2σ. A similar cut for the central radii would show little effect as the CO signal there is strong. 
4.4. Noise effects
In order to evaluate the behavior of the LeroySandstrom (LS) method in the presence of noise, we took typical values of the CO intensity, the H i column, and noise for both, in order to test how the method was affected by noise. We also allow for the presence of CO dark gas, where dark means gas not observed in CO or H i but detected via the emission of the associated dust. Thus, we start with a single value for each of I_{CO}, N(H i) (optically thin assumption), and K_{dark} (dark gas, assumed constant). Assuming a X_{CO} conversion factor, we calculate the gas column density (N(H) = 2 × X_{CO} × I_{CO} + N(H i) + K_{dark}) which we divide by an assumed gastodust ratio (GDR) to obtain a dust surface density Σ_{dust}, similar to what is estimated from analyses of Herschel photometric data (Kramer et al. 2010; Xilouris et al. 2012; Tabatabaei et al. 2014). We then assume a noise level in the same units for each of these quantities and generate 1000 samples (value + Gaussian noise) of each of I_{CO}, N(H i), K_{dark}, and Σ_{dust}. Σ_{dust} after addition of noise is then converted back into a gas surface density using the same GDR. The final step is to test a grid of X_{CO} and K_{dark} values, minimizing the sum of (3)where the quantities are after addition of noise and the sum is over the 1000 samples.
The fiducial model has I_{CO} = 1 ± 0.25 K km s^{1}, N(H i)=8 ± 1 × 10^{20} cm^{2}, and K_{dark} = 1 ± 0.25 × 10^{20} cm^{2} and we assume the uncertainty in the dust surface density is 25%. We inject X_{CO} = 4 × 10^{20} cm^{2}/ (K km s^{1}) in order to calculate Σ_{gas} – this, along with K_{dark}, is what we try to get out of the simulations. The GDR is transparent in that it is used to convert Σ_{gas} into Σ_{dust} but then back into Σ_{gas} after addition of noise so it disappears.
Figure 8 shows the typical degeneracy between the X_{CO} and K_{dark} parameters. The color scale shows the quality of the fit (the lower the better) and contours show the acceptable regions. The black dotted lines indicate the average gastodust ratio for the pixel (i.e., averaged over the 1000 samples for the (I_{CO},K_{dark}) combination). The dotted lines indicate, from left to right, GDRs of 100, 150, 200, and 250. For this example, with I_{CO} = 0.5 ± 0.25 K km s^{1}, the apparently optimal fit is quite far from the input parameters. These values are quite typical of a large number of the pixels in M 33.
Fig. 8
Quality of fit for model with I_{CO} = 1 ± 0.25 K km s^{1}, N(H i) = 8 ± 1 × 10^{20} cm^{2}, and K_{dark} = 1 ± 0.25 × 10^{20} cm^{2}, assuming that the uncertainty in the dust surface density is 25%. Dotted lines represent, from left to right, constant GDR values of 100, 150, 200, 250. The star at X_{CO} = 4 × 10^{20} cm^{2}/ (K km s^{1}) is the input value but the best fit is far from that. 
Figures 9a−f show how the retrieved values of X_{CO} and K_{dark} vary with the CO intensity (before adding noise) and the noise level of the CO observations. The first two figures show the results for N(H i) = 8 ± 1 × 10^{20} cm^{2} and a 25% uncertainty in the dust surface density. The second set of figures shows how the recovered X_{CO} and K_{dark} values depend on the CO intensity and uncertainty in the case where N(H i) = 4 ± 1 × 10^{20} cm^{2}. In the third set, N(H i) = 8 ± 1 × 10^{20} cm^{2} but the uncertainty in the dust (and thus gas) surface density has been reduced to 10%.
Fig. 9
Optimal retrieved values of X_{CO} (left column) and K (right column) as a function of the CO intensity (before adding noise) and the noise level of the CO observations. Top figures: fiducial model. Middle row: fiducial except H i column density reduced to 4 × 10^{20} cm^{2}. Bottom: fiducial except dust uncertainties reduced to 10%. 
The result is striking: in all cases, the X_{CO} conversion factor and the K_{dark} surface density are well recovered for the high CO intensities and small errors but where the intensity or the S/N is lower the recovered X_{CO} decreases systematically and the amount of dark gas increases rapidly. A general tendency is seen towards high K_{dark} and low X_{CO} as the S/N decreases, similar to Fig. 8.
5. Bayesian method
5.1. Principles
This method enables us to take into account the uncertainties in all of the observed quantities and recover the best estimates of the GDR, X_{CO}, and K_{dark} values. This is done in the Bayesian framework of errors in variables.
The generative model is defined as: The above notation means that the quantity observed at pixel i has a Gaussian distribution centered on the true integrated intensity with a dispersion equal to the observational uncertainty σ_{H i,i}. Same for the CO in Eq. (5). The third line states that the true dust surface density Σ_{dust,i} is a function of the true I_{H i,i} and I_{CO,i} and the three model parameters α_{CO}, GDR and K_{dark}. We assume that the H i emission is optically thin such that X_{H i} = 1.823 × 10^{18} cm^{2}/ (K km s^{1}) which converted into units of solar masses per square pc gives α_{H i} = 0.0146 M_{⊙}/ pc^{2}/ ( K km s^{1}). The fourth equation states that the observed dust surface density (left) has a gaussian distribution centered on the true Σ_{dust,i} with dispersion of σ_{Σdust}. We note that the only equality is for the true quantities, not the observations. This method provides an estimate for the true values of Σ_{dust}, I_{CO}, and I_{H i}, as well as the parameters α_{CO}, GDR, K_{dark}.
Because the observations are independent, we can express the likelihood of the parameters knowing the full dataset as the product of the likelihoods of the parameters knowing each individual datapoint. For N observations, (8)where D is the observed dataset , a = α_{H i}/ GDR, b = α_{CO}/ GDR, c = K_{dark}/ GDR. The likelihood is thus the probablility of having an observed set of (i.e., the observed map of Σ_{dust}, I_{CO} and I_{H i}) given a set of values for α_{H i}/ GDR, α_{CO}/ GDR, K_{dark}/ GDR, {}, , and σ_{Σdust}. We know the uncertainty in the I_{CO} and I_{H i} observations (σ_{IH i},σ_{ICO}) and the values are input to the calculation. On the other hand, we do not have a good estimate of the uncertainty in the dust surface density σ_{Σdust} so this is left as a free parameter and becomes an output of the calculation. This σ_{Σdust} will also parameterize Gaussian scatter around the true relationship so σ_{Σdust} may be larger than the measurement error, but accounts for additional scatter in the data (Hogg et al. 2010).
Thus, there are 4 + 2N parameters (α_{H i}/ GDR, α_{CO}/ GDR, K_{dark}/ GDR, σ_{Σdust} and the and for each of the N pixels) to the model and a total of 3N observations ( for each pixel).
Since we are interested in the distribution of the parameters and the likelihood is a probability distribution for the observations, we use the Bayes formula to convert from one to the other. (9)The left hand side of Eq. (9) is the posterior distribution – the distribution function of the parameters given the observations. The first term on the right is the prior distribution of the parameters. In our case, very little information is injected because only unreasonable values are not tested. GDR is varied either uniformly from 0 to 500 or uniformly from 0 to 50 000, the first case enables to check the influence of using a physically justifies prior, namely that the GDR values cannot be higher than 500. α_{CO} is varied from 0 to 30 times M_{⊙}/ pc^{2}/ K km s^{1} and K_{dark} from −10 to 30 M_{⊙}/ pc^{2}. The and parameters are varied between the minimum and the maximum of the observations. The last term is the probability defined above. The posterior distribution is explored using an Monte Carlo Markov Chain (MCMC) code, specifically the EMCEE Python implementation (ForemanMackey et al. 2013) of Affine Invariant Ensemble Sampler described in Goodman & Weare (2010).
The priors can be summarized as: (10)where stands for a uniform distribution between values x_{min} and x_{max}.
5.2. Validation of the Bayesian method
To test the Bayesian method, we simulated a dataset using α_{CO} = 2 × 3.2 M_{⊙}/ pc^{2}/ K km s^{1} (twice the galactic value), GDR = 150, K_{dark} = 10 M_{⊙}/ pc^{2} and σ_{Σdust} = 0.01 M_{⊙}/ pc^{2}. Since we need to input “true” values of I_{CO} and I_{H i} in order to see if we can recover the parameters we inject the observed values of I_{CO} and I_{H i}. These values are then used to create the “true” dust map as per Eq. (6). Since the method starts with observations, we take the simulated observed values to be the real observed values plus noise. Thus, the calculation uses somewhat noisier values than the real data. The tests use datapoints (I_{H i},I_{CO}) characteristic of the inner disk of M 33. Noise is also added to the “true” dust surface density map (created via Eq. (6)).
This model dataset is then used as input into the Bayesian method described in Sect. 5.1. Figure 10 shows the number density of points in the six planes mixing the four parameters α_{CO}, GDR, K_{dark}, and σ_{dust} in grayscale. The orientation of the contours illustrates any degeneracies in the relationship between the parameters. The input parameters to the simulation are shown as solid blue lines. The 4 histograms show the entire set of values for each parameter and the dashed lines show the median and the ± 1σ and ± 2σ. The results contain no obvious bias and are very close to the input parameters. Furthermore, the confidence intervals (± 1σ and ± 2σ) are determined in a selfconsistent way.
Figure 10 is the result of a simulation of the inner kpc of M 33. In the outer parts, the CO emission is very weak and the gas (and thus dust) surface density is dominated by the H i. The Bayesian method as proposed here is not always able to measure the X_{CO} factor where there is little CO emission but an upper limit comes out naturally. On the other hand, GDR can be measured because H i is present in many pixels.
Fig. 10
Test of the Bayesian method. The input values for the simulation are shown as blue lines and these correspond rather well to the peaks of probability distributions determined by the method. The dashed lines indicate the median and ± 1σ and ± 2σ intervals. 
5.3. Application to M 33
M 33 was divided into 324 macropixels measuring 500 pc × 500 pc, each containing 225 independent pixels of H i, CO, and dust data. This size is large enough that the parameters are welldefined but small enough not to be affected by largescale gradients. From the results for the macropixels, it is possible to estimate the radial variation of each parameter. The large number of pixels and macropixels results in an extremely high computation time – about six months CPU using a machine with 12 processors and 128 Gb of memory.
Nearly all (99%) of this time is taken up by the “error in variables” approach (using the full model consisting of all four Eqs. (4) to (7)). Thus, given the prohibitive CPU time, we tested the Bayesian estimation without the error in variables (using a restricted model consisting of only Eqs. (6) and (7)), which runs in a day so we can test different hypotheses. The cases we would like to test are: using the two different dust maps, with different cuts in CO intensity, and with or without limits on the value that GDR can take.
The “error in variables” approach produces slightly lower uncertainties but essentially the same values for the parameters GDR, X_{CO}, and K_{dark}. This can be seen in Fig. 11 which shows the values of K_{dark}, α_{CO}, and GDR for the errorinvariables and the rapid versions. In these simulations, the dust surface density for the variableβ was used, only pixels with CO intensities above 3σ were included, and GDR was allowed to take values between 0 and 500 (5 times the Galactic value).
Therefore, we use the rapid (1 CPUday) computations in the following.
Even with the Bayesian approach, some degeneracy is present. In Fig. 12 (result) and 13 (radial), we show the results for variableβ dust with a 3σ CO cut but without placing a limit on GDR. Both K_{dark} (upper panel) and GDR (lower panel) diverge at large radii, where the CO becomes less of a constraint. This is due to some pixels reaching arbitrarily high GDR values (thousands). If the CO cut is reduced to 0σ, then K_{dark} and GDR diverge at lower radii. The hydrogen mass to dust mass ratio in the Milky Way is about 100, close to 140 if He is included. We thus decided to limit GDR, not allowing it to go above 500 (close to 700 if He is included). Presumably this is well above any true GDR value for a halfsolar metallicity galaxy. The X_{CO} factor is not very affected by the divergence of K_{dark} and GDR although it is difficult to be confident of its value where there is little CO.
Fig. 11
Comparison of rapid and full errorsinvariables Bayesian simulations. The solid line represents the equality of the two quantities and the dashed (resp. dotted) lines are constant ratios of 0.25 (resp. 0.5). 
Fig. 12
Results of Bayesian analysis with a 3σ cut in CO and no cap on GDR. Top row is K_{dark} (left) and uncertainty in K_{dark} (right), both with the color scale to the right and in units of M_{⊙}/ pc^{2}. The second row is X_{CO} (left) and uncertainty in X_{CO} (right), both with the color scale to the right and in units of M_{⊙}/ pc^{2} per K km s^{1}. Bottom row is GDR (left) and uncertainty in GDR (right), both with the color scale to the right. As with the other figures, we have adopted the variablebeta dust surface density shown in Fig. 2. 
Fig. 13
Radial variation of K_{dark}, α_{CO}, and GDR for the simulation with a cut at 3σ for the CO but no cap on GDR. The value is computed using the maps in Fig. 12 and weighting each macropixel according to is area within a given radial annulus. The error bars indicate the dispersion within this ring. A divergence of K_{dark} and GDR can be seen in the outer part. 
Figure 14 shows the maps of the number of measurements used for each of the macropixels for the 0σ and 3σ CO cuts.
Fig. 14
Maps of the number of pixels in each macropixel for the 0σ (left) and 3σ (right) CO cuts. 
Figure 15 is similar to Fig. 4 in that it shows the influence of the choice of the dust emissivity index β on the derived parameters. For the Bayesian method, as for the LS method, the results are consistent for α_{CO} and K_{dark} but differ for GDR with smaller values found for the β = 2 dust maps. This is expected as the β = 2 maps has hight dust surface densities, particularly at higher radii.
Fig. 15
Average values for K_{dark}, X_{CO}, and GDR derived for 0.5 kpc radial bins using the Bayesian method. The black histogram shows results derived with the variable beta of Tabatabaei et al. (2014) and the red uses the standard β = 2 to determine dust temperatures. 
Fig. 16
As for Fig. 12 but with a 500 cap on GDR. The differences can be seen in the outer parts where some of the high GDR pixels from Fig. 12 which were white because they had values over 500. 
Fig. 17
Radial variation of K_{dark}, α_{CO}, and GDR for the simulation with a cut at 3σ for the CO and GDR capped at 500. The value is computed using the maps in Fig. 16 and weighting each macropixel according to is area within a given radial annulus. The error bars indicate the dispersion within this ring. 
Fig. 18
Results of Bayesian analysis with a 0σ cut in CO and GDR limited to 500. Top row is K_{dark} (left) and uncertainty in K_{dark} (right), both with the color scale to the right and in units of M_{⊙}/ pc^{2}. The second row is X_{CO} (left) and uncertainty in X_{CO} (right), both with the color scale to the right and in units of M_{⊙}/ pc^{2} per K km s^{1}. Bottom row is GDR (left) and uncertainty in GDR (right), both with the color scale to the right. As with the other figures, we have adopted the variablebeta dust surface density shown in Fig. 2. 
Fig. 19
Radial variation of K_{dark} , α_{CO}, and GDR for the simulation with a cut at 0σ for the CO and GDR capped at 500. The value is computed using the maps in Fig. 18 and weighting each macropixel according to is area within a given radial annulus. The error bars indicate the dispersion within this ring. The comparison with the preceding figures shows that the cap on GDR is critical to avoid diverging values of GDR and K_{dark}. 
Figures 16 (result) and 17 (radial) show the same as Figs. 12 and 13 but when GDR cannot take values above 500. This essentially avoids finding an optimal result at extremely high GDR and K_{dark}. Where the CO is present in a significant number of pixels (Fig. 14), the limitation (of GDR) is unnecessary but when the equation really only equates GDR and K_{dark} then they are highly degenerate.
Figures 18 (result) and 19 (radial) show the radial variation of K_{dark}, X_{CO}, and GDR for the 0σ and 3σ CO cuts. The similarity shows that when GDR is not allowed to take unphysical values, the CO cut is not critical.
The values of GDR we find in the outer regions using the variableβ approach are actually consistent with the GDR found by Gordon et al. (2014) in the Large Magellanic Cloud. The LMC is a useful comparison as it is only slightly smaller, less massive, and less metallic than M 33 but the LMC is much more irregular.
Several interesting features are present. First of all, even though GDR increases with radius, K_{dark} decreases. This shows that the increase in K_{dark} seen without the limit on GDR was only due to the divergent pixels. The X_{CO} shows no clear radial trend. This is probably unlike large spirals like our own, where a number of works have suggested the X_{CO} increases with radius (Sodroski et al. 1995; Braine et al. 1997), with a particularly low value in the central regions. However, large spirals also show systematic decreases in the CO(2–1)/CO(1–0) ratio whereas M 33 does not (Druard et al. 2014). The value of X_{CO} is only 10% greater than the Galactic value, indicated by a horizontal line in Figs. 13, 17, and 19. This may appear surprising as the X_{CO} factor is expected to increase as the metallicity decreases.
The X_{CO} factor derived here is not directly comparable to the values for the Galactic X_{CO} derived using dust and/or gammaray observations because these calculations did not allow for dark gas and thus attributed all gas (including any CO dark gas) not identified as H i to H_{2} in order to calculate X_{CO}. In order to calculate a comparable ratio, we can add the CO dark gas to the H_{2} column computed as I_{CO} × X_{CO}. While typically modeled as a constant, K_{dark} is not physically expected to be constant as (a) H i is expected to be optically thick only over very small areas and (b) GMC edges, where H is molecular but CO photodissociated, are only expected to be associated with GMCs, which occupy a very small fraction of the disk (Druard et al. 2014). Thus, we can either take the value of K_{dark} derived for the CO detected (0 or 3σ) positions in the macropixel as representative of all positions, or we can assume that the value of K_{dark} derived for the CO detected pixels are only valid for those pixels and assume zero elsewhere. In this way, we may be able to place upper and lower limits to the total X_{CO} values in M 33, including dark gas.
We thus consider Fig. 10 from Druard et al. (2014) and uncorrect for inclination, uncorrect for He, and rescale to a X_{CO} value of 1.1 Galactic – this is equivalent to dividing their values by 1.24. To this, we can add the K_{dark} as computed either in (a) or (b) above.
Expressing the COemitting H_{2} and K_{dark} as surface densities in Fig. 20, it is interesting to note that they are very comparable for a X_{CO} = 1.1X_{gal} where X_{gal} is taken to be 2 × 10^{20} cm^{2}/ K km s^{1}. If we assume that the dark gas is actually molecular gas, then the two columns should be added in order to compare with the Galactic X_{CO} factors based on dust or gammarays. Depending on whether K_{dark} is assumed to be present everywhere at the level derived from the positions respecting the CO threshold or only for those positions, the total X_{CO} (dark H_{2} + COemitting H_{2} divided by I_{CO}) is about twice Galactic with very little radial variation (except for the case where the only pixels with K_{dark} are those above 3σ in CO). The uncertainties increase dramatically beyond 4.5 kpc so we have not been able to derive constraints for the very outer disk of M 33.
Fig. 20
H_{2} surface density derived from CO and K_{dark} derived from the Bayesian analysis. The continuous curve shows Σ_{H2} based on Fig. 10 of Druard et al. (2014) corrected to a X_{CO} factor of 1.1 Galactic and uncorrected for inclination and helium content. The histograms show the CO dark gas surface density. The solid histogram shows K_{dark} as derived assuming that all positions have the same dark column as the positions where CO is detected above the threshold. The dashed and dotted histograms represent K_{dark} assuming that the dark column only is present where CO is detected above the 0σ and 3σ thresholds respectively. 
Although we initially expected K_{dark} to increase (at least with respect to CO) with galactocentric distance as in Pineda et al. (2013, Fig. 15), is not surprising the K_{dark} decreases with radius because the UV field decreases much more quickly than the metallicity. As for the expected increase of X_{CO} with galactocentric radius as is observed in large spirals (Sodroski et al. 1995; Braine et al. 1997), it is not seen in M 33. This was initially a surprise but the constant CO line ratios (2–1/1–0 and 3–2/1–0) support this. In large spirals we see clear decreases in these line ratios (and increases in X_{CO}), but this is not the case in M 33. We did not initially expect K_{dark} to follow the CO column density variation – that came out of the analysis. However, it is natural if the CO dark gas is in the outskirts of GMCs. This implies that there is no large population of diffuse H_{2} clouds (unrelated to GMCs) without CO emission.
Our findings are in apparent disagreement with Pineda et al. (2013). However, Pineda et al. (2013) computed H column densities assuming a constant ratio. Introducing a radially decreasing would at least reduce the difference in our findings. Our findings are in agreement with Mookerjea et al. (2016) who find more CO dark gas near the center than in the BCLMP302 region, although it is very difficult to generalize from a small number of regions. While we describe K_{dark} as decreasing with radius, that is only true in an absolute sense, just like many other quantities decrease with radius (galactocentric distance). Assuming that K_{dark} is not attributable to optically thick H i, a roughly constant mass fraction of molecular material is CO dark, independent of radius. This is in agreement with the findings of Wolfire et al. (2010) where they model the dark gas as the region surrounding molecular clouds where the CO is photodissociated but not the H_{2}. This is in excellent agreement with our observations.
It is worth noting that there is no reason to think that the amount of gas not traced by CO or H i should be constant. Figure 21 shows the dust surface density as a function of the H i column density for 3 macropixels near the center and 3 macropixels between 4 and 5 kpc from the center. Examining the central pixels, it is immediately apparent that the intercept (K_{dark}), varies significantly from one pixel to another, even for neighboring regions. Comparing with the lower panel, we see that K_{dark} tends to be lower in the outskirts although for example, for the brown dots the distribution is rather flat (moderate K_{dark}, infinite GDR) at least when only the H i is plotted. Assuming no CO is present at low H i column density, it is also immediately apparent that there is more dust per unit gas near the center, which is the equivalent of a radially increasing GDR. The low (high) GDR is a factor common to all three pixels at small (large) radii.
Fig. 21
Top: link between H i column density and dust surface density for 3 macropixels near the center of M 33. Each color represents the pixel values of N(H i) and Σ_{dust} for a single macropixel. Bottom: same as above but for 3 macropixels between 4 and 5 kpc from the center. 
6. Conclusions
In order to investigate how GDR, X_{CO}, and K_{dark} vary in M 33, the first step was to take a published estimate of the gas column density N(H)_{dust} based on the Herschel dust observations and plot N(H)_{dust}−N(H i) versus I_{CO}. The systematically positive intercept (Fig. 2) suggests that there is lowcolumn density gas traced by dust but not CO or H i, which we refer to as K_{dark} (Tielens & Hollenbach 1985; Planck Collaboration XIX 2011).
The next step is to construct a map of the dust surface density. Two methods were used – the classical β = 2 dust emissivity (Fig. 1, left panel) and the variableβ (same figure, right panel) developed by Tabatabaei et al. (2014). We adopt the second method because in other subsolar metallicity galaxies (Galliano et al. 2011) the classical approach yields too large a dust mass, presumably due to a change in grain properties with respect to Milky Way dust. Using β = 2 for M 33 also yields a very high dust mass and Tabatabaei et al. (2014) show that β = 2 is a poor approximation for M 33.
We then look for optimal values of GDR, X_{CO}, and K_{dark} to relate the dust surface density to the H i and CO intensities. Except where the S/N is high, major degeneracies are present between these parameters (Fig. 3) such that they all increase (or decrease) simultaneously with similar scatter in log (GDR).
Using simulated data with noise, a similar effect is seen in that the deduced solutions generally have lower GDR, X_{CO}, and K_{dark} than the input values (Figs. 5–7). Setting GDR to the correct (input) value yields reasonably accurate results. Solving only for GDR and X_{CO}, implicitly assuming K_{dark} = 0 when the input value was K_{dark} = 5M_{⊙}/ pc^{2}, yields results for GDR and X_{CO} that strongly depend on the amount of CO with respect to H i. The degeneracies are illustrated by Figs. 8 and 9.
An extremely computationintensive simulation using the Bayesian errorsinvariables approach was used to obtain “true” values of the parameters. Fortunately, a very similar result can be obtained using the Bayesian formalism but without the errorsinvariables approach, as shown from the comparison in Fig. 11. The main difference is the slightly lower uncertainty with the errorsinvariables approach. The degeneracies present using the other methods are (almost) no longer an issue (Fig. 22).
Fig. 22
Search for degeneracies between the GDR, X_{CO}, and K_{dark} in the Bayesian approach. Top panel: X_{CO} (α_{CO}) and K_{dark} as a function of GDR. Bottom panel: link (or absence) between X_{CO} (α_{co}) and K_{dark}. Each point represents a pixel in the maps shown in Fig. 16. 
There is a radial increase in GDR from ~200 near the center to nearly 400 in the outer disk. The X_{CO} ratio remains constant with galactocentric distance, as does the CO(2−1)/CO(1−0) line ratio (Druard et al. 2014) and CO(3−2)/CO(2−1) line ratio (in prep.), unlike what is observed in large spirals. The surface density of dark gas, K_{dark}, decreases from the center (10M_{⊙}/ pc^{2}) to the outer parts (roughly zero) in the same way as the CO emission such that the dark gas represents close to half of the H_{2} assuming that the dark gas is in fact H_{2}. As a result, the ratio of all H_{2} (dark gas plus the H_{2} traced directly by CO), is about twice the local value of 2 × 10^{20} cm^{2}/ K km s^{1}.
Some traces of the degeneracies between K_{dark} and GDR are still present in that some macropixels with little CO find optimal values that are physically unrealistic (typically GDR ~ 5000 with a corresponding divergence of K_{dark}). Limiting the GDR to values less than 500 (5 times the Milky Way value) avoids the problem.
Overall, our results argue for a fairly high GDR in M 33 (GDR ≥ 200), a radially decreasing K_{dark} roughly proportional to the amount of CO emission, and a fairly constant X_{CO} conversion both of the H_{2} directly traced by CO and the total H_{2} content including the dark gas (whose radial distribution is similar to that of the CO).
The results presented here on the link between CO and total molecular gas mass (and/or any optically thick H i) confirm the earlier estimates of the H_{2} mass of M 33. As a result, either the H_{2} is converted into stars more quickly than in large spirals or the starformation rate is overestimated due to for example a change in IMF in this environment.
Acknowledgments
P.G. thanks ERC starting grant (3DICE, grant agreement 336474) for funding during this work. P.G.’s current postdoctoral position is funded by INSU/CNRS.
References
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Boquien, M., Calzetti, D., Combes, F., et al. 2011, AJ, 142, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Braine, J., Brouillet, N., & Baudry, A. 1997, A&A, 318, 19 [NASA ADS] [Google Scholar]
 Braine, J., Gratier, P., Kramer, C., et al. 2010, A&A, 518, L69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Combes, F. 2013, in New Trends in Radio Astronomy in the ALMA Era: The 30th Anniversary of Nobeyama Radio Observatory, eds. R. Kawabe, N. Kuno, & S. Yamamoto, ASP Conf. Ser., 476, 23 [Google Scholar]
 Dickman, R. L., Snell, R. L., & Schloerb, F. P. 1986, ApJ, 309, 326 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T., & Li, A. 2007, ApJ, 657, 810 [NASA ADS] [CrossRef] [Google Scholar]
 Druard, C., Braine, J., Schuster, K. F., et al. 2014, A&A, 567, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Galleti, S., Bellazzini, M., & Ferraro, F. R. 2004, A&A, 423, 925 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Galliano, F., Hony, S., Bernard, J.P., et al. 2011, A&A, 536, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gardan, E., Braine, J., Schuster, K. F., Brouillet, N., & Sievers, A. 2007, A&A, 473, 91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gielen, C., van Winckel, H., Min, M., et al. 2010, A&A, 515, C2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5 [Google Scholar]
 Gordon, K. D., RomanDuval, J., Bot, C., et al. 2014, ApJ, 797, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Gratier, P., Braine, J., RodriguezFernandez, N. J., et al. 2010a, A&A, 512, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gratier, P., Braine, J., RodriguezFernandez, N. J., et al. 2010b, A&A, 522, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gratier, P., Braine, J., RodriguezFernandez, N. J., et al. 2012, A&A, 542, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hogg, D. W., Bovy, J., & Lang, D. 2010, ArXiv eprints [arXiv:1008.4686] [Google Scholar]
 Hunt, L. K., GarcíaBurillo, S., Casasola, V., et al. 2015, A&A, 583, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, C., Buchbender, C., Xilouris, E. M., et al. 2010, A&A, 518, L67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kruegel, E., & Siebenmorgen, R. 1994, A&A, 288, 929 [NASA ADS] [Google Scholar]
 Lada, C. J., Forbrich, J., Lombardi, M., & Alves, J. F. 2012, ApJ, 745, 190 [NASA ADS] [CrossRef] [Google Scholar]
 Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Magrini, L., Stanghellini, L., & Villaver, E. 2009, ApJ, 696, 729 [NASA ADS] [CrossRef] [Google Scholar]
 Matsuura, M., Dwek, E., Meixner, M., et al. 2011, Science, 333, 1258 [NASA ADS] [CrossRef] [Google Scholar]
 Mookerjea, B., Israel, F., Kramer, C., et al. 2016, A&A, 586, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, A&A, 554, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2011, A&A, 536, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 RémyRuyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Sodroski, T. J., Odegard, N., Dwek, E., et al. 1995, ApJ, 452, 262 [NASA ADS] [CrossRef] [Google Scholar]
 Tabatabaei, F. S., Braine, J., Xilouris, E. M., et al. 2014, A&A, 561, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tielens, A. G. G. M., & Hollenbach, D. 1985, ApJ, 291, 722 [NASA ADS] [CrossRef] [Google Scholar]
 Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191 [NASA ADS] [CrossRef] [Google Scholar]
 Xilouris, E. M., Tabatabaei, F. S., Boquien, M., et al. 2012, A&A, 543, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1
Dust surface density [g/cm^{2}] maps of M 33 at 25′′ resolution: (left) for a constant β = 2, (right) radially variable β (2−1.3) as derived in Tabatabaei et al. (2014). The ellipsed correspond to a galactocentric radius of 7 kpc. 

In the text 
Fig. 2
Fit of dustderived N(H_{2}) as a function of I_{CO} for data in radial intervals between 0 and 1 kpc (top), 1−2 kpc (middle), and 4−5 kpc (bottom). No cut in intensity has been applied. The color scale indicates the density of points and the thick red histogram shows the N(H_{2}) data averaged in bins of 0.5 K km s^{1}. The thin green line shows an affine fit between N(H_{2}) and I_{CO}; the corresponding fit results are printed in green. The thin red line is a linear fit without an offset; the corresponding fit results are printed in red. Blue cross: average value of the plotted data. 

In the text 
Fig. 3
Scatter in log (GDR) as a function of X_{CO} and K_{dark}. The color scale and solid white contours indicate the amplitude of the scatter in log (GDR) as measured by the standard deviation for varying X_{CO} and K_{dark} offsets. Radii between 0 and 1 kpc (top), 1−2 kpc (middle), and 4−5 kpc (bottom). The white cross corresponds to the minimum scatter (i.e., best fit). The contours correspond to constant scatter values and give an indication of the uncertainties and degeneracies. The white lines correspond to constant GDR values of 100 (solid), 150 (dashed), 200 (dotted), 250 (dashdotted). 

In the text 
Fig. 4
Average values for K_{dark}, X_{CO}, and GDR derived for 1 kpc radial bins using the LeroySandstrom method. The black histogram shows results derived with the variable beta (Tabatabaei et al. 2014) prescription and the red used the β = 2 to determine dust temperatures. 

In the text 
Fig. 5
Histogram of recovered values for the generative model including noise in all three observables X_{CO}, GDR, and K_{dark}. Bottom row: recovering the 3 parameters. Middle row: recovering only α_{CO}and GDR even though K_{dark} is present in the data. Top row: same as bottom row but imposing the correct value of GDR. This figure is for the central kpc of M 33. Input values are in red. 

In the text 
Fig. 6
Same as Fig. 5 but for 4 kpc <R< 5 kpc. 

In the text 
Fig. 7
Same as Fig. 6 but only considering pixels where I_{CO} > 2σ. A similar cut for the central radii would show little effect as the CO signal there is strong. 

In the text 
Fig. 8
Quality of fit for model with I_{CO} = 1 ± 0.25 K km s^{1}, N(H i) = 8 ± 1 × 10^{20} cm^{2}, and K_{dark} = 1 ± 0.25 × 10^{20} cm^{2}, assuming that the uncertainty in the dust surface density is 25%. Dotted lines represent, from left to right, constant GDR values of 100, 150, 200, 250. The star at X_{CO} = 4 × 10^{20} cm^{2}/ (K km s^{1}) is the input value but the best fit is far from that. 

In the text 
Fig. 9
Optimal retrieved values of X_{CO} (left column) and K (right column) as a function of the CO intensity (before adding noise) and the noise level of the CO observations. Top figures: fiducial model. Middle row: fiducial except H i column density reduced to 4 × 10^{20} cm^{2}. Bottom: fiducial except dust uncertainties reduced to 10%. 

In the text 
Fig. 10
Test of the Bayesian method. The input values for the simulation are shown as blue lines and these correspond rather well to the peaks of probability distributions determined by the method. The dashed lines indicate the median and ± 1σ and ± 2σ intervals. 

In the text 
Fig. 11
Comparison of rapid and full errorsinvariables Bayesian simulations. The solid line represents the equality of the two quantities and the dashed (resp. dotted) lines are constant ratios of 0.25 (resp. 0.5). 

In the text 
Fig. 12
Results of Bayesian analysis with a 3σ cut in CO and no cap on GDR. Top row is K_{dark} (left) and uncertainty in K_{dark} (right), both with the color scale to the right and in units of M_{⊙}/ pc^{2}. The second row is X_{CO} (left) and uncertainty in X_{CO} (right), both with the color scale to the right and in units of M_{⊙}/ pc^{2} per K km s^{1}. Bottom row is GDR (left) and uncertainty in GDR (right), both with the color scale to the right. As with the other figures, we have adopted the variablebeta dust surface density shown in Fig. 2. 

In the text 
Fig. 13
Radial variation of K_{dark}, α_{CO}, and GDR for the simulation with a cut at 3σ for the CO but no cap on GDR. The value is computed using the maps in Fig. 12 and weighting each macropixel according to is area within a given radial annulus. The error bars indicate the dispersion within this ring. A divergence of K_{dark} and GDR can be seen in the outer part. 

In the text 
Fig. 14
Maps of the number of pixels in each macropixel for the 0σ (left) and 3σ (right) CO cuts. 

In the text 
Fig. 15
Average values for K_{dark}, X_{CO}, and GDR derived for 0.5 kpc radial bins using the Bayesian method. The black histogram shows results derived with the variable beta of Tabatabaei et al. (2014) and the red uses the standard β = 2 to determine dust temperatures. 

In the text 
Fig. 16
As for Fig. 12 but with a 500 cap on GDR. The differences can be seen in the outer parts where some of the high GDR pixels from Fig. 12 which were white because they had values over 500. 

In the text 
Fig. 17
Radial variation of K_{dark}, α_{CO}, and GDR for the simulation with a cut at 3σ for the CO and GDR capped at 500. The value is computed using the maps in Fig. 16 and weighting each macropixel according to is area within a given radial annulus. The error bars indicate the dispersion within this ring. 

In the text 
Fig. 18
Results of Bayesian analysis with a 0σ cut in CO and GDR limited to 500. Top row is K_{dark} (left) and uncertainty in K_{dark} (right), both with the color scale to the right and in units of M_{⊙}/ pc^{2}. The second row is X_{CO} (left) and uncertainty in X_{CO} (right), both with the color scale to the right and in units of M_{⊙}/ pc^{2} per K km s^{1}. Bottom row is GDR (left) and uncertainty in GDR (right), both with the color scale to the right. As with the other figures, we have adopted the variablebeta dust surface density shown in Fig. 2. 

In the text 
Fig. 19
Radial variation of K_{dark} , α_{CO}, and GDR for the simulation with a cut at 0σ for the CO and GDR capped at 500. The value is computed using the maps in Fig. 18 and weighting each macropixel according to is area within a given radial annulus. The error bars indicate the dispersion within this ring. The comparison with the preceding figures shows that the cap on GDR is critical to avoid diverging values of GDR and K_{dark}. 

In the text 
Fig. 20
H_{2} surface density derived from CO and K_{dark} derived from the Bayesian analysis. The continuous curve shows Σ_{H2} based on Fig. 10 of Druard et al. (2014) corrected to a X_{CO} factor of 1.1 Galactic and uncorrected for inclination and helium content. The histograms show the CO dark gas surface density. The solid histogram shows K_{dark} as derived assuming that all positions have the same dark column as the positions where CO is detected above the threshold. The dashed and dotted histograms represent K_{dark} assuming that the dark column only is present where CO is detected above the 0σ and 3σ thresholds respectively. 

In the text 
Fig. 21
Top: link between H i column density and dust surface density for 3 macropixels near the center of M 33. Each color represents the pixel values of N(H i) and Σ_{dust} for a single macropixel. Bottom: same as above but for 3 macropixels between 4 and 5 kpc from the center. 

In the text 
Fig. 22
Search for degeneracies between the GDR, X_{CO}, and K_{dark} in the Bayesian approach. Top panel: X_{CO} (α_{CO}) and K_{dark} as a function of GDR. Bottom panel: link (or absence) between X_{CO} (α_{co}) and K_{dark}. Each point represents a pixel in the maps shown in Fig. 16. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.