© ESO, 2014

1. Introduction

The process of phase transition from the important mass reservoir of atomic hydrogen in galactic discs and the dense star-forming molecular phase is a matter of intense research. The Local Group spiral galaxy M 33 is ideally suited for a study of its molecular gas content and of the dependence of the star formation characteristics on the physical conditions across the spatially resolved galactic disk. At a distance of 840 kpc, the Triangulum galaxy is near enough to resolve individual giant molecular clouds (GMCs) with large radio telescopes such as the IRAM 30 meter antenna. M 33 is a small but classical spiral disk with a mass roughly 10% (Corbelli 2003) that of the Milky Way and, unlike Andromeda and our galaxy (the other two local group spirals), it is inclined such that positions and velocities can be determined with no ambiguity. M 33 is chemically young, 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 roughly a factor 2, M 33 represents a stepping stone toward much lower metallicity objects and is the nearest late-type galaxy with a well observed metallicity gradient (Magrini et al. 2007, 2010; Rosolowsky & Simon 2008).

In this work we present a high-sensitivity and high-resolution survey of the CO(2–1) emission from M 33, covering the full optical disk. M 33 has an inclination of i = 56° that makes the position of the clouds in the disk well defined, in contrast to, say, M 31. The observations of the CO(2–1) line presented here have an angular resolution of 12″, corresponding to 49 pc at a distance of 840 kpc for M 33. Single-dish maps do not suffer from missing flux problems often present in interferometric data, an important asset for understanding of the entire molecular phase in galactic disks. Relevant comparison studies are Engargiola et al. (2003), Rosolowsky et al. (2007), and Tosaki et al. (2011) for M 33 at 13″, 15″ and 19.3″ resolutions, PAWS Schinnerer et al. (2013) for M 51, and Kawamura et al. (2009) and Wong et al. (2011) for the Large Magellanic Cloud. None of these works reach a brightness sensitivity or a coverage comparable to the observations presented here.

The data reduction, leading to a data cube, and the different methods of obtaining integrated intensity maps are discussed in detail. These are the main data products of the survey (see, e.g., the integrated intensity map in Fig. 1) and will be made publicly available on the IRAM large program archival webpage¹. This naturally leads to an investigation of the error-beam pickup, something that has not been done in other studies of M 33. The error beam pickup is important to assess the reality of low-level diffuse CO emission; in particular, since we use the H i line information in the reduction process, we would like to determine whether H i-free regions are really CO-free. Finally, we compare the CO and H i line profiles and velocities over the disk, showing that M 33 is a dynamically very cool disk.

Extensive ancillary data are available for M 33, particularly the H i (Gratier et al. 2010b), multiband Spitzer data (Tabatabaei et al. 2007; Verley et al. 2007), and Herschel observations from the HerM33es project (Kramer et al. 2010; Boquien et al. 2011).

Companion papers are in preparation on (i) the N(H₂)/ I_CO(1−0) factor, using the Herschel and H i data and including a possible hidden H₂ component from PDR (photon-dominated regions) layers that are not traced by CO (Gratier et al., in prep.); (ii) a detailed study of the molecular cloud population (Druard et al., in prep.); (iii) a comparison with models (e.g., Blitz & Rosolowsky 2006; Krumholz et al. 2008; Gnedin et al. 2009) that attempt to explain where molecular clouds form (Braine et al., in prep.); (iv) a study of the probability distribution functions of various tracers (dust, CO, H i) (Schneider et al., in prep.). (v) and a comparison of GMC and young stellar cluster positions to deduce a time scale for the GMC life cycle (Corbelli et al., in prep.).

Fig. 1 CO(2–1) integrated intensity map in K km s-1, expressed in the main beam temperature scale and computed as described in Sect. 3.1.3. The contours show H i-poor regions where the H i line does not reach 10 K. The beam size is shown in the lower left corner of the figure. The white ellipse represents a 7.2 kpc radius from the center.

2. Data and observations

2.1. The complete CO(2–1) dataset

M 33 was observed in the CO(2–1) line at 230.538 GHz with the HEterodyne Receiver Array (HERA, Schuster et al. 2004) on the 30 m telescope of the Institut de RadioAstronomie Millimétrique (IRAM) on Pico Veleta in southern Spain². These observations are a follow-up to Gardan et al. (2007) and Gratier et al. (2010b), which covered only part of M 33 (see Fig. 2). Observing parameters (see Table 1) and procedures have been maintained throughout the observing period (2005–2012).

The multibeam 230 GHz receiver HERA consists of a 3 × 3 array of dual polarization pixels (Schuster et al. 2004). HERA was used in the on-the-fly scanning mode in which spectra are taken rapidly while the telescope is moving continuously with a constant angular velocity. Integration times for individual spectra, or dumps, are typically 0.5 s with a scanning speed of 6″ per second. Reference spectra toward positions where no emission is expected are taken regularly. These “OFF” spectra are considerably longer than the individual dumps so the subtraction of the reference adds essentially no noise. We used the WILMA autocorrelator, which measures spectra with a 2.6 km s^-1 (2 MHz) channel spacing. A second autocorrelator (VESPA) with 1.25 MHz channel spacing was also used starting in 2007 (Gardan et al. 2007), but these data have not been used since the original northern region was not observed at that spectral resolution.

M 33 was divided up into fields that could be covered in a convenient time, typically about 30 to 60 min, which is about half the time between pointings. These areas were mostly observed in strips parallel to right-ascension or declination, while some edges of the map were done with slanted strips (see Fig. 2). Between two coverages of the same region, the array was rotated by 90°, 180° or 270° by means of the K-mirror derotator. This keeps the geometry of the 3 × 3 pixel array unchanged but avoids having the same pixels cover exactly the same positions; differences in receiver temperatures from one pixel to another lead to strips of increased noise.

All the data presented here are on main beam temperature scale, converted with a forward efficiency of F_eff = 0.92 and a main beam efficiency of B_eff = 0.56. These are not the same as used by Gardan et al. (2007) and Gratier et al. (2010b) because updated telescope efficiency measurements covering our observing period have recently become available³. The HERA efficiencies published by Schuster et al. (2004) were deduced from measurements made in 2001, prior to the major surface improvements in 2002. The efficiencies and error beam are essentially a function of the telescope but the average coupling of the multibeam receivers is somewhat poorer. Therefore, we used efficiency measurements of the HERA multibeam with measurements of the efficiencies of the single-beam receivers to estimate the difference in coupling and the post-2002 measurements of the single-beam (EMIR) beam efficiency. The forward efficiency is assumed to be the same for the single and multibeam receivers. All values are for a 230 GHz observing frequency.

Thus, in 2001, the HERA B_eff was 0.49 and, for the A230 and B230 receivers respectively 0.52 and 0.50. We then estimated the HERA B_eff = 0.49 / 0.51 ∗ 0.59 ≈ 0.56 where 0.51 is the average for the A230 and B230 receivers and 0.59 is the current value for EMIR at 230 GHz. Our adopted B_eff and F_eff values for the post-2002 period are then different for the observations presented here compared to earlier works (Gardan et al. 2007; Gratier et al. 2010b): it is important to take this change into account.

