Properties and rotation of molecular clouds in M 33

The sample of 566 molecular clouds identified in the CO(2--1) IRAM survey covering the disk of M~33 is explored in detail.The clouds were found using CPROPS and were subsequently catalogued in terms of their star-forming properties as non-star-forming (A), with embedded star formation (B), or with exposed star formation C.We find that the size-linewidth relation among the M~33 clouds is quite weak but, when comparing with clouds in other nearby galaxies, the linewidth scales with average metallicity.The linewidth and particularly the line brightness decrease with galactocentric distance.The large number of clouds makes it possible to calculate well-sampled cloud mass spectra and mass spectra of subsamples.As noted earlier, but considerably better defined here, the mass spectrum steepens (i.e. higher fraction of small clouds) with galactocentric distance.A new finding is that the mass spectrum of A clouds is much steeper than that of the star-forming clouds.Further dividing the sample, this difference is strong at both large and small galactocentric distances and the A vs C difference is a stronger effect than the inner/outer disk difference in mass spectra.Velocity gradients are identified in the clouds using standard techniques.The gradients are weak and are dominated by prograde rotation; the effect is stronger for the high signal-to-noise clouds.A discussion of the uncertainties is presented.The angular momenta are low but compatible with at least some simulations.The cloud and galactic gradients are similar; the cloud rotation periods are much longer than cloud lifetimes and comparable to the galactic rotation period.The rotational kinetic energy is 1-2\% of the gravitational potential energy and the cloud edge velocity is well below the escape velocity, such that cloud-scale rotation probably has little influence on the evolution of molecular clouds.


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. 2010aGratier et al. , 2012Donovan 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 M33 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 M51. 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 −1 pc −1 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 Section 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.

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 2 )/I CO(1−0) factor of 4 × 10 20 cm −2 /(K km s −1 ) and a CO( 2−1 1−0 ) line ratio of 0.8. These values have been validated by Druard et al. (2014), Gratier et al. (2017), andBraine 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 (black) and the values calculated for the others (red symbols). 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). Figure 3 shows the size-linewidth relation for the M33 clouds (red symbols) 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).

Size-linewidth relation
Two things are apparent from Figure 3. First, the sizelinewidth 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 2 )/I CO(1−0) conversion factor, the molecular gas con-sumption 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 −1 pc −1/2 in this plot. The M33 clouds are seen as small black triangles and they are binned by radial intervals (large black triangles) 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 large error bars indicate the dispersion and the small red error bars indicate the uncertainty in the mean value.
The radial distances in M 33 cover a range of about 7kpc, 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 are shown as red pentagons and the distance is given on the red upper scale, covering nearly 12kpc from the solar circle to 20 kpc. As for the black triangles, the large error bars give the dispersion and the small error bars give the uncertainty in the mean value. The decrease in the turbulent line width on a 1 pc scale is very similar to that in M33, confirming that indeed cloud linewidths decrease with increasing galactocentric distance.

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 for molecular lines of individual clouds. Braine et al. (2010a) reported a decrease in their sample of a few clouds in M33 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 Figure 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 M33 results respectively). Notable counter-examples are M 31 and M81. Gratier et al. (2012) showed that for the then available cloud sample in M 33 the clouds became considerably less CObright 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 2 decreases much more quickly, with a scale length of approximately 2 kpc ).

Mass spectrum of M 33 clouds
The sample of molecular clouds analyzed here, whose individual cloud properties have been given in the on-line 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 radial bins 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. The cause 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-forming clouds are on average closer to the center than Link between CPROPS undeconvolved cloud radii (xaxis) and deconvolved radii for the 80% of the clouds with welldefined 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. 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 ) attribute this to continued gas accretion. Since we measure CO emission, and not mass directly, a systematic variation in the N(H 2 )/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 2 )/I CO factor (and "dark gas") by Gratier et al. (2017) found no identifiable variation in the N(H 2 )/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 M33, 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 nonstar-forming clouds being very similar to the outer disk popula- tion 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.

