© The Authors 2024

1. Introduction

Molecular clouds (MCs) are the birthplaces of stars in galaxies and their formation is a complex process influenced by various physical mechanisms. One key process is the gravitational collapse of dense regions within the interstellar medium (ISM) of galaxies. These regions often arise from instabilities within the ISM, which are triggered by processes such as spiral density waves and stellar feedback (McKee & Ostriker 2007; Dobbs et al. 2014; Renaud et al. 2015). Spiral density waves in galaxies like the Milky Way are mediated by gravitational interactions between stars, gas, and dark matter in the galactic disk. As these waves propagate through the disk, they compress and shock the gas along the trailing edge of a spiral arm (Fujimoto 1968; Roberts 1969), leading to the formation of dense MCs. These clouds serve as sites for star formation (SF) due to their high density and low temperatures. Consequently, spiral arms are expected to exhibit a higher star formation efficiency (SFE) than less dense galaxy regions (Lord & Young 1990; Silva-Villa & Larsen 2012; Yu et al. 2021). Indeed, a greater number of young stars are found in the spiral arms, suggesting a higher star formation rate (SFR) in these areas (Bigiel et al. 2008; Schinnerer et al. 2013). However, the rise in the SFR in spiral arms may simply result from higher surface densities, with the spiral’s gravitational potential restructuring and concentrating the gas rather than influencing SF directly. Elmegreen & Elmegreen (1986) and Querejeta et al. (2024) found no increase in the SFR in galaxies with strong spiral patterns. If this is true in general, the SFE should remain constant regardless of the galaxy, as various studies have noted (Moore et al. 2012; Ragan et al. 2016; Urquhart et al. 2021; Querejeta et al. 2021).

Stellar feedback –in the form of stellar winds, supernova explosions, and radiation pressure from massive stars– also plays a significant role in the formation, evolution, and lifetimes of MCs. These processes inject energy and momentum into the ISM, creating turbulence and disrupting the equilibrium of the gas. The compression and turbulence induced by stellar feedback, as well as the shear induced by galactic dynamics (Chevance et al. 2020), can trigger the collapse of MCs, initiating the formation of new stars. However, stellar feedback can also lead to the destruction of MCs (e.g., Bonne et al. 2023). Chevance et al. (2022) suggested that the main causes of cloud destruction in galaxies are early stellar feedback mechanisms, which take place prior to supernova explosions. It is still a matter of debate as to whether SF is more influenced by the environment – for example, central regions versus spiral arms – or by stellar feedback mechanisms (Corbelli et al. 2017; Rey-Raposo et al. 2017; Kruijssen et al. 2019; Chevance et al. 2022; Liu et al. 2022; Choi et al. 2023). Additionally, other factors such as magnetic fields and turbulence within the ISM can influence the formation and evolution of MCs; magnetic fields provide support against gravitational collapse and can regulate the dynamics of the gas, while turbulence contributes to the fragmentation and structure of MCs.

Linking the physical properties of giant molecular clouds (GMCs) in different large-scale galactic environments, such as spiral arms versus central regions, allows the systematic exploration of how the morphology of a galaxy impacts initial SF conditions. Observations on “cloud scales” (Schinnerer & Leroy 2024) of 50 − 100 pc match the sizes of GMCs or giant molecular associations (GMAs), which are up to a few hundred parsecs in size (Nguyen-Luong et al. 2016, NL16 hereafter). The present study focuses mainly on what these latter authors refer to as molecular cloud associations (MCCs) and “mini-starburst” GMCs with an elevated SFR.

Molecular clouds consist mostly of molecular hydrogen H₂; however, it is difficult to detect cool H₂ directly due to its symmetry and small moment of inertia. One approach to determining the H₂ distribution in a galaxy is to employ observations of dust in the far-infrared (FIR) – for example using the Herschel satellite – and to derive a total hydrogen column density map from a spectral energy distribution (SED) fit to the fluxes. The H₂ distribution is then obtained by subtracting the H I component. H₂ maps can also be obtained using carbon monoxide (CO), the second-most abundant molecule, and applying the CO-to-H₂ conversion factor, X_CO (Bolatto et al. 2013). The low-J rotational transitions of CO are established as a good tracer of the cold regions of MCs because these lines have low excitation temperatures (up to a few 10 K) and low critical densities (a few 100 cm⁻³ up to a few 10³ cm⁻³) for collisional excitation.

However, metallicity significantly affects MCs. Lower metallicity leads to less dust and therefore less dust shielding and self-shielding of ultraviolet (UV) radiation. Thus, the far-UV field can penetrate deeper into MCs, photo-dissociating CO, and leaving a larger envelope of emitting ionized carbon (C⁺) around a smaller CO-rich core (Poglitsch et al. 1995; Stark et al. 1997; Wilson 1997; Bolatto et al. 2013). This is intensified by reduced CO self-shielding due to lower CO abundance. H₂ also photo-dissociates via absorption of Lyman-Werner band photons but achieves sufficient column densities to become self-shielded at moderate extinctions (A_V) within C⁺-emitting gas and thus can become a significant mass reservoir. Hence, there exists a substantial molecular hydrogen reservoir outside the CO-emitting area that is referred to as CO-dark H₂ gas (Röllig et al. 2006; Wolfire et al. 2010). This must be taken into account when comparing the mass estimates of dust and CO, as done in this study.

M33, classified as a flocculent Sc-type spiral galaxy, lies at a distance of 847 kpc (Karachentsev et al. 2004), has an inclination of ∼56° (Regan & Vogel 1994), and contains numerous massive SF regions. Our 18.2″ angular resolution corresponds to ∼75 pc, which is the size of GMCs or GMAs in the Milky Way (NL16), allowing us to resolve individual GMCs on a 2D image. Recent interferometric observations (Peltonen et al. 2023; Muraoka et al. 2023) have even higher resolution, resolving < 50 pc scales, but are not discussed here because we focus on large GMCs. The metallicity of M33 was measured using neon and oxygen abundances in H II regions (Willner & Nelson-Patel 2002; Crockett et al. 2006; Rosolowsky & Simon 2008; Magrini et al. 2010) and varies widely, ranging from values comparable to the Milky Way to lower ones (see discussion and references in Magrini et al. 2010). However, the average metallicity is approximately half solar, which is frequently cited in the literature (Gratier et al. 2012; Druard et al. 2014; Corbelli et al. 2019; Kramer et al. 2020), and the total mass (gas and stars) of M33 is ∼10¹¹ M_⊙, roughly 10% of the Milky Way mass (van der Marel et al. 2012; Patel et al. 2018).

The objective of the present paper is to analyze and compare GMC properties based on dust- and CO-derived N_H2 maps, considering the galactocentric radius and environment within M33. We also aim to establish a GMC mass spectrum and compare our findings with Milky Way studies from NL16, Rice et al. (2016) and PHANGS (Leroy et al. 2021), a survey studying galaxy formation and evolution.

The paper is organized as follows. Section 2 summarizes the data and methods for producing H₂ maps at 18.2″ resolution presented by Keilmann et al. (2024, Paper I hereafter). Section 3 introduces the Dendrograms algorithm (Rosolowsky et al. 2008) and presents the extraction results. The equations for determining MC properties (mass, size, density, pressure, etc.) are given in Sect. 4. These results are presented in Sect. 5, where we also compare with Milky Way studies. Section 6 discusses cloud mass distributions and properties with respect to the galactic radius and environment. Section 7 summarizes the paper.

2. Data and methods

In Paper I, we presented two techniques using FIR Herschel data to produce a total hydrogen column density map of M33. The first procedure (Method I) consists of a multiwavelength SED fit and is briefly summarized in the following. A more detailed description of this method can be found in Paper I. The second approach (Method II) mainly served as a consistency check because it uses only the SPIRE 250 μm data to obtain the N_H2 map and does not account for the variable emissivity index β of the dust. Therefore, we only use the N_H2 map derived from Method I for the current study.

We make use of H I data acquired by Koch et al. (2018). The primary benefit of these data lies in the short-spacing corrections, which are absent in the H I data from Gratier et al. (2010) and which were utilized to generate the κ₀ map and final column density maps in Paper I. The incorporation of the H I data from Koch et al. (2018) for this present study has not resulted in any significant changes to the generated maps, especially in the molecular phase of the maps. In these areas, both H I maps are equal to within ∼10%. The deviation increases in the diffuse region of M33 beyond the molecular regions, where GMCs are not detected anyway. All findings and conclusions of the initial Paper I remain unchanged. The updated data products are available at the CDS.

2.1. Herschel dust data

For Method I, we use the level 2.5 archive data at 160, 250, 350 and 500 μm from the Herschel Key Project HerM33es¹ (Kramer et al. 2010). All maps are in units of MJy sr⁻¹ and reprojected to a grid of 6″.

Method I involved several steps. First, we performed spatial decomposition using a modified Planck function to fit the dust temperature and surface density, Σ, to the SED derived from flux densities within the 160 to 500 μm range. We then applied Eq. (17) from Paper I,

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{high}}={\mathrm{\Sigma }}_{500}+({\mathrm{\Sigma }}_{350}-{\mathrm{\Sigma }}_{350}\ast G_{500\_350})+({\mathrm{\Sigma }}_{250}-{\mathrm{\Sigma }}_{250}\ast G_{350\_250}),\end{array}$$(1)

to generate a high-resolution map of gas surface density Σ_high at a resolution of 18.2″. This equation combines surface density distributions at 250 μm, 350 μm and 500 μm, convolved with a Gaussian kernel to the respective resolutions² in Eq. (1). G_{λc_λo} are the Gaussian kernels of width $\sqrt{{{\theta }_c}^2-{{\theta }_o}^2}$ for the convolution, commonly denoted as *. The beam at the required resolution is specified by θ_c and the beam at the original resolution by θ_o, while the index λ_c_λ_o denotes the corresponding wavelengths. SED fitting was conducted for each map using a pixel-by-pixel modified blackbody function with the specific intensity given by Eq. (10) in Paper I with

$$\begin{array}{c}\hfill I_{\nu }={\kappa }_{0,\mathrm{DGR}}{(\lambda /250\mathrm{\mu }\mathrm{m})}^{-\beta }{\mu }_m{m}_{\mathrm{H}}{N}_{\mathrm{H}}{B}_{\nu }(T_d),\end{array}$$(2)

assuming optically thin emission. The mean molecular weight is indicated as μ_m, m_H is the mass of the hydrogen atom, N_H the total hydrogen column density and B_ν(T_d) the Planck function. A dust opacity law similar to Krügel & Siebenmorgen (1994) has been used with

$$\begin{array}{c}\hfill {\kappa }_g(\nu )={\kappa }_{0,\mathrm{DGR}}\times {(\lambda /250\mathrm{\mu }\mathrm{m})}^{-\beta }.\end{array}$$(3)

The dust-to-gas ratio (DGR) is inherently included in our definition of the dust absorption coefficient in units of (cm²/g), which we denote hereafter simply as κ₀, and β is the emissivity index determined by Tabatabaei et al. (2014). We use this β map alongside the dust temperature of the cold component that are given in Figs. 8 and 9, respectively, in Tabatabaei et al. (2014). We derived the dust absorption coefficient κ₀ pixel-by-pixel using Eq. (15) in Paper I. Further details on the determination of κ₀ and β, including interpolation techniques, are described in Section 3.2 of Paper I. Our approach avoids assumptions regarding the X_CO factor or DGR, allowing for a more accurate evaluation of these factors.

Following our application of this technique to M33, we obtained a total column density map of the galaxy with a spatial resolution of ∼75 pc or an angular resolution of 18.2″ (Fig. 1, left). From this map, we derived a molecular gas column density map (Fig. 1, right) by subtracting H I data from the VLA (Koch et al. 2018). The H I data at 12″ angular resolution have been smoothed to 18.2″ and then transformed to column density using Eq. (1) in Paper I, based on Rohlfs & Wilson (1996). The total H₂ gas mass (Fig. 1, right) is 1.6 × 10⁸ M_⊙.

Fig. 1. Total and molecular hydrogen column density maps. Left: High-resolution N(H) total gas column density map obtained from the Herschel flux maps of M33 with 18.2″ angular resolution (indicated by the circle in the lower left corner) using the β map from Tabatabaei et al. (2014). Values below and above a minimum and maximum threshold (1018 cm−2 and 1022 cm−2, respectively) are blanked. Right: High-resolution H2 column density map derived from the total N(H) map by subtracting the H I component. In both maps, the CO contour levels 3 and 6σ of the CO map (Fig. A.1 in Appendix A) are shown.

2.2. IRAM 30 m telescope ¹²CO(2 − 1) data