Fig. 2 Triangulum galaxy M 33. This image shows in plain green the edges of the coverage of the CO(2–1) survey on top of a Herschel 250 μm map (Xilouris et al. 2012) that traces cold dust. The green dotted limits correspond to the previous processed area in Gratier et al. (2010b). The 2′×40′ wide HerM33es strip, in white dots, is also the region where the CO(1–0) transition has been observed with the 30 m telescope. The purple limited areas (white lines in Fig. 1) are the H i-poor regions where we detect no H i above 10 K and inside which we can see that there is no strong 250 μm emission. The white ellipse represents a 7.2 kpc radius from the center.

More than 2 × 10⁷ spectra were acquired in over 400 hours covering a field of view of 55′×40′ amounting to about 400 Gb of data when including both backends.

At the beginning of each observing session, we pointed toward a strong CO(2–1) source located in M 33 at α = 13409.4, δ = +30°49′06′′(J2000), in order to check that the system was correctly tuned and to check the reliability of the incoming spectra before acquiring more data. A spectrum of this position has already been given in Fig. 17 of Gratier et al. (2010b).

2.2. Reduction process

The data reduction was carried out with the GILDAS packages CLASS and GREG⁴. The pipeline used to reduce the CO(2–1) data was adapted from the one already used in the reduction of the first subsets of data (Gardan et al. 2007; Gratier et al. 2010b, see).

The WILMA backend was attached to the receiver so as to center the spectrum in one half of the spectrometer, which is continuous, therefore excluding platforming (a difference in continuum level) within the velocity range of M 33 (roughly −270 to −90 km s^-1). We kept velocities from − 400 km s^-1 to 0 km s^-1 (LSR reference frame) with a channel width of 2.6 km s^-1. Because spectra are taken with a 0.5 sec integration time, the CO(2–1) line is invisible in individual spectra. The first step is to fit a constant continuum level (zeroth-order polynomial baseline) to each spectrum with no baseline window. Thus the root mean square (rms) noise of each spectrum is computed around this fit and then compared to the “theoretical noise” given by the following relation: $\begin{array}{ccc}{\mathit{\sigma }}_{\mathrm{theo}}\mathrm{\approx }\frac{{\mathit{T}}_{\mathrm{sys}}}{\sqrt{{\mathrm{\Delta }}_{\mathit{\nu }}\mathit{t}}}& & \end{array}$(1)with the system temperature T_sys in kelvins, the channel width Δ_ν in Hz, t the integration time in seconds. The spectra presenting a noise level higher than 1.1σ_theo are filtered out; this corresponds to nearly 11% of our dataset.

The disk of M 33 was divided into fields that can be observed in the time between two pointings. Each field was observed multiple times to reduce the rms noise level and thus enhance the signal-to-noise ratio. As a consequence of the many passes, each position in the sky was associated to multiple spectra observed at similar but typically not identical positions. The typical spacing between spectra is about 3″ so the dataset is oversampled. A large table of the individual spectra was made and the spectra were combined to obtain a regularly gridded position-position-velocity data cube, setting the resolution to 12″ and the pixel size to 3″. We used the XY_MAP procedure of GILDAS to convolve the spectra into a data cube. The convolution kernel used in the gridding process is a Gaussian three times the size of the beam FWHM. We can then convert the cube to a sample of regularly gridded spectra with a much lower noise level than the individual initial 0.5-second integration-time spectra.

Up to this stage, no baseline has been subtracted, apart from a constant continuum level. Although most data were taken under good conditions and severe baseline problems were eliminated through a comparison with σ_theo, the data can be improved by fitting a baseline, typically a polynomial. We compared the noise levels obtained from subtracting polynomial baselines of order 1 to 5. When fitting a baseline, a line window is defined and a polynomial is fit to the remaining channels. A low-order polynomial guarantees that no major oscillation will occur within the line window. However, some spectra may require a higher order polynomial to fit the baseline fluctuations. The goal is to subtract the lowest order polynomial that fits the baselines well.

Fitting polynomial baselines of order 1 to 5, we find that the rms noise decreases with increasing order but only very slowly, about 0.2 mK/channel (~1% of the overall noise) for each increase in baseline order. Our preferred baseline is third-order for two reasons: (i) the improvement between orders two and three is greater than between the other increments in baseline order; and (ii) several regions of the cube remain noisy when a first- or second-order polynomial is subtracted, but the problem is solved by using a third-order polynomial and little further improvement appears when going to a fourth- or fifth-order baseline. Thus, we subtracted a third-order baseline from all spectra as the lowest order polynomial allowing the baselines to be consistently well fit. Other tests compared the flux in the moment zero (integrated intensity) maps associated with the cubes reduced with different baseline order fits (see Sect. 3.1.1). For a baseline order below 3, the fit appears to miss some of the CO emission but increasing the order does not increase the total emission.

To fit the baseline, the line is excluded from the fit by using velocity windows based on the H i emission maps (Gratier et al. 2010b). These windows are computed making the assumption that molecular gas mostly forms from atomic gas in M 33 so that H₂ will not be present in M 33 at velocities where there is no H i. Extensive tests were done by using windows based on the M 33 synthetic rotation curve (built from a tilted ring model combined with a deconvolution of the H i Arecibo data and given by Eq. (4) in Corbelli & Schneider 1997) or on the H i peak velocity, each of them with variations in the width of the window. We concluded that the H i-based windows are the ones that yield the lowest mean rms noise level. A further indication that this choice is appropriate can be found in Sect. 4.

The masking method we used required two types of data: the analytic rotation curve (Corbelli & Schneider 1997, Eq. (4)) and the H i data of M 33 already presented and described in Gratier et al. (2010b). We used the VLA H i at 25″×25″ resolution and channel sampling of 1.27 km s^-1 because the noise level is very low – only 2.6 K on average – while keeping an angular resolution comparable to that of GMCs or H i clouds. The H i data are also available via the Centre de Données de Strasbourg (CDS). For each position in the CO data cube, the velocity limits are calculated by locating the peak H i line temperature and then going down to the first ≤0 K channel on each side of this peak. The lower and upper limits of the CO window are the velocities of these H i-free channels. To avoid false detections, only the region within ±30 km s^-1 of the analytical rotation curve is searched for H i. When no emission above 10 K (roughly a 4σ detection) is present, the line window is taken to be 60 km s^-1 centered on the Corbelli & Schneider (1997) rotation curve. The line windows defined in this way are used to subtract baselines from the CO spectra. The peak H i velocities are shown in Fig. 3. This map is not as smooth as an analytical velocity field, but it follows the H i velocity variations much more closely. A large velocity shift (≥10 km s^-1) between two neighboring pixels usually coincides with double-peaked spectra where nearby pixels are dominated by different velocity components. For example, the H i velocity detected around α = 13403, δ = +30°39′35.80′′(J2000) can be −134 km s^-1 for one pixel and −160 km s^-1 for its neighbor. Figure 4 shows the width of the line window over the disk of M 33.

