Free Access
Issue
A&A
Volume 594, October 2016
Article Number A118
Number of page(s) 13
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201527062
Published online 20 October 2016

© ESO, 2016

1. Introduction

A complete understanding of the relationship between the distribution of density peaks in regions of high-mass star formation (i.e. the core mass function, CMF) and the distribution of stellar masses in clusters, the initial mass function, IMF (e.g. Könyves et al. 2010; Alves et al. 2007) may provide an insight in to the mechanisms responsible for the formation of massive stars (e.g. Goodwin et al. 2007).

In order to compare the CMF and IMF in a star forming region, we must be able to accurately characterise its substructure i.e. the level of fragmentation it contains. Studying fragmentation in high mass star forming regions requires high resolution observations of the earliest stages of massive star formation (e.g. Beuther & Schilke 2004; Zhang et al. 2009; Swift 2009; Bontemps et al. 2010; Wang et al. 2011, 2014). The earliest evolutionary stages of high mass protostars are difficult to identify, but have characteristics consistent with those seen in massive infrared dark clouds (IRDCs). These IRDCs are cold (1020 K; Pillai et al. 2006; Ragan et al. 2011), dense regions within giant molecular clouds (with column densities >1022 cm-2; Peretto & Fuller 2009), manifesting themselves as regions of extinction against the mid-IR emission from the galactic plane (e.g. Peretto & Fuller 2009). Here we investigate sub-fragmentation in SDC13, a region comprising 3 Spitzer IRDCs from the Peretto & Fuller (2009) catalogue (SDC13.174-0.07, SDC13.158-0.073, SDC13.194-0.073) at a distance of 3.6 kpc (Peretto et al. 2014).

Peretto & Fuller (2009) use the term fragment to describe local peaks in column density between contours of 8 μm opacity for the IRDCs in their catalogue (see Appendix A in Peretto & Fuller 2009). They find 18 fragments in extinction in SDC13. IRAM 30 m MAMBO (Max-Planck-Millimeter-Bolometer) 1.2 mm dust continuum observations towards SDC13 indicate fragments of a similar size and position to those seen in extinction, however, the map is dominated by two comparatively large fragments, MM1 and MM2 (see Table 2 below and Peretto et al. 2014). Based on the size of MM1 and MM2 (~a few times 0.1 pc), we consider the term fragment to be analogous to the term core, and will refer to MM1 and MM2 as such. MM1 is not seen in extinction and is associated with 8 μm emission, indicative of active star formation. MM2 shows no evidence of star formation activity, lacking both 8 μm emission and 24 μm emission (associated with warm dust), making it a good candidate for a massive prestellar core (Fig. 1).

thumbnail Fig. 1

A three-colour Spitzer image of SDC13 (24 μm in red, 8 μm in green and 3.6 μm in blue) overlayed with IRAM MAMBO 1.2 mm continuum contours (white) at 5, 15, 25, 35 and 45 mJy (left). The highlighted region (right) shows the two largest fragments (i.e. cores) in SDC13, MM1 and MM2, whose positions are marked with white crosses. Four sub-fragments are seen in the SMA 1.3 mm continuum observations of the region. The positions of these sub-fragments, A, B, C and D, are marked with yellow circles.

Open with DEXTER

In this paper we investigate the substructure in MM1 and MM2 using high angular resolution (<3′′) observations at millimetre wavelengths (Sect. 2). This is achieved using the Submillimeter Array (SMA1; Ho et al. 2004), an 8-element radio interferometer located at the summit of Mauna Kea in Hawaii, to obtain 1.3 mm continuum observations of MM1 and MM2 (Sect. 3). The aim is to gain an insight into the mechanisms responsible for massive star formation. As MM2 appears to be in an earlier stage of evolution than MM1, a comparison of the substructure of MM1 and MM2 may also provide an insight into the early evolution of massive stars and star clusters.

2. SMA observations and data reduction

Observations were performed using 6 antennas of the SMA at 230 GHz. At this frequency the FWHM of the primary beam is ~55′′. Both the extended array configuration and compact array configuration were utilised on the 8th March 2012 and the 30th June 2012 respectively. An overview of the observing parameters and the maximum spatial scales that each configuration is sensitive to are given in Table 1. We observed two overlapping fields, centred on MM1 and MM2.

The calibration of the visibility data was performed using MIR, a software package written in IDL for the purpose of reducing SMA data. A time-dependent phase and gain calibration was carried out using quasars 1733130 and 1743038 for observations in both the compact and extended configurations. Observations in the lower sideband (LSB) cover the frequency range 216.8220.8 GHz, and in the upper sideband (USB) 228.8232.8 GHz. There is a uniform spectral resolution of 0.84 MHz (~1.1 km s-1).

For observations in the extended configuration, the quasar 3c279 was used as the band-pass calibrator; Mars was used to calibrate the flux, and the observed system temperature (Tsys) varies from 150200 K. For observations in the compact configuration, Uranus was used to calibrate both the band-pass and the flux, and 120 K <Tsys< 220 K.

After calibration, the visibility data were exported to the radio interferometry data reduction package MIRIAD (Multichannel Image Reconstruction, Image Analysis and Display, Sault et al. 1995) for further processing and image production. Spectral line and continuum data were separated, and line-free continuum data from both antenna configurations in both sidebands were combined for each pointing. In this paper we present results from continuum observations only.

The continuum map shown in Fig. 2 is a linear mosaic of the pointings towards MM1 and MM2, corrected for primary beam attenuation. The mosaic was produced using the linmos command in MIRIAD with the taper option selected. This attempts to counteract excessive noise amplification at the edge of the mosaic to produce uniform noise across the whole image (see Sault et al. 1996). The data used combines observations in the upper and lower side band for both the compact and extended configurations for each pointing. With natural weighting, the map has a synthesised beam with dimensions 3.73′′ × 2.52′′, position angle (PA) 41.55°, and a 1σ rms of ~1 mJy.

