An ALMA archival study of the clump mass function in the Large Magellanic Cloud

ABSTRACT

1 INTRODUCTION

Millimetre observations beyond the Milky Way are beginning to probe the scales of clustered and individual star-forming regions known as molecular gas clumps and cores. While filamentary structure has previously been shown to be ubiquitous in Galactic star-forming molecular gas (André et al. 2010; Arzoumanian et al. 2011; Könyves et al. 2015) and in other galaxies (the Large Magellanic Cloud in carbon monoxide (CO): Fukui et al. 2008; Wong et al. 2011), it is the higher-density clumps and cores within the filaments that directly set the initial conditions for the formation of stars and clusters.

The choice of these functional forms is motivated by the tendency for the mass functions to have a long tail at high masses and a flatter portion or turnover at low masses. Physically, the power-law tail has been attributed to collapse in molecular clouds in whichgravity is dominating the motions (Ballesteros-Paredes et al. 2011), although analysis of simulated and analytic cloud column-densities has shown the power-law tail occurs as a result of the appearance of strong density peaks in molecular clouds, not the transition to gravity-dominated dynamics (Tassis et al. 2010). Lognormal forms have been argued for as following from the density distribution of supersonic turbulence (Padoan, Nordlund & Jones 1997), although the presence of supersonic turbulence is not necessary for the density distribution to exhibit a lognormal shape (Tassis et al. 2010). The combination of various independent stochastic processes impacting the distribution of molecular gas in clouds has been shown to lead to a lognormal distribution through the central limit theorem (Larson 1973; Adams & Fatuzzo 1996). However, observational evidence for the ubiquity of lognormal distributions may be more difficult to confirm than previously thought (Alves, Lombardi & Lada 2017).

In addition to comparing mass distributions between different regions, core mass functions have also been compared to the similarly shaped stellar IMF. Chabrier (2003, 2005) showed that the stellar IMF follows the original Salpeter (1955) power-law index of −2.35 ± 0.3 above 1 |$\mathrm{M}_{\odot}$|⁠, with a lognormal turnover at low masses peaking at about 0.2 |$\mathrm{M}_{\odot}$|⁠. This mass distribution has been observed to be universal across environments such as the Galactic disc, young and globular star clusters, and the spheroid or stellar halo (Chabrier 2003, 2005; Krumholz 2014). A similar power-law slope has been predicted for the stellar IMF and core mass functions (Chabrier & Hennebelle 2010; Guszejnov & Hopkins 2015). The observed similarity in slopes has led to the idea that there exists a star formation efficiency that is independent of mass and is acting to transform the core mass function to the IMF (Alves, Lombardi & Lada 2007; André et al. 2010; Könyves et al. 2015).

At higher masses, the comparison is between the mass distributions of entire star clusters and the clump mass functions within their natal GMCs (Lada 1992). Lada & Lada (2003) compiled a list of 76 young star clusters still embedded within a GMC and from this derived a cluster mass function with a power-law index of −2 between ∼50 and 1000 |$\mathrm{M}_{\odot}$|⁠. Fall & Chandar (2012) found that the mass functions of star clusters from six galaxies were all well fitted by power laws with indices close to −1.9.

Given the wealth of information gathered on the mass distributions of cores and clumps in the Milky Way, the next step is to extend these studies to other galaxies. This provides new physical environments in which to test models, as well as larger samples from which to draw statistical conclusions. Two of our nearest neighbours, the Magellanic Clouds, are now well within the reach of mm/submm observations of molecular and dust emission for studies of the clump mass function. For example, Indebetouw et al. (2013) derived a CO mass function of 103 clumps at 0.46-pc resolution in the Large Magellanic Cloud (LMC) with the Atacama Large Millimeter/Submillimeter Array (ALMA). While extension beyond the Magellanic Clouds promises a further variety of environments and larger sample sizes, such studies push even the most advanced observatories to their limits. Rubio et al. (2015) reported on 10 CO clouds in the Wolf–Lundmark–Melotte (WLM) galaxy with masses of 5900 to 7.3 × 10⁴ |$\mathrm{M}_{\odot}$|⁠. At a distance of 1 Mpc, they achieved ∼5-pc resolution with ALMA. More recently, Schruba et al. (2017) observed NGC 6822 with ALMA, reaching 2-pc resolution at a distance of 470 kpc. A total of 156 CO clumps were extracted with masses of 9 to 3500 |$\mathrm{M}_{\odot}$|⁠. Thus, if we wish to sample the full range of clump sizes and even to start to probe molecular core scales, we are limited to the nearest Local Group members.

The LMC is the ideal next step after our Galaxy for the study of resolved star formation, given its proximity (49.97 ± 1.11 kpc, Pietrzyński et al. 2013) and nearly face-on orientation. It is also a system that is significantly different from the Milky Way in which to study how stars form, with a lower average metallicity of ∼(1/3)–(1/4) |$\mathrm{Z}_{\odot}$| (Rolleston, Trundle & Dufton 2002; Dufour 1984). The lower metal abundance results in smaller quantities of dust and thus in less shielding for star-forming regions from ultraviolet (UV) radiation. This means that molecular gas reservoirs may be lower in mass (Fukui et al. 1999, 2001), which limits the star-forming fuel throughout the galaxy. The presence of fewer metals also means that cooling through line emission is reduced, which can change the energy balance in molecular clouds. Finally, the characteristic peak mass of the stellar IMF has been predicted to shift to higher masses for low metallicity (Bromm 2005, albeit for near-zero metallicity conditions).

