ACA CO($J=2-1$) Mapping of the Nearest Spiral Galaxy M33. I. Initial Results and Identification of Molecular Clouds

The interstellar medium (ISM) is one of the crucial components in galaxies because stars are formed by the contraction of molecular ISM. In the Milky Way (MW), a large fraction of molecular ISM is in the form of giant molecular clouds (GMCs: Sanders et al., 1985), whose typical size and mass are a few $\times$ $10-100$ pc and $10^{4}-10^{6}\,M_{\odot}$, respectively. It is essential to investigate the properties and formation/evolution processes of GMCs because they are known to be major sites of high-mass star formation, which eventually drives the evolution of galaxies.

So far, a lot of studies have investigated various GMC properties and their relationships. In the MW, Larson (1981) found that the internal velocity dispersions of the molecular clouds are well correlated with their sizes and masses, and also reported that these correlations (i.e., scaling relations) can be expressed as the power-law form. Solomon et al. (1987) measured the velocity dispersions, sizes, virial masses, and CO luminosities for 273 GMCs in the Galactic disk, and found that the velocity dispersion is proportional to the 0.5 power of the size. They also found a tight relationship, over four orders of magnitude, between the virial mass and the CO luminosity with a power-law slope of $\sim$0.8.

Such GMC studies were expanded to the Local Group galaxies outside the MW. Fukui et al. (2008) performed a CO survey toward the Large Magellanic Cloud (LMC) at a spatial resolution of $\sim$40 pc, and identified 272 GMCs with a mass range from $2\times 10^{4}\,M_{\odot}$ to $7\times 10^{6}\,M_{\odot}$ (see also Fukui et al., 1999). In addition, Kawamura et al. (2009) examined spatial comparisons of these GMCs with young star clusters (YSCs) and H ii regions and found that the GMCs can be classified into three types: (1) GMCs associated with no H ii regions nor YSCs, (2) GMCs associated only with small H ii regions, but with no YSCs, and (3) GMCs associated with both YSCs and large H ii regions. Such a classification of GMCs according to the activities of high-mass star formation likely reflects their evolutionary sequence. In addition, GMC surveys have been often conducted toward M33, which is one of the nearest spiral galaxies (e.g., Engargiola et al., 2003; Rosolowsky et al., 2007; Gratier et al., 2012; Miura et al., 2012; Corbelli et al., 2017). These studies identified more than 100 GMCs (in particular, more than 500 GMCs by Corbelli et al., 2017) over the M33 disk at $\sim$50 pc resolution, and discussed timescales and the evolutionary stages of GMCs based on the comparison with H ii regions and YSCs as well as the LMC studies.

High-angular resolution observations by millimeter-wave interferometers enabled to perform the unbiased GMC surveys even toward external spiral galaxies. Colombo et al. (2014) reported the GMC catalog, which contains $\sim$1500 individual objects in the grand-design spiral galaxy M51 at $\sim$40 pc resolution using data from the PdBI Arcsecond Whirlpool Survey (Schinnerer et al., 2013). They proposed that large-scale dynamical processes and feedback from high-mass star formation cause environmental variations in the GMC properties and mass distributions, and also suggested that $\sim$30% of GMCs in M51 are unbound. More recently, PHANGS-ALMA survey mapped CO($J=2-1$) line emission at $\sim$1$\arcsec$ resolution toward 90 nearby star-forming galaxies (Leroy et al., 2021). In particular, Rosolowsky et al. (2021) identified 4986 molecular clouds at a common 90 pc resolution and measured their properties for ten subsamples. They found that the physical properties of clouds vary among galaxies, both as a function of galactocentric radius and as a function of the dynamical environment (e.g., bar, spiral arm, and inter-arm).

However, these earlier studies for external spiral galaxies are likely biased toward the massive ($\geq 10^{5}\,M_{\odot}$) population of molecular clouds except for the case of M33 (e.g., a small GMC down to $2.4\times 10^{4}M_{\odot}$ is recovered by Corbelli et al., 2017). To understand the complex hierarchical structures of molecular gas and also to understand the evolution of molecular clouds in galaxies, smaller and less massive ($<10^{5}\,M_{\odot}$) molecular clouds should be investigated (e.g., the slope of molecular cloud mass functions changes with evolution processes; Kobayashi et al., 2017, 2018). Thus, we need further molecular-cloud surveys covering such less massive clouds in nearby galaxies as a complementary study to PHANGS-ALMA survey.