Table 1

An overview of observing parameters and the largest spatial scales that the SMA is sensitive to in the configurations used.

3. SMA 1.3 mm continuum image

The 1.3 mm continuum emission towards MM1 and MM2 is shown in Fig. 2. We use the term sub-fragment to describe the substructures we find within MM1 and MM2 (which is analogous to the term condensation, often used to describe the 0.01 pc-scale structures within cores). If we allow sub-fragment boundaries to be defined by 3σ contours, we find that MM1 and MM2 consist of a total of four sub-fragments, which we have labelled A, B, C, D. Table 2 gives the J2000 coordinates and some of the physical properties calculated for A, B, C and D.

Sub-fragment dimensions a × b of subfragments A and D were estimated based on the 3σ contours in the 1.3 mm continuum map (Fig. 2). The position angles (PA) of A and D were determined by fitting each with a 2D gaussian using the Pick Object function in GAIA (Graphical Astronomy and Image Analysis Tool, part of the Starlink astronomical software package). For the weaker subfragments C and D, we estimate diameters to be 4 ± 1′′.

The integrated and maximum flux value for each sub-fragment was determined by performing aperture photometry on the 1.3 mm continuum image. This was carried out using the imstat function in MIRIAD. Assuming a typical dust temperature of 15 K (e.g. see Rathborne et al. 2007; Peretto et al. 2014), the mass for each sub-fragment can then be estimated from its integrated flux using: (1)where Fν = the integrated flux at frequency ν, d is the distance to the source (3.6 kpc) and Bν = the Planck function at a dust temperature Tdust. The opacity at ν is calculated using (2)(Hildebrand 1983), with β = 1.5 (see Wang et al. 2011), such that κ1.3 mm = 0.8 cm2 g-1. The masses of MM1 and MM2 given in Table 4 of Peretto et al. (2014) were calculated from the 1.2 mm MAMBO integrated flux, assuming a temperature of 15 K for MM1 and 12 K for MM2, and an opacity κ1.2 mm = 0.5 cm2 g-1. Thus we apply conversion factors of 0.6 and 0.5 respectively, to obtain the masses for MM1 and MM2 as quoted in Table 2, allowing a more accurate comparison to the masses we calculate for A, B, C and D.

We obtain masses of 46.8 M and 40.6 M for MM1 and MM2 respectively. These are similar to average core masses found in previous studies of massive star forming regions (e.g. Motte et al. 2007; Rathborne et al. 2006; Zhang et al. 2009). The masses of sub-fragments A, B, C and D in SDC13 range from ~212 M. These should be considered lower limits as they have not been rescaled to take in to account the filtering of extended emission at 1.3 mm (e.g. see Bontemps et al. 2010; Duarte-Cabral et al. 2013), although this is not likely to be a large effect as the sources are close to unresolved. The sub-fragment masses we obtain are in agreement with those obtained by Wang et al. (2011) for sub-fragments in the star forming region G28.34 (1.410.6 M).

In Fig. 3, the contours from the SMA 1.3 mm continuum observations are overlayed on single dish MAMBO 1.2 mm continuum data, obtained using the IRAM 30 m telescope (Peretto et al. 2014).

Three of the sub-fragments A, C, and D appear to be associated with MM1 and MM2. Sub-fragment B appears to be associated with the slight extension between MM1 and MM2. However, whilst it has a mass approaching that of the brightest sub-fragment A (Table 2), it is much less prominent in the MAMBO data, indicating that it is not associated with a significant extended envelope of emission. In addition, while it appears that sub-fragment B is close to the position of source MM18 in Fig. 1 of Peretto et al. (2014), the distance between the peak of B and MM18 is larger than its 2′′ position uncertainty. Neither does it correspond to any peak in JVLA NH3 data (Williams et al., in prep.). In Sect. 4 we investigate possible candidates to explain the appearance of sub-fragment B.

Two of the sub-fragments A and C coincide with MM1, with the brightest of the two, sub-fragment A, at its peak. This could be indicative of fragmentation in MM1. We further investigate this possibility using images from the Herschel Infrared Galactic Plane Survey (Hi-GAL, Molinari et al. 2010, Sect. 5).

Sub-fragment D coincides with the peak of MM2, however, whereas MM1 and MM2 are of similar size and flux at 1.2 mm (Table 2), sub-fragment D is ~5 times fainter than sub-fragment A at 1.3 mm. An object similar to MM2 is seen in the molecular cloud complex Cygnus X. In their MAMBO survey of Cygnus X, Motte et al. (2007) find a 1.2 mm mass of 100 M for the massive cloud core CygX-N40, yet in further studies using the higher resolution of the PdBI (Bontemps et al. 2010; Duarte-Cabral et al. 2013) CygX-N40 is barely detectable. This result is attributed to the filtering of extended emission at 1.3 mm. The absence of a star in MM2 indicates that it is in an earlier stage of evolution than MM1 which might suggest a more diffuse physical structure, providing a possible explanation for the missing flux at 1.3 mm.

We aim to further understand the physical differences that give rise to such observational differences at higher resolutions. We use RADMC-3D, a software package designed for astrophysical radiative transfer calculations in 1D, 2D and 3D geometries (Dullemond 2012) to model MAMBO 1.2 mm observations of MM1 and MM2 (Sect. 6) and use the CASA data reduction software (McMullin et al. 2007) to process the resulting images and simulate our SMA 1.3 mm observations (Sect. 7).

thumbnail Fig. 2

The primary beam corrected linear mosaic of the 1.3 mm continuum data obtained in two SMA pointings towards MM1 and MM2 (Jy/beam). The central portion, outlined in blue, indicates the region shown in Fig. 3. The synthesised beam is shown in the bottom left hand corner of the image. Contours are at − 3σ,3σ,5σ,7σ,9σ where σ ≈ 1 mJy/beam. Positive contours are indicated by the solid lines and negative contours by the dashed lines. Letters A, B, C and D label the 4 largest sub-fragments seen in the data, defined by the 3σ contours.

