Dense molecular cocoons in the massive protocluster W3 IRS5: a test case for models of massive star formation

Dense molecular cocoons in the massive protocluster W3 IRS5: a test case for models of massive star formation

K.-S. Wang Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands; email:    T. L. Bourke Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA    M. R. Hogerheijde Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands; email:    F. F. S. van der Tak SRON Netherlands Institute for Space Research, Landleven 12, 9747 AD, Groningen, The Netherlands Kapteyn Astronomical Institute, University of Groningen, The Netherlands    A. O. Benz Institute of Astronomy, ETH Zürich, 8093 Zürich, Switzerland    S. T. Megeath Ritter Observatory, MS-113, University of Toledo, 2801 W. Bancroft St., Toledo, OH 43606, USA    T. L. Wilson Naval Research Laboratory, Code 7210, Washington, DC 20375, USA
Key Words.:
Stars: massive – Stars: formation – ISM: individual objects: W3 IRS5 – ISM: kinematics and dynamics
offprints: K.-S. Wang,

Context:Two competing models describe the formation of massive stars in objects like the Orion Trapezium. In the turbulent core accretion model, the resulting stellar masses are directly related to the mass distribution of the cloud condensations. In the competitive accretion model, the gravitational potential of the protocluster captures gas from the surrounding cloud for which the individual cluster members compete.

Aims:With high resolution submillimeter observations of the structure, kinematics, and chemistry of the proto-Trapezium cluster W3 IRS5, we aim to determine which mode of star formation dominates.

Methods:We present 354 GHz Submillimeter Array observations at resolutions of (1800–5400 AU) of W3 IRS5. The dust continuum traces the compact source structure and masses of the individual cores, while molecular lines of CS, SO, SO, HCN, HCS, HNCO, and CHOH (and isotopologues) reveal the gas kinematics, density, and temperature.

Results:The observations show five emission peaks (SMM1–5). SMM1 and SMM2 contain massive embedded stars (20 ); SMM3–5 are starless or contain low-mass stars (8 ). The inferred densities are high, cm, but the core masses are small, . The detected molecular emission reveals four different chemical zones. Abundant ( few to ) SO and SO are associated with SMM1 and SMM2, indicating active sulfur chemistry. A low abundance () of CHOH concentrated on SMM3/4 suggest the presence of a hot core that is only just turning on, possibly by external feedback from SMM1/2. The gas kinematics are complex with contributions from a near pole-on outflow traced by CS, SO, and HCN; rotation in SO, and a jet in vibrationally excited HCN.

Conclusions:The proto-Trapezium cluster W3 IRS5 is an ideal test case to discriminate between models of massive star formation. Either the massive stars accrete locally from their local cores; in this case the small core masses imply that W3 IRS5 is at the very end stages (1000 yr) of infall and accretion, or the stars are accreting from the global collapse of a massive, cluster forming core. We find that the observed masses, densities and line widths observed toward W3 IRS 5 and the surrounding cluster forming core are consistent with the competitive accretion of gas at rates of yr by the massive young forming stars. Future mapping of the gas kinematics from large to small scales will determine whether large-scale gas inflow occurs and how the cluster members compete to accrete this material.

1 Introduction

In contrast to low-mass star formation, there is no single agreed upon scenario for the formation of massive stars (Zinnecker & Yorke 2007). Current models either predict that a turbulent core collapses into a cluster of stars of different mass (e.g. McKee & Tan 2003), or that multiple low-mass stars compete for the same mass reservoir, with a few winning and growing to become massive stars while most others remain low-mass stars (e.g. Bonnell et al. 2001; Bonnell & Bate 2006). Both theories are supported by observations (e.g. Cesaroni et al. 1997; Pillai et al. 2011), suggesting that there may be two different modes of massive star formation (Krumholz & Bonnell 2009). However, it is challenging to observe massive star-forming regions since they are at large distances (a few kpc), requiring high-angular resolution to resolve the highly embedded complex structures in which gravitational fragmentation, powerful outflows and stellar winds, and ionizing radiation fields play influence the star-formation processes (Beuther et al. 2007).

Although observationally challenging, progress in testing the models of massive star formation has been made recently with at centimeter, and (sub)millimeter wavelengths. Based on observations of the protocluster NGC 2264 with the IRAM 30m telescope, Peretto et al. (2006) proposed a mixed type of star formation as proposed by Bonnell et al. (2001) and Bonnell & Bate (2006), and McKee & Tan (2003), where the turbulent, massive star-forming core is formed by the gravitational merger of lower mass cores at the center of a collapsing protocluster. Observations of G8.68–0.37 with the Submillimeter Array (SMA) and the Australia Telescope Compact Array (ATCA) imply that to form an O star, the star-forming core must continuously gain mass by accretion from a larger mass reservoir since protostellar heating from the low-mass stars is not sufficient to halt fragmentation (Longmore et al. 2011). Multi-wavelength observations of G29.96–0.02 and G35.20–1.74 suggest that the mass of massive star-forming cores is not limited to the natal cores formed by fragmentation and turbulence is not strong enough to prevent collapse, which favors the competitive accretion model (Pillai et al. 2011). On the other hand, Very Large Array (VLA) observations of IRAS 05345+3157 (Fontani et al. 2012) suggest that turbulence is an important factor in the initial fragmentation of the parental clump and is strong enough to provide support against further fragmentation down to thermal Jeans masses.

W3 IRS5 is a massive star-forming region located in the Perseus arm at a distance of kpc (Imai et al. 2000) with a total luminosity of (Campbell et al. 1995). A cluster of hypercompact ( AU) H II regions are found from high-resolution cm-wavelength observations (Claussen et al. 1994; Wilson et al. 2003; van der Tak et al. 2005). Based on observations at near-IR wavelengths, Megeath et al. (1996, 2005, 2008) found a high stellar surface density of pc and proposed that W3 IRS5 is a Trapezium cluster (Abt & Corbally 2000) in the making. A high protostellar number density exceeding protostars pc is also concluded by Rodón et al. (2008) from their 1.4 mm observations with the Plateau de Bure Interferometer (PdBI). The proximity and the dense protostellar environment make W3 IRS5 an excellent target to study massive star formation in a highly clustered mode.