In this paper, we present the results of a new CO($J=2-1$) survey toward almost the whole molecular-gas disk of M33 conducted with the Atacama Compact Array (ACA) stand-alone mode of ALMA. The distance to M33 is estimated to be 840 kpc (Freedman et al., 1991; Galleti et al., 2004); thus, 1$\arcsec$ corresponds to 4 pc. The inclination of M33 is 55^∘ (Koch et al., 2018). Its proximity and relatively small inclination angle have enabled many researchers to study the ISM and high-mass star formation over the wide area of the M33 disk at a few $\times$ 10 pc scale (Engargiola et al., 2003; Rosolowsky et al., 2007; Onodera et al., 2010; Gratier et al., 2010; Tosaki et al., 2011; Gratier et al., 2012; Miura et al., 2012; Onodera et al., 2012; Miura et al., 2014; Druard et al., 2014; Gratier et al., 2017; Corbelli et al., 2017). In addition, recent studies based on ALMA 12-m array observations revealed complicated internal molecular-gas structures within some especially massive ($\sim 10^{6}\,M_{\odot}$) GMCs of M33 at 1 – 2 pc scale (e.g., Tokuda et al., 2020; Muraoka et al., 2020; Kondo et al., 2021; Sano et al., 2021). Thus, M33 is a unique target to investigate the hierarchical structure of molecular gas in face-on spiral galaxies from parsec to kiloparsec scales. The basic properties of M33 are summarized in Table 1. The main purposes of the new ACA observations are to obtain the spatial distribution in CO($J=2-1$) emission with the higher sensitivity and higher angular resolution compared to earlier studies in M33 and to identify low-mass ($<10^{5}\,M_{\odot}$) clouds as well as high-mass ($\geq 10^{5}\,M_{\odot}$) clouds. This surely becomes an important step to understand the hierarchical structures of molecular gas and the evolution of molecular clouds in galaxies.

The structure of this paper is as follows. In Section 2, we describe the detail of the ACA observations and data reduction. Then, we present the overall molecular-gas structures in CO($J=2-1$) emission at $\sim$30 pc resolution in M33 in Section 3. In Section 4, we describe the procedure of cloud decomposition based on ¹²CO($J=2-1$) cube data, and summarize the basic properties of cataloged molecular clouds. In Section 5, we examine the scaling relations for the molecular clouds in M33. We compare the cataloged molecular clouds in this study with the earlier GMC catalog in M33 summarized by Corbelli et al. (2017) in Section 6. Finally, we discuss the relationship between the properties of molecular clouds and the high-mass star formation in M33 in Section 7.

Observations toward M33 were carried out in Band 6 (211 – 275 GHz) with the ACA 7 m antennas between 2019 August and 2021 August (project code 2018.A.00058.S). The target molecular lines were ¹²CO($J=2-1$), ¹³CO($J=2-1$), and C¹⁸O($J=2-1$). The bandwidths of the correlator settings were 117.19 MHz with 1920 channels for the ¹²CO line and 960 channels for ¹³CO and C¹⁸O lines. The target field was the rectangle with the size of 1100$\arcsec$ $\times$ 1180$\arcsec$ (4.5 kpc $\times$ 4.8 kpc), covering most of the molecular-gas disk of M33. The total number of mosaic fields is 3129. In addition to this, we retrieved the ALMA archival data (project code 2017.1.00901.S and 2019.1.01182.S), which also observed the molecular-gas disk of M33 by ACA 7 m antennas with almost the same spectral settings as our observations. Prior to the imaging process, we concatenated all visibilities obtained in each science goal with a total number of 36. This data reduction strategy is the same as the previously published large-scale ACA mapping project on the Small Magellanic Cloud (SMC) (Tokuda et al., 2021). Figure 1 shows the eventually observed field.

We used Common Astronomy Software Application (CASA) package (McMullin et al., 2007) version 5.4.0 in the data reduction. We applied the standard calibration scheme provided by the ALMA observatory while we performed the imaging process. We used the tclean task with multi-scale deconvolver (Kepley et al., 2020) to recover extended emission as much as possible. In tclean task, we applied the natural weighting and used the auto-multithresh procedure to identify automatically regions containing emission in the dirty and residual images. We continued the deconvolution process until the intensity of the residual image attains the $\sim$1 $\sigma$ noise level. The beam size and the rms noise level for each emission are summarized in Table 2.