In this paper, we calculate and analyse the dust-derived clump mass function in the LMC in order to facilitate comparison with mass functions of star-forming regions in the Milky Way. In Section 2 we summarize the observations from each ALMA project we used and we describe our continuum map-making process. In Section 3 we describe our method of isolating the dust-only emission in each region, our clump identification procedure, and the steps we took to fit and characterize our final mass function. In Section 4 we discuss our results from Section 3 and place them in the broader context of the study of molecular clump mass functions. The paper is summarized in Section 5. Maps of each field included in our mass function are presented in Appendix  A as observed emission as well as being decomposed into dust-only and free–free-only emission.

2 OBSERVATIONS AND DATA REDUCTION

2.1 Spatial and spectral setups

We retrieved six publicly available projects from the ALMA archive containing observations of seven fields in the LMC covering the star-forming regions 30 Doradus-10 (30 Dor-10), N159W, N159E, N113, N166, GMC 225 and PCC 11546. In addition, we obtained new data in Cycle 4 towards 30 Dor-10 with the 7-m Atacama Compact Array (ACA). The projects were observed between 2011 December 31 and 2017 April 21, spanning Cycles 0, 1, 2 and 4.

For all but one region, multiple pointings were observed across each field in both bands to produce mosaic maps. Mosaic pointings were roughly Nyquist-spaced in a given band in order to achieve relatively uniform coverage across the inner portion of the maps when imaged together. 12-m (main array) plus ACA data were obtained for all fields except N113, which was observed as a single pointing and only with the main array. The number of pointings on a field and in a single band ranged from 1 (N113) to 170 (N166). Mapped areas cover between 0.21 arcmin² (N113) and 7.5 arcmin² (N166; refer to Table 3 for the areas of each field in both bands).

In order to measure dust masses in molecular clumps we focused on continuum observations from 86  to 100 GHz (∼3.2 mm) and from 217  to 233 GHz (∼1.3 mm) in ALMA bands 3 and 6, respectively. Observations at these frequencies are likely to contain significant emission from sources other than thermal dust, such as free–free emission. The effective bandwidths used in making the continuum maps ranged between 1.6  and 4 GHz. The locations of the regions analysed here are marked on the infrared image of the LMC in Fig. 1, and a summary of the observational details is given in Table 1.

Figure 1.

SAGE Spitzer IRAC three-colour image of the Large Magellanic Cloud (Meixner et al. 2006). Red, green and blue are 8, 4.5 and 3.6 |$\mu$|m maps, respectively. The positions of the seven regions discussed in this paper are labelled.

Open in new tabDownload slide

2.2 Background on individual fields

The Tarantula Nebula, or 30 Doradus, is one of the most active star-forming regions in the LMC. It harbours the R136 star cluster, which boasts stellar densities of between 10⁴ and 10⁷ |$\mathrm{M}_{\odot}$| pc⁻³ (Selman & Melnick 2013). Multiple generations of star formation have occurred in 30 Doradus over the course of ∼20 Myr (De Marchi et al. 2011; Walborn & Blades 1987). 30 Doradus has been observed as part of several LMC-wide surveys targeting various emission sources. In CO, it has been observed with theNANTEN telescope in ¹²CO (1–0) (Fukui et al. 2008) and with the Mopra telescope in ¹²CO and ¹³CO (1–0) (Hughes et al. 2010; Wong et al. 2011). The HERITAGE Key Project Survey (Meixner et al. 2013) observed it in dust emission along with the rest of the LMC.

The GMC 30 Dor-10 (Johansson et al. 1998) is part of the Tarantula Nebula, and is within a projected distance of 11 pc from R136. Indebetouw et al. (2013) observed 30 Dor-10 with ALMA in ¹²CO and ¹³CO (2–1) emission as well as the 1.3-mm dust continuum at ∼1.9-arcsec resolution. The dust map was used to derive a total H₂ mass for the GMC of (6 ± 1) × 10⁴ |$\mathrm{M}_{\odot}$|⁠. Clump apertures identified in their ¹²CO and ¹³CO cubes were used to calculate independent H₂ masses from both the dust and the CO data. Indebetouw et al. (2013) calculated a CO-derived mass function for their 103 clumps and fitted it with a power law of α = −1.9 ± 0.2.

N159 was originally identified by Henize (1956) as an H ii region and has since been extensively studied. Noted as the strongest CO-intensity cloud in the initial NANTEN LMC survey (Fukui et al. 2008), N159 has been resolved into three major clumps, namely N159W, N159E and N159S (Johansson et al. 1994). Across the region, multiline excitation analyses on ∼10-pc scales estimated molecular gas kinetic temperatures ranging from 10 to 25 K (Heikkilä, Johansson & Olofsson 1999; Bolatto et al. 2000). ALMA observations of N159W suggest that its constituent group of compact H ii regions may have formed through a cloud–cloud collision (Fukui et al. 2015). By adding a dynamical estimate of one of the young stellar objects (YSOs) of 10⁴ yr, Fukui et al. (2015) showed that the collisional triggering of star formation probably occurred very recently. N159E contains multiple developed H ii regions, the most prominent of which is the Papillon Nebula (Mizuno et al. 2010). Saigo et al. (2017) suggested that a three-cloud collision is occurring in N159E, with the Papillon Nebula protostar in the overlap region. They also observed a molecular hole around the protostar filled with 98-GHz free–free emission, indicating that the protostar has recently begun disrupting the molecular cloud from which it formed. N159S does not exhibit any current star formation.

N113 was also identified by Henize (1956) as an H ii region. Observed as part of a sample of four molecular clouds (N159W, N113, N44BC and N214DE) with the Swedish-ESO Submillimetre Telescope (SEST), it exhibited significantly lower gas-phase C¹⁸O/C¹⁷O abundances compared with molecular clouds in the Milky Way and the centres of starbursts, with a mean across clouds of 1.6 ± 0.3 (Heikkila, Johansson & Olofsson 1998). Subsequent observations have shown it to contain the most intense maser in the Magellanic Clouds (Imai et al. 2013), clumpy molecular gas currently forming stars (Seale et al. 2012), and a host of YSOs as identified by Herschel and Spitzer (Ward et al. 2016). N113 is related to three young stellar clusters (NGC 1874, NGC 1876 and NGC 1877; Bica, Claria & Dottori 1992). It also contains a rich assortment of molecular species (Paron et al. 2014).

