Gas, Dust, Stars, Star Formation and their Evolution in M33 at Giant Molecular Cloud Scales

We report on a multi parameter analysis of giant molecular clouds (GMCs) in the nearby spiral galaxy M33. A catalog of GMCs identifed in 12CO(J=3-2) was used to compile associated 12CO(J=1-0), dust, stellar mass and star formation rate. Each of the 58 GMCs are categorized by their evolutionary stage. Applying the principal component analysis on these parameters, we construct two principal components PC1 and PC2 which retain 75% of the information in the original dataset. PC1 is interpreted as expressing the total interstellar matter content, and PC2 as the total activity of star formation. Young (<10Myr) GMCs occupy a distinct region in the PC1-PC2 plane, with lower ISM content and star formation activity compared to intermediate age and older clouds. Comparison of average cloud properties in different evolutionary stages imply that GMCs may be heated or grow denser and more massive via aggregation of diffuse material in their first ~10 Myr. The PCA also objectively identified a set of tight relations between ISM and star formation. The ratio of the two CO lines is nearly constant, but weakly modulated by massive star formation. Dust is more strongly correlated with the star formation rate than the CO lines, supporting recent findings that dust may trace molecular gas better than CO. Stellar mass contributes weakly to the star formation rate, reminiscent of an extended form of the Schmidt Kennicutt relation with the molecular gas term substituted by dust.


Introduction
The transition from gas to stars is a complex process involving various spatial and temporal scales. Observationally, much of the star formation take place within giant molecular clouds (GMCs) which are typically few×10 to ∼ 100 parsecs in size.
However, current knowledge of the general law of star formation is limited to semi-global scales over kilo-parsecs, where the relation between molecular gas and star formation is characterized by a simple power law, the Schmidt-Kennicutt (SK) relation (Kennicutt 1998;Komugi et al. 2005 and references therein).
The SK relation is known to break down at the scale of individual GMCs (Onodera et al. 2010) when using diffuse molecular gas tracers like the 12 CO(J = 1 − 0) emission. This has been attributed to the different evolutionary stage of GMCs (e.g., Kawamura et al. 2009) which become evident at the ∼100 parsec scale (Onodera et al. 2010). Komugi et al. (2012) support this by showing that the dispersion of the SK relation becomes small when star forming regions at the same age are used. When using denser and/or warmer gas such as that traced by the 12 CO(J = 3 − 2) line, the relation is generally tighter (Komugi et al. 2007;Muraoka et al. 2007;Iono et al. 2009;Muraoka et al. 2016) and this has been explained by the higher CO transition tracing molecular gas that is spatially and temporally closer to the actual site of star formation (c.f., Wu et al. 2005 using HCN emission). Onodera et al. (2012) reveal a connection between the two CO tracers and star formation, finding that larger GMCs tend to have higher 12 CO(J = 3 − 2)/ 12 CO(J = 1 − 0) ratios (r31), along with a higher star formation rate (SF R).
The connection of dust and star formation is also important, as dust grain surfaces act as a catalysis for molecular gas production. Recent studies (Wolfire et al. 2010;Paradis et al. 2012;Scoville et al. 2016) also indicate that dust may trace molecular gas better than CO in some regions where gas is not dense enough to shield itself against UV photon dissociation (Hollenbach & Tielens 1999). Stellar mass may be another parameter which controls star formation at various scales. Shi et al. (2011) and Rahmani et al. (2016) show that when the SK relation is extended to include the stellar mass, the dispersion becomes smaller both globally and within galaxies at kilo-parsec scales. Stellar contribution is also theoretically expected from stability analysis (Dib et al. 2017) and pressure regulated star formation (Blitz & Rosolowsky 2006 and others).
It is becoming increasingly clear that at the scale of GMCs, star formation is not explained by a single parameter equation such as the SK relation; it is better described as a time-evolving equilibrium state between ISM and star formation. In this paper, we attempt to formalize this state using a multi-parameter analysis on a sample of GMCs in the nearby spiral galaxy M33.
The catalog by Miura et al. (2012) also includes a measure of the evolutionary stage of the GMCs, Type A through D, based on stellar group identifications using the Massey et al. (2006) star catalog and Padova stellar synthesis tracks (Marigo et al. 2008;Girardi et al. 2010). In this paper, we use only Type B, C, and D GMCs, because only 1 Type A GMC was identified in the catalog, and further it was not detected at 1.1 mm. Type B are young GMCs which are associated with relatively small HII regions but not with young stellar groups (YSGs), and are at a stage approximately 3-7 Myr after the first formation of massive OB stars. Type C GMCs are associated with both HII regions and YSGs less than 10 Myr old, indicating that they are 10-20 Myr in age. Type D GMCs have HII regions and relatively old (10-30 Myr) stellar groups, and are 20-40 Myr in age (see Miura et al. 2012).
We identified cold dust clumps associated with each of these catalogued GMCs using the 1.1 mm continuum map from Komugi et al. (2011a). A box of 54 ′′ on the side, which is two times the detector FWHM of the observation, was centered on the HII region to search for a 1.1 mm peak. The angular size of the 1.1 mm beam covers the entire GMC, so its peak value can be used to obtain the total dust mass. We assume a dust emissivity of κ = 0.114 m 2 kg −1 , and corrected for the dust temperature using the color temperature between Spitzer MIPS 160 µm. Details of the 1.1 mm data and temperature derivation are found in Komugi et al. (2011a).
The stellar mass M * of the corresponding regions were derived from KS (2.1µm) image from the 2MASS Large Galaxy Survey (Jarrett et al. 2003). The image was background subtracted using planar interpolation from regions well away from the galaxy, and then convolved to 40 ′′ resolution to match the 1.1mm image. A Galactic extinction correction of 0.022 mag. was applied, based on the NASA Extragalactic Database (NED). The KS brightness was then converted to stellar mass using relations in Blitz & Rosolowsky (2006) assuming galactic inclination of 55 Peaks identified in 1.1mm and KS were both inspected by eye to ensure that they are discrete clumps and are associated to the HII region and GMCs identified in both CO lines. Out of the 71 regions catalogued in Miura et al. (2012), 58 GMCs were identified in all of the two CO lines, M d , SF R and M * . The final sample consists of 10 Type B, 30 Type C, and 18 Type D GMCs. All following analyses in this paper use this final sample.