(Sub)millimeter observations show that the molecular structure of W3 IRS5 is physically and chemically complex. Multiple bipolar outflows with various orientations are reported via different tracers (e.g. Mitchell et al. 1992; Ridge & Moore 2001; Gibb et al. 2007; Rodón et al. 2008; Wang et al. 2012). Although the outflow driving sources are not determined unambiguously, the infrared pair NIR1 and NIR2 (Megeath et al. 2005) are likely candidates. The detection of near-IR and X-ray emission toward W3 IRS5 (Megeath et al. 2005; Hofner et al. 2002) likely benefits from an outflow oriented near the line of sight. A velocity gradient in the NW–SE direction on scales of a few arcseconds is seen in SO, which may indicate rotation of the common envelope of NIR1 and NIR2 (Rodón et al. 2008; Wang et al. 2012). At resolution, many molecular species are detected toward W3 IRS5, especially sulfur-bearing molecules (Helmich et al. 1994; Helmich & van Dishoeck 1997), implying active hot-core chemistry (e.g. Charnley 1997) or shocks (e.g. Pineau des Forets et al. 1993). The sulfur-bearing species SO and SO peak near NIR1/NIR2, while complex organic molecules such as CHCN and CHOH peak at offset positions (Rodón et al. 2008; Wang et al. 2012).

This paper presents Submillimeter Array (SMA) observations of W3 IRS5 at 354 GHz with an angular resolution of , revealing the complex molecular environment of the proto-Trapezium cluster and tracing the star formation processes in the dense cluster-forming core. The observations and data calibration are summarized in Sect. 2. We present the observational results in Sect. 3, followed by data analysis in Sect. 4. We discuss our findings in Sect. 5, and conclude this paper in Sect. 6.

Figure 1: (a) Continuum emission of W3 IRS5 region at 353.6 GHz (black contours) overplotted with the NICMOS 2.22 emission (color scale). The synthesized beam size is , P.A. . The contour levels are -3, 3, 5, 7, … , 19, 23, 27, … with 1 of 10 mJy beam. The red crosses mark the positions of the five point sources derived from the visibility fit. The red open squares are the positions of near-infrared sources identified by Megeath et al. (2005) while the red filled circles represent the positions of cm-wave sources from van der Tak et al. (2005). (b) Vector-averaged 353.6 GHz continuum visibilities. The red filled circles with error bars are the observed data. The expected zero amplitude is shown as a dashed histogram. The model consisting of five point sources is plotted with open blue squares.
Source aafootnotemark: bbfootnotemark: ccfootnotemark: ddfootnotemark: n(H)eefootnotemark: Associated cm-wave,
(hh mm ss) (dd mm ss) (Jy beam) (Jy) ( cm) () ( cm) IR, or mm source
SMM1 02 25 40.779 +62 05 52.55 0.41 0.43 4.4 0.6 3.2 D2fffootnotemark: , Q5ggfootnotemark: , K7ggfootnotemark: , MIR1ggfootnotemark: , NIR1hhfootnotemark: ,
MM-1iifootnotemark: , SMS1-MM1jjfootnotemark:
SMM2 02 25 40.679 +62 05 51.89 0.30 0.39 3.2 0.5 2.9 Bfffootnotemark: , Q1-3ggfootnotemark: , K2-4ggfootnotemark: , MIR2ggfootnotemark: , NIR2hhfootnotemark: ,
MM-2/3iifootnotemark: , MM-5iifootnotemark:
SMM3 02 25 40.509 +62 05 50.41 0.16 0.24 1.7 0.3 1.8 SMS1-MM2jjfootnotemark:
SMM4 02 25 40.431 +62 05 51.36 0.14 0.22 1.5 0.3 1.6 SMS1-MM2jjfootnotemark:
SMM5 02 25 40.504 +62 05 52.99 0.09 0.13 1.0 0.2 1.0
111 $a$$a$footnotetext: Peak intensity measured from the image. $b$$b$footnotetext: Point source flux density measured from the visibility fit. $c$$c$footnotetext: H column density. We assume K for all submillimeter sources. $d$$d$footnotetext: H mass. We assume K for all submillimeter sources. $e$$e$footnotetext: H volume density. We assume K for all submillimeter sources. A spherical source with a diameter of or 1830 AU is assumed. $f$$f$footnotetext: Claussen et al. (1994) and Wilson et al. (2003). $g$$g$footnotetext: van der Tak et al. (2005). $h$$h$footnotetext: Megeath et al. (2005). $i$$i$footnotetext: Rodón et al. (2008). $j$$j$footnotetext: Wang et al. (2012).
Table 1: Characteristics of the continuum emission in W3 IRS5 at 353.6 GHz

2 Observations and Data reduction

W3 IRS5 was observed with the SMA222The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. (Ho et al. 2004) on 2008 January 7 in the compact configuration and on 2006 January 15 in the extended configuration. Both observations were conducted with 7 antennas and were centered on = 02254078, = 62055250. The weather conditions were good with . The reference was set to km s. The receiver was tuned to 353.741 GHz in the upper sideband, covering frequencies from 342.6 to 344.6 GHz and from 352.6 to 354.6 GHz. The correlator was configured to sample each spectral window by 256 channels, resulting a velocity resolution of km s per channel. For the compact configuration dataset, 3C454.3 was observed as bandpass calibrator. Gain calibration was performed by frequent observations of two quasars, 0136+478 ( away from the source) and 3C84 ( away from the source). Uranus was adopted as absolute flux calibrator. For the extended configuration dataset, 3C273, 3C84 and 3C273 were observed, respectively, as bandpass, gain and flux calibrators. The adopted total flux density of 3C273 is 9.6 Jy, which is the mean value reported for five observations333 during 2006 January 13–30. We estimate that the uncertainty in absolute flux density is about . The -range sampled by the combined dataset is k () to k ().

Data reduction was conducted by using the MIR package (Scoville et al. 1993) adapted for the SMA, while imaging and analysis were performed in MIRIAD (Sault et al. 1995). To avoid line contamination, line-free channels were selected to separate the continuum and line emission in the visibility domain. The continuum visibilities for each dataset were self-calibrated using the brightest clean components to enhance the signal-to-noise ratio. Self-calibrated datasets were combined to image the continuum. To optimize sensitivity and angular resolution, a robust weighting444 parameter of 0.25 was adopted for the continuum imaging, resulting in a rms noise of 10 mJy per beam. Line identification was conducted by using the compact configuration dataset only with the aid of spectroscopy databases of JPL555 (Pickett et al. 1998) and CDMS666 (Müller et al. 2005). For those lines that were also clearly detected in the extended-configuration dataset, we imaged the combined compact+extended dataset.