M33 was observed using the HERA multibeam dual-polarization receiver in the On-The-Fly mapping mode, targeting the ¹²CO(2 − 1) line. The observations were conducted as part of the IRAM 30 m Large Program titled “The complete CO(2 − 1) map of M33” (Gratier et al. 2010; Druard et al. 2014). The CO line-integrated map has been acquired from the https://www.iram.fr/ILPA/LP006/. The archive data have been converted from antenna temperature scale to main beam brightness temperature using a forward efficiency of F_eff = 0.92 and a beam efficiency of B_eff = 0.56 (Druard et al. 2014). We calculated an rms noise level of 0.35 K km s⁻¹, equivalent to 3σ = 1.046 K km s⁻¹, for the smoothed map with a resolution of 18.2″ (Fig. A.1 in Appendix A). We only utilize the IRAM 30 m line-integrated intensity map of CO, since we do not have velocity information from the dust-derived data. Moreover, due to the inclination of 56°, M33 appears mostly face-on, resulting in small line of sight effects. Hence, employing only the line-integrated CO data is justified and serves as a meaningful comparison with the dust-derived data.

2.3. The X_CO conversion factor map

Using dust-derived N_H2 data and IRAM CO data, we computed the X_CO conversion factor map by dividing the N_H2 map by the CO(2 − 1) map, scaled to CO(1 − 0) with a line ratio of CO(2 − 1)/CO(1 − 0) = 0.8 (Druard et al. 2014). This generates a pixel-by-pixel X_CO map, avoiding a uniform value across the galaxy, as often used in the literature (Gratier et al. 2012; Druard et al. 2014; Clark et al. 2023; Muraoka et al. 2023). The X_CO map, based on the dust-derived N_H2 map, is thus influenced by assumptions in the dust-derived results. Figure B.1 shows the X_CO map outlining the main spiral arms of M33. Minimum values in the outer area are about 10¹⁹ cm⁻²/(K km s⁻¹), with maxima in the northern and southern spiral arms exceeding 10²¹ cm⁻²/(K km s⁻¹).

To compare the X_CO values to those reported in the literature for M33, we calculated a single X_CO value in Paper I using a simple mean and a binned histogram with a log-normal fit. We also conducted scatter plots and a radial line profile to investigate radial dependency. All methods show considerable variability in X_CO with no significant correlation with galactocentric radius. However, the simple mean, log-normal fit and scatter plot indicate values approximately at or below the Galactic X_CO, between 1.6 × 10²⁰ cm⁻²/(K km s⁻¹) and 2.1 × 10²⁰ cm⁻²/(K km s⁻¹). Despite the diversity, this spread is expected due to natural fluctuations in the CO-to-H₂ ratio across the interstellar medium, as noted by Bolatto et al. (2013) and recent studies (Ramambason et al. 2024; Chiang et al. 2024).

3. Identification of coherent structures

To derive physical quantities such as sizes, densities, and masses and to compare structures in H₂ column density maps derived from Herschel with CO data from the IRAM 30 m, we need to identify coherent structures in the 2D maps.

The ISM of a galaxy is a complex multi-phase environment, from hot, tenuous atomic gas to cold, dense molecular gas. The simplest approach defines MCs as entities with well-defined borders (Elmegreen 1985) but complex substructure consisting of clumps and filaments (e.g., Stutzki & Guesten 1990; Schneider et al. 2011; Pineda et al. 2023). Extraction algorithms separate dense gas into distinct clouds/clumps, often using velocity information from spectral line observations for statistical analysis. Various methods identify point sources, clumps, clouds and filaments in the Milky Way (Stutzki & Guesten 1990; Williams et al. 1995; Rosolowsky & Leroy 2006; Rosolowsky et al. 2008; Henshaw et al. 2019; Li et al. 2020; Men’shchikov 2021). Li et al. (2020) concluded that GAUSSCLUMPS (Stutzki & Guesten 1990) and Dendrograms (Rosolowsky et al. 2008) perform best in extracting clumps in synthetic data cubes, including noise.

Molecular clouds were extracted from the CO spectral data cube of the IRAM 30 m telescope at original angular resolution by Gratier et al. (2012) and Corbelli et al. (2017) using the CPROPS algorithm (Rosolowsky & Leroy 2006). CPROPS assumes contiguous, bordered emission of clouds by an isosurface in brightness temperature above a threshold, applying moment measurements to derive size, line width and flux from the position-position-velocity data cube. In our study, we employed Dendrograms for extracting GMCs.

3.1. The Dendrograms method and initial values for M33

The Dendrogram algorithm works intrinsically in two and three dimensions. The algorithm searches for the highest value in the map and systematically collects all other data points with lower values as long as three conditions are fulfilled. The first condition is the minimum difference (min_delta) in intensity between two identified peaks, which must be satisfied to consider those as two different structures. This retaining level of two structures is set to 1σ of the rms noise of the maps. We also have tested different values of min_delta, ranging from 1σ to 5σ. The second parameter is the minimum level (min_value) that a structure must have in order to be considered as a coherent clump/cloud. We explored levels from 3σ to 5σ and ultimately chose to set the threshold above the 3σ level of the rms noise, which is ∼6.5 × 10¹⁹ cm⁻² for the dust-derived N_H2 and ∼0.35 K km s⁻¹ for the line-integrated CO map (Gratier et al. 2010; Druard et al. 2014; Keilmann et al. 2024). Our investigations did not reveal notable variations in the final identification of structures, indicating that the Dendrogram results are not substantially influenced by the selection of these two input parameters when applied to our data. The last parameter defines the minimum number of pixels to be considered as a structure (min_npix), which is related to the width of the structure when we assume circular geometry. This parameter is given by the actual number of pixels, which fit into the full width at half maximum (FWHM) beam width multiplied by 1.2 and corresponds to ∼11 pixels to ensure that the MCs are well resolved. We also experimented with factors of 1, 1.5 and 3 (see below). Obviously, a factor of 1 tends to find more smaller structures around the beam size and a factor of 3 “blurs” clouds into fewer but larger structures. We choose a factor of 1.2 as the best compromise to obtain reliable cloud statistics (see also Schneider & Brooks 2004; Kramer et al. 1998). The extraction of sources has then been applied to the dust-derived N_H2 map and to the line-integrated CO(2 − 1) map, for which the final properties of the detected GMCs have been derived using the X_CO factor map from Paper I. We investigate the influence of varying Dendrogram parameters in Appendix C.

Dendrograms distinguish between so called leaves and branches. Branches contain substructures like other branches or leaves, while leaves only consist of themselves. As the distance to M33 is 847 kpc (Karachentsev et al. 2004) and the angular resolution is 18.2″, the minimum resolved convolved structure that we can identify has a size of 75 pc. This corresponds roughly to the size of GMCs (and small GMAs) in the Milky Way (Roman-Duval et al. 2010; Hughes et al. 2013; Nguyen-Luong et al. 2016; Spilker et al. 2022). We mostly concentrate on leaves in our analysis since they best represent a single GMC/GMA. However, leaves do not capture all the emission. Branches, which contain leaves and comprise larger areas, represent the more diffuse, extended H₂ emission that we define as “inter-cloud” medium. An example is the crowded center of M33, where separating individual clouds can become challenging. Branches also include emission just around individual, well-separated GMCs. We refer to this surrounding structure as the “envelope”. We note that this does not constitute an H I envelope. An example of this is the GMC NGC 604 in the northeast spiral arm, where the smaller leave structure is embedded in the larger branch. Unless stated otherwise, branches refer to dust, as CO branches show a relationship similar to dust branches as leaves do.

3.2. Dendrogram extraction of GMCs in the CO and dust-derived N_H2 maps

In the following sections, we discuss the morphology of detected GMCs, NGC 604 and the differing numbers of detected GMCs in dust and CO. We focus on leaf structures only because they represent mostly the GMC/GMA population. Tables 1 and 2 list the main cloud parameters for the first 10 clouds ordered by their surface mass density. See Sect. 4 for details on the calculation of the quantities listed in the tables. The tables show no one-to-one correlation between dust- and CO-derived GMCs. For example, NGC 604 is one single GMC in the dust-derived map, but in CO it splits into two smaller structures. This discrepancy arises because structures identified in both tracers differ in size and mass.

We focus on the distribution statistics in Sect. 5. A summary of the mean values is given in Table 3. The similarity in this table, especially in masses, is because the structures are well identified as leaves in both tracers (see Sect. 3.2). The CO-derived H₂ column density is lower, evident in the central branch masses, where the CO-derived mass is ∼70% of the dust-derived mass. This similarity extends to many parameters, such as radius and densities. Their average values are comparable, varying by less than a factor of 2.

3.2.1. Morphological description

Figure 2 displays the 326 leaf (red) and 142 branch (green) structures identified in the dust-derived N_H2 map (left panel) and the 199 leaf (blue) and 94 branch (green) structures in the CO(2 − 1) map (right panel). To provide a clearer overview, we only show the lowest level of branch extraction, which may include other branches (and leaves). The overlay of CO and dust-derived structures on the X_CO map (Fig. 3) shows a similar morphology for both tracers. However, the dust emission identifies more structures beyond the spiral arms and central region. 199 structures in the CO map are more locally concentrated, with fewer structures between spiral arms or in M33’s outer regions compared to the dust-derived map.

Fig. 2. Dendrogram source extraction from the NH2 and CO maps. Left: GMC detections outlined by the lowest level branches in the dust-derived NH2 map (green contours) and 326 leaf structures (red contours). NGC 604 is marked with a thick pink contour. Right: GMC detections identifying 199 leaf structures in the line-integrated CO(2 − 1) map (blue contours) and similar to the dust map, the branches (green contours). The small circle in the lower left corner of both figures shows the beam size of 18.2″.

Fig. 3. Structures identified using Dendrograms on CO and dust-derived maps superimposed on the XCO map. The structures found in the dust-derived NH2 (red) and in the CO map (blue) are mapped onto the XCO factor (ratio) map, which represents the dust-derived H2 column density over CO intensity from Paper I.

Furthermore, especially the central region in the dust-derived map shows substantial H₂ column density between the identified leaves. The emission distribution contained in branches focuses in the crowded center region of M33, where many GMCs/GMAs potentially overlap along the line of sight, leading to a rather homogeneous plateau of emission. This line of sight effect can be one reason why it is not possible to separate dust emission into smaller leaf structures. However, Koch et al. (2019) employed a Gaussian decomposition on the full spectral CO line cube of the IRAM 30 m data and found that only ∼10% of the CO spectra show multiple components. This finding does not conclusively rule out the possibility that there are several clouds along the line of sight, but overall this effect is probably less important than for other galaxies. It is also plausible that the crowded emission in the center seen in dust constitutes an inter-cloud H₂ medium, similar to what is found in the Milky Way. We note that the flocculent morphology of M33 already points toward an important gas reservoir between the GMCs. However, another explanation is that the inter-cloud gas is warmer and tends to decrease the CO brightness of the low-J lines, which requires future observations of CO(3 − 2) or CO(4 − 3) line emission. In any case, this more widespread gas reservoir in the center contains a significant mass. While dust-derived leaves collectively hold 8.3 × 10⁷ M_⊙ in total, branches excluding leaves contain 3.1 × 10⁷ M_⊙ of the H₂ gas mass. Especially in the central region of M33, the leaves contain 2.6 × 10⁷ M_⊙, while the branches in the center comprise 1.5 × 10⁷ M_⊙.

The CO emission map (Fig. 2, right) shows less homogeneous material in branches than the dust-derived map and reveals that the CO leaves are typically surrounded by a more extended envelope. This finding supports the one of Rosolowsky & Simon (2008), who propose that around 90% of the diffuse emission to within 100 pc of a GMC arises from a population of small, unresolved MCs. However, the CO sensitivity might be too low to detect CO-dark gas or CO might be easily dissociated in the center. Additionally, the H₂ emission from dust can be overestimated due to the complex map production process and the subtraction of a VLA H I map, which has its own detection limits. The dust map might still contain H I, as CO-faint column densities are low (0.5 to 1 × 10²¹ cm⁻²), close to the atomic-to-molecular hydrogen transition level. The CO-identified leaf structures have a total mass of 4.2 × 10⁷ M_⊙, with branches holding 2.1 × 10⁷ M_⊙ (50% of the mass compared to leaves; 37% in the dust-derived map). In the center, the leaves contain 1.7 × 10⁷ M_⊙ and branches 1.1 × 10⁷ M_⊙ (64% of leaves’ mass; 57% in dust-derived map).

3.2.2. The GMA NGC 604