Dust temperatures with and without star formation
Herschel SPIRE data are available for M33 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 M33, 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 starforming clouds have demonstrably higher dust temperatures (see Figure 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. 6. Galactocentric radius (x-axis) versus 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.

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 −1 ), centered on the central velocity of the cloud, to measure the first moment 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 where a = ∂v ∂x = ∂v ∂RA and b = ∂v ∂y = ∂v ∂Dec because x and y are the pixel numbers following the RA and Dec directions.
The process is illustrated in Figure 1, which shows the CO(2-1) map of M33 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 catalog) 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.
M33 is an inclined spiral whose North(-east)ern side is approaching, i.e. a higher negative velocity. The near side of the 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 3 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. disk is the western side. The geometry is most easily visualized if one initially thinks of M33 as oriented N-S. Thus, prograde follows the rotation of the disk, in which velocities become more negative to the North. Rotating M33 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.   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. Fig. 9. Representative dust temperatures for the big grain component in M33. 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. Figure 10 shows the cloud velocity gradients -for the entire sample in black and for only the high signal-to-noise (S/N) clouds in red. 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 (red), not only is the distribution of velocity gradients narrower but it is more skewed towards negative (prograde) values. 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 100mK T * a . Prograde rotation is given a negative sign here because the galaxy rotation velocity increases with decreasing declination.
It may be worth noting that no difference in linewidth or in velocity gradient between star-forming (C clouds, cf. Sections 2.4 and 2.5) and non-star-forming clouds is found in our sample.

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 24pc beam halfpower radius. The goal of this section is to create mock clouds and 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 Section 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. Fig. 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.

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. (2017) 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 −1 pc. Fig. 11. Comparison of the gradients measured pre-and postconvolution with the telescope beam. On average, a correction factor of 1/0.59 should be applied to the observed gradients.
The angular momentum of a rotating disk with a surface density declining as Σ(r) = Σ(r 0 )(r/r 0 ) −1 is L = R 0 r 2 Ω Σ(r) 2 π dr. Dividing by the mass of the disk (2πΣ 0 r 0 R) yields a specific angular momentum of j = ΩR 2 /3 where Ω is the angular velocity and R the outer radius of the disk. A finite sphere with a density decreasing as r −2 yields a similar result. If the disk has a constant surface density, then the specific angular momentum is j = ΩR 2 /2 (Blitz 1993) but this is less likely. We thus consider the specific angular momentum of our clouds to be j = (Ω/0.59) R 2 /3 where the 0.59 corrects for the underestimate in the velocity gradient due to beam smearing. Fig. 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 Figure 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. (2017) Run VI, the distribution of the angular momenta in the M33 clouds appears narrower. The Li et al. (2017) 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 if the spatial resolutions were closer). Without more measurements of angular momenta of clouds, including in higher surface density galaxies than M33, it is difficult to be conclusive about the comparison with simulations.

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 Fig. 12. Specific angular momenta of the M33 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.
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 ) 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 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 −1 . We therefore inject a velocity dispersion 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.
∆V = 3 km s −1 , corresponding to a half-power linewidth of 7.05 km s −1 . 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 −1 ). 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.

No velocity gradient with noise
We first examine what we obtain from mock clouds with no velocity gradient, i.e. with v xy constant. Fig. 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 distrbution (divided by 2.35 to give the equivalent for a gaussian), are approximately 0.008 and 0.015 km s −1 pc −1 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 −1 pc −1 for respectively random and cube noise. For T max = 400 mK, the dispersion decreases to 0.005 and 0.0085 km s −1 pc −1 for respectively random and cube noise.

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. Fig. 14 shows the distribution of the retrieved gradients. The clouds were created with the observed gradients and shapes  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, T max = 200 mK. 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.

Galactic gradient
M 33 itself has a velocity gradient due to rotation. Since the position angle of M33 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 because at 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): where V 0 = 139.2 km s −1 , r 0 = 1.3 kpc, and d = 0.12. Figure 15 shows the local velocity gradients a = ∂v ∂RA and b = ∂v ∂Dec derived from the axisymmetric rotation curve. The plot of ∂v ∂Dec is negative everywhere with fairly high (absolute) values but ∂v ∂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. ∂v ∂RA = 0) is along the minor axis. There are no vertical (i.e. ∂v ∂Dec = 0) iso-velocity curves, such that for velocities decreasing towards the North, ∂v ∂Dec < 0 everywhere. Now let us rotate the diagram counterclockwise slightly.
In the northern half, we will have a locus of ∂v ∂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 ∂v ∂Dec = 0 points. The same is true by symmetry to the South. No such region where ∂v ∂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 ∂v ∂Dec are close to the minor axis. Similarly, ∂v ∂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.

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, i.e. they rotate with the galaxy. Our results are close to this "null" hypothesis (cf. Figures 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 −1 pc −1 . 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 −1 pc −1 . The rotation period is thus about T = 2π Ω ≈ 120 Myr, which is similar to the rotation period of the inner disk of M33.
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 Fig. 15 and Fig. 16 is real: the average observed ∂v ∂RA = 2.7 ± 1.2 m s −1 pc −1 where ∂v ∂RA > 0 in Fig. 15 but the average observed ∂v ∂RA = −2.0 ± 1.5 m s −1 pc −1 where ∂v ∂RA < 0 in Fig. 15 and all averages are negative for ∂v ∂Dec .

Magnitude of velocity gradients
Is rotation a significant hindrance to cloud collapse? Adopting Ω ≈ 0.05 km s −1 pc −1 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 5 M ⊙ and R = 30 pc as representative values, the rotational kinetic energy is E rot ≈ 10 48 ergs whereas the gravitational potential energy is nearly E grav ≈ 10 50 ergs. Similarly, the rotation velocity at the cloud edge could be expected to be v ≈ 1.5 km s −1 but the escape velocity is much higher, v esc ≈ 7.5 km s −1 . This large difference shows also that rotation contributes little to the overall support and line width of the cloud.

Conclusions
In this work, we have shown, for the first time to our knowledge, that molecular clouds rotate and that their rotation is very Fig. 15. Coefficients ∂v ∂RA (left) and ∂v ∂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 . 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 .03 km s −1 pc −1 (.05 km s −1 pc −1 after correcting for beamsmearing), yielding a rotation period greatly superior to the cloud lifetimes of about 15 Myr in M 33 ). 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.