Fig. 3 H i peak velocity map in km s-1. Areas limited in black are regions where there is no 4σ detection of H i. In these regions the map is completed with the analytical rotation curve (see Sect. 2.2). The 25″×25″ beam size for H i data is shown in the lower left corner.

Fig. 4 Map of the window width used in the reduction process, in km s-1, corresponding to the H i based mask used in the baselining process and in the integrated intensity maps computing. Orange regions show areas where no H i is detected, and where the window is thus equal to 60 km s-1. The 25″×25″ beam size for H i data is shown in the lower left corner.

Figure 5 shows the CO velocity field, built in the same way, on top of the H i peak velocity map. The CO velocity, shown within the contoured regions, is only determined where CO is detected above 4σ (contoured in black). The absence of shifts between the H i velocity (outside of contours) and the CO velocity shows that CO and H i are closely linked and are not separated by more than a few kilometers per second. This confirms that we can use the the atomic gas to trace the CO velocity.

Subtracting the baseline yields the spectra from which a final data cube is created, as well as the integrated intensity map. At the end of the data reduction process, we obtain a CO(2–1) cube at 12″ resolution with a 3″ pixel size and with velocity resolution of 2.6 km s^-1 as well as the CO(2–1) integrated line intensity map (see Fig. 1). We also produce a CO(2–1) cube and an integrated intensity map at 25″ resolution following the same procedure for a direct comparison with lower resolution observations such as the CO(1–0) data presented in Sect. 2.4.

2.3. Noise map

Noise maps are computed using emission free channels of the spectra. For this we used the same mask as the one presented in Sect. 2.2, based on H i emission and the rotation curve. The fields with the lower noise correspond to a larger number of coverages and/or better weather conditions during the observations. The radial strip observed with higher sensitivity is the HerM33es strip partially observed in [C ii], [O i], and [N ii] by Herschel. The variations we see in the noise level are emphasized by the color scale but are very weak. The mean rms noise level over the disk up to 7 kpc is 20.33 mK (${\mathit{T}}_{\mathrm{a}}^{\mathrm{\ast }}$) with a fairly homogeneous distribution, which can be seen in Fig. 6. The median of the noise distribution is 20.37 mK.

2.4. CO(1–0) data

The CO(1–0) transition at 115.271 GHz was observed during some poor but not terrible weather periods when CO(2–1) observing was not feasible, for a total of 40 h. The Earth’s atmosphere is less of a problem for observations of CO(1–0). The map is confined to the HerM33es strip, and at a resolution of 0.81 km s^-1, with an average system temperature of 300 K (${\mathit{T}}_{\mathrm{a}}^{\mathrm{\ast }}$). Data reduction was essentially the same as for the CO(2–1) – filtering of poor baselines and construction of a cube at 25″ resolution. The CO(1–0) integrated intensity for maps at 25″ were calculated using H i intensities down to 0 K to define the line windows.

3. Mass and distribution of the molecular gas

3.1. Integrated intensity maps

3.1.1. Determination of line windows

Three methods were tested to compute the integrated CO(2–1) line intensity. The most basic is to use the rotation curve (Corbelli & Schneider 1997) plus a predefined line window to determine the channels to be summed to make the integrated intensity map. The other two methods use the H i line data (Gratier et al. 2010b) as described in Sect. 2.2 but testing two threshold signal levels: 2.6 K (1σ noise level) and 0 K. We used the H i cube to locate the H i peak velocity for emission ≥10 K. The velocity limits are then given by the first channel below a certain threshold (2.6 or 0 K) on each side of the peak. We adopted the latter method (H i line down to 0 K) for the reasons described below.

Figure 3 shows the velocity field of M 33 based on H i peak temperature with the regions where the H i does not reach 10 K (~4σ) indicated by contours. Within these contours, the rotation velocity is assumed to be defined by the Corbelli & Schneider (1997) rotation curve because we consider that the velocity at a threshold below 10 K is less reliable than that of the rotation curve. Figure 4 shows the width of the line window at all positions, determined with the 0 K H i threshold where the peak H i line temperature exceeds 10 K and chosen to be 60 km s^-1 centered on the rotation curve velocity elsewhere. As shown in Fig. 4, there are some regions where a 60 km s^-1 width mask window would not be enough and many others where it would probably be too much.

Fig. 5 Left panel. CO velocity field (black-contoured) on top of H i velocity field map. Inside the black contours the colors represent the CO velocity, while outside they correspond to the H i velocity. The CO velocity is determined by the velocity of the line peak for 4σ detections. This map shows that the CO and the H i peak line temperature are detected at very similar velocities. Right panel. H i-CO velocity shift for the regions where CO velocity is determined. The shift is maximal in NGC604 (α = 13432, δ =+30°47′00′′(J2000)), where the dynamics of the gas is important.

The line windows generated from the H i emission, based on the assumption that CO emission could be present at all velocities at which H i emission was detected, were sometimes very broad, reaching ~100 km s^-1 in a few regions actively forming stars (see Fig. 4). Some of these regions are far from the center so that a simple radial decrease in the assumed linewidth could not be applied to the rotation curve. Using a line window encompassing all of the H i velocities but based on the rotation curve results in large uncertainties in the integrated intensities. Windows larger than necessary make baseline subtraction much more error-prone, allowing baseline fluctuations within the line window. Thus, either very broad windows were used or some velocities at which H i is detected were not included in the line windows if we base the windows on the analytical rotation curve.

Gratier et al. (2010a) used this technique on NGC 6822 with a threshold H i temperature of 10 K. CO emission in M 33 is much stronger than in NGC 6822 so the H i threshold temperatures tested were lower: 2.6 K and 0 K, where 2.6 K is the 1σ noise level in the H i cube. For most of the parts, the H i and rotation velocities agree well since the differences are small (Fig. 3). Double-peaked spectra (see, e.g., spectra in the Appendix of Gratier et al. 2012) and regions where the H i peak was significantly different from the analytical rotation curve were checked by eye to ensure that between two peaks, for example, the emission did not reach the 0 K level. The main potential advantage of using the more restrictive 2.6 K threshold is to frame the CO line more closely, reducing the noise in the integrated intensity map. However, the reduction in window width was small and we felt the potential risk of missing emission due to a dip between peaks or due to negative noise spikes outweighed the slight reduction in noise.

Comparisons between the integrated intensity maps produced with only the rotation curve mask and with the H i mask with thresholds of 0 K and 2.6 K (plus the rotation curve mask for the H i holes) are presented in Fig. 7. For H i based masks, we can see that the 2.6 K-threshold mask apparently misses part of the CO emission as a difference can be seen with the 0 K curve. However, the highest intensities are given when computed with only the rotation curve mask. The telescope error beam picks up emission over a broad range in velocities so increasing the velocity range of the window results in an apparent increase in flux (see Sect. 3.1.2 paragraphs 2 and 3).

