A High-Resolution, Dust-Selected Molecular Cloud Catalogue of M33, the Triangulum Galaxy

TL;DR

This paper presents a detailed catalogue of 165 Giant Molecular Clouds in M33, analyzing their properties, distribution, and mass function, revealing differences from the Milky Way and insights into cloud formation.

Contribution

The study provides the first high-resolution dust-selected GMC catalogue for M33, including measurements of dust properties, gas-to-dust ratios, and mass functions, with comparisons to previous CO studies.

Findings

GMCs range from 10^4 to 10^7 solar masses.

GDR decreases with galactocentric radius, lower than in the Milky Way.

Mass function slope is steeper than in other local galaxies.

Abstract

We present a catalogue of Giant Molecular Clouds (GMCs) in M33, extracted from cold dust continuum emission. Our GMCs are identified by computing dendrograms. We measure the spatial distribution of these clouds, and characterise their dust properties. Combining these measured properties with CO(J=2-1) and 21cm HI data, we calculate the gas-to-dust ratio (GDR) of these clouds, and from this compute a total cloud mass. In total, we find 165 GMCs with cloud masses in the range of 10$^4$-10$^7$ M$_\odot$. We find that radially, $\log_{10}(\mathrm{GDR}) = -0.043(\pm0.038) \,\mathrm{R [kpc]} + 1.88(\pm0.15)$, a much lower GDR than found in the Milky Way, and a correspondingly higher $\alpha_{\rm CO}$ factor. The mass function of these clouds follows a slope proportional to M$^{-2.84}$, steeper than many previous studies of GMCs in local galaxies, implying that M33 is poorer at forming massive clouds than other nearby spirals. Whilst we can rule out interstellar pressure as the major contributing factor, we are unable to disentangle the relative effects of metallicity and HI velocity dispersion. We find a reasonably featureless number density profile with galactocentric radius, and weak correlations between galactocentric radius and dust temperature/mass. These clouds are reasonably consistent with Larson's scaling relationships, and many of our sources are co-spatial with earlier CO studies. Massive clouds are identified at large galactocentric radius, unlike in these earlier studies, perhaps indicating a population of CO-dark gas dominated clouds at these larger distances.

Tables

Table 1. Table 1: SCUBA-2 data reduction parameters for both the main data reduction, and the data calibration.

Parameter	Value	Description
Data reduction
maptol	0.005	Defines when the map has ‘converged’.
com.perarray	0	Calculate a single common-mode signal for all subarrays.
flt.filt_edge_largescale	120 (450µm), 320 (850µm)	Specifies the largest scale structure to be recoverable in the reduction.
ast.zero_mask	1	Provide an external astronomical signal mask, based on the Herschel 500µm image.
ast,flt,com.zero_freeze	0	Calculate these masks every iteration.
com.sig_limit	5	Remove high-frequency ‘blobs’ from the map.
flt.filt_order	4	Reduce ringing around bright sources.
flt.ring_box1	0.5	Reduce ringing around bright sources.
flagslow	300	Flag data where sources are obscured by 1/f noise.
Calibration
ast.zero_mask	0	Do not use an external mask.
ast,flt.zero_circle	0.033	Use a circular mask of 120 arcsec radius.
ast.mapspike	10	Ensure very bright pixels are included in the final map.
dcthresh	10000	Ensure very bright pixels are included in the final map.

Table 2. Table 2: Leaf node parameters for the GMCs. Generally, errors given are 1 σ 𝜎 \sigma errors. However, in the case of an 3 σ 𝜎 \sigma upper limit, the given flux is 0 and the error given is the 3 σ 𝜎 \sigma upper limit.

GMC ID	R.A. (J2000)	Dec (J2000)	R (kpc)	FWHM (pc)	$S_{100}$ (Jy)	${\sigma }_{100}$	$S_{160}$ (Jy)	${\sigma }_{160}$	$S_{250}$ (Jy)	${\sigma }_{250}$	$S_{450}$ (Jy)	${\sigma }_{450}$	$S_{850}$ (Jy)	${\sigma }_{850}$
0	1 33 56.0	30 22 22.4	6.9	152	0.43	0.01	0.43	0.02	0.55	0.01	0.135	0.009	0.036	0.004
1	1 33 38.5	30 23 04.7	6.6	172	0.83	0.03	1.24	0.03	0.84	0.01	0.151	0.008	0.046	0.002
2	1 33 56.9	30 25 14.2	5.8	87	0.06	0.01	0.18	0.01	0.15	0.01	0.020	0.002	0.005	0.001
3	1 33 22.7	30 26 02.6	5.9	236	1.97	0.04	2.83	0.05	1.97	0.02	0.495	0.012	0.142	0.004
4	1 34 12.8	30 25 50.0	5.8	155	0.37	0.02	0.64	0.02	0.54	0.01	0.087	0.009	0.032	0.002
5	1 32 56.3	30 26 02.0	7.1	105	0.17	0.01	0.21	0.01	0.19	0.01	0.021	0.005	0.010	0.002
6	1 32 57.5	30 27 22.9	6.6	178	0.34	0.02	0.54	0.02	0.38	0.01	0.077	0.009	0.028	0.002
7	1 33 11.6	30 27 27.5	5.8	195	3.42	0.04	3.40	0.05	2.11	0.02	0.444	0.010	0.130	0.003
8	1 33 56.8	30 27 17.2	4.9	134	0.25	0.02	0.38	0.02	0.28	0.01	0.038	0.003	0	0.004
9	1 34 14.5	30 27 39.2	5.2	176	0.47	0.03	0.69	0.02	0.52	0.01	0.076	0.007	0.009	0.001
10	1 33 47.3	30 27 28.6	4.8	85	0.22	0.01	0.22	0.02	0.14	0.01	0	0.012	0.004	0.001
11	1 33 38.6	30 28 14.9	4.6	140	0.51	0.02	0.68	0.02	0.43	0.01	0.096	0.006	0.013	0.001
12	1 33 57.5	30 29 12.3	4.2	144	0.18	0.02	0.30	0.02	0.22	0.01	0.040	0.005	0.010	0.002
13	1 33 39.2	30 29 56.4	3.9	146	0.61	0.03	0.96	0.03	0.73	0.01	0.152	0.006	0.044	0.001
14	1 33 09.0	30 29 42.3	5.3	180	2.07	0.07	2.09	0.06	1.34	0.03	0.211	0.012	0.050	0.004
15	1 33 35.8	30 29 39.5	4.1	99	0.24	0.02	0.25	0.02	0.23	0.01	0.039	0.004	0.007	0.001
16	1 33 17.2	30 31 13.0	4.4	177	0.61	0.03	1.02	0.02	0.72	0.02	0.142	0.007	0.029	0.001
17	1 34 01.8	30 30 53.6	3.6	150	0.23	0.01	0.44	0.02	0.30	0.01	0.041	0.005	0.010	0.001
18	1 33 00.7	30 30 49.4	5.5	139	1.49	0.03	1.49	0.04	0.90	0.02	0.201	0.008	0.043	0.002
19	1 33 06.5	30 30 49.2	5.1	81	0.12	0.02	0.16	0.02	0.15	0.01	0.051	0.004	0.006	0.001
20	1 33 45.0	30 31 09.9	3.4	64	0.10	0.01	0.12	0.02	0.10	0.01	0.011	0.003	0.003	0.000
21	1 34 40.4	30 31 13.2	5.4	79	0.10	0.01	0.18	0.01	0.12	0.01	0.015	0.003	0	0.003
22	1 34 42.4	30 32 01.4	5.3	118	0.30	0.02	0.35	0.02	0.32	0.01	0.046	0.005	0.017	0.001
23	1 33 58.2	30 32 04.1	3.1	103	0.41	0.04	0.55	0.04	0.42	0.01	0.077	0.007	0.011	0.001
24	1 33 29.7	30 31 52.5	3.5	94	2.87	0.08	2.43	0.10	1.13	0.05	0.176	0.013	0.038	0.003
25	1 32 58.7	30 31 52.1	5.4	76	0.10	0.01	0.12	0.02	0.13	0.01	0.024	0.003	0.004	0.001
26	1 33 23.7	30 31 55.4	3.8	95	0.14	0.03	0.25	0.03	0.21	0.02	0.056	0.005	0.010	0.001
27	1 34 10.0	30 31 58.9	3.5	70	0.14	0.01	0.17	0.01	0.12	0.01	0.019	0.004	0.005	0.001
28	1 33 44.9	30 32 04.9	3.0	62	0.24	0.02	0.25	0.03	0.19	0.01	0.018	0.005	0.004	0.001
29	1 33 33.9	30 32 07.7	3.3	97	2.13	0.06	2.25	0.08	1.24	0.03	0.202	0.007	0.040	0.002
30	1 33 37.4	30 32 04.4	3.2	53	0.20	0.03	0.23	0.04	0.17	0.02	0.029	0.004	0.005	0.001
31	1 32 51.0	30 33 02.6	5.7	160	0.26	0.02	0.29	0.02	0.37	0.01	0.108	0.007	0.020	0.002
32	1 33 00.7	30 32 33.0	5.1	112	0.22	0.01	0.34	0.02	0.26	0.01	0.033	0.006	0.015	0.001
33	1 34 23.2	30 32 39.5	3.9	93	0.17	0.03	0.27	0.02	0.21	0.01	0.034	0.004	0.010	0.001
34	1 33 47.6	30 32 51.4	2.7	109	1.75	0.08	1.74	0.09	1.03	0.03	0.168	0.006	0.039	0.002
35	1 34 08.6	30 32 50.7	3.1	100	0.23	0.01	0.33	0.03	0.27	0.01	0.069	0.006	0.008	0.001
36	1 33 59.9	30 32 43.4	2.9	47	0.17	0.01	0.16	0.02	0.09	0.01	0	0.008	0.003	0.001
37	1 33 18.6	30 33 13.6	3.7	144	0.33	0.02	0.46	0.02	0.33	0.02	0.061	0.006	0.011	0.002
38	1 33 53.3	30 33 12.5	2.6	90	0.86	0.06	0.76	0.08	0.43	0.04	0.059	0.008	0.013	0.002
39	1 33 43.9	30 33 10.5	2.6	123	1.10	0.06	1.19	0.07	0.82	0.03	0.136	0.008	0.029	0.002
40	1 34 17.7	30 33 43.2	3.3	75	0.90	0.05	0.90	0.06	0.56	0.03	0.115	0.006	0.021	0.001
41	1 34 13.9	30 33 45.1	3.1	86	2.67	0.07	2.53	0.08	1.48	0.03	0.284	0.011	0.048	0.002
42	1 34 09.2	30 34 21.1	2.6	130	0.24	0.04	0.42	0.04	0.29	0.02	0.077	0.007	0	0.006
43	1 33 50.5	30 33 47.6	2.3	52	0.19	0.03	0.20	0.03	0.15	0.02	0.027	0.004	0.004	0.001
44	1 33 01.0	30 34 33.2	4.7	180	1.76	0.03	2.12	0.04	1.28	0.02	0.222	0.008	0.040	0.002