3 Observational results

3.1 Continuum emission

Figure 1 shows the continuum image and visibilities at 353.6 GHz. At a resolution of and P.A. , five major emission peaks, SMM1 to SMM5, were identified based on their relative intensities and a cutoff of 9 . We note that irregular, extended emission is also observed around the identified peaks. More emission peaks may be present in the diffuse surroundings, which requires better sensitivity, angular resolution, and image fidelity for confirmation. The total flux density recorded by the SMA is about 2.5 Jy which is about 6% of the flux contained in the SCUBA image at 850 m (Wang et al. 2012). To derive the positions and flux densities of the emission peaks, we fitted the data with five point sources in the visibility domain (Fig. 1 (b)). The properties of the five sources are summarized in Table 1. With this simplified model, about half of the observed flux originates from the unresolved sources, while the extended emission detected in baselines of 10–25 k contributes the other half. The fit deviates significantly from the data at long baselines ( k) and on short baselines ( k). The former corresponds to fine image detail that we could fit by adding more point sources, but we find that the S/N does not warrant additional “sources”. The latter indicates that on scales there is significant emission that our 5-point model does not represent properly. An extended, Gaussian envelope could match these observations. However, on even larger spatial scales our sampling misses even larger amounts of emission. Therefore, we do not aim to model this extended emission, but instead refer to the SCUBA 850 m emission to trace this component.

Table 1 lists the cm-wave, (sub)mm and IR sources associated with SMM1–5 (Claussen et al. 1994; Wilson et al. 2003; van der Tak et al. 2005; Megeath et al. 2005; Rodón et al. 2008; Wang et al. 2012). In Fig. 1 (a), we see that the emission decreases between SMM1/SMM2 and SMM3/SMM4/SMM5, dividing the cloud into two regions with the eastern part containing clusters of IR and cm-wave emission as well as X-ray emission (Hofner et al. 2002), while the western part shows no IR and cm-wave sources. This dichotomy implies that the eastern part of the source is more evolved and less embedded than the western part. We note that the projected positions of SMM1, SMM2 and SMM3 are aligned nearly in a straight line which could be the result of fragmentation at the scale of 1″–2″( 1800 AU–3600 AU), although this may be a projection effect.

Because SMM1 and SMM2 contain cm-wave sources, the observed submillimeter flux may contain non-thermal contributions. We estimate the free-free contribution from VLA observations at 5, 15 and 22.5 GHz (Tieftrunk et al. 1997). For SMM1 and SMM2, the extrapolated free-free contributions at 353.6 GHz are 8 mJy and 0.4 mJy, respectively, much smaller than the observed fluxes (SMM1: 0.43 Jy, SMM2: 0.39 Jy). We therefore conclude that the 353.6 GHz emission is dominated by thermal emission from dust and ignore any contribution of free-free emission toward SMM1 and SMM2.

Figure 2: Moment 0 maps of molecular emission toward W3 IRS5. The 353.6-GHz continuum emission is also presented for comparison. In panels (a) to (l), the compact dataset is imaged using natural weighting, resulting in a resolution of , (P.A. ) shown as the filled ellipse in the bottom-left corner of each panel. In panels (m) to (r), the combined dataset is imaged using uniform weighting, resulting in a resolution of (P.A. ). The five crosses represent the positions of SMM1 to SMM5 (c.f. Fig. 1). Open squares are the positions of IR sources taken from Megeath et al. (2005), while filled circles are the positions of cm-wave sources reported by van der Tak et al. (2005). The (0, 0) position corresponds to the reference phase center (Sect. 2). In all panels, contour levels start at 3 (solid line) and 3 (dashed line). The contour units are Jy beam km s for panel (a)–(k) and (m)–(r), and Jy beam for panel (l). The 1- noise levels from (a) to (r) are 1.0, 0.9, 0.35, 0.39, 0.72, 0.46, 1.0, 0.45, 0.22, 0.45, 0.36, 0.04, 0.47, 0.94, 0.66, 0.68, 0.84 and 0.68 respectively. The contour steps from (a) to (r) are 10, 10, 2, 2, 10, 2, 4, 2, 1, 1, 2, 4, 4, 5, 20, 5, 6 and 5, respectively.
Moleculeaafootnotemark: Frequency Transition bbfootnotemark: ccfootnotemark:
(MHz) (K) (cm)
CS 342883.0 =7–6 66
SO 344310.6 = 87
SO 343086.1 =, =15/2–13/2 78
343087.3 =, =17/2–15/2 78
343088.1 =, =19/2–17/2 78
343088.3 =, =21/2–19/2 78
SO 342761.6 = 582 ddfootnotemark:
343923.8 = =1 1038 ddfootnotemark:
SO 353741.0 =, =39/2–39/2 213
353741.1 =, =37/2–37/2 213
353741.6 =, =41/2–41/2 213
353741.6 =, =35/2–35/2 213
SO 344245.3 = 88
344581.0 = 167
353002.4 = 212
354277.6 = 580 ddfootnotemark:
354397.8 = 326
HCN 354505.5 =4–3 43
354460.5 =4–3 =1 1067 ddfootnotemark:
HCN 344200.3 =4–3 43
HCS 342944.4 = 91 ddfootnotemark:
HNCO 352897.9 =, =17–16 187
352897.9 =, =16–15 187
352897.9 =, =15–14 187
CHOH 342729.8 = 227
777$a$$a$footnotetext: Spectroscopy data taken from JPL molecular spectroscopy and CDMS. $b$$b$footnotetext: Upper level energy. $c$$c$footnotetext: Critical density at 100 K derived from LAMDA (Schöier et al. 2005). The critical densities for SO, SO and SO are quoted from their main isotopologues. $d$$d$footnotetext: Collisional rate coefficient is not available for this transition.
Table 2: Detected molecular transitions in W3 IRS5

3.2 Line emission

3.2.1 Molecular distribution