The integrated intensity (moment zero) map produced with the H i mask with a 0 K threshold is presented in Fig. 1 and is computed using main beam temperature. The uncertainty in integrated intensity ΔI_CO, Eq. (2), is shown in Fig. 8, based on Figs. 6 and 4. $\mathrm{\Delta }{\mathit{I}}_{\mathrm{CO}}\mathrm{=}\mathrm{rms}\mathrm{\times }\mathit{\delta }{\mathit{\nu }}_{\mathrm{ch}}\sqrt{{\mathit{N}}_{\mathrm{ch}}}\mathrm{=}\mathrm{rms}\mathrm{\times }\sqrt{\mathrm{\Delta }{\mathit{\nu }}_{\mathrm{win}}\mathit{.}\mathit{\delta }{\mathit{\nu }}_{\mathrm{ch}}}$(2)where rms, in kelvin, is the noise value at a given pixel, δν_ch is the channel width (2.6 km s^-1), Δν_win the window width in km s^-1, and ${\mathit{N}}_{\mathrm{ch}}\mathrm{=}\frac{\mathrm{\Delta }{\mathit{\nu }}_{\mathrm{win}}}{\mathit{\delta }{\mathit{\nu }}_{\mathrm{ch}}}$ the number of channels in the window. The mean ΔI_CO value on the disk is 0.22 K km s^-1 and 0.20 K km s^-1 for the regions where H i is detected above 4σ.

Fig. 6 Noise map of the M 33 CO(2–1) data, in kelvin per 2.6 km s-1 channel in antenna temperature. The average rms noise per channel (inside the black contour) is 20.33 mK with a maximum of 50 mK at the map edge.

Fig. 7 Intensities and the corresponding derived mass from the integrated CO(2–1) flux. This figure shows the total cumulative values contained in a given radius disk. The three symbols correspond to different H i windows used in the computation process. The conversion from CO luminosity to H2 mass assumes a N(H2) /ICO(1−0) = 4 × 1020 cm-2/ (K km s-1) and a line ratio of 0.8 (see Sect. 3.3). The error bars are based on the mean ΔICO over the disk as described in the last paragraph of Sect. 3.1.1.

The error bars of Fig. 7 account for statistical uncertainties. The uncertainty is taken to be the average of ΔI_CO (Fig. 8, within the relevant annulus) divided by the square root of the number of lobes $\sqrt{{\mathit{N}}_{\mathrm{pix}}\mathit{/}{\mathit{N}}_{\mathrm{pix}\mathrm{/}\mathrm{lobe}}}$ and expressed as a mass of H₂.

3.1.2. Question of error-beam pickup

The difference between the forward efficiency and the beam efficiency corresponds to power received in the error beams. M 33 is an extended source so some of the flux at any given position will come from error-beam power toward other regions in M 33. The error beams of the IRAM 30 m dish have been significantly reduced in strength since the Greve et al. (1998) publication due to major surface improvement in 2002. Calculations of the post-2002 error beam pattern can be found on the IRAM 30 m webpage in a report by Kramer et al. (2013)⁵. We computed an estimate of the emission from the error beams by adopting a three-Gaussians error beam structure from the report and calculating the pick-up out to the 6% level of the broadest of the error beams (FWHM ~ 800″).

Because the error beams pick up emission up to ~800″ from the pointing center, the velocity of that emission does not necessarily correspond to that of the pointing center and thus to the line window. Among the questions one might ask is whether the difference in flux between the different line windows used in Fig. 7 is due to emission picked up by the error beams at other velocities. This is important because more emission is picked up each time the window is broadened from the 2.6 K to 0 K threshold to a 60 km s^-1 window. If this increase is not due to error beam pickup, then some CO emission must be present at velocities where the H i emission is extremely weak (below 2.6 K) or not detected at all.

To estimate the possibility of significant error beam pickup between our line masks, we calculated the error beam emission cube. We measured the error beam flux as a function of the type of velocity windows we chose. With the analytical rotation curve parameters (a 60 km s^-1 wide window), the error beam flux corresponds to about 2.5 × 10⁷M_⊙ more than the 2.6 K H i level windows and about 1.5 × 10⁷M_⊙ more than the 0 K H i level windows. This is very close to the differences observed in Fig. 7, suggesting that the differences are indeed due to error beam pick-up. Therefore, it is not necessary to invoke CO emission at velocities outside the H i range to explain the differences observed in Fig. 7.

Previous articles on M 33 have not addressed this issue so we do not attempt to subtract this emission in order to keep our maps comparable to those of earlier work. Since we estimate the error-beam pickup from the observed emission, which contains the error-beam emission, this would be ideally accounted for by an iterative process. Such a detailed analysis is beyond the scope of this article. The uncertainty in flux due to error beam pickup is not constant because in regions with a high velocity gradient, the error beam flux is more likely to fall outside the line window and be eliminated with the baseline subtraction.

Fig. 8 rms noise of the integrated intensity map with H i = 0 K mask, in K km s-1.

3.1.3. H i-poor regions

A question that comes naturally is for the regions where the H i emission is weak, not reaching 10 K (~4σ). These areas are marked by contours in Fig. 1 and the line window is set to 60 km s^-1 (Fig. 4). The H i-poor regions of the disk represent about 7% of the CO(2–1) coverage (1% for the first kiloparsec). Is CO emission detected in these regions?

In general, it is supposed that where H i is not present, and gas density/pressure and metallicity are not high enough to cause complete conversion to H₂ (as can be the case in galactic nuclei), CO (H₂) is not expected to be present. This assumption can to some degree be tested as the error pattern should generate weak CO “emission” in the H i holes. In this way, it is possible to estimate the amount of CO by summing the integrated intensity maps corresponding to these H i holes. For these positions, the mask used is based on the analytical rotation curve of M 33. The total integrated signal over these areas corresponds to a value of 4.8 × 10⁶M_⊙. We used the error beam cubes computed in the previous section to compare with this value. The integrated signal due to error beam on the very same regions is 5 × 10⁶M_⊙, which we consider equivalent.

The apparent CO emission we see in H i-poor regions could be entirely due to error beam pickup, such that we have no evidence for CO detection where the peak H i line temperature does not reach 10 K.

3.2. CO($\frac{\mathrm{2}\mathrm{-}\mathrm{1}}{\mathrm{1}\mathrm{-}\mathrm{0}}$) line ratio

Fig. 9 Top: line ratios along the HerM33es strip (triangles) and for annuli of 0.5 kpc width (stars). The dashed line represents the 0.8 value we assume for the ratio in this paper. Middle: radial evolution of the CO(2–1) and CO(1–0) line intensity along the strip. Bottom: radial evolution of the CO(2–1) and CO(1–0) line intensity for the full disk. For all panels, strip CO(1–0) data is described in Sect. 2.4, and disk CO(1–0) data is taken from Rosolowsky et al. (2007). Error bars are derived from statistical and calibration uncertainties and are described in the 2nd paragraph of Sect. 3.2.