N166 is yet another region originally identified by Henize (1956), and Fukui et al. (2008) determined that it contained only H ii regions. Two molecular clouds were observed overlapping its position in the second NANTEN survey, and a follow-up observation with ALMA targeted a position between these two NANTEN-observed clouds. N166 was also observed with the Atacama Submillimetre Telescope Experiment (ASTE) in ¹²CO (3–2), revealing densities of 10² to 10³ cm⁻³ and kinetic temperatures between 25 and 150 K in a 22-arcsec beam (Paron et al. 2016).

GMC 225 is one of 272 GMCs identified by Fukui et al. (2008) in the second NANTEN survey of the LMC in ¹²CO (1–0) and shows no sign of massive star formation (Kawamura et al. 2009). With a CO-derived total mass of 10⁶ |$\mathrm{M}_{\odot}$| and a radius of ∼73 pc, it was observed as part of a follow-up ALMA project. It is ∼1500 pc south of 30 Dor-10 and east of the molecular ridge.

PCC 11546 is an extremely cold (≲15 K) dust source identified in the southern limits of the LMC as part of the Planck Galactic Cold Cloud catalogue (Planck Collaboration et al. 2016). It also exhibits strong CO (1–0) emission in the Planck integrated CO map and the MAGMA LMC CO survey (Wong et al. 2011). There appears to be a lack of massive star formation within PCC 11546, and it contains lower-density gas than clouds closer to the centre of the LMC (Wong et al. 2017).

2.3 Calibration

The visibility data were retrieved from the archive and calibrated. All fields except 30 Dor-10 had only raw data available in the archive. 30 Dor-10 had calibrated visibilities available, and we used those data after recalculating the weights based on the scatter in the visibilities using the Common Astronomy Software Applications (casa) ‘statwt’ command.

The observations with raw data available were either ‘manually’ or pipeline-calibrated at the observatory using casa (McMullin et al. 2007). For manually calibrated data, we ran the full calibration procedure using the observatory scripts in the latest version of casa available at the time, 4.7.2-REL (r39762). Minor editing of the calibration scripts was necessary in order to account for task and parameter changes from the older versions of casa. Edits were also made to ensure that the visibility weights were properly calculated throughout the entire calibration procedure.

For pipeline-calibrated data, we used the casa and pipeline versions closest to those used by the observatory for the original calibration.¹ There was no concern regarding visibility weights because all publicly available pipeline releases were after casa 4.2.2, which contained a major correction to the handling of the weights.

Once the calibrated data were obtained we inspected the visibilities from the calibrator sources to ensure that the results of the calibration were as expected and that all seriously problematic data were flagged. In most cases this inspection did not reveal anything that needed to be done beyond the observatory-provided calibration process. A few cases, however, did expose situations where marginal antennas or edge channels should have been flagged, and we did so before imaging.

2.4 Imaging

All imaging steps, from inspecting the visibilities up to and including cleaning the maps, were carried out in casa 4.7.2² (r39762) for all fields. Strong emission lines were identified with the plotms task and flagged before the visibilities were imaged to produce continuum maps. Because we combined maps in the two bands to create a dust-only emission map, we needed to match the spatial scales to which each pair of maps is sensitive as closely as possible. Matching the spatial-scale sensitivities was achieved through a combination of trimming the shortest uv-spacings and tapering the weighting of the longest uv-spacings. In order to choose which inner uv-trimming to use (the minimum uv-distance to include in imaging), we plotted the source visibility amplitudes versus their radial uv-distance in wavelengths for each band. The longer minimum baseline length of the two was used as the minimum uv-distance to include. For our fields, this minimum uv-distance was always set by the Band 6 data, so the uv-distance trimming was always applied to the Band 3 data.

The outer uv-tapering directly changes the synthesized beam and sets the smallest spatial scale in the cleaned maps. The Band 6 maps always had smaller initial synthesized beams than the Band 3 maps. Once the uv-distance trimming was applied to the Band 3 data, we used the Band 3 synthesized beam as the target beam shape for the Band 6 uv-tapering.

Matching the beam shapes precisely through tapering alone was not possible given the differences in intrinsic uv-coverage between the observations. Therefore, we used the ‘restoringbeam’ parameter in clean to force the Band 6 Gaussian beam shape to be the same as that of the Band 3 map. Note that this means that all noise and any emission left in the residual maps is still at the tapered resolution, and only the emission that was cleaned is exactly matched to the Band 3 resolution. However, we were always able to bring the beams into fairly close agreement with the uv-tapering. The differences between the Band 6 and Band 3 beam axes were <10 per cent (and usually much less), and position angle differences were <5 per cent. Only 30 Dor-10 had a larger position angle difference (∼30° versus ∼60°).