Within the passbands (342.6–344.6 GHz and 352.6–354.6 GHz), we identify 17 spectral features from 11 molecules (Table 2). Most emission features come from sulfur-bearing molecules such as CS, SO, SO, SO, SO, SO and HCS. Others come from HCN, HCN, HNCO and CHOH. This set of lines covers upper level energies from 43 K to 1067 K and high critical densities ( cm). Figure 2 shows a sample of molecular emission images toward W3 IRS5 using the compact-configuration data and natural weighting (panels (a) to (l), with a resolution of , P.A. ) and with uniform weighting (panels (m) to (r), with a resolution of , P.A. ).

Most of the emission from SO and SO and their isotopologues peaks toward the continuum peaks SMM1 and SMM2, while weaker emission is seen toward SMM3 and SMM4. However, other sulfur-bearing molecules such as CS and HCS show a very different spatial distribution, and mainly peak to the north-west of SMM1–5 (Fig. 2 (i) and (m)). Most of the emission from HCN and its isotopologues comes from SMM1 and SMM2 (Fig. 2 (g), (h) and (n)). Additional emission can be seen toward SMM3 and SMM4, as well as east and north-west of the SMM sources. The typical hot-core molecule CHOH exclusively peaks toward SMM3 and SMM4 (Fig. 2 (k)). HNCO is detected toward SMM1 and SMM2 only (Fig. 2 (j)). Combining all these observational facts, the overall molecular emission toward W3 IRS5 can be classified into four zones (A: SMM1/SMM2; B: SMM3/SMM4; C :north-west region; may be associated with SMM5; and D: east of SMM1–5), and summarized in Figure 3.

Figure 3: The main molecular zones identified toward W3 IRS5. The representative molecular species including their isotopologues are listed. The crosses, filled circles and open squares represent the same sources as in Fig. 2.

3.2.2 Velocity channel maps

Among the detected spectral features, most of the emission is compact (in zones A and B). Only CS =7–6, SO = and HCN =4–3 are extended and trace the large scale gas kinematics. Figure 4 shows the channel maps of these molecules from km s to km s and from km s to km s, using robust weighting (robust = 0). These three molecules show emission over a broad velocity range. The emission at extreme velocities is mainly concentrated toward zone A. Weaker emission is seen toward zone B. Zone C is best traced by CS with some emission from SO, HCN, and HCS as well. At velocities ranging from to km s, the emission from all three molecules can be seen toward zone D.

Of these three molecules, the image quality near the systemic velocity ( km s) is poor due to the missing flux of extended emission. However, some additional emission features can be seen near the systemic velocity. Figure 4 shows the channel maps of CS, SO and HCN at velocities ranging from km s to km s. The emission of HCN is much weaker than the other species near km s suggesting that HCN is much more extended, and therefore more heavily filtered out by the interferometer, or that the emission is self-absorbed. Additional extended emission in the N–S and NW–SE directions can be seen in the CS channel maps from km s to km s, which might delineate two collimated outflows or a single wide-angle outflow in the NW–SE direction. The small velocity range suggests that the outflow may be in the plane of sky. In the SO channel maps, a bar-like structure can be seen from km s to km s. Our data suggest that CS, SO and HCN trace different parts of the molecular condensation in W3 IRS5 and show a complex morphology in velocity space.

Figure 4: Velocity channel maps of CS =7–6, SO = and HCN =4–3 toward W3 IRS5. Robust weighting (robust = 0) is used for the imaging, resulting resolutions of , P.A. , , P.A. , and , P.A. , for CS, HCN and SO, respectively. The and -axes are R.A. offset and DEC offset relative to the phase center = 02254078, = 62055250 in arcseconds, respectively. The upper three rows shows the velocity range between km s to km s. The lower three rows show a narrower velocity range from km s to km s. The markers are identical to Fig. 2. Solid and dashed contours represent positive and negative intensities, respectively. The absolute contour levels for CS in the upper panel are 3, 8, 13,… (1 = 0.11 Jy beam). The absolute contour levels for HCN in the upper panel are 3, 8, 13,… (1 = 0.12 Jy beam). The absolute contour levels for SO in the upper panel are 3, 13, 23, 43, 63… (1 = 0.11 Jy beam). The absolute contour levels for CS in the lower panel are 3, 8, 13,… (1 = 0.18 Jy beam). The absolute contour levels for HCN in the lower panel are 3, 7, 11,… (1 = 0.16 Jy beam). The absolute contour levels for SO in the lower panel are 3, 13, 23, 43, 63,… (1 = 0.14 Jy beam).

3.2.3 Line profiles

Figure 5 presents the hanning-smoothed spectra toward zone A (thick line; offset position ,) and B (thin line; offset position ,) where most of the detected molecules peak. The compact-configuration dataset and natural weighting is used for all molecules except CS, SO and HCN, which are imaged with the combined dataset in natural weighting and convolved to the same resolution as the other images.

SO shows a flat-top line profile at both positions, implying the line may be optically thick or is a blend of multiple velocity components. If SO is fully thermalized, optically thick, and fills the beam, the peak brightness temperatures indicate kinetic temperatures of 70 K and 45 K for zone A and B, respectively. These values are lower limits to the kinetic temperature if the beam filling factor is not unity. In addition, a blue-shifted spectral feature is seen at zone B.

CS and HCN show broad line profiles with a dip near the systemic velocity due to missing flux of extended emission. Single-dish JCMT observations of the same transitions (Helmich & van Dishoeck 1997) show single-peaked profiles and rule out self absorption. The line profiles at zone A are comparable with more emission in the blue-shifted line wing, suggesting both species trace similar cloud components. At zone B, both lines also show broad line wings. Near systemic velocity, HCN shows a wider dip than CS, suggesting that the spatial distributions of both lines are different. Vibrationally excited HCN peaks mainly toward zone A and shows an asymmetric line profile. HCN is present toward both zones; its smaller line width toward zone B implying that this zone is more quiescent than zone A.

The line profiles of all the other molecules presented in Fig. 5 are Gaussian toward both zones. There is no significant velocity shift between the two zones. We note that the apparent velocity shift seen in SO and SO is due to the blending of multiple hyperfine transitions; we have defined the velocity axis with respect to the hyperfine component with the lowest frequency (Table 2).