The SF region NGC 604 stands out with the highest mass and largest area, forming a single structure on the dust-derived N_H2 map but several GMCs on the CO map (Fig. 2), similar to the findings of Williams et al. (2019). The discrepancy may arise from the greater extent of the GMC in dust compared to CO, as the 3σ CO signal shows a narrower north-south ridge (see Fig. 2). Another explanation could be CO-dark gas in the enveloping gas, with dust emission reaching 2 × 10²⁰ cm⁻², conducive to CO formation. This and the limited spatial resolution probably explain the divergence of NGC 604 from the majority of the GMC population in this study.

However, Relaño et al. (2013) (and references therein) reported that NGC 604 is not a single H II region, but comprised of a few compact knots and filamentary structures joining the different knots. The whole complex has a radius of 280 pc and forms the second most luminous H II regions association in the Local Group, surpassed only by 30 Doradus in the LMC.

Observations of the Atacama Large Millimeter/submillimeter Array (ALMA) in ¹²CO(2 − 1) and ¹³CO(1 − 0) (Muraoka et al. 2020; Phiri et al. 2021; Peltonen et al. 2023) at an angular resolution of 0.44″ × 0.27″ (1.8 pc × 1.1 pc) (Muraoka et al. 2020) and 3.2″ × 2.4″ (13 pc × 10 pc) (Phiri et al. 2021) confirm that NGC 604 constitutes multiple individual molecular clouds.

3.2.3. Caveats regarding dust- and CO-derived GMCs

Some GMCs are identified only in the CO dataset, while others appear only in the dust-derived dataset. Dust-only detections may indicate CO-dark H₂ gas or may be due to smaller CO(2 − 1) envelopes given its higher critical density compared to CO(1 − 0) (Schinnerer & Leroy 2024). Regions seen only in CO may reflect underestimated molecular hydrogen column densities, possibly from overestimated κ₀ values derived from molecular hydrogen regions, casting doubts on the assumption of a constant κ₀ between atomic and molecular phases.

Furthermore, κ₀ may be overestimated as it requires regions with CO emission below 2σ, hence assuming no CO emission. This leads to a bias due to generally low CO emission in the disk (see Eq. (16) in Paper I). The IRAM CO map might show structures from noise fluctuations. Raising the detection threshold to 5σ can address this potential issue, but this approach also leads to similar detections when consistently increasing the threshold for dust-derived data. The uncertainty in H₂ detection and the prevalent use of the 3σ threshold for CO and dust data make it challenging to conclusively ascertain the origin of these structures.

Increased noise in the dust-derived N_H2 map may cause structures experience a quasi “beam-diluted” effect and to blend into the background due to too low N_H2 (not fulfilling min_value) or failing the minimum size condition (min_npix) to be identified. This aligns with the observation of low column densities in both dust-derived and CO-derived N_H2 maps. Figure 3 supports this, showing that structures detected only in the CO-derived N_H2 data have the lowest X_CO factor. The presence of CO-dark gas and a non-changing κ₀ in both the atomic and molecular phases may explain the greater number of GMCs in the dust-derived data.

4. Determination of physical cloud properties

For each identified structure, we compute several parameters such as the area A in pc², the equivalent radius R in pc, the column density N in cm⁻², the (beam-averaged) number density n in cm⁻³ and the mass M in M_⊙ along with pressure P/k_B in cm⁻³ K. The shape of the identified structures is described by the aspect ratio (AR), that is, the ratio between major and minor axes of the GMC. The following section outlines the methodologies employed to compute these quantities.

To calculate A, the area of each pixel of the identified GMC is summed and scaled by the squared distance, D², to M33. This pixel size is denoted A_pixel = dθ_ra ⋅ dθ_dec ⋅ D², where dθ_ra and dθ_dec represent the angular size of a pixel in radians. Dendrograms provides information on the location of the identified structure, which serves as a mask for the original dataset. This allows for the calculation of the number of pixels associated with a structure.

The radius R is determined as the equivalent radius of a circle with the area A of the Dendrogram structure, A = π R². The radius of each structure is de-convolved by $R{\prime }_i=\sqrt{R_{\mathrm{GMC},i}^2-R_{\mathrm{beam}}^2}$, where R_GMC, i represents the radius of the i-th structure and R_beam corresponds to the beam size.

To calculate the column density of H₂ of a structure using CO(2 − 1) data, scaled to CO(1 − 0) using the line ratio of 0.8 (Druard et al. 2014), we consider all pixels from the X_CO factor map that belong to a detected structure in the line-integrated CO intensity map. We then multiply the corresponding X_CO values with the line-integrated intensities of this map on a pixel-by-pixel basis. This approach provides a more precise estimate compared to using a constant X_CO factor for the entire galaxy, and allows us to uncover intriguing variations in the distribution of GMCs within M33, which is further explored and discussed in Sect. 5. Additionally, to obtain an average X_CO value for each structure, we divided the dust-derived H₂ column density by the corresponding CO line-integrated intensity on a pixel-by-pixel basis. These X_CO values for each structure are presented in Table 2 (as discussed in Sect. 5).

To determine the masses of GMCs using the molecular hydrogen column densities obtained from Dendrograms (both from CO and dust), the pixel size A_pixel = dθ_ra ⋅ dθ_dec ⋅ D² is multiplied by N(H₂)_j. Here, j represents the index of a pixel within an identified structure. This is finally multiplied by the hydrogen mass, m_H, and the mean molecular weight, μ,

$$\begin{array}{c}\hfill M_{\mathrm{GMC}}=A_{\mathrm{pixel}}·{\displaystyle \sum _i}N{({\mathrm{H}}_2)}_j·m_{\mathrm{H}}·\mu ,\end{array}$$(4)

with μ = 2.8 (Kauffmann et al. 2008) to account for Helium (He) and metals.

To calculate the average number density n of a GMC consisting of H₂, we assume a spherical configuration and use the mass and the equivalent radius obtained above. The average density n is then determined as

$$\begin{array}{c}\hfill \frac{n({\mathrm{H}}_2+He)}{{\mathrm{cm}}^{-3}}=14.6\frac{M}{{\mathrm{M}}_{\odot }}·{\left(\frac{4\pi }{3}\frac{R^{\prime 3}}{{\mathrm{pc}}^3}\right)}^{-1},\end{array}$$(5)

where M represents the mass of the cloud in solar masses and R′ denotes the de-convolved equivalent radius of the cloud in parsecs³. Since our spatial resolution is 75 pc, the density can only reflect a beam-averaged density derived by dividing the (beam-averaged) column density by the (beam) size. The detected GMC will be composed of smaller substructures with much higher local densities. We note that, due to the critical density, the density of the clouds should be in the order of 10³ cm⁻³ for the low-J CO transitions to be sufficiently excited.

The surface mass density Σ of the GMCs is determined via

$$\begin{array}{c}\hfill \frac{{\mathrm{\Sigma }}_{\mathrm{GMC}}}{{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}}=\left(\frac{M}{{\mathrm{M}}_{\odot }}\right){\left(\frac{A}{{\mathrm{pc}}^2}\right)}^{-1}.\end{array}$$(6)

The elongation of detected GMCs is quantified by their column density-weighted aspect ratio of semi-major to semi-minor axis.

Finally, the CO luminosity, L_CO, is calculated as

$$\begin{array}{cc}\hfill \frac{L_{\mathrm{CO}}}{\mathrm{K}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{pc}}^2}& =\frac{D^2}{{\mathrm{pc}}^2}{\left(\frac{\pi }{180·3600}\right)}^2{\displaystyle \sum _i}{I}_{\mathrm{int}.}\mathrm{d}{\theta }_{\mathrm{ra}}\mathrm{d}{\theta }_{\mathrm{dec}},\hfill \end{array}$$(7)

where D represents the distance to M33 in parsecs, dθ_ra and dθ_dec are the spatial dimensions of a pixel given in radian and the integrated intensity I_int. in units of K km s⁻¹. Hence, we sum over all pixels of a GMC, multiply each pixel by the size of a pixel, scale it by the distance and convert it to units of radians. The resulting L_CO therefore describes the integrated emission inside a GMC summed over its entire area.

It is formally possible to calculate the gas pressure within the GMCs, P, using the equation

$$\begin{array}{c}\hfill \frac{P/k_{\mathrm{B}}}{{\mathrm{cm}}^{-3}\mathrm{K}}=\frac{n\mu }{{\mathrm{cm}}^{-3}}\frac{T}{\mathrm{K}},\end{array}$$(8)

where k_B is the Boltzmann constant, n the number density of Eq. (5) and T the “mass-weighted” dust temperature. We note, however, that the latter is only the cold component of the dust temperature (Fig. 9 (left) of Tabatabaei et al. 2014), which has also been used to produce the column density maps in Paper I. The dust temperature is around 20 K, which corresponds to a similar gas temperature only when gas and dust are thermally well coupled by collisions (Goldsmith 2001; Goicoechea et al. 2016), which is the case in cool, dense regions. In less dense regions, gas and dust temperatures can differ significantly due to the inefficiency of collisional energy transfer. We thus expect a difference between the inner regions within a GMC, with typical temperatures of around 10 − 20 K, and the inter-cloud gas, which can be significantly higher, corresponding to gas temperatures of > 100 K. In addition, there is possible unresolved substructure in the beam and the density is rather low because of the beam-averaging. Thus, it is not surprising that our pressure results in lower values compared to Hughes et al. (2013), Sun et al. (2020a) and Sun et al. (2020b). The latter two use velocity-resolved CO data from which they obtain higher pressures. The pressure values we derive are therefore only valid for the cool, molecular GMCs and we do not go into great detail in our interpretation.

5. Dendrogram analysis and comparison with the Milky Way

In this section, we discuss the distributions of the key physical cloud properties from the Dendrogram leaves extraction of the dust-derived N_H2 and CO maps individually (Sects. 5.1 and 5.2). We compare our results with the Milky Way GMC statistics from Rice et al. (2016) and NL16 that rely on the Columbia (CfA) ¹²CO(1 − 0) survey (Cohen et al. 1986; Dame et al. 2001), which provides the most comprehensive Milky Way GMC catalog. NL16 derived the cloud properties by eye inspection of line-integrated CO maps and focuses on large (R ⪆ 50 pc) and massive (M ⪆ 10⁶ M_⊙) MCCs. Those MCCs with a SFR larger than 1 M_⊙ yr⁻¹ kpc⁻² are called “mini-starbursts”, an example is the W43 region. However, what we find in M33 are more MCCs without a high SFR (Corbelli et al. 2017), though the SFR densities of MCCs are comparable with the SFR of super giant H II regions in M33 (Miura et al. 2014). Rice et al. (2016) performed a Dendrogram analysis on the velocity-resolved CO data and also included smaller (< 50 pc) and lighter MCs, down to a limit of a few 10⁴ M_⊙, which are beyond our resolutions. Since there are other CO surveys of the Milky Way with extensive datasets, we also partly compare our findings with those. Nevertheless, these studies mainly detect smaller molecular clouds, posing a challenge in making meaningful comparisons with the GMCs we can resolve. For a comprehensive overview of the current CO surveys of regions in the Milky Way, see Park et al. (2023). The most relevant studies utilize data from the Galactic Ring Survey (GRS); see Simon et al. (2001) and Roman-Duval et al. (2010) and cloud compilations presented in (e.g., Kramer et al. 1998; Schneider & Brooks 2004; Su et al. 2019).

We also compare our results with Dobbs et al. (2019), who studied molecular clouds in a simulation of a M33-type galaxy and from the same IRAM CO data of M33 we use. Their models, based on Smooth Particle Hydrodynamic (SPH) codes SPHNG (Bate et al. 1995) and GASOLINE2 (Wadsley et al. 2017), are detailed in Dobbs et al. (2018). They used Friends-of-Friends (FoF) and CPROPS algorithms to determine cloud properties of the simulations.

For completeness and to compare with other studies, Appendix D shows and discusses the ¹²CO(1 − 0) luminosity of M33.

5.1. GMC radii

The calculated mean of the beam-deconvolved cloud equivalent radius reveal a similar overall distribution and mean values of around 58 ± 13 pc and 68 ± 21 pc for dust- and CO-derived GMCs, respectively. The largest structure observed from dust data (NGC 604) exhibits the most notable difference, featuring a radius of approximately 200 pc. For completeness, we report that the branches have a mean radius of 354 ± 152 pc.