We cleaned all the dirty maps that showed obvious emission, which meant that N166, GMC 225 and PCC 11546 were not cleaned. Motivated by the complexity of the emission in the 30 Dor-10 maps and by the desire to have a reproducible method of producing cleaned maps, we implemented an auto-masking algorithm for all cleaning. This algorithm is heavily based on the auto-masking code given in the M100 casa Guide as it appeared in 2016 August. We modified it to work as an automated casa script:³ when given a set of clean parameters, a minimum threshold to clean down to, a casa region text file specifying an emission-free region of the dirty map, and a minimum spatial size to mask, the script iteratively cleaned down to the desired threshold. All auto-masking was carried out with stopping thresholds between 1.5 and 3 times the noise in the map and using minimum mask areas of 0.5 times the map synthesized beam. All fields were cleaned with Briggs weighting (Briggs 1995),⁴ with the ‘robust’ parameter set to 0.5, and in multifrequency synthesis mode, with ‘psfmode’ set to ‘psfclark’ and using a maximum of 10⁴ iterations. Table 2 lists the values used for important clean parameters for each field. The average synthesized beam size was ∼2.3 arcsec, corresponding to ∼0.6 pc at the distance of the LMC. Mapped areas cover between 44 (N113) and 1600 pc².

2.5 Fields without detections

We tried to improve the signal-to-noise ratio (S/N) in the maps with no continuum emission by setting ‘robust’ to 2.0 to favour S/N over resolution. While the noise did drop, it was not enough to enable emission to be detected in any of the non-detection maps. The rms noise and beam shapes listed in Table 3 for GMC 225, N166 and PCC 11546 refer to the maps made with ‘robust’ set to 2.0.

Almost half of the fields we investigated were not detected in continuum emission. This raises the question whether the data for these fields are inherently of poorer sensitivity than the data with abundant emission. We tested this by matching the synthesized beams across all fields to the largest beam in each band. At 3.2 mm, GMC 225 and N166 have the lowest rms noise in our sample by a factor of 2–6. The 1.3-mm observations of PCC 11546 have a smaller rms noise than the fields with detections, and the sensitivities of GMC 225 and N166 at this wavelength are only a factor of 2–3 worse. This comparison suggests that there are intrinsic differences in the physical densities and structures of the non-detection fields compared with those of the four detected regions.

3 DATA ANALYSIS

3.1 Free–free correction

Given the frequencies of these observations, we expect the main contributions to continuum emission to be thermal emission from the ∼30 K dust and free–free emission from ionized gas. Although synchrotron emission is another possible source of continuum emission at these frequencies, H ii regions, protostars and young stars do not produce much synchrotron emission (Ginsburg et al. 2016). Colliding-wind binaries have been observed with spectral indices indicative of synchrotron radiation from particles accelerating in the wind collision zone (De Becker & Raucq 2013). However, following the method of Ginsburg et al. (2016), we estimate a 3.2-mm flux density of ∼6 µJy per binary. These flux densities are 30 to 100 times smaller than the noise in our fields (Table 3).

Fig. 2 compares the 1.3- and 3.2-mm fluxes for clumps identified in the 1.3-mm map. The solid line shows the flux scaling for free–free emission with exponent −0.1, and the two dashed lines show dust emission with exponents (2 + β) for β = 1 and β = 2. Because the points are clustered between the free–free and dust scaling relationships, we conclude that there is significant free–free emission even at 1.3 mm. In order to use the dust emission to estimate the total mass of the clumps in these fields we need to correct for contamination from free–free emission at 1.3 mm.

Figure 2.

Comparison of 3.2-mm (Band 3) and 1.3-mm (Band 6) flux densities for clumps identified at 1.3 mm. The dashed lines show two different thermal dust emission scaling relationships covering the range of realistic dust emissivity spectral indices, and the solid line shows a typical free–free scaling relationship. Upper limits at 3.2 mm are for clumps found at 1.3 mm without significant co-spatial emission at 3.2 mm. Error bars correspond to a 10 per cent calibration uncertainty.

Open in new tabDownload slide

Figure 3.

Dust-only maps of (a) 30 Dor-10, (b) N159W and (c) N159E made from the combination of 3.2-  and 1.3-mm ALMA observations. Contamination from free-free emission has been removed from these maps using equation (6).

Open in new tabDownload slide

We tested our algebraic decomposition using lower-frequency radio observations that are often used to measure free–free and synchrotron emission directly. Observations at 4.8 GHz (6 cm) and 8.4 GHz (3 cm) were made of 30 Doradus (Lazendic, Dickel & Jones 2003) and N159 (Indebetouw, Johnson & Conti 2004) with the Australia Telescope Compact Array (ATCA). We applied spatial filtering to the radio images, to match the ACA uv-coverage of the ALMA observations, through the casa observation simulation functionality. This allowed us to better match the spatial scales measured at the two frequencies (especially for 30 Dor-10, where the ACA uv-coverage filtered out about 50 per cent of the radio flux).

We first measured the total fluxes in each field at 4.8  and 8.4 GHz. Adopting a synchrotron spectral index of −0.83 and a free–free spectral index of −0.1, we iterated in order to estimate the balance of the two emission sources at the two frequencies. We then extrapolated the free–free component to the ALMA Band 3 frequency to obtain an estimate of the expected free–free emission at 3.2 mm.

To compare this estimate with observations, we solved for the free–free-only emission in Band 3 (see Appendix  A for specific decomposition expressions and maps of free–free-only emission). Assuming a 10 per cent absolute flux calibration uncertainty on the radio continuum, the total free–free emission in 30 Dor-10 matches within 15 per cent between extrapolated and decomposed estimates (0.82 and 0.97 Jy, respectively). N159W is consistent within the radio uncertainties, with flux densities of 75  and 78.5 mJy in the ATCA and ALMA maps, respectively.

N159E appears to have about four times too much free–free flux in Band 3 when compared with the extrapolated flux (it has a near-synchrotron-only spectral index in the radio maps). If the free–free emission becomes optically thick between 8.6  and 98 GHz then the free–free spectrum would turn over (e.g. Turner, Ho & Beck 1998), leaving just the synchrotron spectral index and a low free–free flux at low frequencies. Galactic ultracompact H ii regions have been observed with free–free spectral turnovers at frequencies between 10 and 20 GHz (Wood & Churchwell 1989). It is possible that the main region of bright radio emission is originating from a population of ultracompact H ii regions averaged together with supernova remnants within the ALMA beam, leading to synchrotron-dominated emission near 10 GHz and the presence of free–free emission at and above 100 GHz.

