© ESO 2018

1 Introduction

Recent years have seen a sharply increasing number of studies of molecular clouds in external galaxies (Wilson & Scoville 1990; Engargiola et al. 2003; Rosolowsky et al. 2003; Kawamura et al. 2009; Gratier et al. 2010a, 2012; Donovan Meyer et al. 2013; Schinnerer et al. 2013; Druard et al. 2014; Corbelli et al. 2017; Freeman et al. 2017). All of these authors wanted to know whether we could apply knowledge of molecular clouds in the Milky Way to other galaxies and other environments. Of particular interest was the metallicity but also the morphology of the host galaxy.

M 33 is a small spiral galaxy in the Local Group with a thin disk and a metallicity about half that of the solar neighborhood. It is at about 840 kpc (Galleti et al. 2004) and rather ideally inclined with an inclination of 56°, such that cloud velocities are straightforwardly separated and the structure of the disk is plainly visible. The full disk of M 33 was observed in the CO(2–1) line with the IRAM 30meter radiotelescope (see details in Gardan et al. 2007; Gratier et al. 2010b; Druard et al. 2014), yielding a resolution of 12″ or 49 pc. The CPROPS software Rosolowsky & Leroy (2006) was used to extract a 566 cloud sample from the CO data cube (see Gratier et al. 2012; Druard 2014, for details) and published in Corbelli et al. (2017). The 566 clouds identified in M 33 were classified in terms of their star formation by (Corbelli et al. 2017; following Gratier et al. 2012) and here we use this classification to compare star-forming clouds from those with no observed star formation. Figure 1 (left) shows the CO(2–1) integrated intensity map with the cloud contours superposed.

We use the cloud sample to investigate trends in linewidths in M 33 and in the broader context of Local Group galaxies and M 51. We specifically look for a connection between galaxy metallicity and the cloud size-linewidth relation. Some of the work based on a partial cloud sample and presented in Gratier et al. (2012) is redone with the full cloud sample. The link between cloud mass spectra and star formation is examined and Herschel SPIRE photometry data are used to estimate dust temperatures for the star-forming and non-star-forming clouds.

In this work we study not only the general properties of the clouds but also look for velocity gradients across the clouds as signs of possible cloud rotation. Assuming that the medium out of which clouds form is larger than the molecular cloud, angular momentum conservation should result in detectable rotation velocities (see e.g. Rosolowsky et al. 2003). The first to measure cloud velocity gradients, in order to trace cloud rotation, were Kutner et al. (1977) who estimated a large-scale velocity gradient of 0.135 km s⁻¹ pc⁻¹ opposite to galactic rotation (retrograde). These authors argue that the velocity gradient is due to rotation. In the review by Blitz (1993), it is concluded that velocity gradients are probably the result of rotation but that individual clouds (they identify W3) could be exceptions. A more detailed discussion of Milky Way cloud velocity gradients can be found in Phillips (1999). Our observations are GMC-scale, as opposed to much smaller scales where rotating disks are apparent, so influences other than rotation may be present in our measures (see Sect. 3).

Rosolowsky et al. (2003) were the first to look at velocity gradients in an external galaxy. Extragalactic work is quite complementary to the Galactic observations as the biases are not at all the same although the spatial resolution is much poorer. They identified 45 molecular clouds in M 33 and found velocity gradients to be very low, with nearly half in the retrograde direction. Imara et al. (2011) looked at the sample and suggested that in fact the clouds may not be rotating. We use our sample of 566 clouds to re-evaluate the question of cloud rotation.

2 Cloud properties

In addition to estimating 3D cloud boundaries, CPROPS generates information such as deconvolved sizes and cloud luminosities. Outside the CPROPS calculations, we estimate cloud linewidths as in Gratier et al. (2012) by fitting a Gaussian line profile to the average cloud profile. Cloud masses are estimated by assuming a constant N(H₂)∕I_CO(1−0) factor of 4 × 10²⁰ cm⁻²∕(K km s⁻¹) and a CO($\frac{2-1}{1-0}$) line ratio of 0.8. These values have been validated by Druard et al. (2014), Gratier et al. (2017), and Braine et al. (2010b).

CPROPS was able to estimate deconvolved radii for 449 out of 566 clouds. In order to provide radii for the remaining 20%, we use the link between non-deconvolved and deconvolved radii for the 449 clouds to extrapolate to determine deconvolved radii for the remaining clouds. This enables the use of the whole sample but does not change the results. Figure 2 shows the relation we fit with the 449 clouds and the values calculated for the others. From now on, we use the whole set of deconvolved cloud radii when radii are used. Over the whole disk, about half the CO emission comes from identified clouds. The fraction decreases with galactocentric distance: 2/3 of the CO emission is in clouds in the inner disk and 1/3 in the outer disk (Druard 2014).

Fig. 2 Link between CPROPS undeconvolved cloud radii (x-axis) and deconvolved radii for the 80% of the clouds with well-defined deconvolved radii (black symbols). The red triangles show the radii attributed to the clouds for which CPROPS could not estimate radii – the radii follow the envelope defined by the CPROPS values down to a constant value below which we have no confidence. The goal here is to avoid generating clouds which are inappropriately small while still attributing reasonable radii.

2.1 Size-linewidth relation

Figure 3 shows the size-linewidth relation for the M 33 clouds as compared to other nearby galaxies for which similar data are available, including our own Galaxy. The Solomon et al. (1987) relation for the Milky Way is shown as a line. M 51 data come from Colombo et al. (2014), the LMC region from Hughes et al. (2010), and the NGC 6822 data from Gratier et al. (2010a).