For comparison, Gratier et al. (2012) and Corbelli et al. (2017) found mean radii of 42 ± 13 pc and 45 ± 12 pc, respectively, from the IRAM CO map using CPROPS. These sizes are smaller compared to our findings, primarily due to the higher spatial resolution of 12″ of the unsmoothed IRAM CO map they used. Williams et al. (2019) report a median GMC size of 105 pc for their identified GMCs in M33, while we derived a mean value for this catalog of 116 ± 29 pc. They identified with Dendrogram the clouds in the SPIRE 250 μm map at 18″ resolution and then performed an SED fit on the averaged flux values of 160, 250, 450 and 850 μm within one identified structure and determined the cloud mass with a fitted DGR and X_CO factor. They find an X_CO value of ∼6 × 10²⁰ cm⁻²/(K km s⁻¹) by fixing the dust absorption coefficient κ₀ and emissivity index β and bin the GMCs at 500 pc. A DGR and X_CO factor are radially determined via scatter analysis fitting both parameters simultaneously, resulting in possibly degenerate values since different combinations can lead to the same result (their Eq. (6)). They subtracted an H I map without short-spacing from Gratier et al. (2012). A source extraction was also performed on the higher resolution (13″) 450 μm map, reporting a similar size distribution of the clouds.

Rice et al. (2016) has the most complete GMC catalog of the Milky Way, with mean radii of 34 ± 6 pc. Since the subset of NL16 only concentrates on large and massive clouds, the mean value is higher around 90 ± 20 pc (see Table 3 and Fig. 4). Despite this, the trend of the distribution closely mirrors the patterns observed in the M33 data derived from dust and CO data for larger GMCs. There appears to be a size limit of around 150 pc for the largest GMCs/GMAs, in the Milky Way as well as in M33, though both galaxies are different in terms of size, mass and age. Interestingly, Dobbs et al. (2019) find a similar threshold in their simulations of M33 and their cloud extractions of the IRAM CO data set (their Fig. 4). The three distributions exhibit a comparable decline in both the shape and the number of structures as they increase in size. We further do not clearly detect Giant Molecular Filaments (GMFs), which can reach lengths of up to 200 pc in the Milky Way (Wang et al. 2020). However, some of our GMCs have an elongated geometry and aspect ratios larger than 3 so that they formally fit to the definition of GMFs. We come back to this point in Sect. 6.2.5.

Fig. 4. Distributions of GMC properties from our study and the literature. The panels show histograms of radius (top left), surface mass density (top right), mass (bottom left) and beam-averaged number density (bottom right) derived from H2 data from dust (solid red) and CO (solid blue) from this study. The dashed purple and golden lines display the distributions obtained for M33 from dust (Williams et al. 2019) and CO (Corbelli et al. 2017), respectively. The dashed black and light red lines give the distributions for the Milky Way studies of NL16 and Rice et al. (2016), respectively. The vertical black dashed line in the mass distribution signifies the minimum mass limit of the selected structures in the Milky Way by NL16.

A potential mechanism that explains the growth of GMCs in alignment with the results can be attributed to supernovae. Kobayashi et al. (2017) and Kobayashi et al. (2018) show that H I gas is an important mass reservoir for growing GMCs and they show that supernovae can accumulate the H I gas to molecular clouds. In this case, the GMC growth is assumed to depend on the maximum potential sizes of supernovae remnants. Hence, most GMCs are predicted to show sizes of ≲100 pc, with a few exceptions of up to 150 to 200 pc. Another potential explanation is that GMC growth depends on the galactic gas disk scale height, h_z. When a GMC is smaller than h_z, it can grow in all three dimensions. Once it reaches the size of h_z, its ability to grow in the vertical direction will drop. Only the two remaining directions allow the GMCs to expand, but this slows their growth, giving time for stellar feedback or other mechanisms to destroy and regulate cloud sizes. The gas scale height of the galactic disk in the Milky Way ranges from 300 to 400 pc (Carroll & Ostlie 2007). M33 shows a comparable scale height of 320 ± 80 pc (Berkhuijsen et al. 2013). Therefore, this rationale could account for the analogous shape and upper size limit of the largest GMCs in both galaxies. We note, however, that Koch et al. (2019) determined a CO/H I line width ratio of around 0.7 and suggest that M33 has a marginal thick molecular disk and not a thin disk dominated by GMCs and a thicker diffuse molecular disk as seen for the Milky Way and other more massive spirals.

However, we caution that the GMCs identified in M33 potentially have line of sight effects due to limited resolution, the inclination and the increased thickness of the central region, which can blend distinct GMCs into a larger structure that is not one coherent GMA. In addition, as mentioned in Rice et al. (2016), the mass obtained for some GMCs can be inaccurate by up to an order of magnitude due to challenges of reliably determining the correct distances.

5.2. GMC masses and densities

Figure 4 (bottom row) shows the mass and average density distributions of H₂ derived from dust and CO. The black dashed line represents the minimum mass selection used in NL16. The average number density for the binned data set (Fig. 4 bottom right) and the individual clouds (Tables 1 and 2) are low, typically below 30 cm⁻³ for both tracers. The mean of the average densities are similar, with values of n = 5.2 ± 1.5 cm⁻³ for dust-derived and n = 3.0 ± 1.1 cm⁻³ for CO-derived GMCs. Our maps have a spatial resolution of 75 pc, and therefore the identified structures are likely composed of smaller substructures with higher local densities. The densities of GMCs in the Milky Way (29.1 ± 8.0 cm⁻³) have been calculated using the same methodology, based on the data presented in Table 1 of NL16. The branches have a low average density of 1.1 ± 0.4 cm⁻³, which is reasonable given that they span larger areas than the leaves and incorporate a significant amount of inter-cloud and envelope material, both of which are expected to have lower densities.

The mass distributions in M33 (Fig. 4, bottom left) show no significant differences between CO and dust for our study. The maximum GMC mass from CO is ≈5 × 10⁶ M_⊙, whereas for dust it is NGC 604 with ≈8 × 10⁶ M_⊙. We note that there are only a few GMCs in M33 above 10⁶ M_⊙ in both tracers. The mean values derived from the dust data are very similar, with M = (2.8 ± 0.9)×10⁵ M_⊙ compared to M = (2.9 ± 0.9)×10⁵ M_⊙ for CO. The branches have a mean mass of M = (1.3 ± 0.2)×10⁷ M_⊙.

NL16 selected only Milky Way GMCs/GMAs with masses of larger than around 10⁶ M_⊙, and thus it is not surprising that the distribution contains only GMCs in this mass range (GMCs with lower masses are not absent but were not included in the survey). Notably, M33 lacks a significant high-mass GMC population. The procedure for mass determination is the same for our study and that for NL16: for CO, the line-integrated intensity was used to derive the CO column density and then finally the mass using an X_CO factor; and for the dust, the N_H2 column density was derived from an SED fit. Interestingly, the lack of significant high-mass GMCs in M33 also becomes evident by comparing with the comprehensive Rice et al. (2016) Milky Way catalog, which arises from a velocity-based identification of GMCs from CO data, which in addition shows a mean mass similar to our results of (2.4 ± 1.0)×10⁵ M_⊙. The difference in mass thus stems from the lower overall CO luminosity and hydrogen column density in M33.

The H₂ gas mass in the center of M33 (see Fig. F.1 for an outline of the center) is ∼25% of the total dust-derived H₂ mass and amounts to 4.3 × 10⁷ M_⊙. This is an order of magnitude lower than the central molecular zone (CMZ) of the Milky Way (∼1.3 × 10⁸ M_⊙). We note that the overall mass of M33 is one order of magnitude lower than that of the Milky Way. The total H₂ mass of the Milky Way is suggested to be 1.4 times higher than the values found in earlier studies (Sun et al. 2021). Applied to the results reported in García et al. (2014), this leads to a total H₂ mass of 4.2 × 10⁸ M_⊙. Consequently, the proportion of the CMZ of the Milky Way to this mass is ∼30%.

The Milky Way also shows higher number densities, with a mean density of about 30 ± 11 cm⁻³ and 18 ± 6 cm⁻³ for the dataset presented in NL16 and Rice et al. (2016), respectively. Mean values from dust and CO data are roughly five times smaller than those in the Milky Way datasets (and not one order of magnitude). According to the mass-size relation discussed in Sect. 5.3, the density decreases with size. This might explain why mean densities do not show the same trend like GMC masses, central region mass or total H₂ mass of M33, all of which are consistently an order of magnitude lower compared to those of the Milky Way.

Williams et al. (2019), using dust data, found in M33 GMC masses shifted to higher values averaging to ranges from (3.7 ± 1.4)×10⁵ M_⊙ and low mean number densities of 1 ± 0.4 cm⁻³, while the average cloud mass in Gratier et al. (2012) from CO data is (2.4 ± 0.9)×10⁵ M_⊙ with a mean density of 30 ± 7 cm⁻³. Corbelli et al. (2017) find similar results using the same data at the same angular resolution with (2.6 ± 1.1)×10⁵ M_⊙ and 17 ± 5 cm⁻³. The masses match our findings, but the higher number densities are due to detecting smaller structures. This may result from the 12″ resolution of the unsmoothed CO data and a different cloud extraction method (CPROPS).

The effects of limited resolution of our data do not cause non-detections of GMCs with similar masses and densities in M33 compared to the Milky Way. Using a larger beam would inaccurately merge smaller structures into fewer larger structures, consequently inflating the overall mass. The dissimilarity in mass between M33 and the Milky Way, with M33 having only around 10% of the mass of the Milky Way, probably originates from variations in the sizes and evolutionary stages of the galaxies. The diameter of M33 is approximately two-thirds the size of the Milky Way. NL16 used a dataset with an angular resolution of $8\stackrel{\prime }{.}8$, translating into a spatial resolution of ∼60 pc for the most distant GMCs in the Milky Way. For these distant GMCs, our spatial resolution is similar.

By considering sweeping the H I medium by supernovae as we discussed in Sect. 5.1, the typical maximum mass is limited by the gas scale height so that n_ISM ⋅ h_z³, where n_ISM is the volume density of the ISM and h_z is the galactic gas disk scale height. The gas disk scale heights of both galaxies are similar, as discussed above. Thus, the remaining factor influencing the mass growth may be attributed to the density of galaxies. Given that the Milky Way has a higher H₂ density and total mass (from which a greater column density and ultimately a higher number density can be expected), we anticipate that the Milky Way will show higher densities. Therefore, we propose this mechanism as a possible driver. Meanwhile, Kobayashi et al. (2017) and Kobayashi et al. (2018) performed a semi-analytic theory to investigate the impact of cloud-cloud collisions. They show that, even in the Milky Way galaxy, cloud-cloud collisions have a minor impact on GMC growth and are only effective to clouds more massive than 10⁶ M_⊙. We therefore suspect that cloud-cloud collisions are mostly ineffective for M33.

5.3. Mass–size relations

Figure 5 illustrates Larson’s (Larson 1981) relationship between mass and size for data derived from dust and CO. For dust-derived GMCs the slope is 1.8 ± 0.1, while it is 2.2 ± 0.2 for CO-derived GMCs. In the simulations of M33, Dobbs et al. (2019) also found a clear mass-size relation in the observations and simulations. However, they did not quantify the slope of this relation so that we extracted the data points from their Fig. 4 and fitted a linear function with a slope of 1.4 ± 0.1. This estimation carries a significant level of uncertainty, attributed to manual data extraction and the overlapping of numerous data points, preventing clear individual identification. As previously noted in Dobbs et al. (2019), the larger GMCs appear too extended, also compared to the GMCs we identified, leading to a less steep slope.

Fig. 5. Mass–size relation of the identified GMCs. Top: Mass–size Larson relation for GMCs derived from dust. Bottom: Mass–size Larson relation for GMCs derived from CO data. The various colors indicate the average density of individual clouds. The panel displays the slope κ and its corresponding error.

A slope of about 1.6 (Lombardi et al. 2010; Kauffmann et al. 2010) in the mass-size relation suggests the presence of substructures within individual clouds, while a slope of around 2.4 was identified for GMCs in the GRS Galactic plane survey (Roman-Duval et al. 2010). NL16 determined a slope of 1.9 for GMCs and a slope of 2.2 for MCCs. Given that the third Larson relation indicates a power-law connection between mass and size, represented as M ∝ R^κ, with a typical power-law exponent usually around 2, it suggests similar gas surface mass densities for all GMCs. Furthermore, in line with assumptions for a spherical object, the mass can be linked to the size by M ∼ n/R³ ∼ N/R², leading to M ∼ R². This finding aligns well with observations that incorporate a column density threshold (see Schneider & Brooks 2004 and references therein).

5.4. GMC surface mass densities

Comparing cloud masses and sizes across studies can be unreliable due to differing GMC definitions and boundary settings. Resolution limits may also cause undetected clouds or beam smearing. The concept of cloud surface densities inherently considers the cloud size per definition, Σ_GMC = M/A, thereby mitigating the impact of varying resolutions across studies. However, complete resolution uniformity is not achieved for instance when the galaxy is not perfectly face-on, as some large clouds may still merge into one single larger cloud when the beam size is large, resulting in a higher surface mass density. Conversely, smaller clouds, if sufficiently spaced from others, may get smeared within the beam, causing dilution and a decrease in surface mass density. This can be mitigated by excluding too small structures, which we do by only accepting structures 1.2 times the beam size. Nonetheless, comparing surface mass densities can facilitate a less biased evaluation of clouds occupying similar spatial areas.

