EDP Sciences
Press Release
Free access
Issue
A&A
Volume 555, July 2013
Article Number A112
Number of page(s) 10
Section Interstellar and circumstellar matter
DOI http://dx.doi.org/10.1051/0004-6361/201321318
Published online 10 July 2013

© ESO, 2013

1. Introduction

The formation of massive stars remains, in many ways, a mystery (Beuther et al. 2007; Zinnecker & Yorke 2007). More specifically, the key question of which physical processes determine their mass accretion history is yet to be answered. On one hand, some theories predict that primordial fragmentation of globally stable molecular clouds may form compact reservoirs of gas, called cores (with sizes up to 0.1 pc), from which a forming star subsequently accumulates its mass (McKee & Tan 2003; Beuther & Schilke 2004). In an alternative scenario, molecular clouds undergo global collapse (Peretto et al. 2006, 2007), gathering matter from large scales to the centre of their gravitational potential well, where cores, and protostars in them, are simultaneously growing in mass (Bonnell et al. 2004; Smith et al. 2009). Even though only statistical studies of large samples of massive star-forming clouds can definitely tell us which of these scenarios, if any, is most relevant, detailed observations of individual massive star-forming clouds can still provide important hints.

thumbnail Fig. 1

a) Mid-infrared Spitzer composite image (red: 8 μm; green: 4.5 μm; blue: 3.6 μm) of SDC335. The 6 filaments identified by eye are indicated with yellow dashed lines, emphasizing their converging pattern. The diffuse 4.5 μm emission associated with the two IR sources in the centre is usually interpreted as a signature of powerful outflow activity. The positions of the two cores are marked with black crosses. b) Herschel column density image of SDC335. The locations of the filaments and cores are marked similarly as in the a) panel. The final angular resolution of this image is 25′′ (yellow circle), that of Herschel at 350 μm (see text). The contours range from 3.5 × 1022 to 9.5 × 1022 cm-2 in steps of 2 × 1022 cm-2, and from 2.15 × 1023 to 4.15 × 1023 cm-2 in steps of 1 × 1023 cm-2. The two yellow contours define the regions in which we calculated the SDC335 and Centre region masses quoted in Table 1. c) ALMA 3.2 mm dust continuum emission of the central region of SDC335 where two cores are identified, MM1 and MM2. The rms noise is 0.4 mJy/beam. The contours range from 2 to 22 in steps of 5 mJy/beam, and from 22 to 62 in steps of 10 mJy/beam. The yellow ellipse represents the ALMA beam size.

Open with DEXTER

The cloud under investigation is the Spitzer dark cloud SDC335.579-0.292 (hereafter SDC335; Peretto & Fuller 2009), a massive infrared dark cloud (IRDC) located at a distance of 3.25 kpc from the Sun (distance obtained using the Reid et al. 2009 model). The low level of radiative feedback from protostars in IRDCs ensures that the initial conditions for star formation are still imprinted in the gas properties (Rathborne et al. 2006; Peretto & Fuller 2010). Massive IRDCs, such as SDC335, are therefore ideal places to study the earliest stages of high-mass star formation (e.g. Kauffmann & Pillai 2010). SDC335 exhibits a remarkable network of filaments seen in extinction at 8 μm (Fig. 1), reminiscent of hub-filament systems (Myers 2009) observed in a number of low-mass (André et al. 2010; Peretto et al. 2012) and high-mass (Schneider et al. 2012; Hennemann et al. 2012) star-forming regions. The SDC335 filaments intersect at the centre of the cloud where two infrared protostars (Lbol > 2 × 104 L; Garay et al. 2002) excite extended 4.5 μm emission, a tracer of powerful outflow activity (Cyganowski et al. 2008). Consistently, class II methanol masers, unique tracers of massive star formation (Xu et al. 2008), have also been reported towards these sources (Caswell et al. 2011). However, despite these signposts of massive star formation, no 6 cm free-free emission has been detected towards SDC335 down to a limit of 0.2 mJy (Garay et al. 2002). This shows that little gas has been ionized in the centre of SDC335 and suggests that we are witnessing the early stages of the formation of, at least, two massive stars.

The goal of this paper is to map the dense gas kinematics of SDC335 and analyse it in the context of massive star formation scenarios. In Sect. 2 we describe the observations. In Sect. 3 we discuss the mass partition in SDC335, and Sect. 4 presents observations of the SDC335 dense gas kinematics. Finally, we discus our results and their implications in Sect. 5, the summary and conclusions are presented in Sect. 6.

2. Observations

2.1. Spitzer and Herschel observations

We used publicly available1Spitzer GLIMPSE data (Churchwell et al. 2009). The angular resolution of the 8 μm data is ~2′′. We also used the PACS (Poglitsch et al. 2010) 160 μm and SPIRE (Griffin et al. 2010) 350 μm Herschel (Pilbratt et al. 2010) data from the Hi-GAL survey (Molinari et al. 2010). These data were reduced as described in Traficante et al. (2011), using the ROMAGAL map making algorithm. The nominal angular resolution at these two wavelengths are 12′′ and 25′′.

2.2. Mopra observations

In May 2010 we observed SDC335 with the ATNF Mopra 22 m single-dish telescope. We observed transitions including HCO+(1−0), H13CO+(1−0) and N2H+(1−0) in a 5′ × 5′ field centred on SDC335. We performed on-the-fly observations, switching to an off-position free of dense gas emission. Pointing was checked every hour and was found to be better than 10′′. We used the zoom mode of the MOPS spectrometer providing a velocity resolution of 0.1 km s-1. The angular resolution of these 3 mm Mopra observation is ~37″ and the rms noise is 0.1 K on the scale (~0.2 K on the Tmb scale because the beam efficiency factor is ~2 at 93 GHz on Mopra − Ladd et al. 2005).

2.3. ALMA observations