Open with DEXTER

Table 2

Physical properties and J2000 coordinates of MM1 and MM2, and sub-fragments A, B, C and D.

thumbnail Fig. 3

SMA 1.3 mm continuum data overlayed on 1.2 mm MAMBO continuum data (mJy/beam). The contours indicate the SMA continuum at 3, 3, 5,7 and 9 mJy/beam. Positive contours are shown in pink and negative contours are shown dashed in green. Letters A, B, C, D indicate the 4 sub-fragments seen in the region which are identified as peaks in the 1.3 mm SMA continuum. The MAMBO beamsize is 10.7′′.

Open with DEXTER

4. Sub-fragment B

Sub-fragment B, seen in the SMA data, has a mass (based on its 1.3 mm flux) approaching that of sub-fragment A, and yet it is much less prominent in the MAMBO data, indicating that it is not associated with a significant extended envelope of emission. We consider the possibility that sub-fragment B is a background source, an AGB star, the free-free emission of a HII region, or a runaway object.

Maloney et al. (2005) performed a fluctuation analysis on data from the 1.1 mm Bolocam Lockman Hole Survey in order to constrain the slope and amplitude of the number count distribution of high redshift galaxies at λ = 1.1 mm. They find the best-fitting power-law model to have an index δ = 2.7, a differential number density at 1 mJy n0 ≈ 1595 mJy-1 deg2, and an integrated number density N(>1 mJy) ≈ 940. These three parameters are related by, (3)where S is the peak flux density. The peak flux density of B is SB ~12 mJy (Table 2). Applying Eq. (3), we find that N(>12 mJy) 13.7. Based on the SMA beamsize, and the overlap of the two pointings, our observations cover an area of 0.0003 deg2. The probability of finding a background source within this region is ~0.004. Therefore, it is unlikely that B is background source.

Object B is not listed in any AGB catalogues and, since its angular diameter (see Table 2) makes it too large to be an AGB star unless it resides within a few parsecs of the Earth (e.g. Villaver & Livio 2007), it is highly unlikely that it is yet to be identified. Using the surface density of dust envelope AGB stars in the solar neighbourhood, calculated by Olivier et al. (2001), we calculate the probability of finding an AGB star within the region observed to be ~5.7 × 10-5.

In order to determine whether B could be the result of emission from a HII region, we looked at a 5GHz radio continuum image of the SDC13 region obtained by the VLA CORNISH survey. The rms noise in the images is <0.4 mJy beam-1, which would allow detection of an unresolved UCHII region around a B0 star at 16 kpc (Purcell et al. 2013). Source B is not associated with any CORNISH radio source, so it can not explained by emission from a HII region.

One further possibility is that B is a runaway object that has been ejected from MM1, and its continuum emission arises from a circumstellar disk. About 40% of O stars and ~10% of B stars are thought to be runaways, and the ejection of stellar embryos from star forming regions, via many body interactions in clustered environments (the dynamical ejection scenario, DES; Poveda et al. 1967), has been proposed as a method of brown dwarf formation (e.g. Bate 2009). Assuming a DES for B, we can estimate its potential ejection velocity based on its distance from the centre of MM1 (i.e. its distance from sub-fragment A) and the approximate time since ejection. A lower limit on the distance between A and B can be obtained from their observed separation on the sky. This is ~28′′, equivalent to ~1.5 × 1018 cm at a distance of 3.6 kpc. An upper limit on the time since ejection is equivalent to their formation timescale (~105 years for massive prestellar/protostellar objects). Based on these values, we thus calculate a lower limit for the hypothetical ejection velocity of B of ~5 km s-1. The expected escape velocities for O/B runaways are ~200400 km s-1 (e.g. Gvaramadze et al. 2009). Based on our calculated lower limit for the velocity of sub-fragment B it is plausible that it B is a runaway object.

5. Herschel Hi-GAL observations

Figure 5 shows Hi-GAL images at 350, 250, 160 and 70 μm, observed towards SDC13. The 160 μm Hi-GAL flux peaks at A, and the 70 μm flux peaks over A and C, indicative of star formation activity in this region. Conversely, sub-fragment D is dark in the Hi-GAL 70 μm images. This would suggest that D does not contain any embedded sources and may therefore be in an earlier stage of evolution than A and C.

The SMA continuum data (Fig. 2) indicates that MM1 consists of two sub-fragments (A and C). However, A and C appear unresolved in the Herschel observations at all wavelengths. At 70 μm, Hi-GAL observations show an elongated source covering the positions of A and C. The resolution of the image is insufficient, however, to confirm whether this elongated source is the result of emission from two objects.

Source extraction was performed at all Hi-GAL wavelengths with the Hyper (Hybrid Photometry and Extraction Routine) algorithm (Traficante et al. 2015). Hyper is an enhanced version of classical aperture photometry, designed to take into account the strong background variability and source crowding typical of Galactic observations. Sources are modelled with 2D Gaussians, allowing the FWHM to vary between 1 and 2 times the instrumental PSF. The background is subtracted automatically by Hyper. It estimates several backgrounds, modelled with polynomials of different orders, and chooses the best background model based on which results in the lowest residuals. The flux is estimated within a 2D Gaussian region with aperture radii equal to the Gaussians FWHM (Traficante et al. 2015). Hyper identified 2 sources at 70 μm, Her1 and Her2. The results of the source extraction and photometry for Her1 and Her2 are shown in Table 3.