In Fig. 4, we compare the GMC surface mass densities. For Milky Way GMCs (NL16), the mean value of 187 ± 51 M_⊙ pc² is approximately one order of magnitude higher than our dust-derived (22 ± 5 M_⊙ pc²) and CO-derived values (16 ± 6 M_⊙ pc²) in M33. Whereas compared to the more complete cloud catalog obtained by Rice et al. (2016), the mean value is 38 ± 14 M_⊙ pc², approaching similar high surface densities at the higher end of the spectrum as the clouds presented in NL16. Branches show consistent mean surface mass densities of 19 ± 5 M_⊙ pc².

For comparison, Hughes et al. (2013) report a gas surface mass density for M33 of 46 ± 20 M_⊙ pc² using CO(1 − 0) data published by Rosolowsky et al. (2007). This value is roughly a factor of 2 higher than our results. Corbelli et al. (2017) similarly find 44 ± 15 M_⊙ pc². Although they identify smaller GMCs with the CO data at 12″ angular resolution, they still find similar masses, resulting in higher surface mass densities. The fact that they find surface mass densities about twice as high as our data are likely attributed to their application of a X_CO value twice that of the Galactic standard value. However, this has been disputed in Paper I, which finds an average value nearly identical to the Galactic one. Gratier et al. (2012) find similar values for these properties for the same reasons. The data of Williams et al. (2019) exhibit the lowest mean surface mass densities of all with 7.5 ± 2.5 M_⊙ pc² which is probably due to the large sizes of the GMCs. Roman-Duval et al. (2010) find for their Milky Way data a median surface mass density of 144 ± 20 M_⊙ pc². Although the mass-size relations indicate a comparable surface mass density, there is an observed dissimilarity in the distribution shapes, with mean values varying by a factor of approximately one order of magnitude between M33 and the Milky Way. It should be noted that this finding aligns with the GMC masses we find in M33 being approximately an order of magnitude lower compared to the Milky Way GMCs and with the total masses of the two galaxies found by other studies mentioned above. We note that in the simulations of Dobbs (2008) the GMCs are more massive in galaxies with stronger spiral shocks or higher surface densities.

Increased SF activity and higher pressures correlate with increased molecular gas surface mass densities (Heyer et al. 2004; Lehnert et al. 2015; Wang et al. 2017; Krumholz et al. 2018). The difference between our results and the subset in NL16 is most likely due to manual selection of GMCs, involving a threshold applied to their masses. However, considering the Rice et al. (2016) Milky Way catalog reveals a similar range of especially high surface mass densities between both galaxies. Given the smaller mean sizes and higher masses of this catalog compared to our results, both mass and size lead to increased surface mass densities by a factor of ∼2. However, Corbelli et al. (2017) find a distribution similar to the GMCs in the Rice et al. (2016) catalog. We attribute this to the higher spatial resolution of the observations by Corbelli et al. (2017), which yield smaller GMC sizes relative to ours, although they still report mean masses comparable to ours probably due to the use of an X_CO value twice the Galactic standard value.

5.5. Power-law mass spectra

We aim to determine the mass spectrum of GMCs in M33 identified via dust and CO, which may relate to cluster and star mass functions (Kennicutt & Evans 2012, and references therein). Differences in the mass spectrum across regions might reflect variations in the processes that govern the formation, evolution and destruction of clouds (Rosolowsky 2005; Colombo et al. 2014). The mass spectrum typically conforms to a power-law probability distribution. To determine the power-law exponent α and its standard error ($\sigma /\sqrt{N}$) we first linearize the function

$$\begin{array}{c}\hfill p(x)\propto x^{-\alpha },\end{array}$$(9)

and then employ a least-squares approach to fit a linear slope to the data.

Figure 6 shows the distributions with an index determined to be α = 2.32 ± 0.10 for the dust-derived and α = 1.87 ± 0.08 for the CO-derived mass spectrum. The steeper slope of the dust-derived data indicates a larger number of less massive structures. This result is somewhat unexpected, because the low angular resolution of our data would have moreover suggested that many smaller molecular clouds along the line of sight would be artificially grouped into larger complexes, resulting in a flatter slope. Dobbs et al. (2019) reported α = 1.59, using the CPROPS algorithm to identify clouds in M33 from the IRAM CO data. For the data of Corbelli et al. (2017), using the same data and extraction method, a curved pattern with an index of 1.6 has been identified. For the results reported in Gratier et al. (2012), using IRAM CO data as well as CPROPS, we determine a single power-law of α = 2.12 ± 0.08. We note that in the simulations of Dobbs et al. (2019), the spread in α is large, between 1.66 and 2.27 (with an uncertainty of 0.2) and depends on the simulation (SPHNG or GASOLINE2) and the identification algorithm (FoF or CPROPS). Williams et al. (2019) find a higher slope of 2.83. Their result suggests a poorer ability of M33 to form massive clouds. Rosolowsky (2005) report a similarly steep mass spectrum slope of 2.9 ± 0.4, which may be biased by only sampling the high-mass end of the mass spectrum where the slope tends to be steeper.

Fig. 6. Power-law mass spectra of detected GMCs from dust-derived data (red) and CO (blue) as well as from NL16 (Milky Way, black). In black, the fit is performed with all listed structures in NL16, while the light gray dotted-dashed line shows the fit, where the CMZ of the Milky Way is excluded.

For the power-law index of Milky Way clouds, including the CMZ, we determine α = 2.35 ± 0.24 using the data presented in NL16. Excluding the CMZ results in an exponent of α = 2.58 ± 0.28. We note that this comparison relies solely on the 44 structures manually selected by NL16 for structures more massive than 0.7 × 10⁶ M_⊙, which could introduce a potential bias and an undetected systematic error. For the results reported in Roman-Duval et al. (2010), we derive α = 1.61 ± 0.03 (not shown in the figure) for the Milky Way CO data from the Galactic Ring Survey. These findings are in alignment with Rice et al. (2016), who found a slope of 1.6 ± 0.1 for their entire catalog. For the outer Galaxy, they reported a higher slope of 2.2 ± 0.1, whereas for the inner Galaxy, the slope remained at 1.6 ± 0.1. Similarly, Fujita et al. (2023) found generally higher slopes, yet they show a consistent pattern with an index of α = 2.30 ± 0.11 derived from ¹²CO data for distances below 8.15 kpc and α = 2.51 ± 0.14 for distances less than 16.3 kpc. The power-law indices found in several other studies of the Milky Way, all using CO data, typically range between 1.6 and 2 (Kramer et al. 1998; Simon et al. 2001; Schneider & Brooks 2004; Roman-Duval et al. 2010).

The efficiency of cloud formation has been associated with various processes. As discussed in Williams et al. (2019), the influence on the GMC population of the spiral density wave amplitude (e.g., Shu et al. 1972) can be excluded to explain the tentatively higher slopes in M33 due to modeling efforts, which indicate that the spiral arms of M33 are likely due to gravitational instabilities (Dobbs et al. 2018). The interstellar gas pressure might also be influential (Elmegreen 1996; Blitz & Rosolowsky 2006). Kasparova & Zasov (2008) report increased interstellar pressure compared to the Milky Way, potentially leading to the formation of more massive clouds. Thus, interstellar pressure may not be the primary cause of a potential inefficient cloud formation. This contrasts with findings by Blitz & Rosolowsky (2006) and Sun et al. (2018), indicating M33 lies within a lower pressure regime. This scenario aligns with the upper cloud mass limit being influenced by interstellar pressure. Another factor could be the role of metallicity in the transformation of H I-to-H₂ (Krumholz et al. 2008; Kobayashi et al. 2023). If M33 indeed has subsolar metallicity, this conversion would be less efficient, resulting in similarly inefficient cloud formation. However, the determined metallicity of M33 shows a very high dispersion (Willner & Nelson-Patel 2002; Crockett et al. 2006; Rosolowsky & Simon 2008; Magrini et al. 2010). Furthermore, it is proposed that merging H I clouds could form H₂ (e.g., Heitsch et al. 2005), suggesting that larger H I velocity dispersions could lead to more massive clouds. In M33, however, the average H I velocity dispersion is around 13 km s⁻¹ with minimal radial variation (Corbelli et al. 2018). This is consistent with the velocity dispersion of 11 nearby galaxies of ∼10 km s⁻¹ presented in Tamburro et al. (2009). Typical velocity dispersions measured for the Milky Way are in the same range (Malhotra 1995; Marasco et al. 2017). Another potential mechanism remains within supernovae. The power-law index may be considered to represent the balance between GMC mass-growth and destruction by massive stars (Kobayashi et al. 2017, 2018). The supernova frequency per unit volume varies across the galactic disk and the expansion of supernovae remnants compresses the ISM initiating the transition of H I-to-H₂ (Kobayashi et al. 2020, 2022). In this case, the power-law slopes of the GMC mass functions are determined by the balance between the transition rate from H I-to-H₂ and the destruction rate by stellar feedback from massive stars, mainly radiative feedback (Kobayashi et al. 2017). Additional mechanisms like shear may also set the maximum mass and lifetimes (Jeffreson & Kruijssen 2018), especially in a region where the shear rate is high and the orbital speed is fast (e.g., the outer regions of the CMZ in case of the Milky Way galaxy). We cannot determine which of these mechanisms primarily drive the potentially inefficient cloud formation in M33, as suggested by some of the findings discussed above.

In summary, the power-law index α shows a large spread for both M33 (1.6 to 2.9) and the Milky Way (1.6 to 2.5), due to differences in datasets and methods. Despite errors, there is no significant variance between M33 and the Milky Way, except for a slight tendency for higher values in M33. Both exhibit self-similarity from molecular clouds (∼50 pc) to larger GMAs, suggesting similar physical mechanisms for massive GMCs in both galaxies and a limit in sizes and masses despite their high difference in mass. Given the values in the existing literature, it is difficult to determine whether the cloud mass distribution in M33 is significantly different from that in other large spirals within our Local Group.

6. Trends with galactocentric radius and galactic environment in M33

Molecular clouds do not possess a perfectly spherical shape. Instead, their morphology is often influenced by complex processes such as merging or turbulent flows (e.g., Vazquez-Semadeni et al. 1995; Heitsch et al. 2006; Clark et al. 2019; Schneider et al. 2023) or cloud-cloud collisions (Casoli & Combes 1982; Fukui et al. 2021), leading to irregular shapes characterized by clumps and filaments. Variations in cloud properties under different environmental conditions within a galaxy offer valuable insight into the factors shaping cloud formation and evolution (e.g., Sun et al. 2020b).

The molecular gas, for example, forms huge associations as a result of the gravitational attraction of the spiral arm. As the gas exits the spiral arms and experiences significant shear forces, it breaks apart and reverts to smaller elongated structures (La Vigne et al. 2006). Numerous observational and computational studies emphasize the presence of filamentary structures in the areas between the arms (Ragan et al. 2014; Duarte-Cabral & Dobbs 2016, 2017) and the presence of high-mass structures within the spiral arms (Dobbs et al. 2011; Miyamoto et al. 2014). Apart from structure variations, metallicity gradients within a galaxy can also lead to variations in the physical properties of the molecular cloud. We thus examine in the following sections the physical properties of the GMCs in M33 as a function of the galactocentric radius and the galactic environment of M33.

6.1. Trends with galactocentric radius

Figures 7 and 8 display the mass, average density, surface mass density, radius as well as aspect ratio and dust temperature as a function of the galactocentric radius. The relationship between GMC properties and galactocentric radius has also been examined by Gratier et al. (2012), Corbelli et al. (2017) and Braine et al. (2018). A comprehensive discussion of specific properties can be found in Appendix E.

Fig. 7. Mass, average density, surface mass density, and radius of the identified GMCs along the galactocentric radius. The dust-derived GMCs are shown in red, whereas the CO-derived GMCs are represented in blue. The bigger dark red and blue data points mark the GMCs that have a surface mass density exceeding 40 M⊙ pc2.

Fig. 8. Aspect ratio and temperature from dust and CO-derived data along the galactocentric radius. The bigger dark red and blue data points mark the GMCs that have a surface mass density exceeding 40 M⊙ pc2.

In summary, the parameters show only a weak (for high-Σ GMCs) or non-existing (for low-Σ GMCs) trend with the distance from the galaxy’s center, raising the question whether SF is influenced by the galactocentric radius. Only GMCs with the highest surface mass densities (above 40 M_⊙ pc²) show a tendency to have higher values for density and Σ in the center of M33. This finding is similar to what is observed in the Milky Way. In both galaxies, self-gravity and cloud-cloud collisions become more important for these high-Σ GMCs in the respective CMZ.