Two things are apparent from Fig. 3. First, the size-linewidth relation is very weak in M 33, NGC 6822, and M 51, but apparently strong in the Galaxy and the LMC, although Fig. 4a of Hughes et al. (2010) shows that there is a very high dispersion in the relation in the LMC. Second, the smaller (and lower metallicity) galaxies have clouds with narrower lines at a given size. To our knowledge, this is the first time this has been noticed. Given the link between metallicity and galaxy size or mass, it is not clear whether the narrower lines are due to the change in metallicity or the change in stellar surface density.

Nonetheless, this shows that molecular clouds have distinctly different properties in different types of galaxies. A number of studies have shown that, even after correcting for a varying N(H₂)∕I_CO(1−0) conversion factor, the molecular gas consumption time is lower in low-metallicity galaxies (Gardan et al. 2007; Gratier et al. 2010b; Braine et al. 2010b; Dib et al. 2011; Druard et al. 2014). This is likely partially due to the weaker stellar winds, slowing cloud dispersal, in subsolar metallicity stars, but also because molecular gas is likely to form at slightly higher densities due to the reduced dust content. Both factors could result in lower cloud line widths for a given size. Clearly, a change in the stellar Initial Mass Function could greatly affect the elements of the above calculation.

Another way of looking at the size-linewidth relation and its variation is shown in Fig. 4. Here we plot the linewidth normalized by the square root of the cloud radius, in order to eliminate the standard size-linewidth relation. The Solomon et al. (1987) relation would be a horizontal line at 1.7 km s⁻¹ pc^−1∕2 in this plot. The M 33 clouds are seen as small black triangles and they are binned by radial intervals to show how this quantity, which can be thought of as the turbulent line width on a 1 pc scale, varies with galactocentric distance. A clear decrease can be seen with distance from the center.

The radial distances in M 33 cover a range of about 7 kpc, with few clouds beyond 6.5 kpc. All of the M 33 averages fall below the Solomon et al. (1987) relation. For comparison, we took the large sample of Heyer et al. (2001) for the outer disk of the Milky Way. Taking the clouds larger than 1 pc, we calculated the linewidth normalized by the square root of the cloud radius and binned by as a function of the distance from the Galactic center. These points cover nearly 12 kpc from the solar circle to 20 kpc and the distance is given on the upper scale. The decrease in the turbulent line width on a 1 pc scale is very similar to that in M 33, confirming that indeed cloud linewidths decrease with increasing galactocentric distance.

Fig. 3 Size-linewidth relation for M 33 clouds and other galaxies observed with similar or better angular resolution. In the top panel, the clouds with extrapolated radii are shown as magenta symbols and the remaining 80% of the M 33 clouds are shown in red. The non-M 33 clouds are shown as black symbols. The region populated by LMC clouds is hatched and the M 51 data (squares) are binned averages. The lower panel has logarithmic axes and has all M 33 clouds in red. The line labeled SRBY is the Solomon et al. (1987) relation for the Galaxy.

Fig. 4 Turbulent line width on a 1 pc scale as a function of galactocentric distance in M 33 (black) and the outer Milky Way (red). Individual M 33 clouds are shown as small black triangles and the large symbols are the binned averages for the M 33 clouds (large black triangles) and for the Heyer et al. (2001) outer Galaxy clouds (red pentagons and red upper distance scale). The large and small error bars on both symbols indicate respectively the dispersion and the uncertainty in the mean on the average values. Both galaxies show a modest decrease with galactocentric radius.

2.2 Linewidth variation with galactocentric distance

It is known that the HI line widths in galaxies decrease with galactocentric distance (Tamburro et al. 2009) but few measures are available formolecular lines of individual clouds. Braine et al. (2010a) reported a decrease in their sample of a few clouds in M 33 and Gratier et al. (2012) reported a weak decline in linewidth. Figure 5 shows the results for the whole sample of 566 clouds. The decline in linewidth is only slight but is significant at the 8σ level, similar to what was found by Gratier et al. (2012). The uncertainty is obtained by bootstrapping, using 5000 iterations each choosing 566 clouds randomly out of the sample (allowing clouds to be chosen more than once or not at all) and examining the dispersion in the correlation coefficient of the fit between galactocentric distance and linewidth (see details in Gratier et al. 2012). It can be seen that the dispersion in line widths is generally lower for the Gaussian fits so these were used in Fig. 3.

2.3 Cloud brightness as a function of galactocentric distance

It is well-known that the large-scale CO brightness of galaxies decreases with distance from the center in most galaxies (see Young et al. 1995; Druard et al. 2014, for survey and M 33 results respectively). Notable counter-examples are M 31 and M 81. Gratier et al. (2012) showed that for the then available cloud sample in M 33 the clouds became considerably less CO-bright with increasing distance to the center. Here we use the entire cloud sample and two different measures of brightness: Fig. 6 shows the peak line temperature reached within the cloud and the peak of the Gaussian fit to the cloud-averaged line profile. The decrease in both quantities is approximately a factor two over 4 kpc, corresponding to a scale length of about 6 kpc. The average CO surface brightness per kpc² decreases much more quickly, with a scale length of approximately 2 kpc (Druard et al. 2014).