In September and November 2011 we observed SDC335 at 3 mm wavelength with the 16 antennas of ALMA (Cycle 0) in its compact configuration. We performed an 11-pointing mosaic covering the entire area seen in extinction with Spitzer (Fig. 1a). We simultaneously observed the 3.2 mm dust continuum, along with the CH3OH(13−12) and N2H+(1−0) transitions at a spectral resolution of ~0.1 km s-1. Flux and phase calibration were performed on Neptune and B1600-445, respectively. The data were reduced using CASA2 (McMullin et al. 2007). The synthesized beam is 5.6″ × 4.0″ with a position angle of + 97°. The rms noise in the continuum is 0.4 mJy/beam, while for the line we reach an rms sensitivity of 14 mJy/beam (~0.08 K).

As with any interferometer, ALMA filters out large-scale emission. To recover this emission, we used the Mopra single-dish data to provide the short-spacing information, for which we used the GILDAS3 software. This combination significantly improved the image quality, in particular in the central region of SDC335. The rms noise on these combined datacubes is 0.14 Jy/beam (~0.8 K), significantly higher than the ALMA-only dataset. This reflects the higher noise of the Mopra dataset per ALMA beam.

3. Mass partition in SDC335

The mid-infrared composite image of SDC335 is displayed in Fig. 1a. In extinction we easily identify a network of six filaments (indicated as yellow dashed lines, and named F1 to F6), while in the centre of SDC335 we observe bright infrared sources exciting diffuse 4.5 μm emission (in green). In the following we provide mass measurements of the entire SDC335 clump, the filaments, and the cores at the centre.

3.1. Clump and filaments

Table 1

Mass partition in SDC335.

Mid-infrared extinction mapping of IRDCs is a powerful method for measuring their column density distribution at high resolution (Peretto & Fuller 2009; Butler & Tan 2009). From the resulting maps one can measure the masses of these clouds. However, this method is limited by definition to mid-infrared absorbing dust, and in cases such as SDC335, where bright 8 μm sources have already formed, the resulting dust extinction masses become more uncertain. Using Herschel data allows us to circumvent this problem by looking at far-infrared dust emission from 70 μm up to 500 μm. A standard way for recovering the column density distribution from Herschel data is to perform a pixel-by-pixel spectral energy distribution fitting after smoothing the data to the Herschel 500 μm resolution (36′′). To obtain a higher angular resolution, we decided here to use the 160 μm/350 μm ratio map4 of SDC335 as an indicator of the dust temperature, and then reconstruct the column density distribution at 25′′ resolution by combining the dust temperature (Td) and 350 μm (S350) maps, assuming that dust radiates as a modified black-body. The column density is therefore written as:(1)where B350 is the Planck function measured at 350 μm, μ = 2.33 is the average molecular weight, mH is the atomic mass of hydrogen, and κ350 is the specific dust opacity at 350 μm. Using the dust opacity law (Hildebrand 1983; Beckwith et al. 1990) cm2 g-1 with β = 2, we constructed the SDC335 column density map presented in Fig. 1b. We can see that, despite the difference in resolution, the column density structure of SDC335 follows the extinction features we see in the mid-infrared. Moreover, the entire central pc-size region lies above a high column density of 1 × 1023 cm-2.

The dust opacity parameters we used to construct the SDC335 column density map are very uncertain. There are indications (e.g. Paradis et al. 2010) that for cold, dense clouds a spectral index β ≃ 2.4 is possibly more appropriate. The net effect of using a lower β is to overestimate the temperature, and therefore underestimate the column density (and the mass calculated from this column density). Also, the dust emission towards IRDCs is composed of the emission from the the cloud itself and emission from the galactic plane background/foreground dust, which is warmer and more diffuse. A consequence of this is that, once again, the IRDC temperature we calculate is overestimated. The measured average dust temperature towards the coldest parts of SDC335 is ~16 K, which appears to be, indeed, 2−3 K warmer compared with other IRDCs (Peretto et al. 2010). Altogether, the column densities presented in Fig. 1b are likely to be underestimated. To obtain an upper limit on the SDC335 column density we decreased the dust temperature map by 2 K, which relocates SDC335 in the typical temperature range observed in IRDCs. We then recalculated the column density map using Eq. (1). Masses and uncertainties quoted in Table 1 for SDC335, the Centre region, and the F1 and F2 filaments were obtained taking the average mass obtained from the two Herschel column density maps. All masses are background subtracted, which corresponds to the mass enclosed in a specific column density contour for a centrally concentrated source. For SDC335 and the Centre region, the sizes correspond to the diameter of the disc with the same areas, for F1 and F2 they correspond to the two dimensions of the rectangle with the same areas. The surface areas of these two filaments correspond to the polygons we drew on the Herschel column density map around the F1 and F2 filaments. The widths of these polygons are constrained by the filament profiles as seen in the 8 μm extinction map (4′′ resolution), which correspond to ~0.3 pc. Note that for the filaments we performed an independent mass measurement directly using the 8 μm extinction map from Peretto & Fuller (2009), confirming the Herschel mass measurements. Densities were then calculated assuming a uniform density and spherical geometry for SDC335 and the Centre region, and cylindrical geometry for the F1 and F2 filaments.

3.2. Dense cores

The ALMA 3.2 mm dust continuum observations presented in Fig. 1c show two bright sources, MM1 and MM2, which have J2000 coordinates (RA: 16h30m58.76s; Dec: −48°43′53.4′′) and (RA: 16h30m57.26s; Dec: −48°43′39.7′′), respectively. Each of these cores is associated with one infrared source and class II methanol maser spots (Caswell et al. 2011), leaving no doubt that MM1 and MM2 are currently forming massive stars.

thumbnail Fig. 2

ALMA CH3OH(13−12) integrated-intensity image of SDC335 in colour scale; overplotted are the contours of the 3.2 mm dust continuum emission as displayed in Fig. 1c. We can see that the two types of emission coincide spatially. The insert in the top left corner shows the methanol spectra observed at the central position of each core. The cyan solid lines are the best Gaussian fit to the data. The yellow ellipse represents the ALMA beam size.

Open with DEXTER