In the following section, we discuss the more significant trends we observe for different regions (center, spiral arms, outskirts) in M33, as a radial dependence on the galactocentric radius does not entirely unveil systematic differences in the galactic environments.

6.2. Trends with galactic environment

It is not yet clear whether SF is more efficient in particular regions of galaxies and to which extent the SFR and SFE are linked to the physical properties of the GMC population. Observations and simulations indicate that GMCs are concentrated in spiral arms, often with regular spacing, which can be explained when GMCs are formed by gravitational instabilities (Elmegreen 1990; Kim & Ostriker 2002). On the other hand, GMCs can also form by agglomeration of smaller clouds or merging of flows (see references above). A higher SFR can then be just a by-product of the higher material reservoir in the spiral arms. While some studies (Koda et al. 2009; Pettitt et al. 2020; Colombo et al. 2022) report variations between their spiral arms and inter-arm populations, others (Duarte-Cabral & Dobbs 2016; Querejeta et al. 2021) find no discernible differences in the overall properties of the cloud population.

In this section, we systematically investigate if there are variations in the physical properties of the GMC population in certain regions of M33. For that, we use our dust-derived column density map and split the galaxy by eye-view into a central region, the two main spiral arms and the outskirts (Fig. F.1). The two main spiral arms are approximated to extend to a galactocentric radius of roughly 4 kpc, whereas the central area of M33 can roughly be described as an equivalent circle with a galactocentric radius of around 1.3 kpc. The outskirts are considered to be the remaining area of M33’s disk⁴. To determine the spiral arm structure more quantitatively, we additionally employ a similar approach as in Querejeta et al. (2021) and model the spiral arms with a log-spiral function and perform a fit to this model. Details of this procedure and the results are given in Appendix F and in Fig. F.1. The visually estimated borders of the spiral arms already capture the fitted log-spirals very well. We therefore continue to use the masks presented in Fig. F.1 to study the spiral arms and outskirts.

6.2.1. Column density complementary cumulative distributions

Querejeta et al. (2021) reported increased gas surface densities closer to the central regions of galaxies by analyzing the CO(2 − 1) data obtained from the PHANGS-ALMA survey (Leroy et al. 2021). We confirm this finding for M33 using our dust-derived high-resolution N_H2 map (Fig. 9), which shows the complementary cumulative distributions of the entire disk of M33 and the three defined environments⁵. The complementary cumulative distribution function provides the likelihood that an observation from a sample exceeds a certain value on the x-axis. It becomes evident that the central region exhibits column densities throughout the spectrum higher than those of the spiral arms and outer regions. The spiral arms and the outskirts display comparable levels of low N_H2 below approximately 2 × 10²⁰ cm⁻². Beyond this threshold, the spiral arms diverge, maintaining higher column density values. This finding aligns with the results reported in Leroy et al. (2021). The median value for the central region (provided in the panel for all distributions) is roughly three to 3.5 times higher than for the other two regions. Furthermore, the central region shows the steepest slope among all distributions, while the spiral arms and outer regions demonstrate a shallower slope towards higher column densities.

Fig. 9. Complementary cumulative NH2 distributions of the entire galaxy (based on dust-derived data) and the three galactic environments. We note that these distributions are solely based on pixels and are not connected to GMCs.

6.2.2. GMC properties in different environments

The distributions of GMC properties (mass, average density, surface mass density, radius and aspect ratio) are shown in Fig. 10 as a function of galactic environment. The global distribution of the entire disk of M33 is shown in red on the left for comparison. GMCs located in the center are represented in violet, those in the two main spiral arms are in brown and those in the outskirts, excluding the center and the two main spiral arms, are depicted in turquoise. The median is displayed as a straight line within the boxes in beige.

Fig. 10. Box plots of the determined dust-derived parameters categorized based on galactic environments. The lower and upper whiskers of the box plot represent the lowest and maximum values of the dataset, respectively. The colored box shows the distribution’s interquartile spread, or the range from the 25th to the 75th percentile; the median is indicated by the solid beige line inside the box. The distributions’ outliers are shown as circles.

Most of the properties show a weak variation for the median values in different environments. Only the central region of M33 exhibits larger masses and surface mass densities of the GMCs compared to the regions in the remaining disk (see also Sect. 6.1, where we have already observed this trend). Overall, the GMCs in the center are denser, those in the spiral arms are larger, while those in the outskirts are more elongated. Generally, the GMC populations in the spiral arm and outer regions do not exhibit large variations in their properties.

6.2.3. GMC masses in different environments

The masses of the GMCs are noticeably higher in the central region of M33. The median and minimum values indicate significantly higher masses compared to the other two regions. Apart from the exceptional case of NGC 604 in the spiral arm, the highest mass values are comparable to those in the central area, while the GMCs with the lowest masses have even lower values. The outer regions exhibit GMCs with similarly low masses as those in the spiral arms but lack GMCs with such high masses.

One hypothesis is that spiral arms, which contain a larger amount of material, increase the occurrence of cloud-cloud collisions, thereby supporting the creation of high-mass entities (Dobbs 2008). This would result in a tendency for the most massive clouds to be situated in spiral arms. However, the spiral arms exhibit lower densities. This is also true for the surface mass density compared to that in the central region. If larger GMCs gather more mass and thus support SF, then this should yield higher surface mass densities. Since the GMCs in spiral arms are merely larger without possessing higher column densities, this results in lower masses and surface mass densities, which correlate with SF, suggesting that SF should be lower. As discussed above, the impact of cloud-cloud collisions in the Milky Way have been investigated by Kobayashi et al. (2017) and Kobayashi et al. (2018), for which an effective impact has only been found for GMCs more massive than 10⁶ M_⊙. Furthermore, while the most massive GMC (NGC 604) is located in the northern spiral arm, the other GMCs in these environments do not support this picture. Both the median and the 75th percentile values are lower than those of the center. Additionally, most outliers, except for NGC 604, have less mass compared to those in the center. This discrepancy may be due to the limited resolution of 75 pc, whereas Dobbs (2008) simulate molecular clouds with higher resolution. Corbelli et al. (2019) suggested that the formation of more massive clouds in the center may occur due to the rapid rotation of the disk relative to the spiral arm pattern, allowing the clouds to grow further as they traverse the arms.

6.2.4. GMC densities, surface mass densities, and radii in different environments

For GMC densities, the environments show minimal variation. The median values are similar across different regions. The main distinction is observed in the outliers at the outskirts, where the densities do not peak as high as those in the GMCs in other areas.

The surface mass densities demonstrate a pattern similar to that of the masses. The central region contains GMCs with the highest masses, whereas the median and the values at the lower ends of the spectrum decrease in the spiral arms and decrease even more in the outer regions. This finding aligns with the radial trends that have been discussed in Sect. 6.1 and confirms Querejeta et al. (2021), reporting increased gas surface densities closer to the central regions of galaxies by analyzing the CO(2–1) data obtained from the PHANGS-ALMA survey (Leroy et al. 2021).

On the other hand, the GMC radii differ most significantly in the spiral arms, with NGC 604 as the outlier. The center and outskirts have smaller GMC radii. The median radius in the spiral arms is ∼90 pc, while in the outskirts it is ∼63 pc.

6.2.5. GMC elongations/aspect ratios and temperatures in different environments

The center has GMCs with the least elongation at the higher end of the spectrum, while the spiral arm and outskirts have slightly stronger elongated GMCs. It is not clear to which extent the GMCs with an AR larger than 3 represent GMFs that were found in the Milky Way (Ragan et al. 2014; Goodman et al. 2014; Zucker et al. 2015; Wang et al. 2020) and in external galaxies (Hughes et al. 2013; Leroy et al. 2016). These studies typically define GMFs as long filamentary structures with lengths exceeding 50 pc and masses above 10⁵ M_⊙ and suggest that they trace the denser spine region of the spiral arms and the mid-plane of the gravitational potential in the galaxy. We note that while some Galactic GMFs exhibit widths down to ∼1 pc, a scale which remains indistinguishable from our current resolution, other GMFs possess notably larger widths (Zucker et al. 2018). In particular, Wang et al. (2020) presented dust and dense gas tracers of one filament in the Milky Way with an AR of about 3 and a length of 68 pc which would fit formally to some of the GMCs we detect. Figure 10 shows that there is only a very weak environmental dependency of the AR. However, the most elongated GMCs are found within the spiral arms and the outskirts and this could indicate (as discussed in Sect. E.5) a stretching effect due to shear forces on the massive GMCs (or molecular clouds) as they transition from the spiral arms to the inter-arm regions (Koda et al. 2009). It could also be the result of disruption caused by feedback from stars (Meidt et al. 2015; Chevance et al. 2020; Bonne et al. 2023). Duarte-Cabral & Dobbs (2016) discovered in a computational simulation of GMCs within a two-armed spiral galaxy that, while the average characteristics of the inter-arm and spiral arm GMCs are comparable in terms of their ARs, the extremely elongated GMCs in their dataset are predominantly associated with the inter-arm regions. Given that the outskirts exhibit highly elongated GMCs it is possible that stellar feedback contributes to the disruption of GMCs in both environments.

Hence, our proposition is that the shear forces in the center do not account for disrupting the GMCs, since the GMCs in the center exhibit the lowest median elongation. Additionally, the least elongated and most massive GMCs are located in the center (excluding NGC 604), indicating that the center is conducive to the formation of high-mass GMCs. As detailed in Sect. E.6, we argue that the strong galactic potential subjects GMCs to an isotropic pressure, which accounts for the observed distribution of GMC elongation. It is unlikely that stellar feedback plays a significant role in elongating GMCs in the center, as the feedback would be uniformly distributed throughout the center, resulting in disrupted GMCs across the region and dynamically altering GMC boundaries. After typical GMC lifetimes, any cloud detection algorithm would identify new segments of an original GMC as a new GMC, incorporating parts of previously disrupted GMCs. With this iterative process and the isotropic galactic potential, extreme elongation tendencies are expected to diminish, resulting in the non-increasing elongation of GMCs, unless specific conditions, such as the presence of a bar, exist in the central area.

In summary, it is observed that while the overall dynamics on a large scale influences cloud properties, there is no clear indication that SFE is notably enhanced in any specific environment.

6.3. Power-law mass spectra with galactic environment

As pointed out in Sect. 6.2.3, there are noticeable differences in the masses of the GMCs between the center and the remaining disk of M33. We therefore also investigate the distribution of GMC masses within each large-scale environment by conducting power-law mass spectra.

Figure 11 illustrates the mass spectra of the GMCs in the different environments of M33. The central area and the spiral arms of M33 exhibit the highest abundance of high-mass clouds, with α_Center = 2.19 ± 0.19 and α_SpiralArms = 2.16 ± 0.14 having relatively shallow slopes. In contrast, there is a significant reduction in high-mass GMCs toward the outer regions, where GMCs have mainly lower to moderate masses with a steeper slope of α_SpiralArms = 2.32 ± 0.11. Bigiel et al. (2010) observed smaller GMCs at larger galactocentric radii of M33, suggesting a steeper slope in the outskirts of M33, which is supported by the slope we determine. This higher slope suggests that high-mass objects in the outskirts may face challenges in their formation or are rapidly destroyed after formation. This aligns with the findings presented in Sect. 6.2.3 and with the results of Rosolowsky et al. (2021) in their examination of GMCs across spiral galaxies within the PHANGS dataset.

Fig. 11. Power-law mass spectra of GMCs in the dust-derived map for the three galactic environments. The power-law mass spectra is split into the three galactic environments center (pink), spiral arms (brown) and outskirts (turquoise). For comparison, the fit for the whole disk is shown again in red.

The spectra shown in Fig. 11 indicate that GMCs have higher masses in areas with lower galactocentric radii, closer to the center. In the outskirts, the most massive GMC reaches a mass of ∼1 × 10⁶ M_⊙, the lowest of the three environments. This is consistent with having the steepest slope of all three. In the central region, GMCs that are only about three times more massive are found, while in the spiral arm, the most massive GMC (NGC 604) has a mass of ∼8 × 10⁶ M_⊙. This observation indicates that cloud growth may be prevented or that large GMCs are being disrupted in the central area, at least to reach such high masses as observed in NGC 604. This phenomenon could be attributed to complex dynamics and shear forces or to the enhanced interstellar radiation field in the central region. In contrast, the mass distribution of GMCs in the spiral arms, excluding NGC 604, consists of less massive clouds than in the center. NGC 604 leads to a flattening of the slope in the spiral arms. Despite the predominance of low-mass objects in the spiral arms, the conditions in this region appear to be conducive to the growth of larger clouds, maybe due to the absence of a strong interstellar radiation field and/or shear forces disrupting the clouds. However, since the center hosts the GMCs with the highest masses, with the exception of NGC 604, this conclusion remains uncertain.