Table 3. Table 6: Calculated dust and gas properties for the GMCs. Generally, errors given are 1 σ 𝜎 \sigma errors. However, in the case of an 3 σ 𝜎 \sigma upper limit, the reported value is 0 and the corresponding error is that 3 σ 𝜎 \sigma upper limit.

GMC ID	T (K)	${\sigma }_T$	log($M_{\mathrm{dust}}$ [$M_{\odot }$])	${\sigma }_{\log \left(M_{\mathrm{dust}}\right)}$	log($L_{\mathrm{TIR}}$ [$L_{\odot }$])	${\sigma }_{\log \left(L_{\mathrm{TIR}}\right)}$	$L_{\mathrm{CO}}$ (K km s-1)	${\sigma }_{L_{\mathrm{CO}}}$	$L_{\mathrm{CO}}/{\sigma }_{L_{\mathrm{CO}}}\times {10}^x$	${\Sigma }_{H\text{i}}$ (M⊙ pc-2)	${\sigma }_{{\Sigma }_{H\text{i}}}$
0	20.95	0.54	3.76	0.05	5.54	0.02	1.95	0.15	-1	5.17	0.18
1	21.04	0.45	4.05	0.04	5.83	0.02	2.72	0.12	-1	2.73	0.18
2	18.13	0.42	3.43	0.04	4.86	0.03	9.47	1.60	-2	1.87	0.21
3	20.98	0.43	4.43	0.04	6.21	0.02	2.94	0.13	-1	3.62	0.17
4	19.58	0.41	3.92	0.04	5.54	0.02	1.25	0.10	-1	4.52	0.34
5	21.14	0.64	3.34	0.05	5.14	0.03	1.93	0.16	-1	2.16	0.35
6	20.44	0.50	3.74	0.04	5.46	0.02	2.74	0.85	-2	1.83	0.10
7	24.54	0.57	4.26	0.04	6.41	0.03	5.18	0.17	-1	2.97	0.18
8	22.76	0.61	3.38	0.04	5.35	0.03	4.03	1.26	-2	2.95	0.16
9	21.97	0.50	3.71	0.04	5.61	0.02	0	4.11	-2	4.22	0.21
10	25.55	0.78	2.98	0.05	5.23	0.03	2.18	0.28	-1	0	0.60
11	22.20	0.49	3.69	0.04	5.60	0.02	2.72	0.22	-1	1.20	0.17
12	20.20	0.61	3.51	0.05	5.20	0.03	9.63	1.75	-2	1.09	0.21
13	20.24	0.46	4.03	0.04	5.73	0.02	2.25	0.10	-1	3.18	0.23
14	24.50	0.61	4.04	0.04	6.20	0.03	3.34	0.33	-1	1.87	0.26
15	22.22	0.74	3.35	0.06	5.27	0.03	1.87	0.17	-1	1.82	0.35
16	20.29	0.43	4.02	0.04	5.72	0.02	1.93	0.12	-1	3.68	0.21
17	20.26	0.42	3.62	0.04	5.31	0.02	1.67	0.21	-1	1.38	0.11
18	24.56	0.60	3.89	0.04	6.05	0.03	2.01	0.21	-1	3.01	0.42
19	19.93	0.89	3.34	0.08	5.00	0.04	2.13	0.46	-1	0	1.42
20	22.57	0.93	2.95	0.07	4.90	0.05	3.17	0.28	-1	1.51	0.28
21	21.81	0.68	3.08	0.06	4.95	0.03	9.01	2.85	-2	0	1.17
22	21.47	0.58	3.54	0.05	5.38	0.03	2.48	0.19	-1	2.04	0.19
23	21.93	0.66	3.65	0.05	5.54	0.03	1.97	0.28	-1	1.82	0.15
24	28.33	0.85	3.84	0.05	6.34	0.03	1.57	0.14	0	0	3.19
25	20.59	0.74	3.19	0.06	4.92	0.03	2.94	0.26	-1	0	1.09
26	19.01	0.86	3.57	0.08	5.11	0.04	6.97	0.93	-1	0	1.53
27	22.70	0.79	3.07	0.06	5.04	0.04	2.63	0.37	-1	1.58	0.48
28	24.69	0.94	3.10	0.06	5.27	0.04	4.20	0.55	-1	2.38	0.44
29	24.96	0.63	4.01	0.04	6.20	0.03	8.75	0.89	-1	2.81	0.42
30	23.39	1.23	3.15	0.08	5.19	0.06	7.19	1.13	-1	0	2.89
31	19.30	0.63	3.75	0.06	5.33	0.03	9.17	1.47	-2	2.61	0.23
32	20.90	0.54	3.50	0.05	5.27	0.03	1.80	0.19	-1	1.72	0.53
33	20.79	0.84	3.43	0.06	5.18	0.04	4.79	0.35	-1	2.58	0.37
34	24.99	0.71	3.92	0.05	6.12	0.03	8.52	0.32	-1	2.97	0.26
35	20.64	0.51	3.55	0.05	5.29	0.02	5.40	0.33	-1	2.62	0.31
36	28.84	1.51	2.58	0.09	5.12	0.05	2.08	0.48	-1	0	1.67
37	21.52	0.56	3.58	0.05	5.42	0.03	1.23	0.33	-1	2.95	0.34
38	27.23	1.19	3.41	0.06	5.81	0.04	3.69	0.59	-1	0	2.62
39	23.45	0.69	3.88	0.05	5.92	0.03	6.14	0.37	-1	1.23	0.34
40	24.33	0.80	3.69	0.06	5.82	0.03	4.91	0.32	-1	2.24	0.40
41	25.41	0.66	4.06	0.04	6.30	0.03	1.51	0.07	0	0	1.53
42	21.65	0.83	3.50	0.06	5.35	0.04	2.96	0.29	-1	0.99	0.20
43	24.06	1.34	3.05	0.08	5.16	0.06	0	4.82	-1	2.12	0.69