The initial parameters for our RADMC-3D models of MM1 and MM2 (Sect. 6) are derived based on an SED fit to Hi-GAL observations towards MM1 at 160, 250 and 350 μm, assuming a distance to SDC13 of d ≃ 3.6 kpc. We use an elliptical aperture with axes equal to the FWHMs derived from a 2D Gaussian fit at 250 μm. The SED fit is performed using a single-temperature greybody model with fixed β = 1.5, and temperature and mass as free parameters. The results of the fit are shown in Fig. 4. We find an MM1 luminosity LHi − GAL = 1080 L, radius RHi − GAL = 1.93 × 1018 cm (~0.6 pc) and core mass MHi − GAL = 243 M (equivalent to a dust mass Mdust = 2.43 M, assuming a gas:dust mass ratio of 100:1). We used the physical properties of MM1 derived from this SED to model MAMBO observations of MM1 and MM2 (Sect. 6).

Table 3

Photometry results and J2000 coordinates for sources Her1 and Her2.

6. Modelling with RADMC-3D

Using RADMC-3D, we model the dust continuum emission of MM1 and MM2 with 1D logarithmically spaced, spherically symmetric models with inner radii Rin = 1 × 1016 cm (based on typical values for best-fit models to high-mass cores given in Table 5 of Williams et al. 2005). This corresponds to ~0.2′′ at the distance of SDC13. Our initial estimates for the outer radius (Rout = 1.93 × 1018 cm) of our model dust clouds and the total dust mass (Mdust = 2.4 M) are based on the SED fit to the Hi-GAL data for MM1 (Sect. 5). We use dust opacities calculated by Ossenkopf & Henning (1994) for coagulated dust grains with thin ice mantles, at a gas density of 105 cm-3 and a Draine & Lee (1984) dust grain size distribution.

Spitzer MIPSGAL observations of the IRDC SDC13 indicate a 24 μm source at the centre of MM1, indicating the presence of a central star (or multiple system). There is no evidence of such a source at the centre of MM2 (Fig. 1). The effect of including a stellar source at the centre of our models is therefore investigated. We cannot rule out the presence of a multiple system at the centre of MM1, however for simplicity we use a single object at the centre of our models. We base the luminosity of the central star (L = 1080 L) on the luminosity of MM1, which is derived from the Hi-GAL data (Sect. 5). We estimate the contribution of the dust emission to the total luminosity (in the absence of a star) using a modified blackbody function (e.g. see Battersby et al. 2011) at a temperature of 15 K, normalised to the MM1 1.2 mm MAMBO flux. We find a contribution of ~14 L from dust emission, which is negligible compared to the total luminosity. The remaining stellar properties (mass M = 7.36 M, radius R = 4.94 R and temperature T = 14 914 K, corresponding to a class B ZAMS star) are calculated using typical relationships between the luminosity, mass, radius and temperature of main sequence stars (Schulz 2012).

For prestellar cores (without a central source), the dominant source of radiation will be external. To account for this, we include an external radiation field in our models based on the ISRF in the solar neighbourhood. Our models cover wavelengths 0.01 ≤ λ ≤ 1000 μm, therefore our model ISRF covers emission in this range, incorporating contributions from the cosmic microwave background (CMB), infrared emission from dust, and photons of starlight (Draine 2011).

thumbnail Fig. 4

SED fit to Herschel Hi-GAL data observed towards MM1 at 300, 250, 160, 70 μm, based on a distance of 3.6 kpc and using β = 1.5.

Open with DEXTER

The earliest models of core collapse describe isolated spherical cores which are pressure bounded, isothermal and non-fragmenting i.e. Bonnor-Ebert (BE) spheres (Bonnor 1956; Ebert 1955). For a critical BE sphere on the verge of gravitational collapse, there are a family of solutions corresponding to spheres with different initial density distributions. Each consist of a uniform density central core (whose size and density varies with temperature) and an outer envelope with density ρr-2 (where r is the radius). At one extreme of this group of models is the Shu (1977) solution, which describes the inside-out collapse of a singular isothermal sphere (SIS) i.e. a sphere with infinite central density and ρr-2. At the other extreme of the family of BE spheres, the Larson-Penston solution (Larson 1969; Penston 1969) describes the collapse of a cloud with uniform density, that acquires large infall velocities at all radii.

During the inside-out collapse of a SIS, described by Shu et al. (1987), the density increases fastest in central regions, with collapse occurring first at the centre and then propagating outwards via an expansion wave (Shu 1977). The result of the collapse is an infalling central region in free-fall collapse (ρr-1.5) surrounded by a static envelope (retaining the isothermal density distribution, ρr-2). Other theoretical models invoke different density profiles at various points during the process of collapse (e.g. Foster & Chevalier 1993; Basu & Mouschovias 1994; Whitworth & Ward-Thompson 2001), however, Foster & Chevalier (1993) found that following the formation of the central protostar, density profiles tend towards those seen for SISs. More recently, based on numerical simulations of core and protostar formation in supersonic flows, Gong & Ostriker (2009, 2011) find the ρr-2 density profile to be an attractor for the collapse of a molecular cloud core at the point of protostar formation, regardless of the mechanism responsible for the collapse. Similar results have been found for models of massive cores (Tilley & Pudritz 2004).

Observations of low mass star forming regions (e.g. Young et al. 2003) show core density profiles ρrα with − 2.0 ≤ α ≤-1.5, in agreement with the predictions of models of collapsing cores (e.g. Shu et al. 1987; Foster & Chevalier 1993). However the density distributions of high mass cores have been more difficult to characterise due to the relative scarcity of massive cores and their more complex environments. Beuther et al. (2002) presented a study on a large sample of massive star forming regions and found a typical value of α = − 1.6 for their density distributions. Additional previous modelling and observations of high mass star forming regions have found a range of density distributions with -2.25≤ α ≤ 0 (e.g. Wolfire & Churchwell 1994; Garay & Rodriguez 1990; Hatchell & van der Tak 2003; Williams et al. 2005). As such, we model our spherical dust clouds with density distributions ρrα (where -3.0<α ≤ 0). The upper limit on α is a consequence of our mass normalization (see Eq. (4)).