The MM1 and MM2 cores are compact but partially resolved. The results of a 2D Gaussian fit to the 3.2 mm continuum emission of these compact sources give integrated fluxes of 101(± 10) mJy and 12(± 2) mJy and deconvolved sizes of 0.054 pc and 0.057 pc for MM1 and MM2, respectively. To estimate how much of this emission could be free-free, we scaled the 3σ non-detection limit at 6cm from Garay et al. (2002) to 3.2 mm using . Upper limits on α have been determined for a set of massive protostellar objects and HCHII regions (Cyganowski et al. 2011; Hoare 2005) consistent with α = 15. This provides an upper limit for free-free contamination at 3.2 mm of mJy. Clearly, this is negligible for MM1, whereas it could contribute up to 33% of the MM2 flux. In the strict upper-case limit of optically thick free-free emission at both 6 cm and 3.2 mm α equals 2, increasing the upper limit for free-free contamination to 70% of the MM1 ALMA 3.2 mm flux. However, the recent detection of MM1 at 7 mm with ATCA (Avison et al., in prep.) shows no excess emission over that expected from the dust. Still, conservatively assuming that the entire 7 mm emission is from optically thick free-free emission implies an upper limit on the free-free contamination of 30% of the observed 3.2 mm flux. Such a level of contamination would not change any of the results presented here, therefore we neglect any potential free-free contamination in the remainder of the paper. The gas mass and the 3.2 mm flux of the MM cores are related through (2)where d is the distance to the source, F3.2 is the 3.2 mm flux, κ3.2 is the specific dust opacity (accounting for the dust-to-gas-mass ratio) and B3.2(Td) is the Planck function measured at 3.2 mm with a dust temperature Td. The main sources of uncertainties on this mass estimate come from the dust properties, temperature, and opacity. The dust temperatures of these two sources are difficult to determine based on these ALMA observations alone. However, both sources have strong mid/far-IR emission seen with Spitzer (Churchwell et al. 2009) and Herschel (Molinari et al. 2010), class II methanol maser emission (Caswell et al. 2011), and are detected in high-excitation thermal lines of methanol (Sect. 4). These are indicative of dust in the centre of the cores with temperatures >100 K, but it is also clear that on larger scales the dust within the dark SDC335 filaments is cold, with temperatures ~15 K as measured in many other IRDCs (Peretto et al. 2010; Wilcock et al. 2012). For the vast majority of massive protostellar cores in the literature (cf. caption of Fig. 6), the assumed or measured dust/gas temperature (via SED or K-ladder fitting of some specific lines) varies between 15 K and 100 K. Here, we adopted an intermediate dust temperature of 50 K for both MM cores, and considered that a factor of 2 uncertainty on this dust temperature is conservative. In the future, radiative transfer modelling of these sources will be necessary to better constrain their temperature profiles.

We took the same dust opacity law as used for the Herschel data, providing κ3.2 = 8.7 × 10-4 cm2 g-1 (assuming a dust-to-gas-mass ratio of 1%). However, this value is sensitive to the dust model used. It is unclear which model is the most appropriate for protostellar cores, but as shown in Fig. 6, most values adopted in the literature for core-mass measurements agree within a factor of 2. With these assumptions we estimated the gas masses and their associated uncertainties as and .

thumbnail Fig. 3

Spitzer 8 μm image of SDC335 (colour scale) over-plotted with the Mopra HCO+(1−0) spectra. The temperature scale and velocity are indicated in the bottom-left corner. The HCO+(1−0) line is self-absorbed and blue-shifted in the bulk of the cloud. This is usually interpreted as a signature of collapse.

Open with DEXTER

thumbnail Fig. 4

a) Same as in Fig. 1a; b) ALMA-only image of the N2H+(1−0) integrated intensity over the 7 hyperfine structure components. The rms noise on the resulting map is ~ 6 mJy/beam km s-1. The contours go from 0.1 to 1.5 in steps of 0.7 Jy/beam km s-1 and 1.5 to 9 in steps of 1.5 Jy/beam km s-1. The crosses mark the positions of the two dense cores. The ALMA beam is represented as a yellow elliptical symbol in the bottom-right corner of the image. We can see the excellent match between the Spitzer dust extinction of the filaments and the N2H+(1−0) emission; c) ALMA N2H+(1−0) velocity field using the first order moment map. The crosses mark the positions of the cores and the contours are the same as in the b) panel.

Open with DEXTER

4. Dense-gas kinematics in SDC335

In this section we discuss the dense-gas kinematics of the cores, filaments, and clump as observed with the Mopra and ALMA telescopes.

4.1. ALMA CH3OH(13−12) core velocities

To determine the systemic velocity of the MM cores, we mapped the thermal methanol CH3OH(13−12) transition at 105.063761 GHz. Because of the high energy levels of this transition (Eu = 223.8 K), CH3OH(13−12) is preferentially observed in dense and warm regions. Figure 2 shows the ALMA integrated-intensity image towards the cores. We see the excellent agreement between the position of the peak of the dust continuum cores and the methanol emission, indicating that methanol is a good tracer of their systemic velocities. We also note that the methanol emission is more compact (unresolved, i.e. <0.01 pc) than the dust continuum emission, which indicates that it arises from the warm, innermost regions of the cores. Gaussian fits to the methanol spectra observed at the central position of the two cores (insert of Fig. 2) provide the systemic velocities of the cores (VMM1 =  −46.6 km s-1, VMM2 =  −46.5 km s-1) and the gas velocity dispersion in the densest parts of MM1 and MM2 (ΔVMM1 = 4.6 km s-1, ΔVMM2 = 4.8 km s-1).

4.2. Mopra HCO+(1−0) self-absorbed lines

HCO+ is a well-known tracer of dense gas in molecular clouds. In these regions, HCO+(1−0) can be optically thick, in which case the line shape can provide information of the global motions of the gas along the line of sight (e.g. Fuller et al. 2005; Smith et al. 2012). The HCO+(1−0) observations towards SDC335 (Fig. 3) show blue-shifted self-absorbed spectra in the bulk of the cloud. Such line profiles are expected for an optically thick tracer of idealized collapsing clouds in which the excitation temperature is rising towards the centre. What is important to note here is the extent (over at least 12 independent Mopra beams) over which this spectral signature is observed, and the absence of any other line asymmetry. For expanding motions we would expect red-shifted self-absorbed spectra, while in the case of rotation blue-shifted and red-shifted spectra on either side of the rotation axis should be produced. Therefore these HCO+(1−0) observations towards SDC335 already rule out the possibility of a rotating or expanding cloud, and strongly suggest that SDC335 is collapsing.