Dobbs et al. (2019) found a decrease in the power-law index after incorporating SF into their simulations. As clouds become dense, the index drops to values between α ≈ 1.8 and ≈2. Considering delayed SF results in an index that agrees better with observations. This suggests, in general, that SF occurs in later stages of GMC formation. These authors also divided the clouds into “star-forming” and “non-star-forming” clouds (SF clouds and non-SF clouds hereafter). While SF clouds inject energy into the clouds, heating them locally, this leads to a flattened slope of α = 1.8, whereas non-SF clouds exhibit a slope of α = 2.68. They identified that the non-SF clouds tend to reside at a larger galactocentric radius, indicating higher SF activity in the center. Increased surface mass densities are associated with this phenomenon and we also observe higher surface mass densities toward the center, in accordance with this. Additionally, a stronger galactic potential towards the center could account for this finding. Compared with our results, this suggests a higher SF activity in the center and spiral arms than in the outskirts. Dobbs et al. (2019) also provide reasons as to why larger clouds tend to host more SF: these larger clouds are statistically more inclined to have dense areas, thus increasing the likelihood of SF. These clouds probably accumulate more mass as they begin to form stars, suggesting that clouds not undergoing SF may just be in an earlier phase of their lifetime when they have lower masses.

This is consistent with what Braine et al. (2018) found by analyzing the IRAM 30 m CO data of M33. Detected GMCs (from the catalog of Corbelli et al. 2017) have been divided into three radial bins: R_gal < 2.2 kpc, 2.2 < R_gal < 3.7 kpc and R_gal > 3.7 kpc. The power-law indices for these bins are α = 1.36, α = 1.68 and α = 1.87, respectively, showing an increase with radius. They also subdivided the GMCs into three SF classes – no obvious SF (A), embedded SF (B) and exposed SF (C) – based on Corbelli et al. (2017). More evolved GMCs accumulate more mass and show shallower slopes. Star-forming GMCs lie closer to the center than non-SF GMCs. However, A-class GMCs consistently show steep slopes regardless of their position, while C-type GMCs also have similar slopes regardless of their position, indicating that SF activity is more important than galactic environment.

As also discussed in Sect. 5.5, Fujita et al. (2023) found distinct power-law indices of α = 2.30 ± 0.11 and α = 2.51 ± 0.14 in the Milky Way within a galactocentric radius of < 8.15 kpc and beyond < 16.3 kpc, respectively. Taking the errors into account, this is also consistent with our data split into the environments, which represent distinct regions along the galactocentric radius.

7. Conclusion and summary

In Paper I, we presented a novel technique to use the Herschel flux maps and CO(2 − 1) data of M33 to produce N_H2 maps at 18.2″ (∼ 75 pc) resolution, resolving GMCs. A complete X_CO map was applied to the CO map to compute the N_H2 map with values in the range of 1.6 − 2 × 10²⁰ cm⁻²/(K km s⁻¹). This X_CO factor is close to the Milky Way value and thus questions the usual approach of applying a single, adopted X_CO factor for the whole galaxy and simply using a two-times-higher value for M33 due to its lower metallicity.

We then employ the Dendrograms algorithm to identify GMCs from these maps, calculate the physical properties, and compare the results between dust and CO and with Milky Way data from the CO(1 − 0) Columbia survey presented in Nguyen-Luong et al. (2016). In addition, an investigation was conducted to explore the potential influences of the galactocentric radius and galactic environment on GMC properties.

1. We find that M33 lacks the more massive (> 10⁶ M_⊙) and denser GMCs that are present in the Milky Way. The mean GMC masses of M33 are about an order of magnitude lower than those of the Milky Way. A power-law fit to the mass spectrum gives values of α = 2.32 ± 0.10 for dust and α = 1.87 ± 0.08 for CO. These indices align with those found in other studies of M33, Milky Way values, and simulations, which all show a large spread.

2. There appears to be a limit to the sizes of GMCs of around 150 pc, as the distributions for the largest GMCs of M33 and Milky Way show similar shapes and a decline above ∼100 pc. We do not find the equivalent of GMFs in the Milky Way but note that there is an inter-cloud medium at column densities of around 10²¹ cm⁻² that contains a significant mass, in particular in the central region of M33. In the outskirts, the lower-column-density material encloses the GMCs, which is particularly evident in dust.

3. The surface mass densities for M33 are 22 ± 5 M_⊙ pc² from dust and 16 ± 6 M_⊙ pc² for CO, which are about an order of magnitude lower than the same values for the Milky Way. The increased surface mass density may suggest an increase in SFR. Finally, M33 shows similar patterns in some alternative characteristics to those observed in other nearby galaxies in the PHANGS survey.

4. We find no or only weak correlations between physical properties and galactocentric radius, but some results indicate a dependence on the larger-scale environment.

5. The central region of M33 displays slightly higher median values for parameters such as mass, average density, surface mass density, and dust temperature, but contains the GMCs with the smallest aspect ratios. The center hosts the most massive GMCs (except for NGC 604), which also exhibit the highest surface mass densities. However, as the center seems to be the region with the highest influence on star formation, the variations in physical parameters across the environments are predominantly minor in nature. The spiral arms mainly host the largest GMCs, while they contain most of the extreme outliers across different parameters, such as mass, surface mass density, size, and elongation. On the contrary, the outskirts generally feature the lowest median values, with the exception of average density and elongation. However, the majority of the GMCs, despite some outliers, do not seem to be significantly affected by the conditions of the galactic environment.

6. The power-law fits to the mass spectra derived from CO and dust vary with the galactic environment (α = 2.19 ± 0.19 for the center, α = 2.16 ± 0.14 for the spiral arms, and α = 2.32 ± 0.11 for the outskirts). These results are consistent with observations in the Milky Way, suggesting similar indices for both the inner and outer disk of our galaxy. However, the slope of the spiral arms decreases due to the high mass of GMC NGC 604. The remaining high-mass GMCs in the spiral arms have lower masses than those in the galaxy center. This complicates the identification of the physical mechanisms at work, as high interstellar radiation fields and shear forces are likely to disrupt more massive GMCs, whereas the absence of these mechanisms would enhance these parameters in the spiral arms.

Overall, we conclude that the center seems to have a slightly greater influence on GMC properties than the other environments, but that mechanisms operating at the cloud scale – notably stellar feedback – may have a similar or greater impact on GMCs than large-scale dynamics inherent to galactic environments.

Data availability

Full Tables 1 and 2 and the updated H₂ column density data are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/692/A226

Acknowledgments