Fig. 6 Galactocentric radius (x-axis) vs. peak line temperature reached within the cloud (black, left y-scale) and the peak of the Gaussian fit to the line profile (red symbols, right scale). The factor two difference between the left and right scales is because the peak line temperature is from a single position and is about twice as strong as the Gaussian fit to the cloud-averaged line profile. Both temperatures decrease in a similar way.

2.4 Massspectrum of M 33 clouds

The sample of molecular clouds analyzed here, whose individual cloud properties have been given in the online table of Corbelli et al. (2017) is one of the largest of any galaxy and the largest for which a classification in terms of star formation has been established. We fit a truncated power law to the distribution of cloud masses above the completeness limit, following the procedure described by Maschberger & Kroupa (2009) and the methods described in Gratier et al. (2012), to determine mass spectra for the sample and subcategories, making the assumption as elsewhere that the CO luminosity reflects the mass of clouds (cf. Corbelli et al. 2017).

For the entire sample of 566 clouds, we obtain an exposant α = 1.65 defining the slope of the mass spectrum n(m)dm ∝ m^−αdm (Fig 7). The upper panel shows, as seen previously by Gratier et al. (2012) and Rosolowsky (2005) that the mass spectrum steepens in the outer disk. With the larger sample, we are able to divide the sample into 3 radialbins with close to 200 clouds per bin, far more than earlier samples. The slope of the mass function steepens from α = 1.4 to α = 1.9 with radius and would probably continue to steepen but the number of clouds available (above the completeness limit) decreases sharply, such that the slope becomes poorly defined. Thecause of such a steepening remains unclear so we also compared the mass spectra of the clouds at different stages in the star formation process.

The lower panel of Fig 7 compares the sub-sample of non-star-forming clouds (A) with clouds showing embedded and exposed star formation (respectively B and C clouds, Corbelli et al. 2017). The spectra are remarkably different: the C clouds are more massive and have a distinctly flatter spectrum than, particularly, the clouds without star formation. The B clouds, less numerous, are intermediate. This behavior mimics the radial variation and indeed there is some degeneracy as the star-formingclouds are on average closer to the center than those without star formation. While the cloud classification system is not the same, Kawamura et al. (2009) also found that the more evolved clouds were more massive. We (Corbelli et al. 2017) attribute this to continued gas accretion. Since we measure CO emission, and not mass directly, a systematic variation in the N(H₂)∕I_CO conversion factor could generate a similar result. There are two reasons for thinking that this is unlikely. A sophisticated Bayesian analysis of the N(H₂)∕I_CO factor (and “dark gas”) by Gratier et al. (2017) found no identifiable variation in the N(H₂)∕I_CO conversion as a function of radius in M 33. Secondly, the dust-based cloud masses from the Herschel data show the same trends with star formation class.

The sample is large enough to be divided further in order to determine whether the steepening is primarily due to lack of star formation or position in disk. We thus selected A and C clouds beyond 2 kpc from the center of M 33, obtaining 98 A clouds above the completeness limit with an average galactocentric distance of 3.5 kpc and respectively 174 clouds and a distance of 3.9 kpc for the C clouds. These are plotted as solid lines in Fig. 8. Despite the higher average galactocentric distance, the outer disk C clouds have significantly higher masses and a shallower mass spectrum than the A clouds. In the figure, the number of clouds given is the total sample population, not the number above the completeness limit (given above).

Although A clouds are not common in the inner disk, we selected A and C clouds within respectively 2.5 and 2.8 kpc from the center, yielding 47 and 135 clouds above the completeness limit with an average galactocentric distance of 1.6 kpc for both samples. Again, the mass spectra are very different, with the non-star-forming clouds being very similar to the outer disk population even when they belong to the inner disk. Clouds, and their mass spectra, change as star formation develops and progresses. This change is more important than their position in the disk although the change in CO luminosity is greater as a function of galactocentric distance than star-formation class.

Fig. 7 Mass spectra of M 33 molecular clouds. The y-axis gives the cumulative number of clouds with mass above the corresponding x-axis mass. The vertical line shows the completeness limit as in Corbelli et al. (2017). The black line shows the whole sample in both panels. The color-coded α values give the slopes of the mass spectra. In the top panel, clouds are segregated by galactocentric radius as indicated – the sample has been divided into three roughly equally populated radial bins. In the lower panel, the types indicate clouds with exposed star formation (e.g., Hα emission), embedded star formation, and no star formation, denoted respectively C, B, and A types. The division into types is discussed in Corbelli et al. (2017) and Gratier et al. (2012). The solid lines are the results of the fits, with the color corresponding to the (sub)sample.

Fig. 8 Similar to Fig. 7 but separating inner and outer disk A and C clouds. It is immediately apparent that galactocentric distance is less important than star formation. The number of clouds and the slope of the mass function are given for each group. See text for further details.

2.5 Dust temperatures with and without star formation

Herschel SPIRE data are available for M 33 and, following Braine et al. (2010b), we use the 250 μm to 350 μm flux ratio (after convolving the 250 μm maps to the 350 μm resolution) to estimate dust temperatures for the clouds. We use a dust emissivity β = 1.8 for M 33, following Tabatabaei et al. (2014) for inner disk clouds, in order to calculate the temperatures. While the absolute temperatures determined will vary, the fact that star-forming clouds have demonstrably higher dust temperatures (see Fig. 9) than non-star-forming clouds does not change with the value of β or wavelengths used.