Equations

$$F(\beta ,T)={\left(\frac{500\mu \mathrm{m}}{450\mu \mathrm{m}}\right)}^{3+\beta }\times \frac{\exp \left(\frac{\mathrm{hc}}{500\mu \mathrm{m}\times \mathrm{kT}}\right)}{\exp \left(\frac{\mathrm{hc}}{450\mu \mathrm{m}\times \mathrm{kT}}\right)}$$

$$C_{\mathrm{Planck}}=\frac{\int R(\nu )(\nu /353{)}^{-1}\mathrm{d}\nu }{\int R(\nu )(\nu /353{)}^{\alpha }\mathrm{d}\nu }$$

$$S_{\nu }=\frac{{\kappa }_{\nu }{\mathrm{M}}_{\text{dust}}B(\nu ,T_{\text{dust}})}{D^2},$$

$${\kappa }_{\nu }={\kappa }_{{\nu }_0}{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta },$$

$${\Sigma }_{\text{Hsci}}=1.8\times 10^{18}{\text{cm}}^{-2}/({\text{K km s}}^{-1}).$$

$$\mathrm{GDR}\times {\Sigma }_{\mathrm{Dust}}={\Sigma }_{\mathrm{HI}}+{\alpha }_{\mathrm{CO}}\times I_{\mathrm{CO}}.$$

$${\log }_{10}(\mathrm{GDR})=-0.043(\pm 0.038)\mathrm{R}[\mathrm{kpc}]+1.88(\pm 0.15).$$

$${\sigma }_{\mathrm{v}}={\left(\frac{\pi G}{5}\right)}^{1/2}{\Sigma }^{1/2}{R}^{1/2}.$$

$${\sigma }_{\mathrm{v}}=\frac{I_{\mathrm{CO}}}{\sqrt{2\pi }{T}_{\mathrm{peak}}},$$

Full text

††pubyear: 2018††pagerange: A High-Resolution, Dust-Selected Molecular Cloud Catalogue of M33, the Triangulum Galaxy–C

1 Introduction

The study of star-formation and the study of molecular clouds are inexorably linked. As stars are believed to form from the dense molecular gas in these clouds (André et al., 2010; Lada et al., 2010), our understanding of star-formation is ultimately limited by our ability to resolve ensembles of these star-forming regions. Within our own galaxy, we are faced with the challenges of distance ambiguity – to overcome this, we can turn to high-resolution mapping of galaxies for studies of large numbers of these molecular clouds.

One option for locating these molecular clouds is to trace the molecular hydrogen that they contain. However, due to the size and symmetry of the H2 molecule, it is impossible to trace the cold component associated with star formation directly and so a proxy must be employed. Generally, the rotational transitions of CO (the next most common molecule) are favoured, as they are believed to trace the cold molecular gas that resides within these clouds. Resolving these molecular clouds poses a great challenge – with the average Milky Way (MW) GMC size being $\sim$40 pc (Solomon et al., 1979), and $\sim$30 pc in the LMC (Hughes et al., 2010), we are limited to studies in our local Universe (e.g. Israel et al. 1993; Rosolowsky 2007; Hughes et al. 2010). Recently, with the advent of the Atacama Large Millimetre/submillimetre Array (ALMA), these studies can be extended beyond our Local Group of galaxies (e.g. Liu et al. in prep.).

Alternatively, an independent method to probe the properties of GMCs uses the cold dust continuum emission of a galaxy. It has long been established that there is a link between the dust content of a galaxy and its molecular gas (e.g. Hildebrand 1983; Magdis et al. 2012; Eales et al. 2012). Thus, the dust continuum allows us an alternative method to CO measurements to probe the properties of GMCs. However, due to the limited resolution of these instruments and the sizes of clouds this method of probing GMCs is only suitable for some of our most nearby galaxies. Using, for example, the Herschel Space Observatory (Pilbratt et al., 2010), we can resolve an average-sized molecular cloud up to a distance of around 200kpc at 500µm wavelengths. With the Submillimetre Common-User Bolometer Array 2 (SCUBA-2; Holland et al. 2013) on the James Clerk Maxwell Telescope (JCMT), we can resolve these objects up to 600 kpc away (850µm), or 1.2 Mpc (450µm). However, with ground-based sub-mm observatories we must overcome noise from the sky varying over small scales at the sub-mm wavelengths we probe – a harsh sky subtraction process must be performed, which has the drawback of also filtering out large-scale structure in these galaxies. Using a Fourier combination technique, we can use space-based observatories operating at similar wavelengths to add this large-scale structure back in to this data, allowing us to retain both the large-scale structure and the much finer structure these ground-based observatories offer.

M33 provides an excellent laboratory for resolved molecular cloud studies. Located at a distance of 840 kpc (Madore & Freedman, 1991), it is the third massive spiral galaxy of our Local Group, behind our own Milky Way (MW), and Andromeda (M31). Unlike M31, however, M33 is more face-on, with an inclination of 56°  (Regan & Vogel, 1994), and so suffers less from projection effects. It is also actively star-forming across its disk (Heyer et al., 2004), and is host to a large number of GMCs. Previous studies of M33 have identified GMCs using line data from 12CO(J=1–0), such as Wilson & Scoville (1990), surveying the inner 2 kpc of M33 at 7 arcsec resolution, finding 38 GMCs. All-disk surveys of M33 have suffered from poorer resolution than this, such as Engargiola et al. (2003), using the J=1–0 line, and Gratier et al. (2012), using the J=2–1 line, finding 148 and 337 GMCs across the disk of M33, respectively. Both of these surveys have resolutions of $\sim$50 pc, and so many of the GMCs are only marginally resolved.