Figure 9 shows CO($\frac{\mathrm{2}\mathrm{-}\mathrm{1}}{\mathrm{1}\mathrm{-}\mathrm{0}}$) line ratio as a function of radius both for the whole disk and for the HerM33es strip. The CO(2–1) data used here are from the cube calculated at an angular resolution of 25″. The CO(1–0) data along the strip were presented in Sect. 2.4, and are at a resolution of 25″. The CO(1–0) values for the entire disk are, however, derived from the masses given by Rosolowsky et al. (2007) in Fig. 7 (although their coverage is not full above 4 kpc). Their data were divided by 4.3 to convert them from M_⊙ pc^-2 to cm^-2 K km s^-1 accounting for helium and corrected from an inclination of 52° to the 56° we assume. The line ratios are then calculated for each radial bin for the whole map and for the strip. The middle plot shows that although the behavior of the line ratio of the strip is not like that of the whole disk, the line intensities of the two lines follow each other very closely, showing that the line ratio varies over the disk and not just with radius. The bottom plot of Fig. 9 shows the decrease in the line intensity with the distance from the center.

We estimate the uncertainty in the calibration of the data to be 15% for the CO(2–1), 10% for the CO(1–0) IRAM data, and 15% for the Rosolowsky et al. (2007) CO(1–0) data. Statistical variations have also been taken into account by including the rms scatter in the integrated intensity map. Uncertainties are derived by dividing this scatter by the square root of the number of lobes in the covered area $\sqrt{{\mathit{N}}_{\mathrm{pix}}\mathit{/}{\mathit{N}}_{\mathrm{pix}\mathrm{/}\mathrm{lobe}}}$. Statistical uncertainties for full disk CO(1–0) data are taken from Fig. 7 of Rosolowsky et al. (2007). The statistical uncertainties on these values are generally larger than in Fig. 7 because they account for the whole scatter in I_CO and not ΔI_CO. These values are thus the upper limit of the uncertainty.

While the average line intensities clearly decrease with radius, the CO($\frac{\mathrm{2}\mathrm{-}\mathrm{1}}{\mathrm{1}\mathrm{-}\mathrm{0}}$) ratio does not vary in a regular way. We assume a constant line ratio of 0.8, consistent with Fig. 9, but higher than the 0.73 used by Gratier et al. (2010b). In large spirals, the CO($\frac{\mathrm{2}\mathrm{-}\mathrm{1}}{\mathrm{1}\mathrm{-}\mathrm{0}}$) line ratio decreases radially away from the center (Sawada et al. 2001; Braine et al. 1997). The lack of a decrease may be due to the lower metallicity, and presumably lower average optical depth of the CO lines, in M 33.

3.3. Molecular gas distribution

As with most work on CO observations, we use the CO emission as a proxy for the H₂ column density, assuming a constant N(H₂) /I_CO factor such that N_H2 = I_CO(2−1) × N(H₂) /I_CO(2−1). The H₂ mass is then ${\mathit{M}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{=}{\mathit{I}}_{\mathrm{CO}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{\times }{\left(\frac{{\mathit{I}}_{\mathrm{CO}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{1}\mathrm{)}}}{{\mathit{I}}_{\mathrm{CO}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}}}\right)}^{-1}{\mathit{X}}_{\mathrm{CO}}\frac{\mathrm{2}{\mathrm{m}}_{\mathrm{p}}}{{\mathit{f}}_{\mathrm{mol}}}\mathrm{\Omega D}2$(3)where I_{CO(2 − 1)} is the CO(2–1) intensity on the main beam scale in K km s^-1, $\frac{{\mathit{I}}_{\mathrm{CO}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{1}\mathrm{)}}}{{\mathit{I}}_{\mathrm{CO}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}}}$ is the line ratio studied in Sect. 3.2 and taken equal to 0.8 throughout the disk. We take the CO to H₂ conversion factor to be ${\mathit{X}}_{\mathrm{CO}}\mathrm{=}\frac{{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}}{{\mathit{I}}_{\mathrm{CO}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}}}\mathrm{=}\mathrm{4}\mathrm{\times }{\mathrm{10}}^{\mathrm{20}}\mathrm{c}{\mathrm{m}}^{-2}\mathit{/}\mathrm{(}\mathrm{K}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}\mathrm{)}$, twice the Milky Way value to be consistent with the previous work of Gratier et al. (2010b), assuming an inverse relation between X_CO and the metallicity (Wilson 1995), and based on the far-IR dust emission(Braine et al. 2010). f_mol accounts for the mass of the helium in the molecular gas, with a correction of 37%, and m_p is the mass of the proton. D = 840 kpc is the distance to M 33 and Ω is the beamsize in steradians.

Fig. 10 Radial distribution of the CO-derived H2 (red) and H i (green) mass surface density in M⊙ pc-2. The mass surface density is corrected for inclination and includes helium. The red dots correspond to an exponential fit to the first 7 kpc, yielding an exponential length of 2.1 kpc. The blue dots represent two partial fits: one for the inner galaxy (out to 3.5 kpc), and one for the outer galaxy (3.5 kpc <R< 7.0 kpc). Statistical error bars are derived from the rms scatter over the rings, as described in the 2nd paragraph of Sect. 3.3.

The azimuthally averaged radial distribution of the CO emission, and thus presumably H₂ mass, is shown in Fig. 7. Fitting an exponential disk (solid line) yields a disk scale length of 2.1 ± 0.1 kpc when fitting the whole disk to 7 kpc (see Fig. 10). Fits can be done on the halves of the disk up to 3.5 kpc and beyond. This gives scale lengths of 2.2 ± 0.3 kpc (0–3.5 kpc) and 1.9 ± 0.2 kpc (3.5–7 kpc), similar to the results shown in Gratier et al. (2010b). Error bars correspond only to statistical variations, to be comparable with H i data, and are equal to the rms scatter in I_CO divided by the square root of the number of beams in the covered area $\sqrt{{\mathit{N}}_{\mathrm{pix}}\mathit{/}{\mathit{N}}_{\mathrm{pix}\mathrm{/}\mathrm{lobe}}}$.

While the H₂ surface density, assuming we have chosen an appropriate N(H₂) /I_CO conversion factor, is slightly higher than that of the atomic gas in the center, the H i is dominant beyond the inner kiloparsec.

4. Relationship between atomic and molecular gas

How closely linked are the atomic and molecular components? If molecular clouds are seen as dense clumps embedded in a much warmer and more diffuse neutral atomic medium, then their velocities are not necessarily linked. On the other hand, if the H₂ forms quiescently from the denser atomic clouds, then one would expect that the dispersion between the two components would be very small. The interstellar medium of M 33 is dominated by the atomic component virtually throughout, unlike e.g. M 51, such that the H₂ in M 33 forms from the H i rather than the H i forming from photodissociated H₂.

Molecular gas is not always found where the H i column density is high – there are high N(H i) regions without CO emission just as molecular clouds are sometimes observed in regions of moderate H i column density (e.g., “lonely cloud”, Gardan et al. 2007). This suggests that other factors play a role in provoking the conversion of H i into H₂. Various schemes have been suggested to explain large-scale atomic-to-molecular gas ratios (Blitz & Rosolowsky 2006; Gnedin et al. 2009; Krumholz et al. 2008) in galaxies but not the formation of individual GMCs. It has been suggested that “colliding flows” (Elmegreen 1993; Audit & Hennebelle 2005; Clark et al. 2012) or “colliding clouds” (Motte et al. 2014) may provide the compression required to create H₂ from H i. Presumably, this process would increase linewidths proportionally to the shock velocity and increase the dispersion between the H i and H₂ velocities.