Figure 6 compares the MAMBO 1.2 mm continuum flux distributions observed towards MM1 and MM2 for a range of power law density distributions. Using the initial Hi-GAL derived estimate for the core radius Rout = RHi − GAL = 1.93 × 1018 cm, we obtain good fit to the slope of the MAMBO 1.2 mm observations of MM1 using a model with a central source and a density distribution ρr-1.5 (model A2a). This is similar to the observed density profiles of protostellar cores in low-mass star forming regions, and is suggestive of free-fall collapse (Shu 1977; Shu et al. 1987), consistent with infall around the central star seen in the Spitzer data (Fig. 1). In order to also fit the peak of the MM1 observations using this density distribution, we require a model dust mass (i.e. the mass of dust only) Mdust = 2.33 M.

We initially modelled MM2 as a dust envelope without a central source, with Rout = RHi − GAL. A reasonable fit to the slope of the MM2 MAMBO data is then obtained for a model with a power law density distribution ρrα, with − 2.9 <α< − 2.5 (models D3a and D4a respectively). This is steeper than expected for an SIS collapse scenario. The best fit models to MM2 have dust masses in the range 6.17 M<Mdust ≤ 6.49 M. Inspection of the MIPSGAL images suggests there may be some heating due to an enhanced radiation field in this region (Fig. 1). In principle this would lower the mass of MM2, but the effect would be subtle. This heating may however account for the spread in the radial profile of MM2 (e.g. see grey lines in Fig. 6). The temperature and density profiles of the best fit models (with Rout = RHi − GAL) to MAMBO observations of MM1 and MM2 are shown in Fig. 7. This figure shows that the temperature profile of the best fit model to MM2 drops down to ~4 − 5 K at its centre. While temperatures as low as 6 K have been seen in some sources (e.g. Harju et al. 2008), massive star forming regions and IRDCs are typically observed to have temperatures of ~10 − 15 K. We therefore investigated the effect of setting a lower temperature profile limit of 10 K on the resulting RADMC-3D model density profiles. This temperature limit could be justified by the effect of cosmic ray heating, for example. We found than while the 10 K lower temperature limit did result in a shallower density distribution for the best fit model to MM2 (α ~ 2.0) this is still steeper than the density profile found for MM1. This fixed minimum temperature did however result in a decreased model mass of ~230 M for MM2.

In an attempt to further improve the fit to the MAMBO observations of MM2 we tried and successfully fitted a number of different models to the MM2 data which correspond to a variety of physical conditions. These included models with truncated density distributions (expected for cold cores embedded in warmer gas; Fischera 2014); and two-part power law density profiles, consistent with observations of central flattening in prestellar cores (e.g. Beuther et al. 2002; Andre et al. 2000). The SMA data provides additional constraints on these models. In Sect. 7, we describe the use of the CASA data reduction software (McMullin et al. 2007) to process our model images and produce SMA 1.3 mm simulations to investigate if we are able to reproduce our 1.3 mm observations (Fig. 2).

thumbnail Fig. 5

SMA 1.3 mm continuum data overlayed on Herschel Hi-GAL images (Jy/beam) at 350 μm (top left), 250 μm (top right), 160 μm (bottom left) and 70 μm (bottom right). The contours indicate the SMA 1.3 mm continuum at 3, 5, 7 and 9 mJy/beam. Letters A and C indicate 2 of the 4 sub-fragments seen in the region which coincide with MM1. The two sources, Her1 and Her2, indicated in blue, were extracted from the 70 μm Herschel data using the Hyper algorithm (Traficante et al. 2015, see Table 3).

Open with DEXTER

thumbnail Fig. 6

Comparison of the observed MAMBO 1.2 mm flux distributions of MM1 (left) and MM2 (bottom) with RADMC-3D model fluxes, obtained using various model parameters. The solid black lines indicate the observed fluxes plotted against the equivalent radius (where Anσ is the area contined within contours , for n = 3, 4, 6, 7, 8, 9 and σ ~ 5 mJy/beam). The grey lines indicate perpendicular slices of observed flux against radius; these slices give an indication of the range of variation in the observed flux distributions in MM1 and MM2. The coloured lines show the model fluxes against radius for RADMC-3D models with various dust distributions, as indicated in the key. The model dust masses Mdust required to match each peak model flux with the peak observed flux are also given. Models are labelled A1a to A3a (for those fit to MM1 data) and D1a to D4a (for those fit to MM2 data). The dashed curves indicate the MAMBO beam (10.7′′).

Open with DEXTER

thumbnail Fig. 7

Density and temperture profiles of the best fit models (with Rout = RHi − GAL) to MAMBO observations of MM1 (left) and MM2 (right; models A2a and D3a respectively).

Open with DEXTER

7. CASA simulation of SMA observations

Initially, the RADMC-3D models of MM1 and MM2 were added linearly into a blank field to create a single sky model using the CASA imaging toolkit functions. The CASA task simobserve was then used to generate visibility data. This created four measurement sets, one for each SMA array configuration at each sideband. The date and hourangle were set such that the UTC and date for each observation matched the true observations. The zenith opacity was set to give Tsys values for the simulations which matched the highest values measured in the SMA observations with each array. The on source time, comprising multiple scans toward two pointings, matched the observed data in length. An artificial phase calibration cycle was included such that the time coverage of the uv-plane was as close as possible to real observations. The measurement sets (MS) for each sideband for each array were then concatenated together, giving a single MS.

Finally, imaging of the visibility data was carried out using the CASA clean task in mosaic mode, with a clean box at each source position and a threshold of ~0.45 mJy, (a mid point value between the ~0.3 mJy for the compact data, and ~0.6 mJy for the extended data). The cell size (0.52′′ × 0.50′′) and the restoring beam (3.78′′ × 2.52′′, PA = 41.55°) were fixed to match the real data image values. Examples of the results of the CASA simulation process are shown in Figs. 8 and 9.