SDC335 is well enough characterised that we can estimate the HCO+ abundance using the 1D non-LTE RADEX radiative transfer code (van der Tak et al. 2007). This code predicts line intensities based on a set of input parameters for which we have strong constraints: the kinetic temperature (20 ± 5 K, estimated from Herschel data), the cosmic background temperature (2.73 K), the central H2 density averaged over the Mopra beam (6 ± 1 × 104 cm-3, estimated from the column density map presented in Fig. 1), and the velocity dispersion (1.3 ± 0.3 km s-1; cf. Sect. 5.4). Then we iterate on the last input parameter, i.e. the molecule column density, to match the model line intensities with the observed line temperature, i.e. K on the Tmb scale. Doing so, we obtain cm-2, corresponding to an abundance . The corresponding excitation temperature is K, confirming that HCO+(1−0) is not thermalised. Using the same set of parameters, we performed the same exercise for the central H13CO+(1−0) line (se Fig. B.1), which has K on the Tmb scale. For this line we obtain cm-2, corresponding to an abundance . The corresponding excitation temperature is K. Therefore, as for HCO+(1−0), H13CO+(1−0) is not thermalised. Note that the lower excitation temperature of H13CO+(1−0) is most likely due to a lower beam-filling factor. Another important point is that given the abundances we calculated for both molecules, we obtain an abundance ratio 15 ≤ [HCO+] / [H13CO+] ≤ 20. The [12C] / [13C] ratio is known to increase as a function of the galactocentric radius, and at the galactocentric distance of SDC335 (i.e. ~5 kpc) the predicted [12C] / [13C] is around 30 (Langer & Penzias 1993; Savage et al. 2002). The value we find is about half this value, which, considering the uncertainties on these kinds of measurements, is in reasonable agreement. We use the latter value of the fractional abundance for the radiative modelling presented in Sect. 5.3.

thumbnail Fig. 5

(Left) Examples of combined ALMA and Mopra N2H+(1−0) spectra observed at some specific positions in SDC335 (upper panel for filaments − lower panel for cores), along with their best fits as red solid lines. The N2H+(1−0) spectra exhibit a hyperfine structure (HFS) composed of 7 components (the positions are displayed as blue vertical ticks for the F2 filament). Some of these components are close enough to be blended when the velocity dispersion of the gas is supersonic, resulting in three groups of lines. For a kinetic temperature of 10 K, the velocity dispersion along the filaments is supersonic by a factor 1.5 to 3, similar to what is observed in other IRDCs (Ragan et al. 2012). (Right) Schematic representation of the systemic velocity and velocity dispersion of the different structures. The length of each box represents the velocity dispersion (FWHM) of the gas, and its central position the systemic velocity (represented as a filled circle). The colour of the boxes codes the line that has been used for the measurements: red for N2H+(1−0), and cyan for CH3OH(13−12). The vertical dashed lines mark the systemic velocities of the cores.

Open with DEXTER

4.3. ALMA N2H+(1−0) cloud velocity field

Figure 4b shows the ALMA N2H+(1−0) integrated-intensity map of SDC335. The visual comparison with the Spitzer image of SDC335 demonstrates how efficient this molecule is in tracing the network of pc-long filaments seen in dust extinction. This justifies our choice of using this line to trace the filaments kinematics. On the other hand, we can also see that N2H+ is a poor tracer of the cores, where the central heating may have partly removed it from the gas phase (Zinchenko et al. 2009; Busquet et al. 2011).

Figures 4c and 5 show that SDC335 velocity field is homogeneous in each filament, with distinct velocities from filament to filament (e.g. ⟨ VF1 ⟩ =  −47.4 ± 0.1 km s-1; ⟨ VF3 ⟩ =  −45.8 ± 0.2 km s-1). It becomes more complex towards the centre of the cloud. In Fig. 5 we see that two separate velocity components are present close to MM2, while the broad asymmetric line profiles around MM1 suggest their blending, consistent with observations of other massive cores (Csengeri et al. 2011). This line shape cannot be the result of high optical depth since the N2H+(1−0) hyperfine line-fitting (performed with GILDAS) gives an opacity lower than 1 everywhere in the cloud. Kinematically, the gas traced by N2H+(1−0) at the centre of the cloud appears to be composed of a mix of the gas originating from the two main filaments, F1 and F2. Figure 5 (right) presents a schematic view of the velocities of the filaments and cores. We can actually see that the two cores lie at an intermediate velocity between the velocities of the different filaments. This configuration suggests that the cores are at least partly fed by the pristine gas flowing along these filaments at a velocity Vinf ≃ 1 km s-1.

5. Discussion

5.1. SDC335: an OB cluster progenitor

In Sect. 3 we inferred core masses of and in deconvolved diameters 0.05 − 0.06 pc. Figure 6 shows a radius-versus-mass diagram for a significant (but not complete) sample of massive protostellar cores published in the literature. In this diagram we can see that SDC335 MM1 stands out, and for cores with similar sizes MM1 is a factor of between 3 and 20 more massive. However, given the uncertainties on the dust properties and density profile of cores, MM1 could match the mass of the most massive cores observed, on smaller scales, in Cygnus X (Bontemps et al. 2010). Nevertheless, MM1 clearly appears to be one of the most massive, compact protostellar cores ever observed in the Galaxy.

Another interesting source that provides an informative comparison is W51 North. This source is believed to contain an already formed ≥65 M star, with a surrounding 3000 AU accreting disc of 40 M (Zapata et al. 2009). Adding this source to Fig. 6 using its total (star+disc) mass shows it to be an extreme object as well. SDC335 MM1 could represent an earlier version of such an O-type star-forming system. For compact cores, the fraction of mass likely to be accreted onto the star is typically 50% of the total core mass (McKee & Tan 2003; Duarte-Cabral et al. 2013). Despite probable unresolved fragmentation on smaller scales, the MM1 core and its large mass have the potential to form at least one star of 50 M to 100 M.