E.K. and N.S. acknowledge support from the FEEDBACK-plus project of the BMWI via DLR, Project Number 50OR2217. E.K. acknowledges support by the BMWI via DLR, project number 50OK2101. S.K. acknowledges support from the Orion-Legacy project that is supported by the BMWI via DLR, project number 50OR2311. We gratefully acknowledge the Collaborative Research Center 1601 (SFB 1601 subprojects A6, B2 and B3) funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 500700252. This research made use of astrodendro, a Python package to compute Dendrograms of Astronomical data (http://www.dendrograms.org/)

References

Appendix A: IRAM ¹²CO(2 − 1) line-integrated intensity map of M33

Figure A.1 shows the ¹²CO(2 − 1) map of M33 obtained with the IRAM 30m telescope (Gratier et al. 2010; Druard et al. 2014).

Fig. A.1. 12CO(2 − 1) line-integrated intensity map of M33 (Druard et al. 2014). The map has been smoothed to the resolution of 18.2″ and re-gridded to the coordinate grid of the SPIRE 250 μm map.

Appendix B: X_CO factor map of M33

Figure B.1 displays the X_CO factor map defined as the dust-derived N_H2 over CO line-integrated intensity at each position in M33 at 18.2″ and scaled with the $\mathrm{CO}(\frac{2-1}{1-0})$ line ratio (Druard et al. 2014) to CO(1 − 0) intensity. See Paper I for more details.

Fig. B.1. XCO factor (ratio) map of Method I of Paper I. The two ellipses represent a circular radius equivalent to 2 and 4 kpc.

Appendix C: Influence of Dendrogram parameters on the GMC statistics

We conducted a Dendrograms parameter study by changing the min_value and the beam factor for cloud selection. Obviously, increasing these parameters results in the identification of larger and more massive GMCs, whereas changing min_delta practically does not alter the results. Nevertheless, the undetected residual emission became more substantial, while the results change non-significantly considering the uncertainties. Consequently, we conclude that the use of 3σ for min_value and a beam factor of 1.2 are the optimal settings for the Dendrogram analysis. Tables C.1 and C.2 list the mean values of the main properties of the dust- and CO-derived GMCs for a min_value of 5σ and a beam factor of 1.5. We further illustrate the results for a subset of the plots discussed above. Figure C.1 shows the distributions of the radii and surface mass densities for both varied parameters. The size is mainly unchanged, while a min_value of 5σ excludes the low-mass GMCs, leading to an overall shift towards higher values in mass and (surface mass) densities. However, the change is still low.

Fig. C.1. Distributions of GMC properties with varying Dendrogram parameters. The upper panels show histograms of radius and surface mass density with a min_value of 5σ. The lower panels display the same properties for a beam factor of 1.5.

The effect on the galactocentric radius dependence as an example for the mass and density is marginally shifted to higher values. Furthermore, there is no notable trend with the galactocentric radius for any of the other properties, similar to the result found with a min_value of 3σ and a beam factor of 1.2.

For the power-law slopes, the increase in the slope for dust-derived data is insignificant. For CO, the slope rises from 1.87 to 2.03. However, considering the uncertainties, the change falls within the margin of error.

This analysis confirms that our selected Dendrogram parameters are robust and produce reliable results.

Appendix D: CO luminosity

For completeness and to enable comparison to other studies, we display in Fig. D.1 the ¹²CO(1 − 0) luminosity L_CO of M33, using the CO(2 − 1)/CO(1 − 0) line ratio of 0.8 from Druard et al. (2014). The mean CO luminosity is (9.2 ± 2.0)×10⁴ K km s⁻¹ pc², which is lower than the value of 3 × 10⁵ K km s⁻¹ pc² for the CO(2 − 1) IRAM data at 12″, which corresponds to 3.75 × 10⁵ K km s⁻¹ pc² for CO(1 − 0) applying the same line ratio of 0.8. Hughes et al. (2013) present M33 data from their observations of CO(1 − 0), indicating values around 1 × 10⁵ K km s⁻¹ pc². A slope of 1.5 indicates that the CO emission of smaller GMCs is faster saturated compared to larger GMCs. In other words, more CO emission per area can be accounted for in the outer regions of GMCs.

Fig. D.1. CO(1 − 0) luminosity-size relation of the structures identified in the CO data. The horizontal line represents the 2σ sensitivity limit.

Appendix E: Trends with galactocentric radius of M33

It is evident across all parameters that the majority of GMCs do not exhibit a significant trend with the galactocentric radius. This observation is supported by Spearman correlation coefficients⁶ ranging from |8.7 × 10⁻⁴| to |0.4|, showing no or only low correlation at best (except for the temperature, which clearly shows a gradient; see Tabatabaei et al. 2014; Keilmann et al. 2024). Nonetheless, upon closer inspection of the GMCs at the extreme ends of the spectra, there is a subtle trend of mostly decreasing values with increasing galactocentric radius, especially for the surface mass density and average number density. Extreme GMCs, that is, GMCs at the tails of the distributions, enhanced in particular large-scale galactic environments, may suggest the presence of physical mechanisms directly enabling the development of particular cloud types in specific galactic regions, potentially influencing SF. We consider GMCs with surface mass densities above 40 M_⊙ pc² to account for an analysis of whether and how extreme clouds may depend on the galactocentric radius. The corresponding GMCs are depicted as thicker and darker pentagons in the following figures.

E.1. GMC masses with galactocentric radius

Specifically, for dust-traced masses, GMCs located within approximately 2 kpc demonstrate an increase toward the center for those with the lowest masses. This pattern is absent in CO-traced GMCs. This could be due to strong interstellar radiation fields, which photo-dissociate CO (Offner et al. 2014), leading to less CO emission in the center. The remaining GMCs beyond 2 kpc do not show a dependence on the galactocentric radius for both tracers. The Spearman correlation coefficients show practically no correlation with galactocentric radius over the whole data range. For the most massive GMCs observed with both tracers, their highest masses decrease with galactocentric radius, supporting Corbelli et al. (2017), who found a similar decrease beyond 4.5 kpc. The GMCs with the most extreme surface mass density values are not those with the highest masses; they appear to be arbitrarily distributed in terms of mass; see Fig. 7 (top left). The data point around 4 kpc is NGC 604. In the case of branches, the median mass is 2 × 10⁵ M_⊙. Additionally, there are structures with increased masses within low galactocentric radii below 1 kpc. The remaining structures show a similar trend to the GMCs (leaves) with a slightly higher correlation coefficient of −0.5, which is still only a moderate correlation.

According to their simulation, Dobbs et al. (2019) report that the masses of GMCs are influenced by the distance from the center of M33. Their findings exhibit a resemblance to our results, particularly in the case of the most massive GMCs, showing a correlation with the galactocentric radius. Nevertheless, Dobbs et al. (2019) did not provide a quantitative assessment of this dependency, making direct comparisons difficult. Due to the similarity in the plots showing that the majority of the clouds do not seem to depend strongly on the galactocentric radius, it is possible that the actual dependency they state is similar in magnitude to what we quantify.

The lack of correlation between cloud mass and R_gal likely results from a balance between the more compact and therefore more luminous clouds at the galaxy center and the diffuse and hence more extended sources in the outskirts. Williams et al. (2019) identified a higher, yet still weak correlation between mass and galactocentric radius with a Kendall rank correlation coefficient⁷ of 0.12.

E.2. GMC average densities with galactocentric radius

While, as in the case of masses, a subtle trend of the least dense GMCs is also noticeable in dust-traced GMCs in terms of average density, it is not observed in those traced by CO. In contrast, GMCs with the highest densities tend to be more concentrated toward the inner disk of M33, as indicated by dust-derived GMCs and to a lesser extent by the CO-derived GMCs. Almost all of these most dense GMCs are also those which have the highest surface mass densities. This is not surprising, as the number and the surface mass density are closely related. However, the overall correlation is absent, as in the case of the masses. The Spearman correlation coefficients practically do not quantify the correlation. Structures identified as branches show a median density of 1.1 cm⁻³ and exhibit a slightly higher correlation coefficient of 0.3, which is still only a weak correlation.

E.3. GMC surface mass densities with galactocentric radius

Since the masses are divided by the area of a GMC, the surface mass density is somewhat less dependent on the resolution, making it comparably easier to compare with other studies.

The surface mass density (Fig. 7, bottom left) exhibits a similar pattern as for the mass among dust-traced GMCs considering the GMCs with the lowest surface mass densities found within 2 kpc. For both tracers, GMCs located at the higher end of the spectrum demonstrate a slight tendency to exhibit a higher surface mass density as the galactocentric radius decreases.

The only strong connections to the other parameters are the radius (see Sect. E.4 and bottom right of Fig. 7) and the averaged number density (upper right of Fig. 7). The number density is related to the surface mass density in a natural way, which does not reveal new surprising insights.

The overall correlation for all data points in CO is practically absent, while for dust-derived GMCs a low correlation of −0.2 is determined. Thus, also the surface mass density seems to have no strong dependence on the galactocentric radius.

In the case of branches, the dependency is similar with a median value of 18 M_⊙ pc⁻² and a doubled but still only lower moderate correlation coefficient of −0.4.

E.4. GMC radii with galactocentric radius

Considering the radius, the data points suggest that there is no clear pattern with the galactocentric radius, which is further supported by the nearly nonexistent or very weak correlation of −0.3. However, there are also a small number of GMCs with the largest radii close to the center, which are not those with the highest surface mass densities. Instead, the GMCs with the highest surface mass density are preferably the smallest ones.

The branches show a median radii of 185 pc with very large structures within a galactocentric radius below 1 kpc. This mainly increases the correlation coefficient to a moderate value of −0.5.

E.5. GMC elongations with galactocentric radius

In terms of AR (Fig. 8, left), dust-derived GMCs that are most elongated tend to be situated farther away from the galaxy’s center, typically beyond approximately 3 kpc, whereas the majority of GMCs do not show a clear trend, showing a weak correlation coefficient of 0.2. CO-derived GMCs exhibit a similar distribution with little to no correlation, with the most elongated GMCs to a low extent found in the mid-range between ∼ 2 kpc and ∼ 5 kpc of the galaxy. However, the most extreme GMCs do not follow a consistent pattern. For both tracers, there seems to be no distinct pattern regarding the GMCs with the highest surface mass densities. The overall elongation of CO-derived GMCs is less compared to dust-derived GMCs with a median of 1.6, which aligns with the contours of the CO-derived structures in Fig. 2.

An almost unchanged elongation is found for the branches with a median value of 2.1 and a reduced correlation coefficient of 0.1.

E.6. GMC temperatures with galactocentric radius

The data points for the temperature show a dependency with galactocentric radius to some extent (Fig. 8, right). However, some data points beyond roughly 2.5 kpc deviate distinctly from the remaining data points. This leads to the most significant correlation with the galactocentric radius among all the parameters analyzed, which is −0.4 and −0.6 for dust- and CO-derived GMCs, respectively. This is a weak to moderate correlation. Interestingly, Williams et al. (2019) find only a weak correlation of −0.26. In the case of branches, a similar trend is found with a slightly higher median temperature of 21.5 K and a correlation coefficient of −0.7. This is the strongest correlation found in this study. A natural explanation for the decreasing dust temperature is the overall decrease in intensity of the interstellar radiation field with increasing galactocentric radius (Rice et al. 1990).

A compelling relationship emerges when higher pressures appear to result in less elongated GMCs (indicated by higher densities and temperatures in the central region). Such a phenomenon could possibly be attributed to the pervasive pressure within the central region (Sun et al. 2020a,b) as a result of the stronger galactic potential, which acts uniformly, resulting in more isotropically shaped GMCs. In contrast, GMCs located in the outer regions appear to be influenced by pressure originating predominantly from a specific direction, causing forces that are not uniformly distributed across all GMCs. Those at mid-range distances (presumably located in the spiral arms) could be elongated when exiting the spiral arm and therefore its gravitational potential. It is also possible that stellar feedback causes those GMCs to elongate. Nevertheless, this phenomenon appears to be relevant solely to the GMCs lying at the tail of the spectra, specifically those with the greatest elongation.

Appendix F: Galactic environments of M33

We constructed an unsharp-masked image of the N_H2 map and identified the densest points along the spiral arms on our N_H2 map. This is different to Querejeta et al. (2021), who use stellar densities and a morphological decomposition based on Spitzer 3.6 μm. However, since we intend to define the spiral arms for GMCs, which are located in H₂ gas, we use our N_H2 map for this decomposition. The coordinates of the densest points were deprojected to the plane of the galaxy (using an inclination of 56° and a position angle of 23°). The fit was then performed in logarithmic polar coordinates. The log-spiral fit was projected back to the plane of the sky and is shown in Fig. F.1. It matches the areas of the spiral arms very well by eye-inspection. From this result, the asymmetry of both spiral arms becomes obvious. The northern spiral arm is wound stronger with a higher pitch angle starting from the center compared to the southern spiral arm. Overall, we distinguish three galactic environments: center, spiral arms and outskirts onto the dust-derived N_H2 map. We have chosen a circle of 1.3 kpc for the center to distinguish from the spiral arms and outskirts. This radius was selected since it encompasses the maximum column density distribution of the central region, while still maintaining reasonable borders for the fitted spiral arms.

Fig. F.1. The dust-derived NH2 map shows the boundaries of three galactic environments: the center (dashed pink lines), spiral arms (dashed brown lines) and outskirts (remaining area). Center coordinates are RA(2000) = 1h33m50s, Dec(2000) = 30°39′37″ (SIMBAD 2024). The result of the log-spiral fit is shown in red.

All Tables

All Figures

Fig. 1. Total and molecular hydrogen column density maps. Left: High-resolution N(H) total gas column density map obtained from the Herschel flux maps of M33 with 18.2″ angular resolution (indicated by the circle in the lower left corner) using the β map from Tabatabaei et al. (2014). Values below and above a minimum and maximum threshold (1018 cm−2 and 1022 cm−2, respectively) are blanked. Right: High-resolution H2 column density map derived from the total N(H) map by subtracting the H I component. In both maps, the CO contour levels 3 and 6σ of the CO map (Fig. A.1 in Appendix A) are shown.

In the text

Fig. 2. Dendrogram source extraction from the NH2 and CO maps. Left: GMC detections outlined by the lowest level branches in the dust-derived NH2 map (green contours) and 326 leaf structures (red contours). NGC 604 is marked with a thick pink contour. Right: GMC detections identifying 199 leaf structures in the line-integrated CO(2 − 1) map (blue contours) and similar to the dust map, the branches (green contours). The small circle in the lower left corner of both figures shows the beam size of 18.2″.

In the text

	Fig. 3. Structures identified using Dendrograms on CO and dust-derived maps superimposed on the XCO map. The structures found in the dust-derived NH2 (red) and in the CO map (blue) are mapped onto the XCO factor (ratio) map, which represents the dust-derived H2 column density over CO intensity from Paper I.
In the text

Fig. 4. Distributions of GMC properties from our study and the literature. The panels show histograms of radius (top left), surface mass density (top right), mass (bottom left) and beam-averaged number density (bottom right) derived from H2 data from dust (solid red) and CO (solid blue) from this study. The dashed purple and golden lines display the distributions obtained for M33 from dust (Williams et al. 2019) and CO (Corbelli et al. 2017), respectively. The dashed black and light red lines give the distributions for the Milky Way studies of NL16 and Rice et al. (2016), respectively. The vertical black dashed line in the mass distribution signifies the minimum mass limit of the selected structures in the Milky Way by NL16.

In the text

	Fig. 5. Mass–size relation of the identified GMCs. Top: Mass–size Larson relation for GMCs derived from dust. Bottom: Mass–size Larson relation for GMCs derived from CO data. The various colors indicate the average density of individual clouds. The panel displays the slope κ and its corresponding error.
In the text

	Fig. 6. Power-law mass spectra of detected GMCs from dust-derived data (red) and CO (blue) as well as from NL16 (Milky Way, black). In black, the fit is performed with all listed structures in NL16, while the light gray dotted-dashed line shows the fit, where the CMZ of the Milky Way is excluded.
In the text

	Fig. 7. Mass, average density, surface mass density, and radius of the identified GMCs along the galactocentric radius. The dust-derived GMCs are shown in red, whereas the CO-derived GMCs are represented in blue. The bigger dark red and blue data points mark the GMCs that have a surface mass density exceeding 40 M⊙ pc2.
In the text

	Fig. 8. Aspect ratio and temperature from dust and CO-derived data along the galactocentric radius. The bigger dark red and blue data points mark the GMCs that have a surface mass density exceeding 40 M⊙ pc2.
In the text

	Fig. 9. Complementary cumulative NH2 distributions of the entire galaxy (based on dust-derived data) and the three galactic environments. We note that these distributions are solely based on pixels and are not connected to GMCs.
In the text

Fig. 10. Box plots of the determined dust-derived parameters categorized based on galactic environments. The lower and upper whiskers of the box plot represent the lowest and maximum values of the dataset, respectively. The colored box shows the distribution’s interquartile spread, or the range from the 25th to the 75th percentile; the median is indicated by the solid beige line inside the box. The distributions’ outliers are shown as circles.

In the text

	Fig. 11. Power-law mass spectra of GMCs in the dust-derived map for the three galactic environments. The power-law mass spectra is split into the three galactic environments center (pink), spiral arms (brown) and outskirts (turquoise). For comparison, the fit for the whole disk is shown again in red.
In the text

	Fig. A.1. 12CO(2 − 1) line-integrated intensity map of M33 (Druard et al. 2014). The map has been smoothed to the resolution of 18.2″ and re-gridded to the coordinate grid of the SPIRE 250 μm map.
In the text

	Fig. B.1. XCO factor (ratio) map of Method I of Paper I. The two ellipses represent a circular radius equivalent to 2 and 4 kpc.
In the text

	Fig. C.1. Distributions of GMC properties with varying Dendrogram parameters. The upper panels show histograms of radius and surface mass density with a min_value of 5σ. The lower panels display the same properties for a beam factor of 1.5.
In the text

	Fig. D.1. CO(1 − 0) luminosity-size relation of the structures identified in the CO data. The horizontal line represents the 2σ sensitivity limit.
In the text

	Fig. F.1. The dust-derived NH2 map shows the boundaries of three galactic environments: the center (dashed pink lines), spiral arms (dashed brown lines) and outskirts (remaining area). Center coordinates are RA(2000) = 1h33m50s, Dec(2000) = 30°39′37″ (SIMBAD 2024). The result of the log-spiral fit is shown in red.
In the text