Figure 8 shows the result of the CASA simulation using models A2a and D3a (see Sect. 6 and Table 4). The resulting sub-fragments A2a and D3a are of similar size and flux, with peak fluxes FA2a = 1.68 mJy/beam and FD3a = 1.52 mJy/beam respectively (i.e FA2a ~ 10% larger than FD3a). In comparison, the peak fluxes of sub-fragments A and D in the SMA 1.3 mm data (coincident with MM1 and MM2 in the MAMBO data) are FA = 13.6 mJy/beam and FD = 4.9 mJy/beam respectively i.e FA ~ 2.8 times larger than FD (see Table 2). We find we were able to better match the observed peak fluxes FA and FD (and therefore their ratio FA/FD) by adjusting the value of Rout in the RADMC-3D models of MM1 and MM2.

We normalize our RADMC-3D models based on the total dust mass Mdust, which we use to calculate a normalisation constant ρ0 (the density at the outer radius Rout). For a density distribution ρ = cRα and . The density distribution can thus be expressed as ρ = ρ0(R/Rout)α. Based on the dust mass (4)we find the normalisation constant (5)A consequence of this method of normalisation is that MdustRout i.e. an increase/decrease in Rout in our models requires an increase/decrease in Mdust to maintain the same density distribution. This increase/decrease in Mdust results in an increase/decrease in flux from the model.

For example, as shown in Fig. 9, we are better able to match the observed fluxes FA and FD with increased outer radii for the RADMC-3D models of MM1 and MM2 (Model A2b with RA2b = 3.5 × 1018 cm and Model D3b with RD3b = 2.9 × 1018 cm respectively; see Table 4). This subsequently increases the model peak fluxes, giving values FA2b = 10.5 mJy and FD3b = 4.2 mJy resulting in a flux ratio FA2b/FD3b ~ 2.5. This is similar to the ratio of the observed fluxes for sub-fragments A and D, FA/FD ~ 2.8. Figure 9 indicates that adjusting the model outer radii as described above does not have a significant effect on the quality of the fit to the MAMBO observations.

Simulations performed using models of MM2 with truncated and two-part power-law density profiles produce model sub-fragments that are too diffuse to be detected. As such, we consider models A2b and D3b to be the best-fits for MM1 and MM2 respectively, as they provide the closest match to both the MAMBO and SMA observations. The best fit model to MM1 is therefore a cloud of total mass ~360 M (assuming gas:dust = 100:1), radius Rout = RA2b = 3.5 × 1018 cm and density profile rρ-1.5, with a central star of temperature T = 14 914 K. The best fit model to MM2 has a total mass ~824 M, radius Rout = RD3b = 2.9 × 1018 cm and density profile rρ-2.5 and does not contain a central star. Peretto et al. (2014) find masses of 74.8 M in 0.26 pc and 81.1 M in 0.21 pc for MM1 and MM2 respectively (see Table 2). Within the same radii, our best-fit models for MM1 and MM2 contain masses 41.3 M and 254 M respectively.

Table 4

Parameters of RADMC-3D models shown in Figs. 8 and 9 and their corresponding CASA simulated peak flux values.

thumbnail Fig. 8

CASA simulated SMA 1.3 mm continuum map of the SDC13 region (bottom) made using RADMC-3D models A2a (top) and D3a (centre; see Table 4). The map contour is at 1.5 mJy. The rms noise in the centre of the map ~0.09 mJy. Model fluxes for A2a (shown in red) and D3a (shown in green) are compared with plots of the observed MAMBO 1.2 mm flux against equivalent radius for MM1 and MM2 respectively (solid black lines). Model parameters are given in the key of each plot. The grey lines in the top two plots indicate perpendicular slices of observed flux against radius; these slices give an indication of the range of variation in the observed flux distributions in MM1 and MM2. The dashed curves indicate the MAMBO beam (FWHM = 10.7′′).

Open with DEXTER

thumbnail Fig. 9

CASA simulated SMA 1.3 mm continuum map of the SDC13 region (bottom) made using RADMC-3D models A2b (top) and D3b (centre; see Table 4). Map contours are at 2, 4, 6, 8 and 10 mJy. The rms noise in the centre of the map ~0.3 mJy. Model fluxes for A2b (shown in red) and D3b (shown in green) are compared with plots of the observed MAMBO 1.2 mm flux against equivalent radius for MM1 and MM2 respectively (solid black lines). Model parameters are given in the key of each plot, with further details given in the text. The grey lines in the top two plots indicate perpendicular slices of observed flux against radius; these slices give an indication of variation in the observed flux distributions in MM1 and MM2. The dashed curves indicate the MAMBO beam (FWHM = 10.7′′).

Open with DEXTER

8. Discussion and summary

High resolution SMA 1.3 mm observations toward MM1 and MM2, the two largest fragments (i.e. cores) in SDC13 (as observed with MAMBO at 1.2 mm with lower resolution), indicate the presence of 4 sub-fragments A, B, C and D (Fig. 2). Three of the sub-fragments A, C and D are associated with the cloud. One of the sub-fragments, B, does not appear to be associated with a significant extended envelope and is much less prominent in the MAMBO 1.2 mm observations (Fig. 3). The nature of B remains unconfirmed, but statistical analysis indicates that B is unlikely to be a background source, an AGB star or free-free emission of a HII region. However B could plausibly be a runaway object (see Sect. 4).

MM1 consists of two sub-fragments: A, at its centre and a fainter sub-fragment C. The fragmentation of MM1 into two objects is confirmed by Herschel Hi-GAL observations at 70 μm (see Fig. 5). MM2 consists of only one object at its centre (sub-fragment D). Although MM1 and MM2 look similar at the lower MAMBO resolution, with a similar size and brightness, at the higher resolution of the SMA, they look very different (Fig. 3), with the object at the centre of MM1 (sub-fragment A) appearing much brighter than the object at the centre of MM2 (sub-fragment D). We are able to reproduce this effect with CASA simulations of SMA observations, using RADMC-3D models with single power law density distributions and extended outer radii (see Fig. 9, Sect. 7).