Principal Component Analysis
We utilize the Principal Component Analysis (PCA) approach to the catalogue of GMCs. The PCA works on an arbitrary n number of parameters x, and constructs n new axes (principal components; PCs) using a linear combination of these n parameters, so that the new axes "better describe" the characteristics of the sample than the original parameters. Formally, where xi are explanatory parameters of the GMC samples, which are log SF R, log L ′ 3−2 , log L ′ 1−0 , log M * and log M d . Since the parameters xi have different dynamical range and units, they are normalized so that their average is zero and variance is unity. The coefficients α are normalized so that The amount of information that is projected onto the new PC can be measured by their variance (eigen value). In the present study, n = 5, and the total variance of 5 is conserved after the projection of the axes. The coefficients α represent the amount of contribution that each parameter xi has on PCi. The first principal component, PC1, is defined so that the variance of data around this axis is maximized. PC2 is then defined orthogonal to PC1, having the next largest variance, and so on. The n-th principle component PCn will have the minimum variance. In short, PC1 is a parameter which best describes the GMC sample, and PC2 follows. Although molecular clouds are characterized by their CO luminosities and their associated dust, stellar mass and SF R, the SK relation is usually described in terms of surface densities of each quantity. We have applied the PCA in two ways, to both luminosity (mass) and surface densities. The surface density of each parameters were calculated by dividing each quantity by the projected area of the molecular cloud as derived in Miura et al. (2012), and the stellar surface density was derived directly from the KS brightness.

Results
Tables 1 and 2 summarize the variances determined for the PCs. For both cases, PC1, PC2 and PC3 have variance adding up to more than 90% of the total information in the original dataset. For PC4 and PC5, the amount of information is significantly reduced from the original axes; the variance around these PCs decrease from 1 to 0.42 (0.29) and 0.07 (0.07), corresponding to only 8% (6%) and 1% (1%) of the total variance, respectively, where values in parentheses are for surface densities. Tables 3  and 4 summarize the coefficients α. The PCA is a purely mathematical operation, so a physical interpretation must be given to the new PCs. For a robust interpretation, we adopt only the terms that contribute significantly to the new PCs, and ig-nore terms with coefficients smaller than 0.45 (corresponding to α 2 = 1/n in equation 2 for n = 5). Coefficients which contribute significantly to the PCs are simliar for luminosities and surface densities, except for the case of PC3. They are shown by boldface in tables 3 and 4.

Robustness of PCA Analysis
The PCA analysis is sensitive to uncertainties in the original dataset. Therefore, we estimate the uncertainty in the variances, PCA scores and coefficients by using a bootstrapping method . We create a resampled dataset where values for each original parameter for each GMC follow a gaussian distribution with 30% standard deviation with respect to the original data. The assumed uncertainty of 30% corresponds to typical error in the L ′ 3−2 measurement, which has the largest observational error in the datasets, and therefore gives a conservative estimate of how the PCA analysis can be affected by intrinsic uncertainties in the dataset. We create 10 4 such resampled datasets, and run the PCA analysis on each of them. The error estimates in figure  1 and tables 1 through 4 correspond to the standard deviation of each values computed from bootstrapping. The errors for the variances, coefficients and PCA scores show that the PCA analyses are stable and results presented here are not severely affected by observational uncertainties.