We derive the line center LSR velocity (), FWHM line width () and integrated line intensity () of the line profiles in zones A and B using Gaussian decomposition (Table 3). Since the CS and HCN line profiles are highly non-Gaussian, we do not attempt to derive the line parameters. For molecules with hyperfine transitions such as SO, SO and HNCO, we assume that each component has the same FWHM line width and LSR velocity. The relative intensity of hyperfine components are calculated assuming LTE and we fit the combined profile to the data. The results are summarized in Table 3. From the Gaussian fit, we do not find significant differences in LSR velocity between zones A and B. Different line widths are derived in both zones depending on the molecular transitions but the general trend is that the lines are broader in zone A, suggesting a more turbulent environment.

We estimate the amount of missing flux of CS 7–6, SO , HCN 4–3, and SO by comparing our SMA observations with the JCMT observations (Helmich & van Dishoeck 1997). Our SMA observations miss , , and flux for CS, SO, and HCN, respectively. Given the excitation of these lines, the derived missing flux is consistent with the large amount missing flux at 850 m (; Sect. 3.1). The SMA observation of SO likely recovers all emission, because this transition has a high upper level energy (corresponding to 582 K), and the emission therefore preferentially traces the dense, warm, and compact part of W3 IRS5.

Figure 5: Hanning-smoothed molecular spectra toward molecular zone A (thick line; offset position (,)) and B (thin line; offset position (,)). The compact-configuration dataset and natural weighting is used for all molecules except CS 7–6, SO = and HCN 4–3 which are taken from the combined dataset using natural weighting but convolved to the beam size of the compact dataset (, P.A. ). The apparent offset in velocity seen in SO and SO is due to the blending of hyperfine transitions for which the reference velocity is set to the transition with lowest frequency (see Table 2).
Zone A:SMM1/SMM2 Zone B:SMM3/SMM4
Molecule Transition aafootnotemark: bbfootnotemark: ccfootnotemark: aafootnotemark: bbfootnotemark: ccfootnotemark:
(Jy beam km s) (km s) (km s) (Jy beam km s) (km s) (km s)
CSddfootnotemark: 7–6
HCNddfootnotemark: 4–3
4–3 =1
HCN 4–3
888The beam size for all transitions listed in this table is about , P.A. . $a$$a$footnotetext: Integrated intensity. $b$$b$footnotetext: FWHM line width. $c$$c$footnotetext: Line center LSR velocity. $d$$d$footnotetext: Not attempted to perform Gaussian decomposition due to missing flux and complex line profile. $e$$e$footnotetext: With hyperfine transitions. Total is reported. $f$$f$footnotetext: Peaked toward zone C.
Table 3: Gaussian decomposition of the line profiles in W3 IRS5

4 Analysis

4.1 Mass and density of the submillimeter sources

We estimate the H column density and mass of each point source based on the continuum at 353.6 GHz. The H column density at the emission peak can be estimated as


where is the peak flux density, is the gas-to-dust ratio (100), is the beam solid angle, is the mass of atomic hydrogen, is the dust opacity per unit mass and is the Planck function at dust temperature . We apply the interpolated value of cm g as suggested by Ossenkopf & Henning (1994) for gas densities of cm and coagulated dust particles with thin ice mantles. The dust temperature of all five submillimeter sources is assumed to be 150 K (see excitation analysis, Sect. 4.2). The gas mass for each continuum peak is estimated from the total flux derived from the visibility fit (Table 1) via


where is the total flux density of dust emission and (1.83 kpc) is the distance to the source. Assuming a spherical source with in size, the H number density is also derived. We summarized the results in Table 1.

From the density and mass estimates, we find beam-averaged H column densities, masses and volume densities toward SMM1 and SMM2 of about cm, 0.5 and cm, respectively. Toward SMM3, SMM4 and SMM5 values smaller by factors 2–3 are found. The derived H column densities and masses are sensitive to the adopted dust opacity and temperature. If a bare grain model is adopted (Ossenkopf & Henning 1994), the derived values are reduced by a factor of 3. For dust temperatures between 100 and 200 K, the derived values change by a few ten% only. If the dust temperature is 30 K, the reported H column densities and masses increase by factors of 5. The excitation analysis of Sects. 4.2.1 and 4.2.2 suggest such low temperatures are unlikely. The estimated H volume density is sensitive to the assumed source size as size. For example, the densities toward SMM1 and SMM2 increase to cm if the source size decreases to . Since the submillimeter continuum peaks are unresolved, we treat the H volume densities of Table 1 as lower limit. We note that the core masses are all small (1 ).

4.2 Excitation analysis: temperatures

The detected molecular transitions, covering from 43 K to 1067 K and critical densities of cm (Table 2), form an useful dataset for the diagnostics of the physical conditions toward W3 IRS5. Among these molecular lines, multiple detections of the transitions from SO and SO toward zones A and B allow us to perform excitation analysis. Multiple detections of CHCN =12–11 transitions reported by Wang et al. (2012) form another useful dataset for additional constraints of the excitation conditions toward zone B. We describe the details in the following two sections.

Figure 6: (a) Rotation diagram of SO (filled squares) and SO (filled triangles) toward zone A. The straight lines are the model fit. The derived rotational temperatures are indicated. (b) Rotation diagram of CHCN =12–11 toward zone B (open circles) taken from Wang et al. (2012). The best-fit model from the RADEX LVG analysis is overplotted for comparison (filled symbols).

4.2.1 Rotational diagram

The molecular excitation conditions can be estimated via rotation diagram analysis assuming all the transitions are optically thin and the emission fills the beam. In a rotation diagram, the column density of the upper state is given by


where is the total degeneracy of the upper state, is the total molecular column density, is the partition function, is the upper level energy, is the Boltzmann constant and is the rotational temperature. The left-hand side of equation 3 can be derived from the observations via


where and are the major and minor axes of the clean beam in arcsec, respectively, is the integrated intensity in Jy beam km s, and are the spin and projected rotational degeneracies, respectively, is the rest frequency in GHz, is the line strength and is the dipole moment of the transition in Debye. Figure 6 plots the logarithm of the column densities from our data (Eq. 4) versus (Eq. 3) and a fitted straight line with and as free parameters. We adopt a uncertainty of the integrated intensity in the analysis.

Figure 7: RADEX excitation analysis toward zone A traced by SO (left four panels) and zone B traced by CHCN (right four panels). The surfaces near the best-fit solution on each parameter axis are displayed.