Finally, we estimate that diffuse synchrotron emission contributes up to ∼12 per cent of the Band 3 flux. Because our method scales between Band 3 and Band 6 with the free–free spectral index, this synchrotron component is oversubtracted from the Band 6 flux by equation (6). This results in, at most, 4 per cent of the total Band 6 flux being oversubtracted when making the dust-only maps as a result of neglecting a synchrotron term in equation (6). Considering the imperfect uv-coverage matching between the ATCA and ALMA observations and the 10 per cent calibration uncertainty on the ALMA and ATCA observations, we believe that the radio free–free emission estimates agree well enough and that the synchrotron is weak enough in Band 6 that we can continue with our algebraic method.

We are not able to include N113 in our final mass function analysis because of the non-trivial free–free emission illustrated in Fig. 2 and our lack of observations at 3.2 mm for this field. Not being able to include N113 was unfortunate because this object contributed about the same number of clumps in Band 6 as the N159E map, even within the much smaller area.

3.2 Clump finding

In order to identify molecular gas clumps in our dust-only maps we used the ClumpFind algorithm (Williams, de Geus & Blitz 1994) available in starlink⁵ through the cupid package (Berry et al. 2007). We chose to use this algorithm to facilitate direct comparison with previous work that also used ClumpFind. A total of 32 dust clumps were identified in three fields (30 Dor-10, N159W and N159E), with the majority found in 30 Dor-10. Clump properties are summarized in Table 4.

The rms noise we input to ClumpFind was measured in each dust-only map from a region where no obvious emission, artifacts or negative bowls were visible. The lowest contour was set to three times the map rms to ensure that sources could be trusted with some confidence and that most of the emission associated with a clump was recovered. Contour spacings were set to two times the rms as this minimized splitting up sources while still recovering most of the obvious features. This combination of lowest contour level and contour spacing was shown to result in the lowest number of missed and “false” clumps through synthetic source-detection testing by Williams et al. (1994). Our own testing with several choices for these two settings showed that our final choice gave results that picked out distinct but somewhat blended sources without extracting spurious-looking sources, while staying well within believable bounds around obvious sources. Background subtraction was not used because the interferometer’s intrinsic spatial filtering removes large-scale emission, often causing a complicated interplay between bright emission and negative regions adjacent to the bright emission.

During completeness testing (Section 3.4) we found that our setting for the beam full width at half-maximum (FWHM) parameter in ClumpFind caused the algorithm to reject sources as too small when it should not have done so. To avoid this problem we instead set the FWHM parameter to zero and specified the ‘minpix’ parameter to reject clumps that were too small. To determine the correct ‘minpix’ value, we ran ClumpFind once with the FWHM set to the geometric mean of the synthesized beam axes and noted the ‘minpix’ value that ClumpFind calculated internally: 56 for 30 Dor-10 and N159W, and 61 for N159E. This method gave the same resulting clumps as when using the FWHM parameter but it prevented clumps near the size of the beam from being rejected as too small as a result of the lowest contour setting (Berry, private communication). Because this ‘minpix’ setting could select clumps that were of the same area as the synthesized beam but smaller than one of the beam dimensions (e.g. long and skinny clumps), we inspected the clump regions identified with ClumpFind visually for each field.

3.3 Converting dust flux to total gas mass

Using the dust mass calculated with equation (7), we converted to total gas mass using the gas-to-dust (G/D) ratio of 500 from Roman-Duval et al. (2017). This G/D ratio was obtained by stacking the dust SEDs across the entire LMC as measured with IRAS and Planck in bins of varying gas surface density. Roman-Duval et al. (2017) derived atomic gas surface densities from 21-cm Parkes observations (Staveley-Smith et al. 2003) and molecular gas surface densities from ¹²CO (1-0) NANTEN observations (Mizuno et al. 2001). The stacked SEDs were fitted with a modified blackbody to estimate the dust surface density, as well as the temperature and spectral emissivity index. We adopted the G/D from the bin with the highest gas surface density from Roman-Duval et al. (2017), as this is likely to be where stars are forming. Applying this conversion to our dust rms noise measurements from Table 3 results in mass sensitivities (1σ) of 45, 85 and 72 |$\mathrm{M}_{\odot}$| for 30 Dor-10, N159W and N159E, respectively.

Fig. 4 shows the cumulative dust-only clump mass function along with the mass function derived from the 1.3-mm map assuming that the 1.3-mm emission is purely from dust. It shows that our free–free correction is significantly larger than the random noise in each clump and systematically shifts the clump mass function to smaller masses.

Figure 4.

Comparison of cumulative clump mass functions from 1.3-mm maps (black squares), observed dust-only maps (orange triangles), and completeness-corrected dust-only maps (blue circles). The parameters used to convert the flux to gas mass in equation (7) are shown in the bottom left. Error bars are random uncertainties on clump masses from noise in the maps; see Section 3.3 for details.

Open in new tabDownload slide

3.4 Completeness corrections

We estimated the completeness of our mass function through Monte Carlo techniques by injecting synthetic sources of known fluxes and sizes into the dust-only maps one source at a time. We ran ClumpFind on the altered map with the same parameters as for the original maps and counted when synthetic clumps were recovered and when they were not. We then took the ratio of the number recovered to the total number injected as the completeness fraction. Synthetic source positions were randomly placed across the maps with a uniform distribution. All sources were elliptical Gaussians with FWHM and position angle identical to those of the synthesized beam for the corresponding map. We injected 600 synthetic sources per mass bin. Synthetic source masses were spaced evenly in log-space from roughly 21 to 1.2 × 10⁴ |$\mathrm{M}_{\odot}$|⁠. Fig. 5 shows the completeness data for 30 Dor-10, N159W and N159E.