To estimate the velocity dispersion between the H i and CO, we subtracted the H i velocity, as measured by the velocity at the line peak, from the CO spectra. The new cube of CO(2–1) spectra “recentered” to the H i velocity can be used to stack spectra. In this way, we directly obtain the H i-CO velocity dispersion, once a typical CO linewidth has been established. This technique, used in Schruba et al. (2011) and Caldú-Primo et al. (2013) is also very useful to see large-scale radial variations or reveal low-level emission through coherent stacking over large areas.

In addition to recentering the CO line with the H i velocity, we applied this technique to the CO cube using the 4σ detections of CO to determine velocities and recentering the CO with the CO. Similarly, the H i cube was recentered using the velocity of the H i peak temperature.

4.1. Recentered cubes: method

For each line of sight (each spectrum in our regularly gridded data cube), we can associate a velocity based on the H i (or on the analytical rotation curve – see below), and the CO spectrum at this position can be shifted from this velocity to a reference velocity set to be 0 km s^-1 by redefining the velocity axis.

The velocity used to recenter the spectra can be defined in different ways: (i) with the analytical rotation curve (Eq. (4) of Corbelli & Schneider 1997)) given for each point in the disk; (ii) the peak H i channel velocity (computed and used in Sects. 2.2 and 3.1.1) completed with the analytical rotation curve for regions where there is no H i over 10 K (see Fig. 3); or (iii) the first moment of the H i emission calculated in a window within 30 km s^-1 of the rotation curve as follows: $\mathrm{⟨}\mathit{V}\mathrm{⟩}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{channels}}\mathit{Tv}\mathrm{d}\mathit{v}}{{\mathrm{\Sigma }}_{\mathrm{channels}}\mathit{T}\mathrm{d}\mathit{v}}$(4)where T is the temperature of the channel and ν the frequency associated with the channel. Over 90% of the disk of M 33 has H i spectra with S/N> 4 so that the velocity of the peak temperature is well defined. The CO emission covers a much lower fraction of the disk because the vast majority of the lines of sight do not show CO above a 3σ level.

Once all the spectra have been recentered, they can be averaged (i.e., stacked). If the CO emission systematically follows the H i, then we expect to find the CO peak of the stacked spectra at zero velocity. Any systematic difference between H i and CO would create a velocity difference that, while invisible in the individual spectra, might appear in the stacked spectra. The width of the stacked CO spectra comes from the sum of the intrinsic width of the CO spectra, the intrinsic dispersion between the atomic and molecular components, the presence of multipeaked CO, asymmetric profiles, and the error in estimating the H i velocity due to noise, in addition to the small broadening due to the finite channel widths.

Figure 11 shows the stacked spectra corresponding to the full disk coverage of M 33. This means that the entire disk of M 33 is included in these stacked spectra. Each spectrum was computed via a different centering method (see above) applied to the H i and CO spectra: H i peak line temperature, H i first moment, and analytical rotation curve. The “H i peak vel (no baseline)” spectrum comes from a baseline-free CO cube in order to make sure that the subtraction of a third-order baseline does not affect the line wings. Because baseline variations depend on frequency, the effect should statistically cancel each other out when stacked. In the end, the average spectra for the two “H i peak vel” CO data cubes are the same, showing that our reduction process does not introduce any “false” signal in our data cube or filter out anything potentially real.

The H i profiles are wider than those of the CO. This can be explained by the fact that the H i clouds are bigger than the corresponding CO ones and that the velocity dispersion is higher.

Fig. 11 Averaged “recentered” spectra corresponding to the whole disk. Spectra are calculated from H i spectra recentered on the H i peak velocity (light blue) and from the CO cube recentered with the H i peak velocity (black and red), the H i emission moment velocity (green), and the analytical rotation curve velocity (dark blue). The H i spectrum is divided by 2500 to be compared with the CO spectra.

4.2. CO-H I : velocity dispersion and kinematics

This powerful technique also has the advantage of showing the CO velocity dispersion compared to the H i. This parameter is related to the width of our stacked spectra. For example, a CO spectrum with a peak at a velocity different from the H i velocity used during the centering process will broaden the stacked line.

Gaussian profiles are fit to the central channels of each line (≈± 13 km s^-1 around the line center) down to the ~6% level. Table 2 shows the parameters of these fits: ΔV is the full width at half maximum and V_max the velocity shift of the peak of the line, compared to 0, and is computed as V_COpeak − V_{H ipeak}. Even though H i line profiles are clearly not Gaussian (as also seen on the right of Fig. 12), the H i linewidth measurements can be reliable. The uncertainties of the fits are given by the MFIT algorithm of GILDAS.

Fig. 12 Left panel: average CO line profiles for each 1 kpc ring computed from the H i peak velocity recentered CO(2–1) cube. Middle panel: average CO line profiles for each 1 kpc ring computed from the CO peak velocity recentered CO(2–1) cube. Right panel: average H i line profiles for each 1 kpc ring computed from the recentered H i cube. Each spectrum is normalized to unity and separated by adding 0.2.

Figure 11 and Table 2 clearly show that CO linewidths are smallest when using the peak H i velocity for the recentering. Compared to the H i first moment velocity, this shows that H₂ is more likely to form at the velocity of the H i line peak, presumably reflecting the velocity of the highest volume density material rather than at the most representative velocity (the first moment velocity). This phenomenon is seen when examining the cloud catalog for M 33 given by Gratier et al. (2012) where the majority of single CO peaks are associated with the stronger H i peak when there are multiple peaks (e.g., clouds 256, 260, 282, 290) or when two H i peaks match with two CO peaks (e.g., clouds 52, 45, 218). Nonetheless, in a few cases, a single CO peak can be located at the weaker H i peak (e.g., clouds 209, 237, 225) or more rarely between two H i peaks (e.g., cloud 163).

The difference in linewidths (12.5 ± 0.4 km s^-1 versus 15.1 ± 0.1 km s^-1) is beyond what might result from the effect of noise in determining the H i velocity. Recentering with respect to the analytical rotation curve, which is symmetric and determined by fitting the velocity field over the disk, yields a considerably broader line. This is not very surprising because approaching and receding halves of galaxies often do not show identical rotation curves. The comparisons here show that the small-scale wiggles in the rotation curve are followed by both the atomic and molecular components. The width of this “recentered” CO line can be seen as representing the sum of the average CO cloud velocity width and the cloud-cloud dispersion in an axisymmetric potential with no perturbations, although this may not be realistic.

It is also interesting to look at how this CO-H i velocity dispersion evolves with the radius as shown in Figs. 12 for CO and H i. We successively masked the emission outside of concentric rings of 1 kpc width and then summed all the individual spectra in the rings. In these figures, the intensities (normalized to unity) for each ring have been plotted, changing color and adding 0.2 between successive rings. The presence of higher noise in the larger rings can be explained by the fact that there is less CO emission in regions farther away from the center. The parameters of the Gaussian fits parameters are given in Table 2 for the stacked CO spectra.