PC1 and PC2
The first two PCs contain the most information in the dataset. PC1 has the largest variance, contributing ∼ 50% to the total value, thus best characterizes the star forming GMCs. From tables 3 and 4, PC1 is contributed almost equally by L ′ 1−0 , L ′

3−2
and M d which all express quantities of the ISM, all with positive coefficients. PC1 may be interpreted as a scaling factor indicating the total ISM content of the regions. PC2 correlates positively with SF R and stellar mass, both of which are related to star formation in the ISM. A higher SF R, and its integration over time (the stellar mass) lead to a higher PC2. It is also anticorrelated with L ′ 1−0 , which indicates that PC2 increases for a higher SF R per molecular gas mass, i.e., the efficiency of star formation. PC2 may be interpreted as a measure of the star formation activity, both instantaneous and over a longer timescale. The first two PCs combined account for ∼75% of total information. They can be written as PC ′ 2 = 0.60 ± 0.02 (0.65 ± 0.03)log M * + 0.48 ± 0.02 (0.49 ± 0.03)log SF R The eigen values of each of the PCs correspond to the amount of information that is projected onto the PC. The propertion of variance is its ratio to the total amount of information, i.e., the variance divided by n = 5.
where the prime on the left hand side indicate that insignificant terms on the right hand side have been omitted, and tilde on the right hand side indicate that the (logarithm of) values are normalized as explained in section 3. Values in parentheses indicate the case for surface densities (same hereafter). Figure  1 shows PC1 versus PC2, categorized by their evolutionary stages. For both luminosity and surface density, young Type B molecular clouds are distributed in a different region from intermediate (Type C) and older (Type D) clouds, such that they have less ISM content and star formation activity.

PC3
The third PC has a variance that is lower than unity. This indicates that the amount of information has decreased from the original axes, i.e., that PC3 does not explain the properties of molecular clouds better than the original parameters. It is dominated by contribution from the stellar mass, with smaller contribution from SFR in the case of luminosities, and dust in the case of surface densities. The physical interpretation of PC3 is unclear, as the contributing parameters are different between luminosities and surface densities. Here we will not attempt to give PC3 a physical interpretation.

PC4 and PC5
The last two PCs have the least variance in the given datasets, and can be used to detect relations between variables. Since the PC scores of individual clouds average to zero, we can write for luminosity, and

Molecular gas and star formation at GMC scales
A correlation between molecular gas and SF R is non-existent when molecular gas is traced by the 12 CO(J = 1 − 0) line (figure 2 top right and bottom right with correlation coefficient 0.18 and 0.20, respectively) but slightly more correlated when using higher density tracers like the 12 CO(J = 3 − 2) (figure 2 top middle and bottom middle with correlation coefficient 0.40 and 0.42, respectively). This is consistent with previous studies which find that molecular lines which trace denser or warmer molecular gas is more strongly correlated with star formation (Komugi et al. 2007;Muraoka et al. 2007;Iono et al. 2009;Muraoka et al. 2016), presumably because they trace gas that is spatially and temporally closer to star forming molecular gas. The correlation between SF R and dust mass is much more pronounced (equations 7 and 9, figure 2 top left and bottom left with correlation coefficient 0.64 and 0.66, respectively) than the relation for CO. If we assume that molecular gas is the direct ingredient of star formation, and that their quantities are intrinsically correlated, the strong dust-SFR correlation supports recent studies which claim that dust emission traces molecular gas better than CO (Wolfire et al. 2010;Paradis et al. 2012;Scoville et al. 2016). The basis behind this may be the existence of "COdark" molecular gas, in the diffuse outskirts of the molecular clouds where CO can be dissociated by UV photons but hydrogen molecules survive due to self-shielding (Hollenbach & Tielens 1999;Wolfire et al. 2010). In this case, the strong correlation between SF R and dust may be specific to M 33 because of its relatively low metallicity (12 + log(O/H) = 8.48 in the central region; Bresolin 2011), since the amount of COdark gas is inferred to increase in low metallicity systems (Israel 1997;Leroy et al. 2007;Komugi et al. 2011b). Alternatively, a strong correlation between dust and SF R could imply that the 1.1 mm observation is predominantly tracing dense gas that is linked directly to star formation, since millimeter dust is optically thin where CO would be optically thick and unable to accurately trace gas quantity.