In this work, we take an alternative approach to map the GMC content of M33. By combining far-infrared and sub-millimetre data, we probe the properties of GMCs via the cold dust continuum emission of M33. The layout of this paper is as follows: we first present an overview of the data used in our study (Sec. 2), and our method of source extraction (Sec. 3). We then move on to measure the properties of these GMCs (Sec. 4) and a comparison to earlier CO surveys (Sec. 5). Finally, we summarise our main results (Sec. 6).

2 Data

2.1 Far Infrared/sub-millimetre

Our first source of FIR/sub-mm data comes from the Herschel Space Observatory. We make use of observations taken as part of the Herschel M33 extended survey (HerM33es, Kramer et al. 2010), which mapped a 70 arcmin2 region around M33. Data at 100 and 160µm was taken with the Photoconductor Array Camera and Spectrometer (PACS, Poglitsch et al. 2010), with beam sizes of 7.7 arcsec and 12 arcsec, respectively. The details of this data reduction are presented in Boquien et al. (2011) and Boquien et al. (2015). This data has a Root Mean Squared (RMS) noise level of 2.6 mJy pixel*-1* (100µm) and 6.9 mJy pixel*-1* (160µm).

HerM33es simultaneously used the Spectral and Photometric Imaging Receiver (SPIRE, Griffin et al. 2010) aboard Herschel, which mapped M33 at 250µm, 350µm, and 500µm  with a resolution of 18 arcsec, 25 arcsec, and 36 arcsec, respectively. This data covers the same region as the PACS maps, to an RMS noise level of 14.1, 9.2, and 8 mJy beam*-1* at 250, 350 and 500µm, respectively.

Archival SCUBA-2 observations of M33 at 450 and 850µm were taken between 2012-07-01 and 2012-07-12, consisting of $\sim$7 hours of pong1800 (which maps a roughly circular, 30 arcmin field) observations of M33, and $\sim$4 hours of smaller, cv daisy (constant velocity, small field-of-view) observations. For more details of these SCUBA-2 observing modes, we refer readers to the JCMT observing mode webpage111https://www.eaobservatory.org/jcmt/instrumentation/continuum/scuba-2/observing-modes/. These observations were taken in mostly Band 2/Band 3 weather (225 GHz opacity, $0.04\leq\tau_{225}\leq 0.12$). Due to our adopted reduction parameters (see the details on flagslow in Section 2.1.1), we cannot use these daisy maps in our reduction, and so for our purposes, this archival data reaches an RMS noise level of $\sim$6 mJy beam*-1* at 850µm, and $\sim$85 mJy beam*-1* at 450µm (with pixel sizes of 4 arcsec and 2 arcsec respectively). As we are particularly interested in the resolution the 450µm data provides, we found that this RMS noise was inadequate and so between 2017-10-17 and 2017-11-21, under Program ID M17BP003 (PI W.K.G.), we obtained a further 12 hours of pong1800 observations of M33, in good Band 1 weather ($\tau_{225}<0.05$). In the following sections, we describe the data reduction process, which allowed us to create 450µm  and 850µm maps of M33 with RMS noise levels of $\sim$35 mJy beam*-1* and $\sim$4 mJy beam*-1*, respectively. An initial reduction of the data was first presented in Williams et al. (2018), but in this work we detail this new reduction. The entire dataset used to measure the dust continuum of our GMCs can be seen in Fig. 1.

The resolution of this SCUBA-2 data is 7.9 arcsec at 450µm and 13 arcsec at 850µm (Dempsey et al., 2013), corresponding to 32 pc and 52 pc at the distance of M33. However, due to atmospheric variations, extended large-scale structure is filtered out in the reduction process. In order to restore this, we make use of complementary Herschel 500µm data for the 450µm data and Planck 353GHz data for the 850µm map. A similar technique has previously been employed with Atacama Pathfinder Experiment Telescope (APEX) Large APEX BOlometer CAmera (LABOCA) data (Csengeri et al., 2016) to recover large-scale, extended structure in the Galactic plane, but we have tailored this technique to SCUBA-2 data.

2.1.1 SCUBA-2 Data Reduction

The SCUBA-2 data reduction pipeline, makemap, is described in detail in Chapin et al. (2013), and we refer readers to this work for a full description. We used a modified version of this algorithm, called skyloop, which performs a single makemap iteration each time, including data from all individual observations simultaneously. This helps to constrain the map, and reduce spurious extended emission, which is particularly important for SCUBA-2 observations of local, extended galaxies.

makemap is invoked with a file containing the parameters for the map maker. We have attempted to recover some large-scale structure in the SCUBA-2 maps, and so have based our reduction strategy on that of the JCMT Plane Survey (JPS, Eden et al. 2017). Our most important, non-default parameters are summarised in Table 1 – for a more detailed description of these parameters, we refer the reader to the SCUBA-2 Data Reduction Cookbook222http://starlink.eao.hawaii.edu/docs/sc21.htx/sc21.html.

2.1.2 SCUBA-2 Calibration

makemap produces an output file in units of pW, so it is necessary to apply a flux conversion factor (FCF) to the data, to convert it into units of Jy beam*-1*. The standard FCFs have been calculated to be 491 Jy beam*-1* pW*-1* at 450µm, and 537 Jy beam*-1* pW*-1*  at 850µm (Dempsey et al., 2013), but can vary during the night due to effects such as variations in seeing. Particularly for observations near the start of the night, dish cooling can have a major impact on the measured FCF. It is also important to note that the standard FCFs are calculated using a standard configuration file tailored for bright, compact sources, and the configuration parameters can also have an effect. We therefore calibrated our data using FCFs calculated from standard calibrators taken on the same night as the observations. These calibration observations are taken from Mars, Uranus, CRL618, CRL2688, or HL Tau. For observations of M33 between calibrator observations, we take a linear interpolation between the nearest calibrator FCF before and after. In the case that we did not have a calibrator observed either before or after, we took the FCF of the nearest calibrator. We reduced these calibrator observations using the same configuration file as our M33 reduction, with some small modifications (see Table 1). Along with these, we also removed the flagslow parameter, as since these calibration observations are daisys, rather than the larger pongs, the telescope was moving slowly enough that all data were flagged.

Using this reduction method, we find an average FCF of 522$\pm$51 Jy beam*-1* pW*-1* at 450µm (6% higher than the standard FCF), and 518$\pm$44 Jy beam*-1* pW*-1*  at 850µm (4% lower than the standard FCF). The scatter in FCF is similar to the 10% at 450µm found by (Dempsey et al., 2013), but higher than the standard 5% scatter at 850µm. Having calculated an FCF for each observation, we then multiplied the raw data by the ratio of the calculated to the standard FCF. After then reducing the data using skyloop, we multiplied the final map by the standard FCF value. We found that calibrating the data in this way led to an increase in flux of $\sim$3% in the 450µm map, and a negligible change in the 850µm map compared to simply using the standard FCF on the final map. We also found a decrease in noise of $\sim$3% in the 450µm map, and $\sim$15% in the 850µm map.

2.1.3 Combination with Herschel and Planck Data

As previously mentioned, the SCUBA-2 data reduction process necessarily removes extended structures in the map. However, using a method similar to interferometric ‘feathering’, we can restore this extended structure. Previous work has shown that this technique can work to combine Planck and LABOCA data (Csengeri et al., 2016), but we have tailored this code for SCUBA-2.