The middle panel of Fig. 12 shows the CO spectra centered on the CO peak velocity detected above 4σ, which recovers less signal than the other CO-centered spectra (left figure). The parameters of the corresponding fits show the intrinsic dispersion of the CO gas. However, since only a small fraction of M 33 has a CO brightness above 4σ, the “CO recentered CO” spectra cover relatively few lines of sight, whereas recentering with the H i velocity provides a virtually complete coverage and includes the regions with weak CO emission. The right panel shows the H i recentered spectra, stacked in 1 kpc wide rings. The half-power widths of all of these lines are plotted in Fig. 13.

Fig. 13 Evolution of the linewidth along the radius for the CO(2–1) and the H i fits.

A clear decrease in the average CO and H i linewidth with radius can be seen in Fig. 13. At all radii, the CO lines are narrower than the corresponding H i lines despite both lines being recentered with respect to the H i. This means that the disk becomes dynamically cooler with radius (van der Kruit & Shostak 1982) and probably that the CO-H i velocity dispersion decreases as well. While the former is expected, the latter is not necessarily because the outer disk is less efficient (because less gas-rich with longer rotation times) in circularizing velocities than the inner disk. The average lifetime of a GMC is only a small fraction of a rotation period.

5. Probability distribution functions

Probability distribution functions (PDFs) of observables, such as column density and temperature, have been largely used for Galactic cloud studies (e.g., Lombardi et al. 2006; Kainulainen et al. 2009; Schneider et al. 2013), and in numerical modeling (e.g., Federrath et al. 2008). They have also been applied as an analytical tool for studying the intensity and temperature distribution from CO observations in galaxies (M 51, Hughes et al. 2013). We here produce a PDF of column density, derived from our CO observations of M 33. For that, we calculate the H₂ column density from the integrated CO intensity using the conversion factor 4 × 10²⁰ cm^-2/ (K km s^-1) (Gratier et al. 2010b). This is a straightforward approach to compare with Galactic observations (e.g., Schneider et al. 2012; Russeil et al. 2013) and simulations (Federrath & Klessen 2012) to interpret the physical origin of the features observed in a PDF.