The non-star-forming (“A”) clouds have significantly lower dust temperatures and FIR fluxes than the star-forming clouds, although the fluxes have a very broad distribution for all cloud types. The “B” clouds, with embedded star formation but little or no Hα emission, are presumably on average younger than those (“C”) with exposed star formation. The dust temperatures in the Herschel bands are not distinguishable between B and C clouds but the 250 μm fluxes are slightly lower for the B clouds.

Fig. 9 Representative dust temperatures for the big grain component in M 33. The decomposition into subsamples in terms of their star-forming properties shows that non-star-forming clouds have lower dust temperatures, although the bands used to estimate dust temperatures are completely independent of those used to classify the star formation in the clouds.

3 Rotation of molecular clouds

As in all previous work on the subject (see references in introduction), we assume that velocity gradients reflect cloud rotation. Rotation clearly results in velocity gradients. However, other processes such as turbulence may be able to create velocity gradients (see e.g. Burkert & Bodenheimer 2000, for small scales). The velocity gradient is our only measurable and the systematic aspect of the gradients (Sect. 3.2) argues against other processes. As in earlier work (e.g. Rosolowsky 2005), the first moment of each spatial pixel in the cloud was calculated. In our implementation, we used 5 velocity channels (13 km s⁻¹), centered on the central velocity of the cloud, to measure the first moment $$v_{(x,y)}={\displaystyle \sum _{i=\text{cen}-2}^{i=\text{cen}+2}{v}_i}{T}_i\text{d}v/{\displaystyle \sum _{i=\text{cen}-2}^{i=\text{cen}+2}{T}_i}\text{d}v.$$(1)

This was found to avoid bringing in too much noise (i.e., as when more channels are used) while still covering the velocities occupied by the cloud.

This yields a velocity for each position in the cloud. A plane is then fit to these velocities $$v_{(x,y)}=ax+by+c,$$(2)

where $a=\frac{\partial v}{\partial x}=\frac{\partial v}{\partial \text{RA}}$ and $b=\frac{\partial v}{\partial y}=\frac{\partial v}{\partial \text{Dec}}$ because x and y are the pixel numbers following the RA and Dec directions.