Figure 5.

Completeness curves for point sources in 30 Dor-10, N159W and N159E (black line). The best-fit logistic curves (equation 8) are shown as orange dashed lines, with the best-fit parameters in the bottom right. 50 per cent completeness is marked by the vertical gray lines and reported as |$M_{50\%}$| in the bottom right.

Open in new tabDownload slide

A surprising feature of the completeness curves shown in Fig. 5 is that none of the maps actually attains 100 per cent completeness. This is related to how crowded each map is: a synthetic source added to a crowded map may land on a real clump, such that ClumpFind does not find the synthetic clump. 30 Dor-10 is the most crowded map and has the lowest maximum completeness, while N159E is the least crowded and has nearly 100 per cent recovery at high masses. This trend shows that it is not only the noise and resolution that limit our ability to find clumps but also the crowdedness of our maps.

It is important to note that we tested source recovery only for sources that were the same size as the synthesized beam, while most of the clumps that we have identified are much larger than the synthesized beam. Choosing an optimal completeness treatment is a general challenge when studying these types of objects, because molecular clumps are intrinsically amorphous in shape and vary significantly in physical size. Because larger clumps of the same integrated flux density as smaller clumps have a lower S/N, our completeness estimates probably overestimate the total fraction of clumps recovered for a given mass. While more sophisticated methods could be employed (e.g. see appendix B of Könyves et al. 2015), having a simplistic estimate of the completeness is better than having no estimate.

We use the fits to the completeness data to correct the mass function numbers by calculating N_corr = N(≥M)/f(M). We used the 30 Dor-10 completeness curve because it contributes the majority of the clumps in our mass function. Application of this completeness correction produces the mass function shown as blue circles in Fig. 4, which can be compared with the observed mass function shown as orange triangles.

4 MASS FUNCTIONS IN THE LMC AND MILKY WAY

4.1 LMC clump mass function

An empirical mass distribution is typically plotted as either a differential mass function or a cumulative mass function. A differential mass function is calculated by defining mass bins spanning the mass range of objects and counting the number of objects within each bin. This form allows for a simple Poissonian counting statistics approach to uncertainties for fitting, but the results can depend on the choice of bin widths and centres. A cumulative mass function avoids complications from binning because it is just an ordered tally of masses at or above a given object’s mass. The trade-off is that there is not a simple statistical approach for handling uncertainties for fitting.

Given our relatively small sample of clumps and the simple model, we carried out our fits with the standard non-linear Levenberg–Marquardt least-squares minimization (More 1977). To estimate fitting weights we followed the approach used by Reid & Wilson (2006a), which specifies the y-data uncertainties as the cumulative number N(≥M) for each clump. The best-fit double-power-law parameters are given in Table 5 for the mass function determined directly from the clump masses, along with two variations that attempt to account for incompleteness. The procedure and inputs for the fits are the same across all three variations.

Given the form of the fitting function, there are uncertainties only in the x-values. To estimate the absolute uncertainties on the best-fit parameters we used Monte Carlo simulation, as in Reid & Wilson (2006a). This process involved generating 10⁵ artificial mass functions and fitting each in the same way as the observed mass function. The artificial mass functions were made by drawing normally distributed random deviates centred on each measured clump mass and with standard deviations equal to the measurement uncertainties on the masses from the dust-only maps. Each newly generated sample of 32 masses was then sorted into descending order to prepare for fitting. Not only does this Monte Carlo approach simplify estimation of the fit parameter uncertainties but it also nicely includes the effects of neighbouring masses swapping places in the sorted order as a result of their uncertainties. From the distributions of fit parameters, we used the inner 68 per cent to obtain the 1σ uncertainties on each fit parameter (Table 5).

Fig. 6 shows the completeness-corrected mass function with the best-fit double power law (Table 5). We used all the clumps in our sample in this fitting, with our lowest mass of 205 |$\mathrm{M}_{\odot}$| at the 41 per cent complete level. We also tested restricting the fit to masses above 500 |$\mathrm{M}_{\odot}$| in the uncorrected mass function, because this is roughly where the 30 Dor-10 completeness curve flattens out at maximum completeness. The best-fit double power laws for the two methods end up being statistically indistinguishable based on our Monte Carlo uncertainty estimates (Table 5), so we focus our analysis on the completeness-corrected mass function.

Figure 6.

Completeness-corrected clump mass function (all points are the same as the blue circles in Fig. 4 but now with the originating source indicated), with best fits plotted as black dashed lines; the double power law is on the left and the lognormal function is on the right. Parameters of the fits (equations 9 and 10) are shown in the bottom left corners. Clumps are colour-coded by the region from which they originate, with 30 Dor-10 as black squares, N159W as orange triangles, and N159E as blue circles. Error bars are the same as in Fig. 4.

Open in new tabDownload slide

It may appear unusual that the numbers of high-mass clumps need to be boosted by our completeness correction. An overlapping pair of sources (merged clumps) would be identified in the data as a more massive clump and thus not be missed. We tested removing the completeness correction to the high-mass clumps by stitching together two sigmoid curves (the fit shown in Fig. 5 and a curve with the same width parameter but an amplitude of 1.0 and shifted slightly to higher masses) to create a smoothly varying sigmoid matching our measured completeness at low masses and asymptoting to 100 per cent complete at high masses. Refitting the mass function corrected with this piecewise completeness curve did steepen both indices and it moved the break to a higher mass, but all within the uncertainties of our original fit.

4.2 Comparison with other studies