First, the units of the two input maps are converted to Jy beam*-1*, if necessary. If a SCUBA-2 map is provided in units of pW, the standard FCF is applied. Generally, SPIRE 500µm maps are in units of MJy sr*-1*, so we convert to Jy beam*-1* using a beam size of 1665 arcsec2 (As reported in the SPIRE Handbook333http://herschel.esac.esa.int/Docs/SPIRE/html/spire_om.html). The Planck maps (which are publicly available in HEALPIX format444https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/) are provided in units of KCMB temperature units, so we convert to Jy beam*-1* using a conversion factor of 287.45 MJy sr*-1* K${}_{\mathrm{CMB}}^{-1}$ (Planck Collaboration et al., 2014), and beam FWHM of 5.19 arcmin and 4.52 arcmin (Planck Collaboration et al., 2013). We also subtract the contribution of the Cosmic Microwave Background (CMB) from this data following Planck Collaboration et al. (2015), as the CMB varies over scales similar to the extent of M33. We then reproject these maps to the image size and pixel scale of the SCUBA-2 data using Python’s reproject package.

There are two corrections that must also be applied to the data, to account for the difference in central wavelengths, and colour corrections due to differences in spectral response. In the case of combining SCUBA-2 850µm and Planck 353GHz, the central wavelength correction is negligible. For the Herschel data, we perform a central frequency correction, assuming a modified blackbody (MBB), so

[TABLE]

where $\beta$ is the dust emissivity index (if not specified, defaults to 2) and T is the dust temperature (with a default value of 20K).

The colour correction to the Planck data is calculated using

[TABLE]

where $R(\nu)$ is the Planck 353GHz passband. $\alpha$ is the index of the source spectrum. In the Rayleigh-Jeans spectral regime, $\alpha=2+\beta$, which gives a default correction factor of 0.854. In the case of the Herschel data, we use a factor 1.0049, the colour correction given in Table 5.2 of the SPIRE Handbook for extended sources.

We perform a background subtraction on the Planck and SPIRE 500µm data, using a 3$\sigma$ clipped median. As the SCUBA-2 reduction pipeline models and subtracts the sky, we perform no further sky subtraction on the SCUBA-2 data. Our code applies a Gaussian filter when combining the data, specified by an inputted FWHM. In our case, we set the FWHM to 36 arcsec for the 450µm data, and 8 arcmin for the 850µm data. If this value is too small, negative bowling will be present around bright sources, and conversely, if set too high the fine detail desired is lost. We perform Fast Fourier Transforms (FFTs) on the data and the filter to transform them into the uv plane, and create parity between the Jy/beam units by multiplying by the volume ratio of the high- and low-resolution beams. The filter is normalised such that its amplitude at the centre of the uv plane is 1.

The FFT of the low-resolution data is then filtered by multiplying by the FFT of the filter, added to the FFT of the high-resolution data and transformed back into the image plane. There is an inherent uncertainty due to errors in $\beta$ and T, but in practice these are negligible. The total flux density should be determined by the low-resolution map, and we find that the flux density of the low-resolution data alone and the combined data are consistent to well within the calibration uncertainty of the SCUBA-2 data.

We homogenised this dataset to a common resolution (that of the SPIRE 250µm  image) and pixel scale. We convolved the data using the method of Aniano et al. (2011), and regrid to pixel sizes of 6 arcsec, to ensure that our maps are Nyquist sampled. This regridding is performed using Python’s reproject routine, which also astrometrically aligns each image.

2.2 Gas Data

We also make use of atomic and molecular gas data in this study. Hi is traced via the 21cm line from archival VLA555https://science.nrao.edu/facilities/vla/archive/index B, C, and D array data (reduced by Gratier et al. 2010). The CO data used in this investigation was taken as part of IRAM’s M33 Survey Large Program666http://www.iram.fr/ILPA/LP006/ (Gratier et al., 2010; Druard et al., 2014), which traces the molecular gas out to a radius of 7 kpc using IRAM’s Heterodyne arRAy (HERA, Schuster et al. 2004) instrument. This data has an angular resolution of 12 arcsec and a spectral resolution of 2.6 km s*-1*. These maps can be seen in Fig. 2.

3 GMC Catalogue

3.1 Identifying GMCs

Disentangling sources from regions of complex emission is a non-trivial task, and several source extraction methods have been developed to achieve this goal (see Men’shchikov et al. 2012 for descriptions of a number of source extraction algorithms). Initial testing using the algorithms clumpfind (Williams et al., 1994) and fellwalker (Berry 2015, an algorithm developed to deal with some issues in clumpfind) revealed shortcomings in these more traditional methods – given that much of the emission at these wavelengths is diffuse, the entire galaxy becomes segmented into unreasonably large “sources”. We also attempted source extraction using SExtractor (Bertin & Arnouts, 1996), which can deblend overlapping sources, but this source extraction software only produces an ellipse, and so fails to take into account the irregular nature of many of the structures we are attempting to recover. The structure of a galaxy is hierarchical and interconnected, and so computing a dendrogram of this structure is one way of identifying sources within the galaxy (Rosolowsky et al., 2008); dendrograms also have the additional benefit of extracting nested structure, which is vital in this study. In this work, we use the astrodendro dendrogram package777http://www.dendrograms.org. We refer the readers to the documentation on the astrodendro website for a more thorough description of the algorithm, but briefly a tree is constructed by arranging the pixels in order of flux. The first structure is centred on the brightest pixel, then the next brightest pixel is checked to see whether it should be considered a new structure or merged into another. The code moves down in flux until neighbouring regions touch, and if the difference between the maxima is significant, these ‘leaf’ structures are merged into a ‘branch’. The code works down to a minimum value and the structure is complete – a series of leaves connected to branches, with a ‘trunk’ at the bottom of each structure. These leaves are analagous to traditional sources, and it is these that we consider as our molecular clouds. For a visual comparison of these various algorithms on this data, see Appendix A. We also note that this data does not include kinematic information. This may lead to unrelated, but co-spatial along the line-of-sight, clouds becoming associated to one source when integrating along that line-of-sight. This is highlighted in Sect. 5.

As we wish to compute dust properties, we require sufficient data across the dust continuum peak, and into longer wavelengths where the bulk of the mass is contained. We therefore choose five wavebands across this peak, as a balance between spectral coverage and spatial resolution. These are the PACS 100 and 160µm  data, the SPIRE 250µm  data and the SCUBA-2 450 and 850µm  maps. We compute our dendrogram on the SPIRE 250µm data, as we found that after regridding and smoothing to the resolution and pixel scale of our lowest resolution data (the 250µm), that this map had the highest S/N. We select only regions with flux greater than 3$\sigma$ in each pixel, and regions must have a difference of greater than 3$\sigma$ to be considered significant and separate. This extraction criteria is selected to be as analogous as possible to Kirk et al. (2015), in order to make our results comparable to this earlier work. We also impose conditions that the region must be bigger than the SPIRE 250µm beam, and that none of these regions can touch the edge of the data. We find 165 leaves (i.e. no resolved substructure) in this dendrogram, which we assume to be our GMCs. The dendrogram for M33 can be seen in Fig. 3, and the positions of these clouds in Fig. 4. The majority of our analysis was performed on these clouds, although we also highlight the effect of performing this extraction on the SCUBA-2 450µm  map (our highest resolution data) in the size distribution (Sect. 4.1).

3.2 Flux Extraction and SED Fitting

Various parameters of these leaf nodes can be seen in Table 2. For each node, we list the mean position of the structure, and the deprojected distance from the centre of M33 (01h33m50.9s, +30$\degr$39′37*′′*; Plucinsky et al. 2008). We also calculate a FWHM of the cloud, based on its intensity-weighted second moment.

We have also computed fluxes in the PACS and SCUBA-2 bands for each of these clouds. astrodendro outputs a mask for each node, and we measured the flux beneath each mask in each waveband for all of the nodes. We estimated a local background from the median of the isocontour surrounding the mask, which we subtracted from the pixels before summation. An estimate of the local RMS error is given by the standard deviation of the pixels in this isocontour. The fluxes are also listed in Table 2 – for fluxes less than 3$\sigma$, we list the flux as an upper limit. The error listed in this table reflects only the RMS error; we have not included any calibration error in this value.

For the clouds with fluxes > $3\sigma$ in 3 or more bands (in our case, every cloud), we fit a single modified blackbody (MBB) of the form

[TABLE]

where $S_{\nu}$ is the flux at frequency $\nu$, $\kappa_{\nu}$ is dust absorption coefficient at frequency $\nu$, i.e.

[TABLE]

M${}_{\text{dust}}$ is the dust mass, $B(\nu,T_{\text{dust}})$ is the Planck function at frequency $\nu$ and dust temperature $T_{\text{dust}}$, and D is the distance to the source. We normalise $\kappa_{\nu}$ using the value calculated by Clark et al. (2016), $\kappa_{500}=0.051^{+0.070}_{-0.026}\,\mathrm{m}^{2}\,\mathrm{kg}^{-1}$. We note that this only holds true for the optically thin case, but as the theoretically expected value for when the optical depth becomes unity is 100µm (Draine, 2006), and experimentally only affects points shorter than 50µm (Casey, 2012), this is a reasonable assumption for our fits. To minimise the number of free parameters in this fit, we assumed a fixed $\beta$ of 1.5, which Tabatabaei et al. (2014) find to be a good fit to M33. We include correlated uncertainties in the PACS and SPIRE bands (as the SCUBA-2 450µm  data includes the SPIRE 500µm  map). This is implemented by employing the full covariance matrix. We performed our fitting within an MCMC framework using emcee888http://dfm.io/emcee/current/, and we quote the errors as the 84th percentile minus the 50th percentile, as we find that our errors are symmetric. Our initial guess for dust mass and temperature were set from a simple least-squares fit. We allowed the dust temperature to freely vary from 0 to 200K, and the dust mass from 0 to 1013 M⊙. An example SED fit is shown in Fig. 5. We also calculated the total infrared (TIR) luminosity of this cloud by integrating the MBB from 3-1100µm. We find that these clouds contribute around 50% of the total dust luminosity of M33, despite only occupying around 20% of the area. This indicates that these clouds are, in general, compact and bright in their dust emission. All of our derived SED parameters are given in Table 6. The dominant error in the dust mass and luminosity is error in $\kappa_{500}$ ($\sim$0.32 dex). As this is a systematic error we do not include this in Table 6. We do, however, include this uncertainty in our analysis.

We also include a measurement of the CO(J=2–1) luminosity (in K km s*-1*) in Table 6. Finally, we calculated Hi surface densities (in M⊙ pc*-2*) of each of our sources. A surface density is calculated, assuming (Rohlfs & Wilson, 1996)

[TABLE]

With this gas data, we performed the same procedure as for the FIR/sub-mm flux extraction – convolution and regridding to the same pixel scale, as well as local background subtraction. Similarly to the FIR/sub-mm fluxes, we list upper limits for intensities less than 3$\sigma$.

4 Cloud Properties

4.1 Size Distribution

For each source, we take the ellipse enclosing the cloud as computed by astrodendro from the half-width at half maximum (HWHM) of the second moments. We calculate a FWHM for each cloud from the average of these HWHM. The size distribution of the clouds can be seen as a Kernel Density Estimator (KDE) plot in Fig. 6. The median size of these clouds is 105 pc, close to the FWHM of the SPIRE 250µm beam, and so may initially be assumed to be complexes of smaller clouds. However, when performing the same extraction on our higher-resolution SCUBA-2 450µm  and 850µm  data, which has a minimum FWHM of 28 and 56 pc as defined in our extraction criteria, very similar trends are seen (Fig. 6). This would indicate that these objects are either (a) genuinely more extended than seen in the MW or (b) complexes of many very small clouds, rather than several larger clouds. Roman-Duval et al. (2010) find cloud sizes of 0.2 to 35pc, with a mean size of $\sim$8pc. More recently, Miville-Deschênes et al. (2017) find MW cloud sizes up to $\sim 500$pc, with a mean size of $\sim$30pc. Given these results from the MW, this would indicate these sources are likely complexes of smaller clouds. Additionally, comparisons to CO surveys (see Sect. 5) show that scenario (b) is more likely the case.

4.2 Dust Temperatures

One of the fitted parameters in the MBB is the dust temperature, and the distribution of this is shown in Fig 7. We find that the clouds have a median temperature of $23\pm 4$K, somewhat warmer than found for clouds in M31, which have a dust temperature of $18\pm 2$K (Kirk et al., 2015). We attribute this to the fact that M33 is much more actively star-forming than M31 (Heyer et al., 2004), and thus this dust is more strongly irradiated by these young stars. The distributions of size (Fig. 6) and dust temperature (Fig. 7) look somewhat similar. However, a calculation of the Kendall rank correlation coefficient (Kendall, 1938), where $\tau=+1$ indicates a perfect correlation, and $\tau=-1$ a perfect anti-correlation gives a weak anti-correlation of -0.17. A two-sample Kolmogorov-Smirnov test gives a p-value $\ll$1%, indicating that these two distributions are significantly different.

4.3 Cloud Masses

4.3.1 Calculating Masses

We simultaneously calculate a gas-to-dust ratio (GDR) and CO conversion factor ($\alpha_{\mathrm{CO}}$) in a fashion similar to that of Sandstrom et al. (2013). A dust mass surface density can be converted to a total gas mass surface density via

[TABLE]

Sandstrom et al. (2013) find the best fit of these two unknown parameters simultaneously by minimising the scatter in the log of the dust to gas ratio (DGR), and we perform this fitting using an MCMC analysis, accounting for errors in the dust mass surface density, Hi surface density and CO intensity.

We performed this fitting by grouping our clouds into bins of increasing galactocentric radius. We split our clouds into radial bins of 0.5 kpc and simultaneously fit $\alpha_{\mathrm{CO}}$ and GDR. The results of this can be seen in Fig. 8, and the slope for the GDR is given by

[TABLE]

Our maximum GDR is somewhat lower than seen in nearby galaxies (Sandstrom et al., 2013), with a value of around 90. However, variation in $\kappa_{\nu}$ can easily lead to huge variations in dust mass. Given that our adopted $\kappa_{\nu}$ is on the low end of literature values (Clark et al., 2016), this is not unexpected. However, we note that whilst our adopted $\kappa_{\nu}$ is lower than many other literature estimates, it is still compatible with the $\kappa_{\nu}$ of Draine & Li (2007), which Sandstrom et al. (2013) use in their work. Thus, this low value for the GDR cannot simply be attributed to our choice of $\kappa_{\nu}$. Using this calculated GDR, we transformed the dust mass into a total gas mass using Equation 7 and then calculated a total cloud mass (the sum of the dust and gas mass). The mass distribution can be seen in the left panel of Fig. 9.

We also calculated $\alpha_{\rm CO}$ simultaneously radially, and we show this in the bottom panel of Fig. 8. This is a CO conversion factor for the J=1–0 line, assuming CO(2–1)/CO(1–0) = 0.7 (Sandstrom et al., 2013). There is little radial variation, unlike the GDR, but the value is much higher than seen in other, nearby galaxies (Sandstrom et al., 2013). Even given variation in $\kappa_{\nu}$ that could decrease these values by a factor $\sim$2, this would indicate a CO conversion factor, $\alpha_{\rm CO}$, that is higher than seen in other, nearby galaxies. This is likely due to the subsolar metallicity of M33, with CO molecules becoming more easily dissociated by UV radiation at low metallicity (Glover & Clark, 2016).

We estimated the point source mass sensitivity by taking a limiting flux of 68.9mJy (a 3$\sigma$ point source as defined by our dendrogram extraction criteria), and sampled the GDR and dust temperature from distributions given by the distributions of our clouds (T = $23\pm 4$K, log(GDR) = $1.74\pm 0.09$). Bootstrapping this 10,000 times, we find a point-source sensitivity of $4.63^{+0.28}_{-0.22}\,\log_{10}(\mathrm{M}_{\odot})$. This is shown in Fig. 9, and cannot fully account for the deviation from a power-law at the low-end of the mass distribution.

We next estimated the completeness by injecting point sources of given cloud mass into a fake map with the same noise properties as the SPIRE 250µm  data, and a background similar to that of M33. We sample the dust temperature and GDR as with the point source sensitivity, and inject 100 sources of each mass into this map. We performed the same extraction criteria as we did with the real data and calculated the completeness for each mass. We find that we are 95% complete above a mass of 105.5 M⊙. This means that the observed downturn is simply due to incompleteness, and is not a genuine turnover. However, we must stress that this is only an approximation of the true completeness limit. We have here assumed only point sources present in a constant background, but given that these sources are extended, and embedded in a complex background, the true completeness limit will be a function of mass, radius, cloud shape and position within the map. Accounting for this complex completeness is beyond the scope of this work.

4.3.2 Power-Law Fitting

We fit a power-law of the form $N$ (M) $\propto$ M${}^{\alpha_{\rm M}}$ to the high-end of the mass distribution. However, in a standard distribution the fit can become biased by small number statistics at high-mass (Maschberger & Kroupa, 2009), and so it is more reliable to fit to the cumulative mass distribution (shown in the right panel of Fig. 9). In this case, the power-law takes the form $N(>{\rm M})\propto{\rm M}^{\alpha_{\rm M}-1}$. To avoid incompleteness, we fit only to values with a cloud mass greater than $10^{5.5}$M⊙. We find a value of $\alpha_{\rm M}$ of $-2.83^{+0.24}_{-0.15}$, steeper than the value of $\alpha_{\rm M}=-2.0\pm 0.1$ found previously in M33 by Gratier et al. (2012) using CO(J=2–1), and $-2.6\pm 0.3$ from the CO(J=1–0) work of Engargiola et al. (2003). Work by Bigiel et al. (2010) has hinted at a steeper slope in the outskirts of M33, and our calculated slope appears to confirm this. We also find that this value is steeper than molecular clouds in the MW, which has an exponent of around -1.5 (e.g. Sanders et al. 1985; Solomon et al. 1979). The slope is also steeper than that found in M31 ($-2.34\pm 0.21$, Kirk et al. 2015, $-2.55\pm 0.2$, Blitz et al. 2007). The steepness of this slope appears to indicate that M33 is more dominated by smaller clouds than in, e.g., the MW. Given that Gratier et al. (2012) use the CO luminosity as a proxy for molecular hydrogen, whilst the dust content should be an independent tracer of total gas content, we can rule out this steep slope being due to a lower CO intensity per H2. It would appear that M33 is intrinsically poorer at cloud assembly than other local spirals.

The efficiency of cloud assembly has been linked to a variety of processes. The amplitude of the spiral density wave can have an effect on the GMC population (e.g. Shu et al. 1972). However, given that recent modelling work has shown that the spiral arms of M33 are most likely driven by gravitational instabilities (Dobbs et al., 2018), we believe this is unlikely to be the case. The interstellar pressure of gas (Elmegreen & Parravano, 1994; Blitz & Rosolowsky, 2006) may also be a factor. However, we see higher interstellar pressure in M33 than in the MW (Kasparova & Zasov, 2008), so given this hypothesis we would expect more massive clouds. We can therefore rule out the interstellar pressure as the main driver of this inefficient cloud formation. Alternatively, metallicity can play a role in the conversion of Hi to H2 (Krumholz et al., 2008). Given the subsolar metallicity of M33, we would expect this conversion to be less efficient, and therefore cloud formation similarly inefficient. Finally, it is believed that H2 can form from merging Hi clouds (e.g. Heitsch et al. 2005), so we may expect from larger Hi velocity dispersions, more massive clouds may form. The average Hi velocity dispersion in M33 is of the order 13 km s*-1*, with little radial variation (Corbelli et al., 2018), whilst the outer MW shows much more turbulent Hi gas, with velocity dispersions of 74 km s*-1* (Kalberla & Dedes, 2008). Our results are unable to distinguish which of these two mechanisms are the main driving force behind this inefficient cloud formation, but it is clear that the cloud mass distribution is significantly different in M33 than the other massive spirals in our Local Group. The exact cause of this is currently unclear, but high-resolution surveys of many galaxies with a wide range of properties will be able to explain the diversity in cloud populations seen even between the galaxies of our Local Group.

4.4 Radial Variation in Cloud Properties

We also investigated any radial variation in the cloud properties. Fig. 10 shows the number density (the number of clouds per annular area) of these clouds with galactocentric radius. We see that up to a radius of 2.5 kpc, our cloud distribution and the GMC distribution of Gratier et al. (2012) agree very well – however, after $\sim$3 kpc, the distribution of GMCs from Gratier et al. (2012) is systematically lower. We believe that this is due to the fact that our data is wider than the CO map on which they perform their extraction. As they perform this analysis on an incomplete map of M33 (only the area covered by the Herschel PACS and Heterodyne Instrument for the Far Infrared (HIFI) spectrometers), they do not map the entire disk of M33. We would suggest our distribution is therefore less biased, and gives a more representative view of the GMC number density. Along with a peak in the distribution at the centre of the galaxy, there is a step in this distribution from around 2 kpc to 4 kpc radius, which corresponds to the positions of the spiral arms in M33. However, the spiral arms are less pronounced in this distribution than M31 (Kirk et al., 2015), where the positions of the spiral arms have clear peaks, and the SFR in the rest of M31 is very low.

We also investigated the radial variation in our two fitted MBB properties – dust temperature and dust mass (Fig. 11. In both cases, we see that the radial correlations are weak – in the case of dust temperature, weakly negative (i.e. dust temperatures are lower at higher galactocentric radii), and weakly positive in the case of dust mass. We find a dust temperature gradient $-0.71\pm 0.01$ K kpc*-1*, and a dust mass gradient of $0.053\pm 0.001$ dex kpc*-1*. The decrease in dust temperature is naturally explained by a general decrease in the strength of the interstellar radiation field (ISRF) at increasing galactocentric radius (Rice et al., 1990). This gradient is also similar to that seen by Tabatabaei et al. (2014), when considering the global properties of M33. The invariance in dust mass is likely due to a balance of generally more compact but brighter clouds in the centre of the galaxy, whereas in the outskirts we tend to find somewhat more diffuse (but extended) sources.

5 Comparison to CO

Using kinematic CO data, we can relate the properties of these clouds to the observed scaling relations of Larson (1981). This can be neatly represented in the plane of the cloud surface density ($\Sigma_{\rm cloud}$) with the ratio of the velocity dispersion ($\sigma_{\rm v}$) to the square root of the cloud radius (R), as demonstrated in Heyer et al. (2009). If these clouds are ideally virialised (Larson’s second law), and the velocity dispersion is related to the cloud’s radius as $\sigma_{\rm v}\propto R^{1/2}$ (Larson’s first law), then it can be shown

[TABLE]

Larson’s third law states that $\Sigma_{\rm cloud}$ is approximately equal for any cloud, so we would expect little dynamic range in this quantity. We calculate the velocity dispersion from the CO data cube of Gratier et al. (2012), following the “equivalent width” as defined in Heyer et al. (2001):

[TABLE]

with $I_{\rm CO}$ the integrated CO in each leaf contour, and $T_{\rm peak}$ the peak line intensity. This relationship is shown in Fig. 12. We find that there is a weak correlation ($\tau$ = 0.21) between $\sigma_{\rm v}/R^{0.5}$ and $\Sigma_{\rm cloud}$ for our clouds. This is somewhat weaker than that found by Heyer et al. (2009) for a selection of clouds in the MW, but still showing a dependence of $\sigma_{\rm v}/R^{0.5}$ on the gas surface density. Given the proximity of these clouds to the line of virial equilibrium (the median deviation below this line is a factor of 1.5), we conclude that these clouds are likely ideally virialised.

Finally, we make comparison to the locations of GMCs identified with earlier CO studies. Fig. 13 shows the positions of our clouds against those of Gratier et al. (2010) and Engargiola et al. (2003). These earlier studies have somewhat better resolution than we have achieved in this investiation (18 arcsec versus 12 arcsec), but clearly the cloud distribution is broadly similar, indicating that these particularly dusty regions are indeed associated with GMCs. However, what often appears to be a single region, even in the 450µm  data, is identified as several clouds in the CO data. Fig. 13 shows that at the 450µm  resolution, larger sources are beginning to break up into smaller clouds, but this is not statistically significant enough for the extraction criteria to define them as seperate leaves. Some of these sources may also be co-spatial along the line-of-sight of the galaxy, but given no kinematic information we cannot separate these as in the CO surveys.

One notable difference between our detected sources and earlier CO surveys is that the works of Engargiola et al. (2003) and Gratier et al. (2012) find a dearth of massive clouds beyond a galactocentric radius of 4kpc. However, we find a nearly flat distribution of dust mass with galactocentric radius. Given that we find our clouds to be co-spatial to these earlier works in the inner region of M33, we would expect these earlier surveys to detect these clouds. We do not believe this is a selection effect due to noise in these CO maps. Gratier et al. (2012) map an area significantly beyond 4kpc, with similar noise as in the centre of the map (see Fig. 3 of Gratier et al. 2010). Engargiola et al. (2003) estimate their map to be complete out to 5.2kpc and more than 50% complete up to 8kpc, so this variation cannot simply be attributed to completeness. Results from Planck (Planck Collaboration et al., 2011) have shown a significant reservoir of molecular hydrogen that is not traced by CO. Gratier et al. (2017) find that this “CO-dark” gas forms around 50% of the total molecular hydrogen mass of M33. The amount of CO dark gas is also expected to increase at lower metallicity (such as in the outskirts of M33), where the CO is more susceptible to photo-dissociation. Given that the dust continuum is not subject to these same caveats, the dust may offer a more representative view of the cloud population than CO surveys in these lower-metallicity environments.

6 Conclusions

In this work, we have combined archival Herschel FIR and sub-mm data with deep SCUBA-2 observations to probe the properties of GMCs using their dust content. Using wavelengths from 100 to 850µm, we have probed the cold dust continuum emission of these sources. We performed source extraction using dendrograms, which found a total of 165 GMCs with sizes (FWHM) of 46-280 pc, and a median size of 105 pc. By fitting a one-temperature MBB, we have calculated the dust mass and temperature for these 165 sources, and compared this to archival CO and Hi data. Using a method similar to that of Sandstrom et al. (2013), we find a weak radial variation in the GDR of these sources, and use this GDR to calculate a total cloud mass.

These cloud masses span the range of $10^{4}-10^{7}{\rm M}_{\odot}$, and the mass function can be fit with a power law slope proportional to ${\rm M}^{-2.84}$, steeper than seen in previous CO studies of M33 and the MW. Whilst we can rule out pressure as the major driver of this inefficient cloud assembly, we are unable to distinguish whether metallicity or turbulent Hi velocities contribute more to this inefficiency. The dust temperatures of these clouds range from 17-32 K, and dust masses from 102-10${}^{5}\,{\rm M}_{\odot}$. In terms of these clouds’ dust properties, we find only weak radial trends with dust mass and dust temperature.

A comparison to CO data shows an $\alpha_{\rm CO}$ factor several times higher than found in nearby galaxies. We attribute this to the subsolar metallicity of M33, where CO is likely a less suitable tracer of molecular hydrogen. We have examined these clouds in the framework of Larson’s scaling relations, and we find a dependence of $\sigma_{\rm v}/R^{0.5}$ with the gas surface density, much like that of Heyer et al. (2009). It would also appear that our clouds are ideally virialised (to within a median factor of 1.5). Finally, a comparison with earlier CO studies shows that the GMCs we are detecting and those found using CO data are generally co-spatial, but due to the limited resolution of the SPIRE 250µm data that we convolve and regrid to, the crowded complexes of clouds seen in this CO data is generally confused into one large cloud in our source extraction. We also find clouds beyond 4kpc in galactocentric radius, unlike earlier CO surveys. This may be due to the CO being a poorer tracer of molecular hydrogen at these larger distances (and lower metallicities), and these clouds being dominated by CO-dark gas.

Acknowledgements

M.W.L.S acknowledges support from the European Research Council (ERC) Forward Progress 7 (FP7) project HELP. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observatories of China and the Chinese Academy of Sciences (Grant No. XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada.

This research made use of montage (http://montage.ipac.caltech.edu/), which is funded by the National Science Foundation under Grant Number ACI-1440620, and was previously funded by the National Aeronautics and Space Administration’s Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology.

This research has made use of Astropy, a community-developed core Python package for Astronomy (http://www.astropy.org/; Astropy Collaboration et al. 2013, 2018). This research has made use of NumPy (http://www.numpy.org/; van der Walt et al. 2011), SciPy (http://www.scipy.org/), and MatPlotLib (http://matplotlib.org/; Hunter 2007). This research made use of APLpy, an open-source plotting package for Python (https://aplpy.github.io/; Robitaille & Bressert 2012).

Appendix A Source Extraction Comparison

Here, we present a brief overview of the various source extraction algorithms used in our initial testing – fellwalker, SExtractor, and astrodendro. We have attempted to homogenise the extraction criteria to make testing these algorithms as fair as possible. In all cases, we detect objects only if they are 3$\sigma$ above the background, with an area larger than the beam. For fellwalker and astrodendro, we also set a minimum significance for the structure to be 3$\sigma$ (else the peaks will be merged into a single peak). For SExtractor, we turn off the deblending threshold. The results of these various algorithms are shown in Fig. 14. There is good correspondence between the three algorithms, and each detect a similar number of sources (169 for fellwalker, 165 for astrodendro, and 188 for SExtractor). However, we can see that fellwalker essentially partitions all of the emission in the image, leading to clearly unreasonably large structures. astrodendro, however, finds more compact sources of emission. SExtractor can deblend some of the most crowded regions, resolving single sources in the other algorithms into several smaller sources, but produces an ellipse of emission, rather than an exact contour. Given that SExtractor can have overlapping ellipses (which is always the case in these crowded areas), whereas astrodendro separates sources by default, we have opted to use astrodendro in this work.

Appendix B Leaf Node Parameters

Appendix C SED Parameters