Assuming now that within SDC335 (M = 5500 ± 800 M) a fully sampled standard initial mass function forms (Kroupa 2002; Chabrier 2003), then, in addition to the early O-type star in MM1, a thousand solar mass cluster of ~320 stars with masses from 1 to 50 M should emerge from SDC335. Including lower mass stars in this calculation, we would reach a star formation efficiency ≥50%, the necessary condition to form a bound open cluster (e.g. Lada & Lada 2003). As a whole, SDC335 could potentially form an OB cluster similar to the Trapezium cluster in Orion (Zinnecker & Yorke 2007).

thumbnail Fig. 6

Mass-radius diagram of massive protostellar cores. The cyan circles correspond to the values as published in the literature (Peretto et al. 2007; Ren et al. 2012; Rathborne et al. 2011; Wang et al. 2011; Zhang & Wang 2011; Bontemps et al. 2010; Rathborne et al. 2007, 2008; Beuther et al. 2002, 2003, 2005, 2012; Molinari et al. 1998; Rodón et al. 2012), while the empty circles correspond to the same set of sources for which we recalculated the mass using the same dust opacity law as in this paper. The MM1 and MM2 sources are indicated as red filled circles. The green filled square marks the position of the W51 North star+disc system from Zapata et al. (2009). The shaded area indicates the region where most sources lie. The cross in the bottom right corner indicates a factor of 2 uncertainty in both the masses and sizes, typical of the results presented here.

Open with DEXTER

5.2. A large mass reservoir for MM1

We can estimate the conditions under which MM1 formed within the context of gravo-turbulent fragmentation models. Using the standard (lognormal) volume-density probability density function (PDF) of non-self-gravitating turbulent clouds (e.g. Padoan et al. 1997; Hennebelle & Chabrier 2008), we calculated (see Appendix A) that less than 0.01% of the gas is expected to be above a density of 107 cm-3, while >10% of the SDC335 mass is above this threshold in the form of cores (see Table 1). Therefore, gravity must have brought together such a large mass in such a small volume. A first possibility is that the material currently lying in MM1 was initially part of a larger volume that then collapsed. To calculate the diameter Dini of this volume we first need to calculate the density ρini above which 10% of the gas is lying, assuming a lognormal density PDF. Using the observed parameters of SDC335 (see Appendix A), we found that Dini of this initial volume must have been ~15 times larger than the current MM1 size, which means Dini ~ 0.8 pc. This size is in fact a lower limit since the calculation implicitly assumes that the entire dense gas above ρini lies within a single dense region. The second possibility is that MM1 initially had the same diameter as observed today. It is then possible to calculate the maximum mass that this volume can contain to match the lognormal PDF. We can show (see Appendix A) that the maximum initial mass of such a core is ~3 M. This low mass means that most of the current MM1 mass must have been accreted from its surroundings. Using the current average density of SDC335, we calculated that the region from which MM1 accreted matter would had to have a diameter of 1.2 pc. Either of these scenarios, therefore, requires large-scale, rather than local, accretion/collapse to form MM1.

5.3. Collapse on large scales

The Mopra HCO+(1−0) spectra presented in Fig. 3 are suggestive of global gravitational collapse. A simple analytical model (Myers et al. 1996) allows one to infer an infall velocity from these spectra. Using this model, we obtained an infall velocity of ~0.4 km s-1. However, as noted by De Vries & Myers (2005), this model underestimates the infall velocity by a factor of ~2. We therefore decided to run a more sophisticated radiative transfer model to better constrain this infall velocity. For this purpose we used the RATRAN 1D Monte Carlo radiative transfer code (Hogerheijde & van der Tak 2000). The input parameters for the calculations are the mass of the cloud, its radius, density profile, kinetic temperature profile, turbulent velocity dispersion, the infall velocity profile, and the abundance profile of the line to be modelled. Obviously, a 1D model cannot describe the detailed kinematics of the filamentary structures observed in SDC335, and for this reason we decided to model only the central HCO+(1−0) and H13CO+(1 − 0) spectra. We used the SDC335 mass and size quoted in Table 1, a cloud-density profile ρ ∝ r-1.5, and a constant temperature profile of 20 K. Based on the discussion in Sect. 4.2, we fixed the HCO+ abundance (relative to H2) to 7 × 10-10 and an abundance ratio [HCO+]/[H13CO+] of 30. We then ran a grid of models varying the last two input parameters, i.e., the infall velocity and the velocity dispersion. Figure 7 shows the results of the HCO+(1−0) modelling of the central pixel for infall velocities ranging from 0.4 km s-1 to 0.9 km s-1, and velocity dispersions from 0.8 km s-1 to 1.2 km s-1. For each spectrum we calculated a reduced χ2 parameter representative of the quality of the fit. This parameter is given in the top right-hand corner of each spectrum. The corresponding H13CO+(1−0) predictions from the model are displayed in Appendix B. From Fig. 7 we consider that 0.5 km   s-1 ≤ Vinf ≤ 0.9 km   s-1 and 0.9 km   s-1 ≤ σturb ≤ 1.1 km   s-1 provide reasonable fits to the central HCO+(1−0) spectrum. We also performed models varying the radius of the collapsing region Rinf. For Rinf < 0.5 pc the modelled HCO+(1−0) spectra remain symmetric, which is inconsistent with the observations. Only for Rinf ≥ 0.8 pc is the asymmetry large enough to resemble the observed one. This shows that the observed HCO+(1−0) self-absorbed spectra towards SDC335 do trace global collapse.

thumbnail Fig. 7