The sample of 11 core/clump mass functions from seven distinct star-forming regions analysed by Reid & Wilson (2006b) is a useful one for comparing our LMC clump mass function with Milky Way mass functions. All mass functions were re-calculated from the original literature sources and fitted with a double power law (equation 9) homogeneously across all regions. While there are various differences in how those data were acquired and processed as well as in how individual cores and clumps were identified, the use of the best-fit parameters reported by Reid & Wilson (2006b) will remove any variability introduced in the fitting step. Fig. 7 compares the high- and low-mass slopes from Reid & Wilson (2006b) with our fit from Section 4.1 for the LMC.

Figure 7.

Comparison of the double-power-law mass function fits from the seven Galactic star-forming regions summarized in Reid & Wilson (2006b) with the low- and high-mass fit indices for our Large Magellanic Cloud (LMC) mass function. Circles and triangles show the Galactic low- and high-mass indices, respectively. The orange and blue dashed lines and shaded regions show the low- and high-mass indices for the LMC, respectively. Milky Way clouds are arranged with their median core/clump masses increasing to the right. Open triangles denote less robust high-mass slope determinations for Galactic clouds M8 (clumps identified by eye) and M17 (mass function fit better with the lognormal function).

Open in new tabDownload slide

The LMC break mass is greater than the individual clump masses for all of the Milky Way objects except for W43 and RCW 106 (and greater than all Milky Way break masses, which range from 0.2 ± 0.1 to 400 ± 300 |$\mathrm{M}_{\odot}$|⁠). As a result, it is probably more reasonable to compare the LMC low-mass slope with the high-mass slopes of the Galactic clouds. Fig. 7 shows that high-mass indices for ρ Oph at 850 |$\mu$|m, NGC 7538 (measured at both 450 and 850 |$\mu$|m) and W43 are consistent with the low-mass index for the LMC. All low-mass Galactic indices, except for Orion B from Motte et al. (2001), are systematically shallower than the LMC low-mass index. The observed consistency between the LMC low-mass index and the high-mass indices from ρ Oph, NGC 7538 and W43 suggests that the shape of the mass function may extend smoothly from the lower-mass clumps of the Milky Way regions to the higher-mass clumps of the LMC. In other words, it appears that the mass function slopes at high masses in the Milky Way extend up to the masses at the bottom of our LMC mass function.

In addition, none of the Galactic mass functions had completeness corrections applied and they are all probably suffering from incompleteness at their low-mass ends. Limiting our comparison to the high-mass end of the Galactic mass functions, where incompleteness should naturalli have a minimal effect, gives a fairer comparison with our LMC mass function because the completeness correction should ensure the effect of incompleteness is also minimized.

We note that there are four measurements of Galactic high-mass indices that are consistent with our high-mass index. Orion B as measured by both Motte et al. (2001) and Johnstone et al. (2001) is consistent within the uncertainties, despite values that are about 20 per cent shallower than ours, because the indices have uncertainties of almost 25 per cent. The remaining fields whose high-mass indices matches ours, M8 and M17, have large enough differences in analysis that we do not believe they present an issue for our conclusions. The cores from M8 were identified by eye rather than algorithmically as well as being fitted with an assumed Gaussian profile, unlike most of the other clump-finding procedures, including our own. For M17, the mass functions were either fitted as well as or better by a lognormal form, so it may not be a fair comparison with our mass function if the M17 mass distributions are not well characterized by a double power law.

The break masses from the Milky Way regions are all significantly smaller than the best-fit break mass for the LMC. Galactic break masses range from 0.2 ± 0.1 to 400 ± 300 |$\mathrm{M}_{\odot}$|⁠. These differences are not surprising, as Reid et al. (2010) showed a correlation between distance and break mass using mass functions derived from synthetic observations of a simulated GMC. This trend can even be seen among the 11 mass functions from Reid & Wilson (2006b), as the break masses increase in step with decreasing spatial resolution. Our results for the break mass in the LMC are reasonable, given the much greater distance to the clumps.

Fig. 8 shows this comparison in another way. All Milky Way mass functions are reproduced from Reid & Wilson (2006b), while the LMC mass function is scaled such that the lowest mass (205 |$\mathrm{M}_{\odot}$|⁠) is equal to the median maximum Milky Way mass (120 |$\mathrm{M}_{\odot}$|⁠) and the lowest mass has N(≥M) = 1. Although this scaling is rather arbitrary, it illustrates the potential connection of the mass function shape between the lower masses in the Milky Way and the higher masses in the LMC.

Figure 8.

All Milky Way mass functions described in Reid & Wilson (2006b) reproduced from their fig. 8 in grey. The LMC mass function, in orange, has been normalized such that the lowest-mass clump is equal to the median of the maximum masses from the Galactic mass functions, and N(≥M) has been normalized such that the lowest mass has N(≥M) = 1.

Open in new tabDownload slide

Other recent observations have also measured the mass function of cores and clumps in the Milky Way. Netterfield et al. (2009) observed 50 square degrees of the Galactic plane in Vela with the balloon-borne large-aperture submillimetre telescope (BLAST) and fitted cold cores (<14 K) with a power law of –3.22 ± 0.14 and warmer cores with an index of −1.95 ± 0.05. A total of 738 cores were extracted from Herschel Gould Belt Survey observations of Aquila by Könyves et al. (2015). These authors fitted their mass function with an index of −2.33 ± 0.06, but after applying an age correction to the distribution they fitted an index of −2.0 ± 0.2. Pattle et al. (2017) observed the Cepheus flare region with SCUBA-2, which included four molecular clouds, namely L1147/L1158, L1174, L1251 and L1228. Individual power-law fits gave indices of −1.8 ± 0.2, −2.0 ± 0.2, −1.8 ± 0.1 and −2.3 ± 0.3, respectively, while fitting all regions together with a double power law gave a low-mass index of −1.9 ± 0.1 and a high-mass index of −2.6 ± 0.3. These results are in line with our conclusion that the high-mass Milky Way mass function slope is similar to the low-mass LMC slope, except when the entire Cepheus flare region is analysed together.