To evaluate the missing flux of the ACA observations, we measured the global ¹²CO($J=2-1$) luminosities over the M33 disk obtained by the ACA 7 m antennas and by the IRAM 30 m telescope (Druard et al., 2014). We found the global ¹²CO($J=2-1$) luminosity with ACA 7 m antennas $L_{\rm CO}^{\rm ACA}=7.6\times 10^{6}$ K km s^-1 pc² over the observed region, and that with the IRAM 30 m telescope $L_{\rm CO}^{\rm IRAM}=2.1\times 10^{7}$ K km s^-1 pc² for the same area. This indicates the global missing flux of ¹²CO($J=2-1$) emission of $60-70$%, which mainly corresponds to diffuse components of molecular gas. To compensate for such diffuse components, we combined the ACA 7 m-array ¹²CO($J=2-1$) data with the IRAM 30 m data using the feather task. Hereafter, we refer to the pre-combined ACA 7 m-array ¹²CO($J=2-1$) data as “stand-alone ACA ¹²CO($J=2-1$)” data, and to the combined ¹²CO($J=2-1$) data as “ACA+IRAM ¹²CO($J=2-1$)” data. The beam size and the rms noise level of ACA+IRAM ¹²CO($J=2-1$) data are the same as those of the stand-alone ACA ¹²CO($J=2-1$) data.

From the reduced three-dimensional cube data, we examine the zeroth moment (i.e., velocity-integrated intensity) in ¹²CO($J=2-1$) and ¹³CO($J=2-1$) emission. To minimize the effect of the noise, we determined the velocity channel in which the CO emission is expected to appear using the atomic (H i) data (Koch et al., 2018) as follows. Firstly, we convolved the H i data whose original angular resolution is 20$\arcsec$ to 40$\arcsec$ in order to reduce the effect of the anomalous H i velocity components. Then, we regridded them to match our CO data and determined the representative H i velocity $V_{\rm rep}$ in each pixel. Finally, we calculated the zeroth moment in ¹²CO($J=2-1$) emission from $V_{\rm rep}$ $-$ 30 km s^-1 to $V_{\rm rep}$ $+$ 30 km s^-1. Although Rosolowsky et al. (2007) reported that 90% of the velocity separation between CO and H i is within 20 km s^-1, each CO line typically has a velocity width of 5 – 10 km s^-1. In fact, we found that some molecular clouds dropped $\sim$30% of CO flux if we apply the velocity range of $V_{\rm rep}$ $\pm$ 20 km s^-1. To correctly measure the CO intensity in M33, we needed the velocity range of $V_{\rm rep}$ $\pm$ 30 km s^-1 for the calculation of the ¹²CO($J=2-1$) zeroth moment. We also calculated the ¹³CO($J=2-1$) zeroth moment from $V_{\rm rep}$ $-$ 30 km s^-1 to $V_{\rm rep}$ $+$ 30 km s^-1.

Figure 2 shows the integrated intensity maps in ¹²CO($J=2-1$) from the stand-alone ACA data and the ACA+IRAM data, respectively. These ¹²CO($J=2-1$) maps clearly depict the molecular-gas structure within M33 at 30 pc resolution. We can easily find a lot of individual molecular clouds over the M33 disk. The ACA+IRAM ¹²CO($J=2-1$) map properly recovers diffuse components of molecular gas, which are missed in the stand-alone ACA map. We show an evident case, ¹²CO($J=2-1$) integrated intensity map for GMCs associated with the giant H ii region NGC 604, in Figure 3.

The integrated intensity map in ¹³CO($J=2-1$) emission over the M33 disk is shown in the left panel of Figure 4. A lot of ¹³CO($J=2-1$) sources are detected. They correspond to moderately dense gas whose density is $\gtrsim$10³ cm^-3 within the ¹²CO cloud. The zoomed-in view of the NGC 604 region is shown in the right panel of Figure 4.

Note that we found no significant C¹⁸O($J=2-1$) emission in the ACA map. The rms noise level of 22 mK yields a 3 $\sigma$ upper limit of 66 mK. To check the validity of the upper limit, we retrieved the ALMA archival data (project code 2017.1.00461.S) and examined the C¹⁸O($J=2-1$) emission in a GMC associated with NGC 604. We found that the peak temperature of the strongest C¹⁸O($J=2-1$) emission is $\sim 1$ K at an angular resolution of 0$\farcs$3 (1.2 pc) and its spatial extent is less than 1$\arcsec$. Then, we convolved the C¹⁸O($J=2-1$) emission to 7$\farcs$5 and found that the peak temperature decreases down to $\sim 30$ mK, which corresponds to 1.4 $\sigma$ in the ACA C¹⁸O($J=2-1$) map. Thus, we consider that the beam smearing effect makes C¹⁸O($J=2-1$) emission undetectable in the ACA 30 pc-resolution map.