MM1 is associated with 8 μm emission indicative of active star formation. MM2 shows no evidence of emission at 8 μm or 24 μm making it a good candidate for a prestellar core (Fig. 1). MM2 requires a steeper density profile and higher mass to model its emission (rρ-2.5; Mdust ~ 824 M) compared to those required for MM1 (rρ-1.5; Mdust ~ 360 M). The relatively steep denisty profile required to model MM2 depends on a significant temperature decrease at its centre (see Sect. 6), which can be justified by the lack of star formation in the core. Increasing the minimum temperature decreases the required model mass, but for temperatures ~10 K, MM2 still requires a steeper model density profile than MM1.

A further example of a higher-mass core at an earlier evolutionary stage than a neighbouring lower-mass core, can be found in Stephens et al. (2015).

8.1. The future stellar population in SDC13

Based on our models, MM1 contains one star with a mass of ~7.36 M. We can use this to estimate the possible future stellar population which MM1 could produce, by assuming that the stars will have an IMF distribution of masses, and conservatively normalising the IMF to have one star in the mass range 6 MM ≤ 9 M. Adopting a power law form of the IMF, dNMαdM (where 2.3 ≤ α ≤ 2.35, Salpeter 1955; Kroupa 2002), we find that MM1 could potentially form ~2526 stars in the mass range 1 MM ≤ 120 M. Kroupa (2002) estimate that stars in this mass range make up ~6.1% of the total number of stars in the galactic field (for an IMF with α = 2.3). This suggests that MM1 has the potential to form at least ~409 stars in total (with masses 0.01 MM ≤ 120 M). Based on the average stellar mass () given by Kroupa (2002), this equates to a total stellar mass ~155 M. Using the total mass of MM1 (~360 M, derived from the best-fit RADMC-3D model; Sect. 7), this is equivalent to a star formation efficiency of ~43%, consistent with the relatively high star formation efficiencies seen in regions of dense core and cluster formation (e.g. Wilking & Lada 1983; Hillenbrand & Hartmann 1998; Koenig et al. 2008). MM2 might be expected to form a similar population of stars to MM1 in the future.

8.2. The evolution of MM1 and MM2

The absence of a star in MM2 suggests that it is in an earlier stage of evolution than MM1. However, based on the proximity of MM1 and MM2, one could assume that they were formed at the same time. This would suggest a shorter evolution timescale for MM1 than for MM2. Differences in star formation timescales could be the result of differences in the size and density of the region in a large scale dynamical collapse scenario (e.g. Peretto et al. 2014) or differences in turbulent and/or magnetic support between regions (Myers & Fuller 1993; McKee & Tan 2003). Peretto et al. (2014) find a slightly higher value for the velocity dispersion in the vicinity of MM2 (σtot ~ 0.84 km s-1) than in the vicinity of MM1 (σtot ~ 0.81 km s-1), however the difference is small and does not translate into enough additional support to prevent the collapse of MM2. The slight difference in the measured velocity dispersion of the turbulent motions around MM1 and MM2 is therefore not a plausible explanation for their different evolutionary stages.

In order to explain the observed distribution of protostellar and prestellar cores in the Perseus molecular cloud, Hatchell & Fuller (2008) suggest that higher-mass prestellar cores will have shorter evolutionary timescales. If this is the case, MM1 may have had a larger initial mass than MM2. Based on our best fit models however, MM2 is currently ~460 M more massive than MM1 (Sect. 7). It is possible that MM1 could have experienced mass loss after the formation of its central protostar via accretion on to the star, but, based on our best fit model, only a small proportion of the mass difference can be accounted for in the star at the centre of MM1 (~7 M).

A further possibility is that the observed mass difference between MM1 and MM2 is the result of outflows from the stars in MM1. The three-colour image of SDC13 shown in Fig. 1 indicates an extended green (8.0 μm) biconical structure in the vicinity of MM1, characteristic of emission from PAHs excited by UV photons. IR emission in the Spitzer IRAC bands, observed towards regions of star formation, has been shown to trace the illumination of outflow cavities by scattered light from central star forming objects (Qiu et al. 2006; Tobin et al. 2007; Velusamy et al. 2011). The presence of 8 μm biconical emission around MM1 provides evidence for the existence of a bipolar outflow cavity (and therefore a bipolar outflow) in this region. However, the removal of 460 M would imply an unusually high mass loss rate ~5 × 10-4 − 5 × 10-3M yr-1 over 105 − 106 yr, compared to a typical value of ~10-4M yr-1 for similarly sized clumps (e.g. Zhang et al. 2005; de Villiers et al. 2014).

An alternative explanation for the mass difference between MM1 and MM2 could be that MM2 may have continued to increase in mass once MM1 started star formation. However, using the mass inflow rate estimated by Peretto et al. (2014) (~2.5 × 10-5M yr-1) this would imply an age difference ~107 yr between MM1 and MM2. This is unrealistically long given the presence of embedded young stars in MM1 and the ~105 yr timescale for the formation of massive stars (e.g. McKee & Tan 2002).

The density profile in MM1 is consistent with that expected for a collapsing region whereas the profile of MM2 is steeper, or at least has a steeper outer envelope. This difference may reflect how a clump density profile evolves as the clump evolves towards forming stars. With its apparent coeval, similar mass, but star forming, neighbouring clump in a near identical environment, MM2 is an important laboratory for further study towards understanding how massive clumps evolve towards the formation of massive stars. The steep density profile of MM2 may also make it more likely to form a single massive central object. This is based on evidence presented by Girichidis et al. (2011), who find that massive stars predominantly form from cores with initial density profiles that are strongly centrally condensed (independent of support from radiation or magnetic fields). They conclude that the initial density profile of star forming cores is perhaps the most important factor in determining the fragmentation, evolution and final mass distribution of massive star forming regions.