Figure 6 (a) shows the rotation diagrams of SO and SO toward zone A, which sample level energies from 88 K to 1038 K (Table 2). SO is detected in five transitions with ranging from 88 to 580 K, and indicates a of K. However, the curvature in the distribution of the points implies that there might be a cooler and a warmer component along the line of sight. A two-component fit to the data gives K and K. The gas component traced by SO ( = 582 and 1038 K) implies a rotational temperature of 233 K, compatible with the temperature of the warm component derived from SO but obviously less well constrained because only two data points are available.

The curvature in the SO rotation diagram (Fig. 6 (a)) may also be due to line opacity (unlikely for the SO isotopologue) or subthermal excitation (see next section). The relatively abundant main-isotope SO may have optically thick lines, leading to an overestimate for . However, we detected the transition of both SO and SO. Equation 4 yields a (SO)/(SO) ratio of at K, consistent with the S/S ratio of 22 in the interstellar medium (Wilson & Rood 1994), suggesting that the transition of SO (and its isotopologues) is optically thin. The transition of SO is also likely optically thin since its Einstein-A coefficient is only about a factor of 2 greater than that of the transition. These considerations suggest the a temperature gradient in zone A explains the rotation diagrams, with a cooler region at and a warmer region of (Table 4).

We do not have enough data to obtain a reliable temperature estimate toward zone B (Table 3). Only the lowest two transitions of SO in level energy are detected here, and imply a rotational temperature of about 133 K. To have a better constraint of the excitation conditions toward zone B, we took the CHCN lines obtained by Wang et al. (2012) with the SMA at a resolution of . This molecule peaks toward zone B exclusively and shows a much narrower FWHM line width ( km s) compared to SO ( km s), suggesting a different molecular condensation in zone B. Figure 6 (b) shows the CHCN =0–4 rotation diagram. The scatter in the data indicates that the lines are not optically thin. Therefore, we can not carry out rotation diagram analysis on the CHCN lines, and instead perform an LVG analysis in Sect. 4.2.2.

4.2.2 RADEX LVG model

To further investigate the physical conditions such as kinetic temperature (), column density () and H volume density () toward zone A and B, we performed statistical equilibrium calculations to solve the level populations inferred from the observations. The purpose here is to have a crude estimate. Detailed analysis of temperature and density profiles in W3 IRS5 is beyond the scope of this work. We used the code RADEX (van der Tak et al. 2007) which solves the level populations with an escape probability method in the large-velocity-gradient (LVG) regime. Limited by the available observed data and collisional rate coefficients (Schöier et al. 2005), we perform calculations only for the SO lines and CHCN lines (Wang et al. 2012). For SO, we adopted the collisional rate coefficients from SO (Green 1995)999Recent calculations by Spielfiedel et al. (2009) and Cernicharo et al. (2011) show that the collisional rate coefficients calculated by Green (1995) are lower by a factor 3–5. The listed critical densities in Table 2 may be 3–5 times smaller. as an approximation since both species have similar molecular properties. The collisional rate coefficients of CHCN are taken from Green (1986). To find the best-fit solution, we applied minimization to the 4-dimensional parameter space (, , and , the beam filling factor).

Figure 7 shows the surfaces near the best-fit solution on each parameter axis, and Table 4 summarizes the results. Toward zone A, we use four SO lines (excluding the transition with K) in the calculations due to the limited table of collisional rate coefficients. The RADEX calculations suggest that the kinetic temperature traced by the SO lines is about 90 K, consistent with the temperature derived from rotation diagram analysis and suggesting the excitation is thermalized. This requires a high H volume density ( cm) toward zone A. Apparently, the cool component toward zone A is very dense. Toward zone B, our RADEX calculations indicate a kinetic temperature of 160 K. Interestingly, the H volume density traced by CHCN is only about 10 cm, lower than the critical density of a few cm, implying subthermal excitation. The beam filling factor is about 0.8. The best-fit model for CHCN is plotted as filled squares in Fig. 6 (b).

Alternatively, the excitation conditions toward zone A and B can be estimated via the line ratio of the over the transition of SO of 1.8. We assume that both transitions have the same filling factor. Using a FWHM line width of 5.0 km s (the average FWHM line widths in zones A and B are 6.0 km s and 4.0 km s, respectively) in the RADEX calculations, we calculate the intensity ratio as function of and for four different values of (Fig. 8). We again assume that the two transitions in the RADEX calculations have the same filling factor. These calculations suggest a high density of cm traced by SO toward zone A and B. The kinetic temperatures are K, depending on the H volume density.

Combining the results from the rotation diagram analysis and the RADEX LVG calculations, we conclude that the gas traced by SO is dense ( cm) in zones A and B. There is a temperature gradient along the line of sight in zone A, characterized by a cool region of K and a warm region of K. Toward zone B, the temperature gradient is less prominent. However, a quiescent (FWHM line width 1.6 km s) and less dense region ( cm) traced by CHCN is also present in zone B.

Zone A:SMM1/SMM2 Zone B:SMM3/SMM4
Molecule (K) (cm) aafootnotemark: (K) (cm) bbfootnotemark: Note
Rotation diagram: =
SO 233 zone A: 2 trans.
SO 133 zone A: lowest 4 trans, zone B: 2 trans.
zone A: highest 3 trans.
RADEX LVG model: =
SO 85 zone A: , cm
CHCNccfootnotemark: 160 zone B: , cm
Column densities at K
SO 150 150 opacities () in both zones
SO 150 150
SO 150 150
SO 150 150
SO 150 150
HCN 150 150
HNCO 150 150
CHOH 150 150
101010$a$$a$footnotetext: Fractional abundance. (H)= cm is derived from our data. $b$$b$footnotetext: Fractional abundance. (H)= cm is derived from our data. $c$$c$footnotetext: Data taken from Wang et al. (2012).
Table 4: Molecular excitations in W3 IRS5

4.2.3 Molecular column density

We estimate the beam averaged () molecular column density of each detected molecule toward zone A and B via Eq. 5 by assuming optically thin emission and a single rotational temperature of 150 K which is a representative value in W3 IRS5 (see previous two sections).


We only use the ground vibrational state for the estimates. If multiple transitions of a given molecule are detected, we adopt the averaged value. The uncertainty of the exact rotational temperature results in an error in the column density of a few 10% (for from 100 K to 200 K). As a consistency check, we estimate the line center opacity of each transition via