As shown in Figure 2, the structure of molecular clouds is highly complex and hierarchical over the M33 disk. To identify individual emission structures in an objective way and to investigate the properties of molecular clouds, we employed PYCPROPS (Rosolowsky et al., 2021), a Python implementation of the algorithm to catalog molecular clouds, CPROPS (Rosolowsky & Leroy, 2006). We use the ACA+IRAM ¹²CO($J=2-1$) data in the following analyses.

Firstly, we convolved the ACA+IRAM ¹²CO($J=2-1$) cube data to 7$\farcs$5 in order to identify molecular clouds by a circular beam. Then, we made masked cube data using cprops_mask provided by Rosolowsky et al. (2021). The criteria of the emission mask are that high significance emission is required to be more than 4 $\sigma$ in continuous three velocity channels and low significance emission which is adjacent to high significance emission is required to be more than 3 $\sigma$ in continuous three velocity channels at least over the size of the ACA 7$\farcs$5 beam. Although the default settings of cprops_mask are 4 $\sigma$ and 2 $\sigma$ for high and low significance masks, respectively, we found that the emission masks with these default settings are not suitable for the ACA+IRAM ¹²CO($J=2-1$) cube data. In particular, the low significance mask does not reject fake emission (i.e., noise) at the cloud edge. Thus, we carefully tuned the rms thresholds, and finally we adopted 3 $\sigma$ for low significance masks.

PYCPROPS firstly searches for all local maxima in the emission-masked cube data and measures the peak temperature in each local maximum, $T_{\rm max}$. When a local maximum has at least one other neighbor whose peak temperature is $T_{\rm merge}$, PYCPROPS compares $T_{\rm max}$ and $T_{\rm merge}$. The neighbor is rejected if $T_{\rm max}-T_{\rm merge}$ is less than 2 $\sigma$, which means that such a local maximum is likely a noise fluctuation. The criterion of 2 $\sigma$ is a default value recommended by Rosolowsky et al. (2008). Then, PYCPROPS determines if the spatial and spectral separations between local maxima are adequate or not. We adopted the beam size (7$\farcs$5) as a minimum spatial separation and also adopted 7 km s^-1 as a minimum spectral separation, which corresponds to a typical velocity width of a GMC with the size of $\sim$30 pc considering the Galactic size-linewidth relation (Solomon et al., 1987). If either spatial separation or spectral separation between local maxima does not satisfy the above threshold, the local maximum which has a smaller $T_{\rm max}$ is rejected. Through these processes, PYCPROPS identifies a set of significant local maxima. We treat these local maxima as seeds to assign all the emission to molecular clouds. To do this, we use a watershed algorithm, which associates all the cube pixels in the emission-masked data with a local maximum. Some pixels are already assigned to a single local maximum, while the remainder (including rejected local maxima in the above processes) are assigned to any of the local maxima by the watershed algorithm. More details on the PYCPROPS algorithm are summarized in Rosolowsky et al. (2021). Finally, PYCPROPS identified 886 molecular clouds.

PYCPROPS gives the basic properties of the identified clouds, including the extrapolated 2nd moment of the emission along the major and minor axes $\sigma_{\rm maj}$ and $\sigma_{\rm min}$ in parsec, the position angle of the major axis $\phi$, the extrapolated velocity dispersion $\sigma_{v,{\rm ext}}$, and the integrated ¹²CO($J=2-1$) flux $S$ within each cloud. The extrapolated cloud properties are calculated to reduce observational bias (Rosolowsky & Leroy, 2006). In Figure 5, we showed frequency distributions of the ratio between the extrapolated value and the observed value for cloud size and velocity dispersion. The extrapolated value is typically 10 – 20% larger than the observed value.

We calculated the intrinsic spherical radius $R$ by the deconvolution of the ACA+IRAM beam, $\sigma_{\rm beam}$, as follows:

Here, $\sigma_{\rm beam}$ is calculated as $7\farcs 5\times 4$ = 30, where the factor of 4 is the spatial size in parsec of 1$\arcsec$ at the distance of M33 (840 kpc). The coefficient 1.91 converts the rms size to the effective spherical radius of the cloud (e.g., Solomon et al., 1987; Rosolowsky & Leroy, 2006). We treat $R$ as the intrinsic spherical radius of each molecular cloud. Note that if $\sigma_{\rm min}$ is smaller than $\sigma_{\rm beam}$, the resultant $R$ is not properly defined. We do not consider such small clouds further in this paper.

We also deconvolved the velocity dispersion $\sigma_{v,{\rm ext}}$ as follows:

where $\sigma_{v,{\rm chan}}$ is the velocity resolution element, which is related to the velocity channel width ($\Delta V_{\rm chan}$ = 0.7 km s^-1) as $\sigma_{v,{\rm chan}}=\Delta V_{\rm chan}/(2\sqrt{2\,{\rm ln}\,2})$.

We do not consider molecular clouds lying at the edge of the ACA field-of-view (FOV) further because the primary beam correction causes larger uncertainties in the obtained properties of molecular clouds. Thus, we consider 848 molecular clouds after excluding small clouds with undefined radius and clouds at the edge of ACA FOV from the originally identified clouds by PYCPROPS. Figure 6 shows the spatial distribution of 848 molecular clouds, whose deconvolved sizes and the measured position angle are represented, in the M33 disk.

Starting with the pioneering work by Larson (1981), a lot of earlier studies suggest that basic properties of molecular clouds are quantitatively related to each other via some kind of scaling relations, which are often referred to as “Larson’s laws”. In this section, we examine such scaling relations based on our molecular-cloud catalog.

Based on the ACA+IRAM ¹²CO($J=2-1$) data of M33, we cataloged 848 molecular clouds. In this section, we compare the ACA+IRAM cloud catalog with the earlier GMC catalog generated by Corbelli et al. (2017). To give a fair comparison between the two catalogs, we extracted 362 GMCs from the Corbelli et al. (2017) catalog which are located within the ACA FOV. Hereafter, we refer to the GMCs in the Corbelli et al. (2017) catalog as “IRAM GMCs”. This comparison between the two catalogs provides new insights into the hierarchical structure of molecular gas.

As shown in Figure 6 (and also described in subsection 4.1), many high-mass clouds are associated with the strong 8 $\mu$m emission in the spiral arm region, while low-mass clouds tend to be apart from such $8\,\mu$m-bright sources and to exist in the inter-arm region. Here, we quantitatively evaluate the relationship between the molecular clouds and the 8 $\mu$m-bright sources. To do this, we regridded the IRAC 8 $\mu$m map (Dale et al., 2009) to match the ACA+IRAM ¹²CO($J=2-1$) map and obtained mean $8\,\mu$m flux by averaging the pixel values included within each molecular cloud. Then, we constructed the cumulative distribution functions (CDFs) both for 423 high-mass clouds and 425 low-mass clouds. Figure 15 clearly shows that the $8\,\mu$m-bright sources are closely associated with high-mass clouds rather than low-mass clouds; the strong ($>2$ MJy sr^-1) $8\,\mu$m emission is found in 72% of high-mass clouds, but only in 36% of low-mass clouds, respectively. Note that this trend does not change even if we exclude the diffuse components of 8 $\mu$m emission (see appendix). Our result indicates that high-mass star formation tends to be associated with high-mass clouds rather than low-mass clouds. Such a trend is consistent with the extensive study by Corbelli et al. (2017); they identified mid-infrared (MIR) emission with GMCs and found that a GMC with bright MIR sources tends to have a large CO luminosity mass.

Since high-mass star formation generally starts from the gravitational instability of molecular gas, the virial parameter $\alpha$, which expresses the degree of gravitational binding, is useful to examine star formation in molecular clouds. Figure 16 shows $\alpha$ as a function of $M_{\rm CO}$ for each ACA+IRAM cloud in M33. $\alpha$ generally decreases (i.e., becomes more unstable against gravitational collapse) with the increase in $M_{\rm CO}$. A similar trend is also observed in the MW (e.g., Miville-Deschênes et al., 2017) and external galaxies (IC 342; Hirota et al., 2011). In M33, high-mass clouds whose mass is larger than $10^{5}$ $M_{\odot}$ seem to be almost virialized; in other words, the self-gravitation is predominant rather than the internal turbulence of the cloud. This indicates that the high-mass star formation likely onsets within such high-mass clouds by the gravitational instability, which is well consistent with the observed feature; many high-mass clouds are associated with $8\,\mu$m-bright sources (Figures 6 and 15). In addition, a large $\alpha$ for the low-S/N clouds with a median of 5.3 can be explained; the low-S/N clouds largely correspond to low-mass clouds, which are gravitationally unbound and not associated with star-forming regions.