Grid of HCO+(1−0) spectra obtained from RATRAN modelling of a collapsing cloud (see text). All input parameters are fixed with the exception of the infall velocity (Vinf) and velocity dispersion (σturb). Each modelled spectrum (in red) has been obtained for the corresponding Vinf − σturb displayed in the top and right-hand sides of the figure. The numbers in the top right-hand corner of each spectrum give a relative idea of the fit quality. The HCO+(1−0) spectrum observed at the centre of SDC335 is plotted in black. Since we aimed to keep the observed spectra displayed in scale, we applied a 0.5 factor to the modelled spectra to take into account the Mopra main-beam efficiency.

Open with DEXTER

5.4. Energy balance

To be collapsing, the gravitational energy of a cloud has to overcome the kinetic energy that supports it. This occurs if the virial parameter is lower than 1 (Bertoldi & McKee 1992). In this equation σturb, R, and M are the 1D velocity dispersion, the cloud radius, and the cloud gas mass. In the case of SDC335, we estimated σturb = 1.3 km s-1 from the averaged N2H+(1−0) spectrum over SDC335 as observed with Mopra. Using Mopra 13CO(1−0) data, which trace less dense gas, we obtained σturb = 1.6 km s-1. Note also that these velocity-dispersion measurements include any systematic motions within the beam, such as infall, which artificially increase the velocity-dispersion estimate (Peretto et al. 2007). Taking this into account, and because the filaments are well traced by N2H+, we estimate σturb = 1.3(±0.3) km s-1. With R = 1.2 pc and M = 5500( ± 800) M we find . Additional support against gravity could be provided by the magnetic field. Following previous studies (Pillai et al. 2011), we estimated that the magnetic field strength | B | vir necessary to virialize SDC335 is | B | vir = 300 μG, which is at least three times higher than observations of clouds at similar densities suggest (Crutcher 2012). Finally, note that the support provided by centrifugal forces can potentially stabilise a cloud against gravity. However, calculating the rotational energy of SDC335 by assuming that it is a homogenous rotating sphere with an angular velocity ω = 1 km s-1/pc, we estimated that it is ~10 times smaller than its kinetic energy as measured from the velocity dispersion. In other words, it is negligible.

5.5. Large-scale velocity field and accretion rates

To illustrate some of the expected signatures of globally collapsing clouds we present, in Fig. 8 a snapshot of a published MHD simulation modelling the evolution of a turbulent and magnetized 10 000 M cloud, that was initially designed to reproduce some of the observational signatures of the DR21 region (Schneider et al. 2010,see Appendix C for more details on the simulation). Overall, this simulation shows some similarities with SDC335, i.e. massive cores in the centre, the formation of filaments converging towards these cores, and a velocity field resembling the one observed in SDC335 (see Fig. 4c). But most importantly, Fig. 8 shows that although a fraction of the gas is indeed collapsing along the filaments, a large portion is collapsing off filaments. In this case the filamentary accretion observed along the filaments represents only the tip of the entire accretion towards the cloud centre.

In the context of a global collapse scenario, the observed velocity field along the filaments is the consequence of the inflowing cold gas. We can therefore estimate the current infall rate of gas running through the filaments using , where Nfil is the number of filaments, Rfil is the filament cross-section radius, Vinf is the gas infall velocity, and ρfil is the gas volume density. With six filaments, an infall velocity of 0.7( ± 0.2) km s-1, a filament section radius of 0.15 pc, and a density of 4( ± 1) × 104 cm-3, we derive an infall rate of 0.7( ± 0.3) × 10-3 M/yr. At this rate, a total mass of 210( ± 90) M would have been gathered in the centre by filamentary accretion within a free-fall time tff ~ 3 × 105 yr. This is slightly less than half of the cumulated core masses. However, less than 20% of the SDC335 mass is lying within the filaments (cf. Sect. 3). Assuming that the remaining gas is collapsing off filaments at a similar infall velocity, as observed in the simulations (see Fig. 8a), the total accretion rate becomes , where Rsph is the radius of the considered spherical volume and ρsph is the density at that radius. At the radius of the Centre region, Rsph = 0.6 pc and ρsph = 1.3( ± 0.2) × 104 cm-3, which leads to inf = 2.5( ± 1.0) × 10-3 M/yr. With this accretion rate 750( ± 300) M of pristine gas is trapped inside the Centre region every cloud free-fall time. This is enough to double the mass of material currently present in the Centre region in cloud free-fall times. This is few free-fall times longer than the typical timescale over which simulations modelling the evolution of massive star-forming clouds evolve (e.g. Smith et al. 2009; Krumholz et al. 2012). However, free-fall times in simulations are estimated at the initial average density while here tff is calculated from the SDC335 current average density . Therefore, a fair comparison with models would imply that we estimate tff from SDC335 , which would then increase tff as . Altogether, evidence indicates that a significant fraction, if not all, of the SDC335 core mass could have been built through the parsec-scale collapse of their parental cloud over a few cloud free-fall times.

thumbnail Fig. 8

Snapshot of a MHD simulation of a 10 000 M collapsing cloud (see Appendix C for more details; Schneider et al. 2010). a) Column density (colour and contours) smoothed to the resolution of the ALMA data (5′′). The arrows show the plane-of-the-sky velocity field. We see that gas flows along filaments and also off the filaments. b) Dense gas line-of-sight velocity field (colour scale) smoothed to the resolution of the ALMA data. We highlighted the filaments by white dashed lines. The contours are the same as in panel a).

Open with DEXTER

6. Summary and conclusion

SDC335 is a massive (5500 ± 800 M) IRDC with two massive star-forming cores located in its centre, one of which is likely to be an early O-type star progenitor. This core has an estimated mass of in a deconvolved diameter of ~0.05 pc, which makes it one of the most massive protostellar cores ever observed in the Galaxy. A theoretical argument based on volume-density PDFs of molecular clouds suggests that such a concentration of mass must occur through the large-scale collapse of the surrounding cloud. This scenario was supported by several observational facts presented in this paper: optically thick molecular line observations showing extended collapse signatures; a virial parameter distinctly lower than 1; a velocity field consistent with models of globally collapsing molecular clouds; accretion rates that are high enough to provide an additional 750( ± 300) M of pristine gas to the central pc-size region of SDC335 per cloud free-fall time. Altogether, these observations strongly suggest that the SDC335 massive star-forming cores managed to build-up their large masses as a result of the supersonic global collapse of their surrounding cloud.