where is the speed of light, the Einstein-A coefficient, and the FWHM line width. All lines are consistent with the optically thin assumption (); only SO has line center opacities as large as –0.5). Therefore, a lower limit of SO column density is derived. A more representative value for SO can be derived from SO by assuming the isotopic ratio S/S of (Wilson & Rood 1994; Chin et al. 1996). To convert the molecular column densities into fractional abundances with respect to H, we adopt H column densities of cm and cm toward zone A and B, respectively (derived from the continuum image in Fig. 2 (l) with the same assumptions described in Sect. 4.1). We summarize the results in Table 4. Toward zone A and zone B, SO and SO are very abundant with fractional abundances up to few , several orders of magnitude higher than found in dark clouds (few , Ohishi et al. 1992). This is indicative of active sulfur chemistry. Among the detected molecules, we do not see significant differences in fractional abundance between zones A and B, except for CHOH which is a factor of 10 more abundant in zone B.

Figure 8: RADEX analysis of the SO line ratios toward zone A and B assuming a FWHM line width 5.0 km s. The intensity ratios / are plotted in color scale. The observed integrated intensities of the transition are plotted in curves (zone A: Jy beam km s, zone B: Jy beam km s). Toward both zones, the line ratio is about 1.8 (Table 3). Our results suggest that a high H volume density ( cm) is needed to reproduce the observed line ratios.
Figure 9: First moment maps of different molecules observed toward W3 IRS5. The upper four maps are derived from the compact-configuration dataset with natural weighting, while the bottom four maps are made from the combined dataset using uniform weighting. The markers are identical to the ones in Fig. 1. The zero position is the phase center. We note that the velocity ranges are the same for all maps except HCN 4–3.

4.3 Kinematics

Complex velocity fields in W3 IRS5 are observed in various molecules. Figure 9 shows several first moment maps derived from different molecules. The upper four maps are made from the compact-configuration dataset using natural weighting, while the bottom ones are derived from the combined dataset with uniform weighting. Different velocity gradients are observed in SO (NW–SE), HCN (E–W) and HCN (NE–SW). CHOH shows an interesting velocity distribution with a slightly blue-shifted emission peak (SMM3 and SMM4). Comparing the velocities near SMM3 and SMM4 derived from SO and CHOH, we suggest that there are two distinct regions along the line of sight. Indeed, the line widths of these two molecules are very different (Table 3). At one-arcsecond resolution (combined dataset with uniform weighting), SO peaks toward SMM1 and SMM2, and shows a velocity gradient consistent with the NW–SE gradient found from lower angular resolution observations. For the three strongest lines, CS, SO and HCN, we adopted uniform weighting to form the images in order to minimize the strong sidelobes which otherwise corrupt the images. CS and HCN show similar velocity patterns with some patchy velocity shifts. The velocity field traced by SO is very different from the ones seen from other molecules in Fig. 9, which is likely due to opacity effects. Near SMM1 and SMM2, HCN and its vibrationally excited transition show similar velocity patterns. The velocity components traced by CS and HCN toward zone C are red-shifted with respect to the systemic velocity ( km s). However, a clear velocity jump is seen near SMM5 if we compare the maps of SO to CS and HCN, implying that the gas component in zone C is not closely related to the ones in zones A and B. Another velocity jump is seen in the HCN map if we compare the velocities in zones A and D, suggesting that the gas component in zone D is also not closely related to zones A and B.

In Fig. 9, we see that the velocity gradients seen in SO and HCN are roughly perpendicular to each other. To confirm this, we fit the emission of HCN 4–3 ( K) and SO ( K) (compact configuration with natural weighting) channel-by-channel with a Gaussians to determine the movement of the peak positions. In Fig. 10, the velocity gradient seen in SO is roughly perpendicular to the line joining SMM1 and SMM2, while the velocity gradient observed in the vibrationally excited HCN is parallel to this line. The vibrationally excited HCN may trace outflow or jet emission since the observed velocity gradient is parallel to the free-free emission knots observed by Wilson et al. (2003). We suggest that SO traces a rotating structure orthogonal to this outflow axis.

In our data, only the CS , HCN and SO = emission show a wide velocity range (about 30–40 km s; Fig. 5), presumably due to outflow activity in W3 IRS5. To study the distribution of the high velocity components (Fig. 11), we integrate the emission over three different velocity ranges (blue: 50 to 43 km s, green: 43 to 35 km s and red: 35 to 28 km s). We took the combined-configuration dataset and used natural weighting for the analysis. The images are plotted with contours starting well above to minimize the impact of sidelobes, especially for the green component. As seen in Fig. 11, the high-velocity components of all three molecules peak near zone A and close to SMM2, suggesting a bipolar outflow with an inclination close to the line of sight. The existence of a line-of-sight outflow also implies a rotating structure in the plane of the sky, in which SMM1 and SMM2 form a binary system. Therefore, the small scale velocity gradients (Fig. 10) seen in HCN 4–3 and SO may highlight the kinematics of the envelope. Additional blue-shifted components are seen toward zone D. The green component of CS highlights the gas in zone C. In the HCN map, the emission in zone D together with the emission peaking near the offset position (, ) may form another bipolar outflow in the E–W direction since the line connecting these peaks passes through zone A which contains star-forming cores. In this case, the one-arcsecond scale velocity gradients seen in Fig. 10 imply a complex motion (rotation plus expansion or contraction) in the common envelope of SMM1 and SMM2. To summarize, the kinematics in W3 IRS5 is very complex and requires additional observations at high resolution with good image fidelity for further analysis.

Figure 10: Channel-by-channel locations of the emission peaks of the SO and HCN 4–3 . The emission peaks are color-coded with respect to their in km s shown in the color wedge. Two distinct, orthogonal velocity gradients are seen.
Figure 11: Integrated intensity maps at different velocity ranges (blue: 50 to 43 km s; green: 43 to 35 km s; red: 35 to 28 km s). We adopted the combined-configuration dataset in natural weighting in order to show the additional weak emission features. The contour levels of each Doppler-shifted component for CS, HCN and SO are 15, 20, 30, 40,…%, 10, 20, 30,…% and 5, 10, 20, 30,…% of the emission peak, respectively. We note that the first contours are all noise level. The markers are identical to the ones in Fig. 1.

5 Discussion

5.1 Star formation in W3 IRS5