The process is illustrated in Fig. 1, which shows the CO(2–1) map of M 33 with the cloud contours superposed and a zoom on a cloud showing the velocities as measured using Eq. (1) in the top right panel, the fit to the velocities using Eq. (2), and to the bottom right the residual (v − v_fit). This cloud (#4 in the catalogue) was chosen to illustrate the process for its rather clear gradient despite its elongated form. It is worth noting that there is emission which CPROPS was not able to decompose into clouds and that the fraction of the CO emission not decomposed into clouds increases with radius.

M 33 is an inclined spiral whose north(-east)ern side is approaching, in other words, this side has a higher negative velocity. The near side of the disk is the western side. The geometry is most easily visualized if one initially thinks of M 33 as oriented N–S. Thus, prograde follows the rotation of the disk, in which velocities become more negative to the north. Rotating M 33 counter-clockwise by 22.5° to its true orientation on the sky, things change little in that northern parts of a cloud have increasingly negative velocities when cloud rotation is prograde. The gradients are measured as $\left|\right|\nabla V\left|\right|=\sqrt{{(\frac{\partial v}{\partial \text{RA}})}^2+{(\frac{\partial v}{\partial \text{Dec}})}^2}$ where $\frac{\partial v}{\partial \text{RA}}$ and $\frac{\partial v}{\partial \text{Dec}}$ come from the fits of a plane to the cloud velocities. The sign of the gradient (negative for prograde, positive for retrograde in order to fit the true orientation of M 33) follows the sign of $\frac{\partial v}{\partial \text{RA}}\text{sin}(22.5)+\frac{\partial v}{\partial \text{Dec}}\text{cos}(22.5)$.

Figure 10 shows the cloud velocity gradients – for both the entire sample and for only the high signal-to-noise (S/N) clouds. The fact that the distribution of velocity gradients extends to higher (absolute) values for the lower luminosity clouds is an immediate suggestion that noise is contributing to the observed gradients. Moving to the more luminous clouds, not only is the distribution of velocity gradients narrower but it is more skewed towards negative (prograde) values. It may be worth noting that no difference in linewidth or in velocity gradient between star-forming (C clouds, cf. Sects. 2.4 and 2.5) and non-star-forming clouds is found in our sample.

Fig. 1 Left: CO(2–1) integrated intensity map of M 33 with the cloud contours superposed in red and a zoom on cloud #4 showing (top right) the velocities as measured using Eq. (1), middle right: the fit to the velocities using Eq. (2), and to the bottom right the velocity residual, all in km s −1.

Fig. 10 Histogram of cloud velocity gradients – for the entire sample in black and for only the stronger clouds in red, where the line temperature averaged over the whole cloud is over 100 mK T${}_a^{*}$. Prograde rotation is given a negative sign here because the galaxy rotation velocity increases with decreasing declination.

3.1 Beam smearing

Molecular clouds are comparable in size to the angular resolution of our observations. A typical (deconvolved) cloud radius is about 45 pc (Fig. 3), to be compared with the 24 pc beam half-power radius. The goal of this section is to create mock cloudsand compare the gradients derived before and after convolving with the telescope beam. Intuitively, one expects gradients to weaken as the angular resolution is degraded.

We create mock clouds with exactly the same gradients and sizes as those we have measured. The gradients are injected as being linear and along the direction of the observed gradient. The peak line temperature decreases linearly with distance to the cloud center until it reaches the noise level. The lines are assumed Gaussian. This process is described in more detail in Sect. 4. No noise is added in this step as we wish to measure the effect of resolution and not noise (done later).

After creating 566 mock clouds, we convolve them with a Gaussian beam of half-power width 12″. We then follow exactly the same procedure to measure the gradients. Figure 11 shows the comparison of the injected and recovered gradients (after convolution). As could be expected, the post-convolution gradients are smaller, with an average ratio of 0.59. We use this to correct the observed gradients in order to estimate rotation periods, angular momenta, and rotational energies. None of the gradients after convolution are higher.

Fig. 11 Comparison of the gradients measured pre- and post-convolution with the telescope beam. On average, a correction factor of 1/0.59 should be applied to the observed gradients.

3.2 Comparison with simulations

Simulations appear to favor prograde rotation but the comparison with our measures is not always straightforward. Dobbs (2008) and Li et al. (2018) use the specific angular momentum, defined as the angular momentum per unit mass averaged over the cloud, instead of the gradient. The units are km s⁻¹ pc.

The angular momentum of a rotating disk with a surface density declining as $\Sigma (r)=\Sigma (r_0){(r/r_0)}^{-1}$ is $L={\displaystyle {\int }_0^R{r}^2}\Omega \Sigma (r)2\pi \text{d}r$. Dividing by the mass of the disk (2πΣ₀r₀R) yields a specific angular momentum of j = ΩR²∕3 where Ω is the angular velocity and R the outer radius of the disk. A finite sphere with a density decreasing as r⁻² yields a similar result. If the disk has a constant surface density, then the specific angular momentum is j = ΩR²∕2 (Blitz 1993) but this is less likely. We thus consider the specific angular momentum of our clouds to be j = (Ω∕0.59) R²∕3 where the 0.59 corrects for the underestimate in the velocity gradient due to beam smearing.

Figure 12 shows the distribution of the angular momenta of our clouds. The top panel is a histogram similar to the velocity gradients. The bottom panel is intended to be useful for comparison with Fig. 8 of Dobbs (2008), giving the angular momenta as a function of cloud mass. In both panels, the strong clouds are in red and retro/prograde rotation are separated either by sign (top) or by symbol (bottom). Prograde rotation clearly dominates, as in the simulations where the self-gravity plays a role in the Dobbs (2008) simulations. Comparing with Li et al. (2018) Run VI, the distribution of the angular momenta in the M 33 clouds appears narrower. The Li et al. (2018) simulations have lower cloud masses and angular momenta for the high resolution Run VI, although the resolution of our observations is worse (suggesting that the difference would be greater ifthe spatial resolutions were closer). Without more measurements of angular momenta of clouds, including in higher surface density galaxies than M 33, it is difficult to be conclusive about the comparison with simulations.

Fig. 12 Specific angular momenta of the M 33 clouds. Red lines or symbols indicate data for the CO-strong clouds. The top panel shows clearly that the asymmetry favoring prograde gradients (Fig. 10) is also present for the angular momenta. The lower panel shows that the variation of angular momentum with cloud mass is very weak and was designed for comparison with simulations.

4 Tests of velocity gradients

We first need to convince ourselves that we measure real gradients. The fact that the pro/retrograde differences are more pronounced for the high S/N sample is certainly a sign that this is the case. In the following subsections, we examine the effect of beam smearing due to the resolution of our observations and create mock clouds with properties very similar to the real clouds in order to test our ability to retrieve cloud rotation in the presence of noise.

The mock clouds are created using the masks of the real clouds, such that the sizes and shapes are perfectly represented. The noise level of the cube is about 20 mK per channel (Druard et al. 2014) and we generate mock clouds with peak CO line temperatures T_max= 100, 200, 400, 800 mK in order to obtain varying S/N levels. The pixels are then given a temperature $$T_{xyj}=\left[20+(T_{\text{max}}-20)\times \left(1-\frac{R}{R_{\text{max}}}\right)\right]\text{exp}\left(-\frac{{(v_j-v_{xy})}^2}{2\Delta V^2}\right),$$(3)

where x, y, and j represent respectively the pixel numbers along the RA, Dec, and velocity axes, and R the distance from the center of the cloud (the position of the center of the cloud is returned by CPROPS). The gradient is injected through the function relating v_xy and the position through Eq. (2). Thus, the line is centered on the velocity v_xy and follows a Gaussian with dispersion ΔV . As can be seen in Fig. 5, the linewidth at half power is roughly 7 km s⁻¹. We therefore inject a velocity dispersion ΔV = 3 km s⁻¹, corresponding to a half-power linewidth of 7.05 km s⁻¹. The central temperature decreases linearly with from T_max to the noise level. The sampling in space and velocity is the same as for the real cube (3″, 2.6 km s⁻¹).

The next step is to add noise. We add random Gaussian noise using the well-tested noise random number generator within GILDAS. However, the data cube has undergone many transformations and the true noise may not be precisely Gaussian. Thus, we extract contiguous channels from signal-free regions of the cube which we use as noise. In fact, as any transformations were designed to preserve particularly the region where signal is present, this is a sort of worst-case noise.

We first examine what we obtain from mock clouds with no velocity gradient, that is to say with v_xy constant. Figure 13 shows the results for T_max = 200 mK using both purely random noise and noise taken from signal-free but unoptimized regions of the cube. The gradients, for equivalent clouds and noise levels, are clearly more dispersed with the “real” cube noise. The dispersions, as measured by the full width at half power of the distribution (divided by 2.35 to give the equivalent for a Gaussian), are approximately 0.008 and 0.015 km s⁻¹ pc⁻¹for the random and cube noise injection. T_max = 200 mK corresponds roughly to the median signal in the cloud sample.

We have also tested with higher and lower S/N levels. For T_max= 100 mK, the dispersion increases to 0.014 and 0.029 km s⁻¹ pc⁻¹ for respectively random and cube noise. For T_max = 400 mK, the dispersion decreases to 0.005 and 0.0085 km s⁻¹ pc⁻¹ for respectively random and cube noise.

Fig. 5 Variation of the cloud linewidth vs. galactocentric distance. Black symbols are the line widths from CPROPS and red symbols are from fitting Gaussians to the cloud spectra. Errorbars indicate the dispersion within each bin, showing that generally the Gaussian fits yield fewer extreme results (as found by Gratier et al. 2012).

Fig. 13 Histogram of gradients from noise. The black line shows the distribution of the gradients for purely random noise. The red line shows the distribution when noise is from signal-free regions of the cube. No velocity gradient has been injected here.

5 Evaluating uncertainties on cloud rotation

The same operations were done with the observed gradients, creating the same clouds but injecting the velocity gradient deduced from the observations for each cloud. Should there be a link between size or shape and the velocity gradient deduced from calculating the first moment, the link would be preserved in these tests. Figure 14 shows the distribution of the retrieved gradients. The clouds were created with the observed gradients and shapes and then noise was added and the gradients remeasured. This can be directly compared with Fig. 10. The distribution in Fig. 14 is of course wider because the gradients from Fig. 10 have been injected and then noise added. The process was repeated with T_max= 100 mK and the distribution is significantly wider but the prograde-retrograde asymmetry is nonetheless preserved. We are thus confident that the gradients retrieved reflect the true distribution of velocity gradients, although the distribution is likely broadened by the presence of the noise in the datacube.

Fig. 14 Like Fig. 13 except that the observed gradients (see Fig. 10) have been injected, noise added, and the gradients recovered. Black line shows the distribution of the gradients for purely random noise. The red line shows the distribution when noise from signal-free regions of the cube. As in Fig. 13, Tmax = 200 mK.

5.1 Galactic gradient

M 33 itself has a velocity gradient due to rotation. Since the position angle of M 33 is close to vertical, there is little gradient expected along the RA (x) axis but there is a negative gradient along the Dec (y) axis becauseat higher Declinations the velocity is more negative. We calculate this for axisymmetric rotation assuming the rotation curve given in Eq. (18) of López Fune et al. (2017): $$V(r)=V_0\frac{(r/r_0)+d}{(r/r_0)+1},$$

where V ₀ = 139.2 km s⁻¹, r₀ = 1.3 kpc, and d = 0.12.

Figure 15 shows the local velocity gradients $a=\frac{\partial v}{\partial \text{RA}}$ and $b=\frac{\partial v}{\partial \text{Dec}}$ derived from the axisymmetric rotation curve. The plot of $\frac{\partial v}{\partial \text{Dec}}$ is negative everywhere with fairly high (absolute) values but $\frac{\partial v}{\partial \text{RA}}$ has both negative and positive regions with a positive average. These values come exclusively from the rotation curve and thus include differential rotation and thus shear.

Figure 15 is not very intuitive. In order to qualitatively understand the negative and positive zones, let us think of isovelocity curves of a differentially rotating spiral disk with a monotonically increasing rotation curve (the so-called spider diagram). When the major axis is vertical (N–S), then the only horizontal iso-velocity curve (i.e., $\frac{\partial v}{\partial \text{RA}}=0$) is along the minor axis. There are no vertical (i.e., $\frac{\partial v}{\partial \text{Dec}}=0$) iso-velocity curves, such that for velocities decreasing towards the north, $\frac{\partial v}{\partial \text{Dec}}<0$ everywhere. Now let us rotate the diagram counterclockwise slightly. In the northern half, we will have a locus of $\frac{\partial v}{\partial \text{Dec}}=0$ points just to the left of the major axis, where the iso-velocity curves are briefly horizontal. Slightly above the minor axis and to the left, the isovelocity curve which went slowly upwards pre-rotation now is approximately flat, leading to another series of $\frac{\partial v}{\partial \text{Dec}}=0$ points. The same is true by symmetry to the south. No such region where $\frac{\partial v}{\partial \text{Dec}}=0$ is present to the upper right or lower left. The magnitude can be understood by imagining how closely spaced (along RA or along Dec) the isovelocity curves are (for equal velocity spacing). This is why the highly negative regions of $\frac{\partial v}{\partial \text{Dec}}$ are close to the minor axis. Similarly, $\frac{\partial v}{\partial \text{RA}}$ is high where isovelocity curves are closely spaced and close to vertical. In all cases, the velocity gradients due to galactic rotation are larger near the center where the rotation rises sharply.

5.2 Comparison of cloud and galactic gradients

Having fit a plane to the velocities of the pixels making up each cloud, we have the gradients along the RA and Dec axes and we can look for patterns. Given Fig. 10 which shows that we can have more confidence in the high-luminosity clouds, we plot the gradients for the stronger clouds in Fig. 16 in a way that can be compared directly with Fig. 15. Let us consider the “null” hypothesis to be that clouds on average are not rotating with respect to their surroundings, that is they rotate with the galaxy. Our results are close to this null hypothesis (cf. Figs. 15 and 16). If correct, an implication is that the cloud formation mechanism has little influence on the velocity gradient.

From Fig. 10, a typical prograde rotation velocity is ≲0.03 km s⁻¹ pc⁻¹. Fig. 16 shows that this is a good representative value for the CO-strong clouds. Including a factor 1∕0.59 to compensate for the beam smearing discussed earlier, this yields Ω≲0.05 km s⁻¹ pc⁻¹. The rotationperiod is thus about $T=\frac{2\pi }{\Omega }\approx 120$ Myr, which is similar to the rotation period of the inner disk of M 33.

Thus, not only are real gradients measured in these clouds but we are able to show that they are dominated by prograde rotation despite the extremely low values. The rotation periods are longer than cloud lifetimes and comparable to the Galactic rotation period. The link between Figs. 15 and 16 is real: the average observed $\frac{\partial v}{\partial \text{RA}}=2.7\pm 1.2$ m s⁻¹ pc⁻¹ where $\frac{\partial v}{\partial \text{RA}}>0$ in Fig. 15 but the average observed $\frac{\partial v}{\partial \text{RA}}=-2.0\pm 1.5$ m s⁻¹ pc⁻¹ where $\frac{\partial v}{\partial \text{RA}}<0$ in Fig. 15 and all averages are negative for $\frac{\partial v}{\partial \text{Dec}}$.

Fig. 15 Coefficients $\frac{\partial v}{\partial \text{RA}}$ (left) and $\frac{\partial v}{\partial \text{Dec}}$ (right) calculated for adjacent pixels using the rotation curve above. Units are meter s−1 pc−1 and the color wedges are shown at the top of each panel. In the left panel, contours are drawn at 0, 10, and 20 m s−1 pc−1 and in the right panel at −10, −20, and −30 m s−1 pc−1.

Fig. 16 Coefficients $a=\frac{\partial v}{\partial \text{RA}}$ and $b=\frac{\partial v}{\partial \text{Dec}}$ for the 222 stronger clouds (TCO > 0.11K), shown only when above 8 m s−1 pc−1 in absolute value in order to reduce the influence of noise and make the numbers legible. Negative values are shown in red and positive in black. The numbers of negative and positive values are given in the panels. Because the panels only show the higher (absolute) values, the disproportion in the right panel is actually greater: 31 values > 8 and 116 values < −8.

5.3 Magnitude of velocity gradients

Is rotation a significant hindrance to cloud collapse? Adopting Ω ≈ 0.05 km s⁻¹ pc⁻¹as typical of a “rotating” cloud, we can compare the rotational kinetic energy with the gravitational potential energy or the edge velocities with escape velocities. Adopting M = 2 × 10⁵ M_⊙ and R = 30 pc as representative values, the rotational kinetic energy is E_rot ≈ 10⁴⁸ ergs whereas the gravitational potential energy is nearly E_grav ≈ 10⁵⁰ ergs. Similarly, the rotation velocity at the cloud edge could be expected to be v ≈ 1.5 km s⁻¹ but the escape velocity is much higher, v_esc ≈ 7.5 km s⁻¹. This large difference shows also that rotation contributes little to the overall support and line width of the cloud.

6 Conclusions

In this work, we have shown, for the first time to our knowledge, that molecular clouds rotate and that their rotation is very slow but measurable from our high-quality data. This relies on the assumption, as in previous work, that rotation can be deduced from velocity gradients. The rotation tends to be prograde. The majority of molecular clouds have an angular velocity below 0.03 km s⁻¹ pc⁻¹ (0.05 km s⁻¹ pc⁻¹ after correcting for beam-smearing), yielding a rotation period greatly superior to the cloud lifetimes of about 15 Myr in M 33 (Corbelli et al. 2017). The rotation contributes (very) little to the support of the cloud against gravity. Simulations as well as classical calculations (e.g. Rosolowsky et al. 2003) tend to find higher angular velocities. At (much) smaller scales rotation is clearly present: stars and proto-stars have disks (which rotate) and rotation was also observed in the massive proto-stellar core W43-mm1 (Jacq et al. 2016).

Not only do molecular cloud mass spectra steepen with galactocentric distance, but the mass spectrum appears to depend even more strongly on whether the clouds host active star formation. At equivalent galactocentric distance, molecular clouds which form stars have considerably flatter mass spectra than those without star formation.

Comparing the molecular clouds in M 33 with those in other nearby galaxies, a displacement in the size-linewidth relation appears in that lower metallicity systems have narrower CO lines for comparable cloud size. There is also a trend for cloud linewidths to become narrower with increasing galactocentric distance. Some degeneracy is present in these measurements as both metallicity and stellar surface density decrease with galactocentric distance and subsolar metallicity galaxies tend to have lower stellar surface densities.

Acknowledgements

The authors would like to gratefully acknowledge several students from the University of Bordeaux who did projects on various aspects of the data: Jimmy Mata, Marc-Robert Antoine, Thomas Goncalves, and Cyril Lenain.

References

All Figures

Fig. 2 Link between CPROPS undeconvolved cloud radii (x-axis) and deconvolved radii for the 80% of the clouds with well-defined deconvolved radii (black symbols). The red triangles show the radii attributed to the clouds for which CPROPS could not estimate radii – the radii follow the envelope defined by the CPROPS values down to a constant value below which we have no confidence. The goal here is to avoid generating clouds which are inappropriately small while still attributing reasonable radii.

In the text

Fig. 3 Size-linewidth relation for M 33 clouds and other galaxies observed with similar or better angular resolution. In the top panel, the clouds with extrapolated radii are shown as magenta symbols and the remaining 80% of the M 33 clouds are shown in red. The non-M 33 clouds are shown as black symbols. The region populated by LMC clouds is hatched and the M 51 data (squares) are binned averages. The lower panel has logarithmic axes and has all M 33 clouds in red. The line labeled SRBY is the Solomon et al. (1987) relation for the Galaxy.

In the text

Fig. 4 Turbulent line width on a 1 pc scale as a function of galactocentric distance in M 33 (black) and the outer Milky Way (red). Individual M 33 clouds are shown as small black triangles and the large symbols are the binned averages for the M 33 clouds (large black triangles) and for the Heyer et al. (2001) outer Galaxy clouds (red pentagons and red upper distance scale). The large and small error bars on both symbols indicate respectively the dispersion and the uncertainty in the mean on the average values. Both galaxies show a modest decrease with galactocentric radius.

In the text

Fig. 6 Galactocentric radius (x-axis) vs. peak line temperature reached within the cloud (black, left y-scale) and the peak of the Gaussian fit to the line profile (red symbols, right scale). The factor two difference between the left and right scales is because the peak line temperature is from a single position and is about twice as strong as the Gaussian fit to the cloud-averaged line profile. Both temperatures decrease in a similar way.

In the text

Fig. 7 Mass spectra of M 33 molecular clouds. The y-axis gives the cumulative number of clouds with mass above the corresponding x-axis mass. The vertical line shows the completeness limit as in Corbelli et al. (2017). The black line shows the whole sample in both panels. The color-coded α values give the slopes of the mass spectra. In the top panel, clouds are segregated by galactocentric radius as indicated – the sample has been divided into three roughly equally populated radial bins. In the lower panel, the types indicate clouds with exposed star formation (e.g., Hα emission), embedded star formation, and no star formation, denoted respectively C, B, and A types. The division into types is discussed in Corbelli et al. (2017) and Gratier et al. (2012). The solid lines are the results of the fits, with the color corresponding to the (sub)sample.

In the text

	Fig. 8 Similar to Fig. 7 but separating inner and outer disk A and C clouds. It is immediately apparent that galactocentric distance is less important than star formation. The number of clouds and the slope of the mass function are given for each group. See text for further details.
In the text

	Fig. 9 Representative dust temperatures for the big grain component in M 33. The decomposition into subsamples in terms of their star-forming properties shows that non-star-forming clouds have lower dust temperatures, although the bands used to estimate dust temperatures are completely independent of those used to classify the star formation in the clouds.
In the text

	Fig. 1 Left: CO(2–1) integrated intensity map of M 33 with the cloud contours superposed in red and a zoom on cloud #4 showing (top right) the velocities as measured using Eq. (1), middle right: the fit to the velocities using Eq. (2), and to the bottom right the velocity residual, all in km s −1.
In the text

	Fig. 10 Histogram of cloud velocity gradients – for the entire sample in black and for only the stronger clouds in red, where the line temperature averaged over the whole cloud is over 100 mK T${}_a^{*}$. Prograde rotation is given a negative sign here because the galaxy rotation velocity increases with decreasing declination.
In the text

	Fig. 11 Comparison of the gradients measured pre- and post-convolution with the telescope beam. On average, a correction factor of 1/0.59 should be applied to the observed gradients.
In the text

	Fig. 12 Specific angular momenta of the M 33 clouds. Red lines or symbols indicate data for the CO-strong clouds. The top panel shows clearly that the asymmetry favoring prograde gradients (Fig. 10) is also present for the angular momenta. The lower panel shows that the variation of angular momentum with cloud mass is very weak and was designed for comparison with simulations.
In the text

	Fig. 5 Variation of the cloud linewidth vs. galactocentric distance. Black symbols are the line widths from CPROPS and red symbols are from fitting Gaussians to the cloud spectra. Errorbars indicate the dispersion within each bin, showing that generally the Gaussian fits yield fewer extreme results (as found by Gratier et al. 2012).
In the text

	Fig. 13 Histogram of gradients from noise. The black line shows the distribution of the gradients for purely random noise. The red line shows the distribution when noise is from signal-free regions of the cube. No velocity gradient has been injected here.
In the text

	Fig. 14 Like Fig. 13 except that the observed gradients (see Fig. 10) have been injected, noise added, and the gradients recovered. Black line shows the distribution of the gradients for purely random noise. The red line shows the distribution when noise from signal-free regions of the cube. As in Fig. 13, Tmax = 200 mK.
In the text

	Fig. 15 Coefficients $\frac{\partial v}{\partial \text{RA}}$ (left) and $\frac{\partial v}{\partial \text{Dec}}$ (right) calculated for adjacent pixels using the rotation curve above. Units are meter s−1 pc−1 and the color wedges are shown at the top of each panel. In the left panel, contours are drawn at 0, 10, and 20 m s−1 pc−1 and in the right panel at −10, −20, and −30 m s−1 pc−1.
In the text

Fig. 16 Coefficients $a=\frac{\partial v}{\partial \text{RA}}$ and $b=\frac{\partial v}{\partial \text{Dec}}$ for the 222 stronger clouds (TCO > 0.11K), shown only when above 8 m s−1 pc−1 in absolute value in order to reduce the influence of noise and make the numbers legible. Negative values are shown in red and positive in black. The numbers of negative and positive values are given in the panels. Because the panels only show the higher (absolute) values, the disproportion in the right panel is actually greater: 31 values > 8 and 116 values < −8.

In the text