Finally, we briefly discuss the evolution of molecular clouds. In M33, many high-mass clouds exist in the spiral arm region, while the inter-arm region is dominated by low-mass clouds. Since the $8\,\mu$m-bright sources are loosely along spiral arms in M33, the stellar potential may play a vital role in the accumulation (and the resultant mass growth) of molecular clouds. The evolution of molecular clouds crossing the spiral arm and the high-mass star formation within them are often discussed for grand-design spiral galaxies such as M51 (Egusa et al., 2011) and IC 342 (Hirota et al., 2011) based on the interferometric CO($J=1-0$) observations at a spatial resolution of a few $\times$ 10 pc. In M51, Egusa et al. (2011) suggested that smaller molecular clouds collide to form smooth giant molecular associations (GMAs) at spiral arm regions and then star formation is triggered in the GMA cores. Hirota et al. (2011) divided the GMCs in the spiral arm of IC 342 into two categories according to whether they are associated with star formation activity or not, and reported that the GMCs with H ii regions are typically more virialized and massive compared to the GMCs without H ii regions. These results are consistent with the picture of molecular clouds and the high-mass star formation in M33 although it is a flocculent galaxy whose spiral arm structures are relatively weak. In a forthcoming paper, we will report a detailed study on the evolutionary stage of GMCs based on the comparison with H ii regions (Konishi al., in preparation). Although the GMC evolution in M33 was investigated in earlier studies (e.g., Miura et al., 2012; Corbelli et al., 2017), the new ACA CO($J=2-1$) data enable us to study the dense-gas formation based on ¹³CO($J=2-1$) emission as well as the evolution of basic properties of clouds (e.g., size, linewidth, mass, and virial parameters).

We have performed ALMA-ACA 7 m-array observations in ¹²CO($J=2-1$), ¹³CO($J=2-1$), and C¹⁸O($J=2-1$) line emission toward the molecular-gas disk in M33 at an angular resolution of 7$\farcs$31 $\times$ 6$\farcs$50 (30 pc $\times$ 26 pc). We combined the ACA 7 m-array ¹²CO($J=2-1$) data with the IRAM 30 m data to compensate for diffuse molecular-gas components. The summary of this work is as follows:

1. 1.

   The ACA+IRAM combined ¹²CO($J=2-1$) map clearly depicts the cloud-scale molecular-gas structure over the M33 disk. In addition, we detected a lot of ¹³CO($J=2-1$) sources which correspond to moderately dense molecular gas.

2. 2.

   We decomposed individual cloud components from the ACA+IRAM ¹²CO($J=2-1$) cube data employing PYCPROPS, and cataloged 848 molecular clouds with a mass range from $10^{3}$ $M_{\odot}$ to $10^{6}$ $M_{\odot}$. We found that high-mass clouds ($M_{\rm CO}\geq 10^{5}\,M_{\odot}$) tend to associate with the $8\,\mu$m-bright sources in the spiral arm region, while low-mass clouds ($M_{\rm CO}<10^{5}\,M_{\odot}$) tend to be apart from such $8\,\mu$m-bright sources and to exist in the inter-arm region.

3. 3.

   We found that most of the molecular clouds in M33 show smaller velocity dispersions than the Galactic $R-\sigma_{v}$ relation at a given radius. This is presumably due to low-surface density molecular clouds, which may be maintained even by the small turbulence.

4. 4.

   We found that a small IRAM GMC is identified as a single molecular cloud even in ACA+IRAM CO data, while a large IRAM GMC (typically its $M_{\rm CO}$ is larger than $3\times 10^{5}M_{\odot}$) can be resolved into multiple ACA+IRAM clouds. The velocity dispersion of a large IRAM GMC is mainly dominated by the line-of-sight velocity difference between small clouds inside the GMC rather than the internal cloud velocity broadening.

5. 5.

   Based on the comparison between $M_{\rm CO}$ and $M_{\rm Vir}$ for ACA+IRAM clouds, we found that high-mass clouds in M33 are almost virialized. This indicates that the high-mass star formation likely onsets within such high-mass clouds by the gravitational instability.