Another interesting comparison we can make is to the clump mass function for 30 Dor-10 derived by Indebetouw et al. (2013) using ¹³CO (2–1) and ¹²CO (2–1) data. Their single power-law fit to the differential logarithmic mass distribution of CO clumps above 500 |$\mathrm{M}_{\odot}$| resulted in a slope of α = −1.9 ± 0.2, which is consistent with the low-mass index measured from dust masses in this work. Fitting the differential form of the mass function resulted in only about five bins, two of which only have one or two clumps, so it is reasonable that the bins below our break mass would dominate their fit. The similar mass range between our dust-estimated masses and their independent CO-estimated masses is evidence for the accuracy of our dust-only map generation and our estimated dust properties.

Reid & Wilson (2006b) also presented lognormal fits to the Galactic mass functions listed in Fig. 7. The Milky Way values for M₀ are all significantly lower than the LMC value, with only RCW 106 coming within ∼2σ, but this is expected because M₀ is correlated with the mass function median mass similarly to the correlation between M₀ and the double-power-law break mass. We find agreement in the LMC σ value and six of the Galactic mass functions. Both the ρ Oph measurements, and M8 and NGC 7538 measured at 850 |$\mu$|m show significantly different best-fit values for σ.

We can also search for a link between the LMC clump mass function and the mass distribution of groups of stars as they are first formed. The slope of −2 for young clusters measured by Lada & Lada (2003) is quite similar to our low-mass slope over a similar mass range. Analogous to how the similar shape of core mass functions and the stellar IMF has been posited as indicating a constant star-formation efficiency, the similarity of the clump and cluster mass functions may result from a constant cluster formation efficiency.

However, dust obscuration challenges observational measurements of the youngest, most deeply embedded star clusters. This hinders our ability to fully constrain the initial mass distribution of clusters when their stars are first born. Cluster formation simulations can bridge the gap between starless molecular clumps and the first star clusters as they emerge from their natal molecular clouds. Fujii & Portegies Zwart (2015) used a hydrodynamical simulation of the parent GMC, and after about a free-fall time instantaneously removed the gas and continued a solely N-body simulation of the stars formed. They measured mass function indices of −1.73 at 2 Myr and of −1.67 at 10 Myr for a range of initial GMC masses from 4 × 10⁴ to 5 × 10⁶ |$\mathrm{M}_{\odot}$|⁠. Howard, Pudritz & Harris (2017) aimed to simulate cluster formation with the inclusion of radiative feedback in their hydrodynamical simulations in order to study its impact on cluster-forming efficiency. Starting from GMCs of similar masses to in Fujii & Portegies Zwart (2015), they found a close agreement between their simulated mass function power-law index of −1.99 ± 0.14 and observed embedded clusters. This was after 5 Myr and for masses from 10 to 2000 |$\mathrm{M}_{\odot}$|⁠. For clusters >2000 |$\mathrm{M}_{\odot}$|⁠, they fitted a steeper power law of index −3.8. These theoretical results are broadly consistent with the LMC clump mass function presented here.

5 CONCLUSIONS

We have presented a clump mass function measured in continuum dust emission from star-forming molecular clouds in the LMC derived from public archival ALMA observations. A total of 32 clumps were observed in the 30 Dor-10, N159W and N159E regions, with masses ranging between 205 |$\mathrm{M}_{\odot}$| and 5740 |$\mathrm{M}_{\odot}$|⁠. In order to derive these total masses from the thermal dust emission we presented a technique of combining 3.2- and 1.3-mm maps to remove free–free emission contamination and produce dust-only emission maps. We implemented a point-source completeness correction to the clump numbers to account for the sensitivity and finite resolution of the observations

From our clump mass function we derived a double-power-law best fit with parameters |$\alpha _{\mathrm{low}} = -1.76^{+0.04}_{-0.05}$|⁠, α_high = −3.3 ± 0.2, and |$M_{\mathrm{break}} = 2500^{+300}_{-200} \ \mathrm{M}_{\odot }$|⁠. Our low-mass power-law index is similar to the high-mass indices of Milky Way regions hosting lower-mass cores and clumps. This may indicate an extension of the Milky Way power law to higher masses. The addition of more clumps from more regions in the LMC, along with mass functions derived from spectral line emission, will be necessary to better characterize the mass distribution of these objects in the LMC and to investigate its consistency with the Milky Way mass functions.

ACKNOWLEDGEMENTS

We thank the anonymous referee for detailed comments that improved the content of this paper. This paper makes use of the following ALMA data:

• ADS/JAO.ALMA#2011.0.00471.S,

• ADS/JAO.ALMA#2012.1.00554.S,

• ADS/JAO.ALMA#2012.1.00603.S,

• ADS/JAO.ALMA#2013.1.00832.S,

• ADS/JAO.ALMA#2013.1.01091.S,

• ADS/JAO.ALMA#2013.1.01136.S, and

• ADS/JAO.ALMA#2016.1.01533.S.

ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration, 2013). This research made use of APLpy, an open-source plotting package for Python (Robitaille & Bressert 2012). This research has made use of the AstroBetter blog and wiki. CDW acknowledges financial support from the Canada Council for the Arts through a Killam Research Fellowship. The research of CDW is supported by grants from the Natural Sciences and Engineering Research Council of Canada and the Canada Research Chairs program.

Footnotes

REFERENCES

APPENDIX A: OBSERVED, DUST-ONLY, AND FREE–FREE-ONLY MAPS