It is always adventurous to draw general conclusions based on a single example. However, there are several sources now for which large-scale infall has been suggested to play a major role in building up the mass of star-forming cores (e.g. Peretto et al. 2006; Schneider et al. 2010; Barnes et al. 2010; Galván-Madrid et al. 2010; Liu et al. 2012). Global infall also naturally explains the mass segregation observed in SDC335 and other regions. Although a number of questions remain to be answered, it is becoming clear that large-scale evolution of molecular clouds is the key to the mass determination process of massive stars.


4

Note that we did not use the 250 μm image because of the significant fraction of saturated pixels at the centre of SDC335.

5

Note that in the same study the authors also determine a lower limit of α > 1.7 for one source, but this measurement is based on a single 4.2σ detection of a very weak source.

Acknowledgments

We thank the anonymous referee for his/her thorough report, which significantly improved the quality of this paper. N.P. was supported by a CEA/Marie Curie Eurotalents fellowship and benefited from the support of the European Research Council advanced grant ORISTARS (Grant Agreement no. 291294). A.D.C. was supported by the PROBeS project funded by the French National Research Agency (ANR). J.E.P. has received funding from the European Community Seventh Framework Programme (/FP7/2007-2013/) under grant agreement No. 229517. N.S. was supported by the ANR-11-BS56-010 STARFICH project. We also acknowledge the support of the European ALMA Regional Centre (ARC) and the UK ARC node. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00474.S. ALMA is a partnership of ESO (representing its member states), NFS (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The Mopra radio telescope is part of the Australia Telescope National Facility, which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. The University of New South Wales Digital Filter Bank used for the observations with the Mopra telescope was provided with the support from Australian Research Council.

References

Appendix A: Volume-density PDF calculations

Volume-density PDFs of turbulent, non-self-gravitating molecular clouds can be described as a lognormal function of the logarithmic density contrast (Padoan et al. 1997; Hennebelle & Chabrier 2008): (A.1)where σ0 is the standard deviation of the distribution, and . Furthermore, the standard deviation of this PDF can be written as , where ℳ is the Mach number and b ≃ 0.25.

Integrating Eq. (A.1), we can estimate the fraction of a cloud that lies above a certain density threshold, ρth, before gravity becomes important. Setting the two free parameters of Eq. (A.1) to the SDC335 observed values ( cm-3 and ℳ = 6), we find that less than 0.01% of the gas should lie above ρth = 1 × 107 cm-3. This is more than three orders of magnitude different from what is observed in SDC335.

Now we can estimate the density ρini (and ) at which the following relation is fulfilled: (A.2)which is equivalent to (A.3)This gives ρini = 3.5 × 104 cm-3. Then we can calculate the volume diameter in which the MM1 mass was initially contained, using Dini = 2 × [3MMM1/(4πρini)] 1/3 = 0.8 pc. Note that here we used ρini as the mean density of the initial volume, while it formally is the minimum density within the volume under consideration. The true mean density is necessarily higher, although it cannot be too centrally concentrated either because this would not satisfy the volume density PDF. It is therefore reasonable to use ρini as the mean density, especially since the dependency of Dini on ρini is weak.

Alternatively, we can estimate the maximum Mini in the current MM1 volume that satisfies the volume density PDF: (A.4)with Mini = ρini   VMM1. We find that Mini ≃ 3 M, which means that in this case nearly the whole MM1 mass must come from its surrounding. We can estimate the volume of this region by taking pc.

Appendix B: H13CO+(1−0) RATRAN modelling

Figure B.1 presents the optically thin H13CO+(1−0) modelled spectra obtain with RATRAN for the cloud-collapse model discussed in Sect. 5.3. We see that the modelled lines overall match the observed spectrum quite well even though, in nearly all cases, the modelled spectrum is slightly too narrow. This can potentially be explained by a more complex infall profile than the one we used for these calculation, resulting in a slightly broader line. The modelled spectra have peak temperatures close to the observed one, which supports our choice of HCO+ and H13CO+ abundances.

thumbnail Fig. B.1

Same as Fig. 7, but for the H13CO+(1−0) line.

Open with DEXTER

Appendix C: Additional details on the MHD simulation

The simulation presented in this study (Fig. 8 of the paper) was initially performed to model the DR21 region (Schneider et al. 2010). It is an MHD simulation of a self-gravitating cloud performed with the AMR RAMSES code. The initial conditions of the simulation consisted of a 10 000 M ellipsoidal cloud with an aspect ratio of 2 and a density profile as ρ(r,z) = ρ0/ [1 + (r/r0)2 + (z/z0)2], where , z0 = 2r0, r0 = 5 pc, and ρ0 = 500 cm-3. The density at the edge of the cloud is ρedge = 50 cm-3, and is maintained in pressure equilibrium with an external medium at lower density. The magnetic field is perpendicular to the main axis of the cloud, with an intensity proportional to the cloud column density and a peak value of 7 μG. By the time of the snapshot presented in this paper, the magnetic field had increased to ~100 μG in the densest regions, with an average value over the dense gas of ~20 μG. The simulation is isothermal at a temperature of 10 K. A turbulent velocity field was seeded to initially obtain 2T + W + M ≃ 0 over the entire cloud, where T is the kinetic energy (thermal and turbulent), W the gravitational energy, and M the magnetic energy. Turbulence was undriven and allowed to decay. These conditions lead to W ≃ 2T ≃ 9M. Although globally, the turbulent and magnetic energy compensate for the gravitational energy of the cloud, it quickly becomes sub-virial due to the compressive nature of turbulence and because its energy quickly dissipates. The consequence of this is the fragmentation and global collapse of the simulated cloud.

All Tables

Table 1

Mass partition in SDC335.

All Figures

thumbnail Fig. 1

a) Mid-infrared Spitzer composite image (red: 8 μm; green: 4.5 μm; blue: 3.6 μm) of SDC335. The 6 filaments identified by eye are indicated with yellow dashed lines, emphasizing their converging pattern. The diffuse 4.5 μm emission associated with the two IR sources in the centre is usually interpreted as a signature of powerful outflow activity. The positions of the two cores are marked with black crosses. b) Herschel column density image of SDC335. The locations of the filaments and cores are marked similarly as in the a) panel. The final angular resolution of this image is 25′′ (yellow circle), that of Herschel at 350 μm (see text). The contours range from 3.5 × 1022 to 9.5 × 1022 cm-2 in steps of 2 × 1022 cm-2, and from 2.15 × 1023 to 4.15 × 1023 cm-2 in steps of 1 × 1023 cm-2. The two yellow contours define the regions in which we calculated the SDC335 and Centre region masses quoted in Table 1. c) ALMA 3.2 mm dust continuum emission of the central region of SDC335 where two cores are identified, MM1 and MM2. The rms noise is 0.4 mJy/beam. The contours range from 2 to 22 in steps of 5 mJy/beam, and from 22 to 62 in steps of 10 mJy/beam. The yellow ellipse represents the ALMA beam size.