All observed pixels from the map shown in Fig. 1 (only excluding the noisy edges seen in Fig. 8) are considered and binned, and normalized to the average column density obtained from the same pixel statistic. The resulting PDF is shown in Fig. 14, expressed as a probability p(η) (see also Federrath et al. (2008) for their definition of a 2D-PDF) with $\begin{array}{ccc}\mathit{\eta }\mathrm{\equiv }\ln \frac{{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}}{\mathrm{⟨}{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{⟩}}\mathrm{·}& & \end{array}$(5)To derive the characteristic properties of the PDF (width, peak, deviations(s) from the log-normal shape), we fit the log-normal function: $\begin{array}{ccc}{\mathit{p}}_{\mathit{\eta }}\hspace{0.17em}\mathrm{d}\mathit{\eta }\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}}\exp \left[\mathrm{-}\frac{\mathrm{(}\mathit{\eta }\mathrm{-}\mathit{\mu }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]\mathrm{d}\mathit{\eta }& & \end{array}$(6)where σ is the dispersion and μ is the mean logarithmic column density. We do this systematically by performing several fits on a grid of parameters for η and μ and then calculate the positive and negative residuals. Because excess is expected to lie above the log-normal form, we select fits with the least negative residuals. We then determine the range of log normality, when the difference between the model and p_η is less than three times the statistical noise in p_η.

Fig. 14 Probability distribution function of H2 column density of M 33. The column density was derived from all pixels in the CO map from Fig. 1, excluding only the noisy edges. The left y-axis gives the normalized probability p(η), the right y-axis the number of pixels per log bin. The upper x-axis is in units of H2 column density, and the lower x-axis is the logarithm of the normalized column density. The green curve indicates the fitted PDF. The dispersion of the fitted PDF is indicated by σ.

The PDF we obtain (Fig. 14) shows a clearly defined log-normal distribution for low column densities with a peak around 0.5 × 10²¹ cm^-2 and excess above ~1.7 × 10²¹ cm^-2. We emphasize that it is not possible to fit a much broader log-normal PDF in order to cover this higher column density range. Assuming that the underlying property indeed has a log-normal distribution⁶, it is important to fit a log-normal PDF to the peak, the column densities left of the peak until the noise limit (indicated in Fig. 14 as a dotted line), and the column densities right of the PDF peak. We then obtain a log-normal PDF with a width of 0.75 (in units of ln(N/ ⟨ N ⟩).

The excess we observe is in the column density range ~1.7 × 10²¹ cm^-2 to ~4 × 10²¹ cm^-2. If a power law is fit to the data in this column density range, then the slope is s ≈ −2.4, typical of the slopes of the power-law tails in Galactic molecular clouds (Russeil et al. 2013; Schneider et al. 2013, 2014) and attributed mainly to self-gravity. Recently, this interpretation was confirmed in analytic collapse models by Girichidis et al. (2014). Other physical processes can affect the shape of the PDF. External compression such as expanding H ii-regions leads to a broader PDF (Schneider et al. 2013) and/or to a “double-peaked” PDF (Schneider et al. 2012; Tremblin et al. 2014). However, these detailed features in the PDF of Galactic clouds are diluted in the galaxy-wide PDF, such as the one for M 33.

The most straightforward explanation is that self-gravity is causing the excess at high column densities in M 33. In their study of M 51 at a comparable spatial resolution, Hughes et al. (2013) did not find such an excess in the CO intensities observed in M 51 (e.g., left panel of their Fig. 2). However, we believe that this is due to their fitting procedure and that a break probably occurs close to log(I_CO) = 1.8. Numerical models (e.g., Girichidis et al. 2014) suggest that the power-law tail for molecular clouds is caused by self-gravity of the clumps and cores located inside the cloud, not necessarily by free-fall contraction of the whole cloud. Numerical hydrodynamical models of galactic disks come to different predictions for the PDF. Wada & Norman (2007) find density PDFs (note that we observe the column density PDF) that resemble more log-normal shapes in their simulations including self-gravity and heating/cooling, while Dobbs (2008) obtain a density PDF with a power-law tail in their SPH simulations that include self-gravity and an adiabatic two-phase gas. They conclude that both agglomeration of small clouds and self-gravity produce GMCs in spiral galaxies.

6. Conclusions

In this paper, the first complete map of the CO(2–1) emission in M 33 up to R_opt (~7 kpc) and the associated integrated intensity map are presented. The average noise level of 20.33 mK per 2.6 km s^-1 channel and the angular resolution is 12″ or 49 pc at the assumed distance of M 33. In addition to the CO(2–1) observations, CO(1–0) has been observed along a radial strip. Our main conclusions are the following.

• 1.

  The total CO(2–1) luminosityis 2.8 × 10⁷ K km s^-1, corresponding to a molecular gas mass of 3.1 × 10⁸M_⊙ assuming a conversion factor of N(H₂) /I_{CO(1 − 0)} = 4 × 10²⁰ cm^-2/ (K km s^-1), twice the classical Milky Way value. The uncertainty in the CO luminosity is dominated by calibration uncertainty of ~15%. The surface density of molecular gas decreases exponentially with radius with a scale length of 2.1 kpc.

• 2.

  Down to a resolution of ≈50 pc (GMC size scale), the velocity dispersion between atomic and molecular gas is very low. The CO(2–1) peak temperature follows the atomic gas peak brightness very closely, suggesting a tight connection between the atomic and molecular components. While shifting the CO(2–1) spectra by the velocity of the H i peak (Sect. 4) and stacking the spectra, the CO line is very narrow (12.4 km s^-1) even when summed over the whole disk. In addition, the linewidths of both components decrease with galactocentric distance, owing to a lower molecular-atomic velocity dispersion and/or intrinsically narrower CO(2–1) lines.

• 3.

  The CO(2–1) emission observed toward H i-poor regions of M 33 is at the level expected from the error beam pickup. Therefore, there is no evidence for molecular gas formation where the H i peak temperature is below 10 K.

• 4.

  The CO($\frac{\mathrm{2}\mathrm{-}\mathrm{1}}{\mathrm{1}\mathrm{-}\mathrm{0}}$) line ratio varies significantly over the disk but not in a regular fashion. The mean value is 0.8, which we apply to the ensemble of our data.

• 5.

  The probability density function of the H₂ column density as traced by the CO emission exhibits a log-normal profile with considerable excess in the high column density regime, presumably owing to the onset of gravitational contraction.

Acknowledgments

We thank the IRAM staff for help provided during the observations and for data reduction. N.S. was supported by the French ANR (Agence Nationale pour la Recherche) project “STARFICH” number ANR-11-BS56-010.

References

All Tables

All Figures

	Fig. 1 CO(2–1) integrated intensity map in K km s-1, expressed in the main beam temperature scale and computed as described in Sect. 3.1.3. The contours show H i-poor regions where the H i line does not reach 10 K. The beam size is shown in the lower left corner of the figure. The white ellipse represents a 7.2 kpc radius from the center.
In the text

Fig. 2 Triangulum galaxy M 33. This image shows in plain green the edges of the coverage of the CO(2–1) survey on top of a Herschel 250 μm map (Xilouris et al. 2012) that traces cold dust. The green dotted limits correspond to the previous processed area in Gratier et al. (2010b). The 2′×40′ wide HerM33es strip, in white dots, is also the region where the CO(1–0) transition has been observed with the 30 m telescope. The purple limited areas (white lines in Fig. 1) are the H i-poor regions where we detect no H i above 10 K and inside which we can see that there is no strong 250 μm emission. The white ellipse represents a 7.2 kpc radius from the center.

In the text

	Fig. 3 H i peak velocity map in km s-1. Areas limited in black are regions where there is no 4σ detection of H i. In these regions the map is completed with the analytical rotation curve (see Sect. 2.2). The 25″×25″ beam size for H i data is shown in the lower left corner.
In the text

	Fig. 4 Map of the window width used in the reduction process, in km s-1, corresponding to the H i based mask used in the baselining process and in the integrated intensity maps computing. Orange regions show areas where no H i is detected, and where the window is thus equal to 60 km s-1. The 25″×25″ beam size for H i data is shown in the lower left corner.
In the text

Fig. 5 Left panel. CO velocity field (black-contoured) on top of H i velocity field map. Inside the black contours the colors represent the CO velocity, while outside they correspond to the H i velocity. The CO velocity is determined by the velocity of the line peak for 4σ detections. This map shows that the CO and the H i peak line temperature are detected at very similar velocities. Right panel. H i-CO velocity shift for the regions where CO velocity is determined. The shift is maximal in NGC604 (α = 13432, δ =+30°47′00′′(J2000)), where the dynamics of the gas is important.

In the text

	Fig. 6 Noise map of the M 33 CO(2–1) data, in kelvin per 2.6 km s-1 channel in antenna temperature. The average rms noise per channel (inside the black contour) is 20.33 mK with a maximum of 50 mK at the map edge.
In the text

Fig. 7 Intensities and the corresponding derived mass from the integrated CO(2–1) flux. This figure shows the total cumulative values contained in a given radius disk. The three symbols correspond to different H i windows used in the computation process. The conversion from CO luminosity to H2 mass assumes a N(H2) /ICO(1−0) = 4 × 1020 cm-2/ (K km s-1) and a line ratio of 0.8 (see Sect. 3.3). The error bars are based on the mean ΔICO over the disk as described in the last paragraph of Sect. 3.1.1.

In the text

	Fig. 8 rms noise of the integrated intensity map with H i = 0 K mask, in K km s-1.
In the text

Fig. 9 Top: line ratios along the HerM33es strip (triangles) and for annuli of 0.5 kpc width (stars). The dashed line represents the 0.8 value we assume for the ratio in this paper. Middle: radial evolution of the CO(2–1) and CO(1–0) line intensity along the strip. Bottom: radial evolution of the CO(2–1) and CO(1–0) line intensity for the full disk. For all panels, strip CO(1–0) data is described in Sect. 2.4, and disk CO(1–0) data is taken from Rosolowsky et al. (2007). Error bars are derived from statistical and calibration uncertainties and are described in the 2nd paragraph of Sect. 3.2.

In the text

Fig. 10 Radial distribution of the CO-derived H2 (red) and H i (green) mass surface density in M⊙ pc-2. The mass surface density is corrected for inclination and includes helium. The red dots correspond to an exponential fit to the first 7 kpc, yielding an exponential length of 2.1 kpc. The blue dots represent two partial fits: one for the inner galaxy (out to 3.5 kpc), and one for the outer galaxy (3.5 kpc <R< 7.0 kpc). Statistical error bars are derived from the rms scatter over the rings, as described in the 2nd paragraph of Sect. 3.3.

In the text

Fig. 11 Averaged “recentered” spectra corresponding to the whole disk. Spectra are calculated from H i spectra recentered on the H i peak velocity (light blue) and from the CO cube recentered with the H i peak velocity (black and red), the H i emission moment velocity (green), and the analytical rotation curve velocity (dark blue). The H i spectrum is divided by 2500 to be compared with the CO spectra.

In the text

Fig. 12 Left panel: average CO line profiles for each 1 kpc ring computed from the H i peak velocity recentered CO(2–1) cube. Middle panel: average CO line profiles for each 1 kpc ring computed from the CO peak velocity recentered CO(2–1) cube. Right panel: average H i line profiles for each 1 kpc ring computed from the recentered H i cube. Each spectrum is normalized to unity and separated by adding 0.2.

In the text

	Fig. 13 Evolution of the linewidth along the radius for the CO(2–1) and the H i fits.
In the text

Fig. 14 Probability distribution function of H2 column density of M 33. The column density was derived from all pixels in the CO map from Fig. 1, excluding only the noisy edges. The left y-axis gives the normalized probability p(η), the right y-axis the number of pixels per log bin. The upper x-axis is in units of H2 column density, and the lower x-axis is the logarithm of the normalized column density. The green curve indicates the fitted PDF. The dispersion of the fitted PDF is indicated by σ.

In the text