1

The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Instituten of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica.

References

All Tables

Table 1

An overview of observing parameters and the largest spatial scales that the SMA is sensitive to in the configurations used.

Table 2

Physical properties and J2000 coordinates of MM1 and MM2, and sub-fragments A, B, C and D.

Table 3

Photometry results and J2000 coordinates for sources Her1 and Her2.

Table 4

Parameters of RADMC-3D models shown in Figs. 8 and 9 and their corresponding CASA simulated peak flux values.

All Figures

thumbnail Fig. 1

A three-colour Spitzer image of SDC13 (24 μm in red, 8 μm in green and 3.6 μm in blue) overlayed with IRAM MAMBO 1.2 mm continuum contours (white) at 5, 15, 25, 35 and 45 mJy (left). The highlighted region (right) shows the two largest fragments (i.e. cores) in SDC13, MM1 and MM2, whose positions are marked with white crosses. Four sub-fragments are seen in the SMA 1.3 mm continuum observations of the region. The positions of these sub-fragments, A, B, C and D, are marked with yellow circles.

Open with DEXTER
In the text
thumbnail Fig. 2

The primary beam corrected linear mosaic of the 1.3 mm continuum data obtained in two SMA pointings towards MM1 and MM2 (Jy/beam). The central portion, outlined in blue, indicates the region shown in Fig. 3. The synthesised beam is shown in the bottom left hand corner of the image. Contours are at − 3σ,3σ,5σ,7σ,9σ where σ ≈ 1 mJy/beam. Positive contours are indicated by the solid lines and negative contours by the dashed lines. Letters A, B, C and D label the 4 largest sub-fragments seen in the data, defined by the 3σ contours.

Open with DEXTER
In the text
thumbnail Fig. 3

SMA 1.3 mm continuum data overlayed on 1.2 mm MAMBO continuum data (mJy/beam). The contours indicate the SMA continuum at 3, 3, 5,7 and 9 mJy/beam. Positive contours are shown in pink and negative contours are shown dashed in green. Letters A, B, C, D indicate the 4 sub-fragments seen in the region which are identified as peaks in the 1.3 mm SMA continuum. The MAMBO beamsize is 10.7′′.

Open with DEXTER
In the text
thumbnail Fig. 4

SED fit to Herschel Hi-GAL data observed towards MM1 at 300, 250, 160, 70 μm, based on a distance of 3.6 kpc and using β = 1.5.

Open with DEXTER
In the text
thumbnail Fig. 5

SMA 1.3 mm continuum data overlayed on Herschel Hi-GAL images (Jy/beam) at 350 μm (top left), 250 μm (top right), 160 μm (bottom left) and 70 μm (bottom right). The contours indicate the SMA 1.3 mm continuum at 3, 5, 7 and 9 mJy/beam. Letters A and C indicate 2 of the 4 sub-fragments seen in the region which coincide with MM1. The two sources, Her1 and Her2, indicated in blue, were extracted from the 70 μm Herschel data using the Hyper algorithm (Traficante et al. 2015, see Table 3).

Open with DEXTER
In the text
thumbnail Fig. 6

Comparison of the observed MAMBO 1.2 mm flux distributions of MM1 (left) and MM2 (bottom) with RADMC-3D model fluxes, obtained using various model parameters. The solid black lines indicate the observed fluxes plotted against the equivalent radius (where Anσ is the area contined within contours , for n = 3, 4, 6, 7, 8, 9 and σ ~ 5 mJy/beam). The grey lines indicate perpendicular slices of observed flux against radius; these slices give an indication of the range of variation in the observed flux distributions in MM1 and MM2. The coloured lines show the model fluxes against radius for RADMC-3D models with various dust distributions, as indicated in the key. The model dust masses Mdust required to match each peak model flux with the peak observed flux are also given. Models are labelled A1a to A3a (for those fit to MM1 data) and D1a to D4a (for those fit to MM2 data). The dashed curves indicate the MAMBO beam (10.7′′).

Open with DEXTER
In the text
thumbnail Fig. 7

Density and temperture profiles of the best fit models (with Rout = RHi − GAL) to MAMBO observations of MM1 (left) and MM2 (right; models A2a and D3a respectively).

Open with DEXTER
In the text
thumbnail Fig. 8

CASA simulated SMA 1.3 mm continuum map of the SDC13 region (bottom) made using RADMC-3D models A2a (top) and D3a (centre; see Table 4). The map contour is at 1.5 mJy. The rms noise in the centre of the map ~0.09 mJy. Model fluxes for A2a (shown in red) and D3a (shown in green) are compared with plots of the observed MAMBO 1.2 mm flux against equivalent radius for MM1 and MM2 respectively (solid black lines). Model parameters are given in the key of each plot. The grey lines in the top two plots indicate perpendicular slices of observed flux against radius; these slices give an indication of the range of variation in the observed flux distributions in MM1 and MM2. The dashed curves indicate the MAMBO beam (FWHM = 10.7′′).

Open with DEXTER
In the text
thumbnail Fig. 9

CASA simulated SMA 1.3 mm continuum map of the SDC13 region (bottom) made using RADMC-3D models A2b (top) and D3b (centre; see Table 4). Map contours are at 2, 4, 6, 8 and 10 mJy. The rms noise in the centre of the map ~0.3 mJy. Model fluxes for A2b (shown in red) and D3b (shown in green) are compared with plots of the observed MAMBO 1.2 mm flux against equivalent radius for MM1 and MM2 respectively (solid black lines). Model parameters are given in the key of each plot, with further details given in the text. The grey lines in the top two plots indicate perpendicular slices of observed flux against radius; these slices give an indication of variation in the observed flux distributions in MM1 and MM2. The dashed curves indicate the MAMBO beam (FWHM = 10.7′′).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.