Open with DEXTER
In the text
thumbnail Fig. 2

ALMA CH3OH(13−12) integrated-intensity image of SDC335 in colour scale; overplotted are the contours of the 3.2 mm dust continuum emission as displayed in Fig. 1c. We can see that the two types of emission coincide spatially. The insert in the top left corner shows the methanol spectra observed at the central position of each core. The cyan solid lines are the best Gaussian fit to the data. The yellow ellipse represents the ALMA beam size.

Open with DEXTER
In the text
thumbnail Fig. 3

Spitzer 8 μm image of SDC335 (colour scale) over-plotted with the Mopra HCO+(1−0) spectra. The temperature scale and velocity are indicated in the bottom-left corner. The HCO+(1−0) line is self-absorbed and blue-shifted in the bulk of the cloud. This is usually interpreted as a signature of collapse.

Open with DEXTER
In the text
thumbnail Fig. 4

a) Same as in Fig. 1a; b) ALMA-only image of the N2H+(1−0) integrated intensity over the 7 hyperfine structure components. The rms noise on the resulting map is ~ 6 mJy/beam km s-1. The contours go from 0.1 to 1.5 in steps of 0.7 Jy/beam km s-1 and 1.5 to 9 in steps of 1.5 Jy/beam km s-1. The crosses mark the positions of the two dense cores. The ALMA beam is represented as a yellow elliptical symbol in the bottom-right corner of the image. We can see the excellent match between the Spitzer dust extinction of the filaments and the N2H+(1−0) emission; c) ALMA N2H+(1−0) velocity field using the first order moment map. The crosses mark the positions of the cores and the contours are the same as in the b) panel.

Open with DEXTER
In the text
thumbnail Fig. 5

(Left) Examples of combined ALMA and Mopra N2H+(1−0) spectra observed at some specific positions in SDC335 (upper panel for filaments − lower panel for cores), along with their best fits as red solid lines. The N2H+(1−0) spectra exhibit a hyperfine structure (HFS) composed of 7 components (the positions are displayed as blue vertical ticks for the F2 filament). Some of these components are close enough to be blended when the velocity dispersion of the gas is supersonic, resulting in three groups of lines. For a kinetic temperature of 10 K, the velocity dispersion along the filaments is supersonic by a factor 1.5 to 3, similar to what is observed in other IRDCs (Ragan et al. 2012). (Right) Schematic representation of the systemic velocity and velocity dispersion of the different structures. The length of each box represents the velocity dispersion (FWHM) of the gas, and its central position the systemic velocity (represented as a filled circle). The colour of the boxes codes the line that has been used for the measurements: red for N2H+(1−0), and cyan for CH3OH(13−12). The vertical dashed lines mark the systemic velocities of the cores.

Open with DEXTER
In the text
thumbnail Fig. 6

Mass-radius diagram of massive protostellar cores. The cyan circles correspond to the values as published in the literature (Peretto et al. 2007; Ren et al. 2012; Rathborne et al. 2011; Wang et al. 2011; Zhang & Wang 2011; Bontemps et al. 2010; Rathborne et al. 2007, 2008; Beuther et al. 2002, 2003, 2005, 2012; Molinari et al. 1998; Rodón et al. 2012), while the empty circles correspond to the same set of sources for which we recalculated the mass using the same dust opacity law as in this paper. The MM1 and MM2 sources are indicated as red filled circles. The green filled square marks the position of the W51 North star+disc system from Zapata et al. (2009). The shaded area indicates the region where most sources lie. The cross in the bottom right corner indicates a factor of 2 uncertainty in both the masses and sizes, typical of the results presented here.

Open with DEXTER
In the text
thumbnail Fig. 7

Grid of HCO+(1−0) spectra obtained from RATRAN modelling of a collapsing cloud (see text). All input parameters are fixed with the exception of the infall velocity (Vinf) and velocity dispersion (σturb). Each modelled spectrum (in red) has been obtained for the corresponding Vinf − σturb displayed in the top and right-hand sides of the figure. The numbers in the top right-hand corner of each spectrum give a relative idea of the fit quality. The HCO+(1−0) spectrum observed at the centre of SDC335 is plotted in black. Since we aimed to keep the observed spectra displayed in scale, we applied a 0.5 factor to the modelled spectra to take into account the Mopra main-beam efficiency.

Open with DEXTER
In the text
thumbnail Fig. 8

Snapshot of a MHD simulation of a 10 000 M collapsing cloud (see Appendix C for more details; Schneider et al. 2010). a) Column density (colour and contours) smoothed to the resolution of the ALMA data (5′′). The arrows show the plane-of-the-sky velocity field. We see that gas flows along filaments and also off the filaments. b) Dense gas line-of-sight velocity field (colour scale) smoothed to the resolution of the ALMA data. We highlighted the filaments by white dashed lines. The contours are the same as in panel a).

Open with DEXTER
In the text
thumbnail Fig. B.1

Same as Fig. 7, but for the H13CO+(1−0) line.

Open with DEXTER
In the text