GMC evolution
The distribution of GMCs in figure 1 has implications for how GMCs evolve after stars form. Directly after the onset of star formation represented by Type B clouds, both the total ISM content (PC1) and SF activity (PC2) increase significantly, both in luminosity and surface density. This is a natural consequence of all parameters contributing to PC1 and PC2 increasing from Type B to Type C/D (see table 5). Although the increase of each parameters contributing to PC1 and PC2 are insignificant or small considering their standard deviations, the PCs with these parameters combined result in the distinct distribution of ISM content (PC1) and SF activity (PC2) between Type B and Type C/D.
Nominally, an increase in log L ′ 3−2 can be attained either from collapse (increase in density) or heating of diffuse gas traced in 12 CO(J = 1 − 0). r31 indeed increases significantly from Type B to C (Miura et al. 2012). The increase in r31 could be due to gas heating by the increased star formation event, or alternatively, could indicate that ISM is being aggregated from the diffuse material. ISM build-up may explain the 60% increase in dust mass from Type B to C, although a part of the increased dust mass may be due to enrichment by red supergiants and/or the death of massive stars occuring during the ∼10 Myr lifetime of Type B clouds (Dunne et al. 2003;Hernández et al. 2006;Eldridge & Relaño 2011).
Clouds with larger mass tend to have higher r31, higher SF R and increased star formation efficiency (Onodera et al.  Figure 1 (right) shows that the density of ISM content and SF activity increase as well. This is because the average size of the clouds determined from 12 CO(J = 3 − 2), given in table 5 (taken from table 6 in Miura et al. (2012)) , indicate no size difference between Type B and C clouds. A slight increase in cloud size for Type D may explain the slight surface density decrease in PC1 and PC2 from Type C to D.

Dynamical Equilibrium of the Interstellar and Stellar Phases
Although parameters which contribute weakly (i.e., with coefficients smaller than 0.45 in tables 3 and 4) are dropped from the analyses above, it is worth noting some relations with these weakly contributing parameters to compare with previous studies. From tables 3 and 4, it can be seen that PC5 is almost completely dominated by contributions from L ′ 3−2 and L ′ 1−0 . However, the effect of log SF R is visible when the three parameters are seen in three dimensional space, where the two CO luminosities and the SF R comprise a plane expressed by for luminosities, and log Σ L ′ 3−2 = 1.0 log Σ L ′ 1−0 − 0.10 log ΣSF R + 2.9 (PC ′ 5 = 0) A number of previous studies have pointed out the possible contribution of stellar surface density to star formation (Blitz & Rosolowsky 2006;Shi et al. 2011;Rahmani et al. 2016;Dib et al. 2017). In this study also, the stellar mass contributes mildly to PC4 in case of surface densities. The plane formed by log ΣSF R, log Σ d and log Σ * are log ΣSF R = 3.1 log Σ d + 0.67 log Σ * − 5.5 (PC ′ 4 = 0) (15) log ΣSF R = (2.1 ± 0.2) log Σ d + (0.5 ± 0.2) log Σ * − (5.9 ± 0.4) least squares fitting (16) with a scatter of 0.4 dex (figure 3 right panel). As discussed in section 5.1, if dust indeed traces molecular gas better than CO lines at scales of individual molecular clouds, the relation given in equation 15 can be interpreted as an extended form of the SK relation that connects SF R with gaseous and stellar content, valid at GMC scales. The power law index of the stellar term (0.67 from PCA, 0.5 from least squares fitting) is roughly consistent with that in recent literature (0.36 in Shi et al. 2011, 0.74 in Rahmani et al. 2016) and theoretical expectations in the case of star formation triggered by two-fluid instabilities (Dib et al. 2017). Blitz & Rosolowsky (2006) show that for pressure regulated star formation, the stellar contribution has an index of roughly 0.5 and molecular gas density has an index of 2.0 (where equations 15 and 16 indicate that the dust term has an  index of 2 to 3) . It should be noted, however, that stellar contribution can be seen only for surface densities ((−0.26) 2 ∼ 6%) but is 0% in case of luminosities. This difference may either indicate that the extended star formation law is indeed sensitive to density, or simply that the appearance of a stellar contribution is due to small number statistics in this study.
An interesting feature of equations 13 and 15 is that the relation holds for clouds regardless of evolutionary stage. GMC types B, C and D mark evolution timescales of 10 Myr, over which the conversion from gas to stars proceed in tandem with feedback processes from stars, resulting in significant changes in distribution of ISM and stars. Thus, one may expect that clouds of different types would lie on a different relation (pro-vided that any relations exist at this scale). The fact that GMCs of different types spanning ∼40 Myr lie on the same planes shown in figure 3 indicate that equations 13 and 15 act as boundary conditions for how physical parameters can change, i.e., can be interpreted as a kind of equation of state describing the dynamical equilibrium between the ISM and stellar phases.

Future Work
This paper presents a case study on a late type galaxy which is relatively metal poor. There are a number of other physical parameters which are expected to relate to star formation. In particular, metallicity and dynamical parameters (velocity dis- persion and larger scale shear) are not included in the PCA analysis presented here. Further studies should include a significant increase in sample size, spanning a broader range of quantified and catalogued characteristics.