5.1.1 Different evolutionary stages within W3 IRS5

Based on the continuum image at 353.6 GHz, we identified five compact sources (SMM1 to SMM5) toward W3 IRS5 (Fig. 1 and Table 1). SMM1 and SMM2 also show compact cm-wave emission (Wilson et al. 2003; van der Tak et al. 2005) indicating that stars sufficiently massive to ionize their surroundings have formed (8 ). The non-detection of compact cm-wave emission toward SMM3, SMM4 and SMM5 implies that currently those cores have not (yet) formed massive stars. By comparing the 353.6-GHz continuum with the NICMOS 2.22 m emission and cm-wave emission (Fig. 1(a), Megeath et al. 2005; Wilson et al. 2003; van der Tak et al. 2005), we suggest that the eastern part (SMM1 and SMM2) of the source is more evolved than the western part (SMM3, SMM4 and SMM5). The small core mass ( ) suggests that SMM3–5 may be forming low-mass stars or are starless cores.

For a virialized star-forming core (), the 1D virial velocity dispersion (; c.f. McKee & Ostriker (2007)) is


where is the radius of the source, is the gravitational constant, and is the virial mass. Assuming a core radius of , core gas masses of 0.3 , and stellar masses of SMM3–5 less than 8 , the 1D virial velocity dispersions should be no more than 1.3 km s. From our observations of CHOH line toward SMM3/4 and CS line toward SMM5, the inferred 1D velocity dispersions ( = /) are 1.0 km s and 1.4 km s. The large observed line widths suggest that either these cores contain stars just under 8 solar masses, or the clumps may have smaller mass stars, and the large velocities may be the influence of outflows from SMM1/2. In the latter, case, it is not clear if the gas in the cores is gravitationally bound. Subarcsecond observations of dust continuum and outflow tracers should tell if SMM3–5 are low-mass starless leftover cores after the formation of SMM1–2 or harbor low-mass protostars under the influence of outflows from SMM1–2.

5.1.2 Fragmentation and Jeans analysis

From the continuum analysis, we found the core masses of the identified SMM sources are small (0.2–0.6 ) and the projected distance between the SMM sources are about (1800 AU – 3600 AU). This observed structure, in the meantime, is still embedded in a large massive core with size about and mass of several hundred (W3-SMS1 in the SCUBA 850 m image, Di Francesco et al. 2008; Wang et al. 2012). At scale, Megeath et al. (1996) found that there are low-mass stars in the core surrounding W3 IRS5 with a stellar density of few 1000 pc. Therefore, it is interesting to investigate the cloud structure from scale to to see what determines the overall core/stellar mass distribution and their spatial distribution. To study this question, we perform Jeans analysis in W3 IRS5 based on our SMA continuum image and the JCMT SCUBA 850m image (Di Francesco et al. 2008; Wang et al. 2012). The Jeans length and Jeans mass are expressed as (c.f. Stahler & Palla 2005, Eq. (9.23) and (9.24)):




where is the kinetic temperature and is the H volume density. These two quantities are more sensitive to the adopted H volume density. To estimate the scale H density, we measured the total flux density in a by box centered toward W3 IRS5 using the SCUBA 850 m image. Using the same prescription outlined in Sect. 4.1, a total flux density of 160 Jy corresponds to a mass of 950 assuming a dust temperature of 40 K. The mean H volume density is cm. Therefore, at we derived the Jeans length and Jeans mass to be 16000 AU (9″) and 2 , respectively. Interestingly, these two characteristic quantities seem to match the observed spatial separations of low-mass stars (18000–12000 AU, Megeath et al. 1996), implying that the spacing of the young low mass stars are consistent with gravitational fragmentation. The Jeans masses are greater than the mean stellar mass for a standard initial mass function (0.5 ), which is expected given a typical efficiency of low-mass star formation (Alves et al. 2007). Toward the central regions of W3 IRS5, the temperature (140 K) and density (few cm) are higher than the values on the scale of the entire core. The corresponding Jeans length and Jeans mass are 2700–3800 AU (–2.1) and 1 , respectively. The Jeans length is consistent with the observed spacing of the SMM cores, implying that gravitational fragmentation may be occurring in the dense, warm center of the cluster. On the other hands, the implied Jeans masses are much smaller. This immediately leads to the question how the massive stars in SMM1/2 accreted their current masses.

5.1.3 Accretion rates

The largest accretion rate that an object can sustain is essentially given by its mass divided by its free-fall time,


where and is the mean density. The masses of the SMM1–5 cores come from Table 1. We use the observed size for the SMM1–5 cores, which are set to in diameter based on Fig. 1 (a). We adopt stellar masses of 20 for SMM1 and SMM2 based on the total luminosity of . We use the combined stellar and core masses and the core radii to determine the average baryonic density; the corresponding free-fall times are 1000 yr for SMM1/2, and yr for SMM3–5. These imply mass accretion rates of  yr for SMM1/2, and  yr for SMM3–5. The upper limits for the latter follow from the upper limit on the embedded stellar masses of 8 . The accretion rates are consistent with those expected for the formation of massive stars (Keto 2003). The combination of small core masses and high accretion rates imply very short time scales of a few thousand years. If the massive stars in SMM1/2 form from accretion of the local cores ( scales), the low core masses imply that they are now in a stage that most of core mass has been accrete onto the stellar components. On the other hand, the SMM cores may be accreting material that have been resolved out by the interferometers on scales . In this case, the reservoir of material from which the cores are accreting is larger than the observed projected spacing of the protostars, implying that the accretion is coming from a global collapse of the molecular core in which the proto-Trapezium is embedded.

To determine whether it is possible that the massive stars in W3 IRS5 are being fed from the collapse of the surrounding core, we estimate the potential accretion rates from this accretion. The free-fall time for the central W3–SMS1 condensation, with a total mass of 950 and an embedded stellar mass of , is 72000 yr and the free-fall mass accretion rate is  yr. We note that the derived accretion rate is the maximum rate of infall assuming free-fall collapse. In the central region, the amount of gas accreted by each star is given by


(Bonnell et al. 2001), where is the infall gas velocity, is the density of the gas surrounding the proto-Trapezium, and is the radius within which the infalling gas is accreted onto a specific star. We assume  cm to compute , which is intermediate between the density of the surrounding diameter core ( cm