© ESO, 2013

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 30 m 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 1.83 ± 0.14 kpc (Imai et al. 2000) with a total luminosity of 2 × 10⁵ L_⊙ (Campbell et al. 1995). A cluster of hypercompact (<240 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 ~10 000 pc^-2 and proposed that W3 IRS5 is a Trapezium cluster (Abt & Corbally 2000) in the making. A high protostellar number density exceeding 10⁶ protostars pc^-3 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 14″ 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 CH₃CN and CH₃OH peak at offset positions (Rodón et al. 2008; Wang et al. 2012).

This paper presents SMA observations of W3 IRS5 at 354 GHz with an angular resolution of 1″–3″, 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.

Fig. 1 a) Continuum emission of W3 IRS5 region at 353.6 GHz (black contours) overplotted with the NICMOS 2.22 μm emission (color scale). The synthesized beam size is , PA − 14°. The contour levels are –3, 3, 5, 7, ..., 19, 23, 27, ... σ with 1σ of 10 mJy beam-1. 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.

2. Observations and data reduction

W3 IRS5 was observed with the SMA¹ (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 , δ(2000) = 62°05′5250. The weather conditions were good with τ_225 GHz ~ 0.1. The reference V_LSR was set to − 39.0 km s^-1. 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 ~0.35 km s^-1 per channel. For the compact configuration dataset, 3C 454.3 was observed as bandpass calibrator. Gain calibration was performed by frequent observations of two quasars, 0136+478 (~ 16° away from the source) and 3C 84 (~22° away from the source). Uranus was adopted as absolute flux calibrator. For the extended configuration dataset, 3C 273, 3C 84 and 3C 273 were observed, respectively, as bandpass, gain and flux calibrators. The adopted total flux density of 3C 273 is 9.6 Jy, which is the mean value reported for five observations² during 2006 January 13–30. We estimate that the uncertainty in absolute flux density is about 20%. The uv-range sampled by the combined dataset is 10 kλ (21″) to 210 kλ (1″).

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 (S/N). Self-calibrated datasets were combined to image the continuum. To optimize sensitivity and angular resolution, a robust weighting³ parameter of 0.25 was adopted for the continuum imaging, resulting in a 1σ rms noise of 10 mJy per 1″ beam. Line identification was conducted by using the compact configuration dataset only with the aid of spectroscopy databases of JPL⁴ (Pickett et al. 1998) and CDMS⁵ (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 PA−14°, 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. 1b). 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 (>150 kλ) and on short baselines (<20 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 >8″ 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 uv 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. 1a, 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–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.

Fig. 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 , (PA − 12°) 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 (PA − 17°). The five crosses represent the positions of SMM1 to SMM5 (cf. 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-1 km s-1 for panel a)–k) and m)–r), and Jy beam-1 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.

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 H₂CS. Others come from HCN, HC¹⁵N, HNCO and CH₃OH. This set of lines covers upper level energies from 43 K to 1067 K and high critical densities (10⁷–10⁸ cm^-3). 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 , PA − 12°) and with uniform weighting (panels (m) to (r), with a resolution of , PA − 17°).

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 H₂CS show a very different spatial distribution, and mainly peak to the north-west of SMM1–5 (Figs. 2i and m). Most of the emission from HCN and its isotopologues comes from SMM1 and SMM2 (Figs. 2g, 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 CH₃OH exclusively peaks toward SMM3 and SMM4 (Fig. 2k). HNCO is detected toward SMM1 and SMM2 only (Fig. 2j). 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 Fig. 3.

Fig. 3 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 J = 7–6, SO N_J = 8₈–7₇ and HCN J = 4–3 are extended and trace the large scale gas kinematics. Figure 4 shows the channel maps of these molecules from − 49 km s^-1 to − 30 km s^-1 and from − 41 km s^-1 to − 35 km s^-1, 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 H₂CS as well. At velocities ranging from − 49 to − 42 km s^-1, the emission from all three molecules can be seen toward zone D.

Of these three molecules, the image quality near the systemic velocity (~− 39.0 km s^-1) 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 − 41 km s^-1 to − 35 km s^-1. The emission of HCN is much weaker than the other species near − 39 km s^-1 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 − 41 km s^-1 to − 38 km s^-1, 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 − 41 km s^-1 to − 37 km s^-1. 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.

Fig. 4 Velocity channel maps of CS J = 7–6, SO NJ = 88–77 and HCN J = 4–3 toward W3 IRS5. Robust weighting (robust = 0) is used for the imaging, resulting resolutions of , PA − 15°, , PA − 13°, and , PA − 14°, for CS, HCN and SO, respectively. The x- and y-axes are RA offset and Dec offset relative to the phase center α(2000) = 02h25m4078, 50 in arcseconds, respectively. The upper three rows shows the velocity range between − 49 km s-1 to − 30 km s-1. The lower three rows show a narrower velocity range from − 41 km s-1 to − 35 km s-1. 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-1). The absolute contour levels for HCN in the upper panel are 3, 8, 13,... σ (1σ = 0.12 Jy beam-1). The absolute contour levels for SO in the upper panel are 3, 13, 23, 43, 63... σ (1σ = 0.11 Jy beam-1). The absolute contour levels for CS in the lower panel are 3, 8, 13,... σ (1σ = 0.18 Jy beam-1). The absolute contour levels for HCN in the lower panel are 3, 7, 11,... σ (1σ = 0.16 Jy beam-1). The absolute contour levels for SO in the lower panel are 3, 13, 23, 43, 63,... σ (1σ = 0.14 Jy beam-1).

Fig. 5 Hanning-smoothed molecular spectra toward molecular zone A (thick line; offset position (− 0.42″, − 0.43″)) and B (thin line; offset position (− 1.89″, − 1.78″)). The compact-configuration dataset and natural weighting is used for all molecules except CS 7–6, SO NJ = 88–77 and HCN 4–3 which are taken from the combined dataset using natural weighting but convolved to the beam size of the compact dataset (, PA − 12°). The apparent offset in velocity seen in 33SO and 33SO2 is due to the blending of hyperfine transitions for which the reference velocity is set to the transition with lowest frequency (see Table 2).

3.2.3. Line profiles

Figure 5 presents the hanning-smoothed spectra toward zone A (thick line; offset position − 0.42″, − 0.43″) and B (thin line; offset position − 1.89″, − 1.78″) 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. HC¹⁵N 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 (V_LSR), FWHM line width (dV) and integrated line intensity (W) 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 8₈–7₇, HCN 4–3, and ³²SO₂ 34_3,31–34_2,32 by comparing our SMA observations with the JCMT observations (Helmich & van Dishoeck 1997). Our SMA observations miss ~70%, ~40%, and ~80% 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 (~94%; Sect. 3.1). The SMA observation of ³²SO₂ 34_3,31–34_2,32 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.

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 (1)where I_ν is the peak flux density, a is the gas-to-dust ratio (100), Ω_b is the beam solid angle, m_H is the mass of atomic hydrogen, κ_ν is the dust opacity per unit mass and B_ν(T_d) is the Planck function at dust temperature T_d. We apply the interpolated value of κ_ν(845μm) ≈ 2.2 cm² g^-1 as suggested by Ossenkopf & Henning (1994) for gas densities of 10⁶–10⁸ cm^-3 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 (2)where S_ν is the total flux density of dust emission and d (1.83 kpc) is the distance to the source. Assuming a spherical source with 1″ 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 3–4 × 10²³ cm^-2, 0.5 M_⊙ and 3 × 10⁷ cm^-3, respectively. Toward SMM3, SMM4, and SMM5 values lower 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^-3. For example, the densities toward SMM1 and SMM2 increase to 2.4 × 10⁸ cm^-3 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 M_⊙).

4.2. Excitation analysis: temperatures

The detected molecular transitions, covering E_u from 43 K to 1067 K and critical densities of ~10⁷–10⁸ cm^-3 (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 CH₃CN J = 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.

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 N_u is given by (3)where g_u is the total degeneracy of the upper state, N_tot is the total molecular column density, Q is the partition function, E_u is the upper level energy, k is the Boltzmann constant and T_rot is the rotational temperature. The left-hand side of Eq. (3) can be derived from the observations via (4)where θ_a and θ_b are the major and minor axes of the clean beam in arcsec, respectively, W is the integrated intensity in Jy beam^-1 km s^-1, g_I and g_K are the spin and projected rotational degeneracies, respectively, ν₀ is the rest frequency in GHz, S 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 E_u/k (Eq. (3)) and a fitted straight line with T_rot and N_tot as free parameters. We adopt a 20% uncertainty of the integrated intensity in the analysis.

Fig. 6 a) Rotation diagram of SO2 (filled squares) and 34SO2 (filled triangles) toward zone A. The straight lines are the model fit. The derived rotational temperatures are indicated. b) Rotation diagram of CH3CN J = 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).

Fig. 7 RADEX excitation analysis toward zone A traced by 34SO2 (left four panels) and zone B traced by CH3CN (right four panels). The χ2 surfaces near the best-fit solution on each parameter axis are displayed.

Figure 6a 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 E_u ranging from 88 to 580 K, and indicates a T_rot of 141 ± 25 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 T_rot = 89 ± 10 K and T_rot = 197 ± 55 K. The gas component traced by SO₂ (E_u = 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. 6a) 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 T_rot. However, we detected the 34_3,31 − 34_2,32 transition of both SO₂ and ³⁴SO₂. Equation (4) yields a N_u(SO₂)/N_u(³⁴SO₂) ratio of ~20 at 50 − 300 K, consistent with the ³²S/³⁴S ratio of 22 in the interstellar medium (Wilson & Rood 1994), suggesting that the 34_3,31 − 34_2,32 transition of SO₂ (and its isotopologues) is optically thin. The 24_2,22 − 23_3,21v₂ = 1 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 34_3,31 − 34_2,32 transition. These considerations suggest the a temperature gradient in zone A explains the rotation diagrams, with a cooler region at ~100 K and a warmer region of ≥200 K (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 CH₃CN J = 12 − 11 lines obtained by Wang et al. (2012) with the SMA at a resolution of 4.0″ × 2.6″. This molecule peaks toward zone B exclusively and shows a much narrower FWHM line width (~1.6 km s^-1) compared to ³⁴SO₂ (~4 km s^-1), suggesting a different molecular condensation in zone B. Figure 6b shows the CH₃CN J = 12 − 11 K = 0 − 4 rotation diagram. The scatter in the data indicates that the lines are not optically thin. Therefore, we cannot carry out rotation diagram analysis on the CH₃CN 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 (T_kin), column density (N_tot) and H₂ volume density (n_H2) 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 CH₃CN J = 12 − 11 lines (Wang et al. 2012). For ³⁴SO₂, we adopted the collisional rate coefficients from SO₂ (Green 1995)⁶ as an approximation since both species have similar molecular properties. The collisional rate coefficients of CH₃CN are taken from Green (1986). To find the best-fit solution, we applied χ² minimization to the 4-dimensional parameter space (T_kin, N_tot, n_H2 and f, 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 E_u = 580 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 (>10^9 − 10 cm^-3) 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 CH₃CN is only about 10⁵ cm^-3, lower than the critical density of a few 10⁶ cm^-3, implying subthermal excitation. The beam filling factor is about 0.8. The best-fit model for CH₃CN is plotted as filled squares in Fig. 6b.

Alternatively, the excitation conditions toward zone A and B can be estimated via the line ratio of the 19_1,19 − 18_0,18 over the 10_4,6 − 10_3,7 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^-1 (the average FWHM line widths in zones A and B are 6.0 km s^-1 and 4.0 km s^-1, respectively) in the RADEX calculations, we calculate the intensity ratio as function of T_kin and N_tot for four different values of n_H2 (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 >10⁹ cm^-3 traced by ³⁴SO₂ toward zone A and B. The kinetic temperatures are ≥150 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 (n_H2 ≥ 10⁹ cm^-3) in zones A and B. There is a temperature gradient along the line of sight in zone A, characterized by a cool region of ~100 K and a warm region of ≥200 K. Toward zone B, the temperature gradient is less prominent. However, a quiescent (FWHM line width 1.6 km s^-1) and less dense region (n_H2 ~ 10⁵ cm^-3) traced by CH₃CN is also present in zone B.

4.2.3. Molecular column density

Fig. 8 RADEX analysis of the 34SO2 line ratios toward zone A and B assuming a FWHM line width 5.0 km s-1. The intensity ratios 191,19 − 180,18/104,6 − 103,7 are plotted in color scale. The observed integrated intensities of the 191,19 − 180,18 transition are plotted in curves (zone A: ~63 Jy beam-1 km s-1, zone B: ~15 Jy beam-1 km s-1). Toward both zones, the line ratio is about 1.8 (Table 3). Our results suggest that a high H2 volume density (> 109 cm-3) is needed to reproduce the observed line ratios.

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). (5)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 T_rot from 100 K to 200 K). As a consistency check, we estimate the line center opacity of each transition via (6)where c is the speed of light, A_ul the Einstein-A coefficient, and ΔV the FWHM line width. All lines are consistent with the optically thin assumption (τ ≪ 1); only SO has line center opacities as large as τ ~ 0.3–0.5). Therefore, a lower limit of SO column density is derived.

Fig. 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.

A more representative value for SO can be derived from ³³SO by assuming the isotopic ratio ³²S/³³S of ~132 (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 1.5 × 10²³ cm^-2 and 8.9 × 10²² cm^-2 toward zone A and B, respectively (derived from the continuum image in Fig. 2l 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 10^-6, several orders of magnitude higher than found in dark clouds (few 10^-9, 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 CH₃OH which is a factor of 10 more abundant in zone B.

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), HC¹⁵N (E–W) and HCN v₂ = 1 (NE–SW). CH₃OH 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 CH₃OH, 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 HC¹⁵N 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 (~− 39 km s^-1). 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.

Fig. 10 Channel-by-channel locations of the emission peaks of the SO2 343,31 − 342,32 and HCN 4–3 v2 = 1. The emission peaks are color-coded with respect to their VLSR in km s-1 shown in the color wedge. Two distinct, orthogonal velocity gradients are seen.

In Fig. 9, we see that the velocity gradients seen in ³⁴SO₂ and HCN v₂ = 1 are roughly perpendicular to each other. To confirm this, we fit the emission of HCN 4–3 v₂ = 1 (E_u = 1067 K) and SO₂ 34_3,31 − 34_2,32 (E_u = 582 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.

Fig. 11 Integrated intensity maps at different velocity ranges (blue: −50 to −43 km s-1; green: −43 to −35 km s-1; red: −35 to −28 km s-1). 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 ≥3σ noise level. The markers are identical to the ones in Fig. 1.

In our data, only the CS 7–6, HCN 4 − 3 and SO N_J = 8₈ − 7₇ emission show a wide velocity range (about 30–40 km s^-1; 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^-1, green: −43 to −35 km s^-1 and red: −35 to −28 km s^-1). We took the combined-configuration dataset and used natural weighting for the analysis. The images are plotted with contours starting well above 3σ 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 v₂ = 1 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 (− 8″, 0″) 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.

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 M_⊙). 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. 1a, 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 (~0.3 M_⊙) suggests that SMM3–5 may be forming low-mass stars or are starless cores.

For a virialized star-forming core (M_core + M_star ≈ M_vir), the 1D virial velocity dispersion (σ_vir; cf. McKee & Ostriker 2007) is (7)where r is the radius of the source, G is the gravitational constant, and M_vir is the virial mass. Assuming a core radius of 0.5″, core gas masses of ~0.3 M_⊙, and stellar masses of SMM3–5 less than 8 M_⊙, the 1D virial velocity dispersions should be no more than 1.3 km s^-1. From our observations of CH₃OH line toward SMM3/4 and CS line toward SMM5, the inferred 1D velocity dispersions (σ = ΔV/) are 1.0 km s^-1 and 1.4 km s^-1. 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 M_⊙) and the projected distance between the SMM sources are about 1″ − 2″ (1800–3600 AU). This observed structure, in the meantime, is still embedded in a large massive core with size about 1′ and mass of several hundred M_⊙ (W3-SMS1 in the SCUBA 850 μm image, Di Francesco et al. 2008; Wang et al. 2012). At 1′ scale, Megeath et al. (1996) found that there are ~300 low-mass stars in the core surrounding W3 IRS5 with a stellar density of few 1000 pc^-3. Therefore, it is interesting to investigate the cloud structure from 1′ scale to 1″ 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 850 μm image (Di Francesco et al. 2008; Wang et al. 2012). The Jeans length and Jeans mass are expressed as (cf. Stahler & Palla 2005, Eqs. (9.23) and (9.24)): (8)and (9)where T is the kinetic temperature and n_H2 is the H₂ volume density. These two quantities are more sensitive to the adopted H₂ volume density. To estimate the 1′ scale H₂ density, we measured the total flux density in a 1′ by 1′ 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 M_⊙ assuming a dust temperature of 40 K. The mean H₂ volume density is 2.4 × 10⁵ cm^-3. Therefore, at 1′ we derived the Jeans length and Jeans mass to be 16 000 AU (~9′′) and 2 M_⊙, respectively. Interestingly, these two characteristic quantities seem to match the observed spatial separations of low-mass stars (18 000–12 000 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 M_⊙), which is expected given a typical efficiency of low-mass star formation ~30% (Alves et al. 2007). Toward the central regions of W3 IRS5, the temperature (~140 K) and density (few 10⁷ cm^-3) are higher than the values on the scale of the entire core. The corresponding Jeans length and Jeans mass are 2700–3800 AU (1.5″–2.1″) and ~1 M_⊙, 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, (10)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 1″ in diameter based on Fig. 1a. We adopt stellar masses of 20 M_⊙ for SMM1 and SMM2 based on the total luminosity of 2 × 10⁵L_⊙. 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 >1700 yr for SMM3–5. These imply mass accretion rates of 5 × 10^-4M_⊙ yr^-1 for SMM1/2, and < 2 × 10^-4M_⊙ yr^-1 for SMM3–5. The upper limits for the latter follow from the upper limit on the embedded stellar masses of 8 M_⊙. 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 (1″ 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 >1″. 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 1′ W3–SMS1 condensation, with a total mass of 950 M_⊙ and an embedded stellar mass of ≲50 M_⊙, is 72 000 yr and the free-fall mass accretion rate is 1 × 10^-2M_⊙ yr^-1. 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 (11)(Bonnell et al. 2001), where v_inf is the infall gas velocity, ρ is the density of the gas surrounding the proto-Trapezium, and R_acc is the radius within which the infalling gas is accreted onto a specific star. We assume n_H2 = 10⁶ cm^-3 to compute ρ, which is intermediate between the density of the surrounding 1′ diameter core (10⁵ cm^-3) and the core of the inner 1″ (10⁷ cm^-3). Given FWHM line widths of 2–5 km s^-1, we adopt a v_inf equal to the 3D gas velocity dispersion implied by the FWHM line width (), which is 1.5–3.7 km s^-1 with a mean value of 2.6 km s^-1. The uncertainty of v_inf gives a factor of <3 in error of accretion rate. The possible values for R_acc are the radius set by the motions of the stars through the gas (Bondi-Hoyle radius) and the radius set by tidal interaction of the stars with the gravitational potential of the cluster core (tidal radius). The actual accretion rate is given by the minimum of these two radii (Bonnell et al. 2001). To estimate R_acc due to tidal interactions, we use the equation from Bonnell et al. (2001) and rewrite in terms of the gas density around the star, (12)In this case, are 4700 AU and 3500 AU for a 20 M_⊙ and 8 M_⊙ central star embedded in a core with gas density 10⁶ cm^-3, respectively. The accretion rates are then 2 × 10^-4M_⊙ yr^-1 and 1 × 10^-4M_⊙ yr^-1 for the 20 M_⊙ and 8 M_⊙ stars, respectively. Alternatively, the motions of the stars through the cloud can limit the gas accretion; only the gas within a radius where the escape velocity is less than the velocity of the star relative to the gas can be accreted. The R_acc in this case can be approximated by the Bondi-Hoyle radius: (13)The resulting are 5200 AU and 2100 AU for the 20 and 8 M_⊙ stars, and the mass infall rates are 3 × 10^-4 M_⊙ yr^-1 for a stellar mass of 20 M_⊙ and 4 × 10^-5M_⊙ yr^-1 for a stellar mass of 8 M_⊙.

As a result, from our basic analysis, global collapse of the core and subsequent accrete of the material from the collapse onto the stars can sustain infall rates in excess of 10^-4 M_⊙ yr^-1 for SMM1/2 and are sufficient to sustain high mass star formation. For SMM3–5, the rates are 10 times lower. However, the local density around SMM3–5 is closer to 10⁷ cm^-3, the actual densities and the accretion rates may be higher. The fact that the stellar accretion rates is two orders of magnitude less than the free-fall accretion rate for the large cores suggests that there is sufficient infall of material to feed the accretion of the individual stars, even if the actual collapse times is ten free fall times.

In summary, we find two alternatives. Either the massive stars accrete locally from their local cores; in this case the small core masses imply W3 IRS5 is at the very end stages (1000 yr) of infall and accretion. Alternatively, the stars are accreting from the global collapse of a massive, cluster forming core. This later scenario is similar to the competitive accretion models of massive star formation (e.g., Bonnell & Bate 2006). For the two massive objects SMM1/2, the observed densities and velocity widths are consistent with those needed for the global collapse to feed accretion onto the massive stars at rates of ~ 10^-4 M_⊙ yr^-1. For the remaining objects, SMM 3–5, a lower infall rate is estimated; perhaps these are forming intermediate mass stars, although the rate of accretion may rise as they increase in mass as precede by the models of competitive accretion. However, a clear detection of the infall of gas has not yet been observed. This makes the W3 IRS5 cluster a perfect region to further test the competitive accretion model, if the kinematics of the gas motions can be traced on scales from 1′ down to <1″. Future single-dish and interferometric observations of this region, e.g., combining SMA or PdBI observations with IRAM 30 m measurements, are required to fully map the gas flow across these scales.

5.2. The molecular environment of W3 IRS5

5.2.1. Distinct molecular zones

A JCMT spectral line survey characterized W3 IRS5 as a remarkable massive star-forming region with rich spectral features from sulfur-bearing molecules (Helmich et al. 1994; Helmich & van Dishoeck 1997). Helmich et al. (1994) also reported that the CH₃OH abundance is remarkably low, even compared to typical dark clouds. From our SMA data, and the results from Wang et al. (2012), we find that sulfur-bearing molecules such as SO, SO₂ and their isotopologues, strongly peak toward the submillimeter sources SMM1 and SMM2 (zone A) with extended emission toward SMM3, SMM4 and SMM5 (Fig. 2). This suggests that SMM1 and SMM2 are the driving center of the sulfur chemistry. Interestingly, typical hot-core molecules such as CH₃OH (Fig. 2k) and CH₃CN (Wang et al. 2012) peak toward SMM3/4 (zone B) exclusively, underlining the chemical differences between zones A and B. Moreover, toward zone B, there are two distinct regions along the line of sight with different physical characteristics. The narrow FWHM line widths seen in CH₃OH (2.3 km s^-1; Table 3) and CH₃CN (1.6 km s^-1; Wang et al. 2012) suggest that this “hot core” source is less turbulent compared to the gas component traced by SO and SO₂, which have larger line widths (~4–7 km s^-1; Table 3). As implied by the excitation analysis (Sect. 4.2.2), this “hot core” is less dense (~10⁵ cm^-3) than the gas traced by ³⁴SO₂ (≥10⁹ cm^-3). We suggest that this “hot core” may created by feedback from the star-forming activity associated with SMM1 and SMM2. It is, however, hard to tell if this core is heated internally or externally since we do not find significant temperature differences between ³⁴SO₂ and CH₃CN (Sect. 4.2.2). In addition, there is another interesting chemical signature in W3 IRS5, where SO and SO₂ are peaked toward zone A, while CS and H₂CS are peaked toward zone C (Fig. 2). Although these species are chemically related (e.g. Charnley 1997), this spatial de-correlation implies that different chemical processes are involved, such as hot core versus shock chemistry (Hatchell & Viti 2002).

5.2.2. Abundances of SO, SO₂ and CH₃OH

It has been proposed that the abundance ratios of sulfur-bearing species can be used to measure the time elapsed since the start of ice-mantle evaporation in star-forming cores if proper physical conditions are applied (e.g. Charnley 1997; van der Tak et al. 2003; Wakelam et al. 2004). Unfortunately, the rough estimates of the SO and SO₂ abundances provided by our data limit us to a qualitative comparison only. SO and SO₂ are abundant toward W3 IRS5. The fractional abundance of SO and SO₂ (Table 4) are, respectively, 8.6 × 10^-6 and 9.4 × 10^-7 toward zone A, and 5.4 × 10^-6 and 2.8 × 10^-7 toward zone B. This corresponds to SO₂/SO ratios of 0.1 and 0.05 for zones A and B, respectively. The uncertainty of these ratios are factors of few, and SO₂ is clearly less abundant than SO toward W3 IRS5.

Charnley (1997), Hatchell et al. (1998) and Wakelam et al. (2004) present chemical models that relate the SO₂/SO ratio to the time since the onset of grain-mantle evaporation. Although these models are very sensitive to the initial mantle composition, cosmic ray ionization rate, density and temperature, these models show that the evolution of the fractional abundances of SO and SO₂ peak after roughly 10⁴–10⁵ yr after grain-mantle evaporation. Such time scales are consistent with the free-fall time for the entire cluster calculated in Sect. 5.1.3 of (0.8–1.4) × 10⁵ yr.

Toward the emission peak in zone B, we derive a CH₃OH abundance of 5.1 × 10^-8, lower by a factor of 10–100 compared to the abundance found toward the Orion Compact Ridge (Menten et al. 1988). This low CH₃OH abundance implies that the “hot core” toward zone B is chemically young since large amounts of CH₃OH are still observed in the solid state (Allamandola et al. 1992). Since zone B does show emission from SO and SO₂, which is indicative of outflow activity (Fig. 11), it is possible that the “hot core” has just been turned on by external heating from zone A.

5.3. Non-detection of CO⁺: an indication of FUV origin

Just as much as detected emission can provide information, so can the non-detection of emission tell us something. Molecular ions and radicals like CN, NO, CO⁺, SO⁺ and SH⁺, are thought to trace high energy (far-ultraviolet, FUV, and X-ray) radiation toward star-forming regions (Sternberg & Dalgarno 1995; Maloney et al. 1996). Toward W3 IRS5, Stäuber et al. (2007) reported detection of two hyperfine transitions of CO⁺N = 3–2, F = 5/2–3/2 (353.7413 GHz) and F = 7/2–5/2 (354.0142 GHz) with the James Clerk Maxwell Telescope (JCMT) at 14″ resolution. They concluded that CO⁺ is present in FUV-irradiated cavity walls which may be part of the outflows along the line of sight, similar to AFGL 2591 (Stäuber et al. 2007; Benz et al. 2007; Bruderer et al. 2009). However, the detection of X-ray emission toward W3 IRS5 (Hofner et al. 2002) with a luminosity of ~ 5 × 10³⁰ erg s^-1 implies that CO⁺ may have an X-ray origin as well. The geometrical dilution of this X-ray flux down to the cosmic ray ionization level suggests that the source size of CO⁺ should be about 0.7″ if X-rays are the dominant source of CO⁺. Therefore interferometric observations would be helpful to determine if CO⁺ can have an X-ray origin.

The same set of CO⁺N = 3–2 hyperfine transitions are also covered in our SMA data. Emission at 353.7413 GHz is clearly detected toward zone A and shows a compact structure (cf. Fig. 2d). However, we do not detect the transition at 354.0142 GHz (Fig. 12). In LTE the intensity ratio of the CO⁺N = 3–2 hyperfine transitions (including F = 3/2–1/2 at 353.0589 MHz) is 1:14:20. A similarly compact source at 354.0142 GHz should therefore have been detected. Also, the CO⁺ abundance we derive from the 353.7413 GHz intensity (~1.3 ×  10^-9) is much larger than the prediction from Stäuber et al. (2007) and Stäuber & Bruderer (2009) (~ 10^-12 − 10^-10). We therefore rule out CO⁺ as the origin of the emission at 353.7413 GHz.

Instead, we assign the line near 353.741 GHz to ³³SO₂ 19_4,16 − 19_3,17 transition, as was suggested by Stäuber et al. (2007). The anomalous CO⁺ line ratio reported by (Stäuber et al. 2007) is explained by blending of the CO⁺ line with ³³SO₂ at 353.741 GHz. If this assignment is correct, our SMA data should contain three other lines from ³³SO₂ (with hyperfine transitions), near 344.031 GHz (E_u = 397 K), 344.504 GHz (E_u = 232 K) and 354.246 GHz (E_u = 297 K). These lines are not detected, but they have either smaller line strengths (344.504 GHz) or decreased intensity due to hyperfine splitting (344.031, 354.246 GHz). In fact, the detected transitions near 353.741 GHz is expected to be the strongest in the passband. The assignment of this line to ³³SO₂ is further supported by the detection of its 11_1,11 − 10_0,10 line by Wang et al. (2012). From our data, we find a column density ratio of ³²SO₂/³³SO₂ of ~147 assuming optically thin lines and LTE conditions (Table 4). This ratio is consistent with the ISM ³²S/³³S ratio of ~132 (Wilson & Rood 1994; Chin et al. 1996).

The fact that CO⁺, detected in the JCMT beam, and undetected in our SMA data, suggests an extended (>21″) distribution of the emission consistent with an origin of the emission in the FUV-irradiated cavity walls, and rules out an X-ray origin. Alternatively, the non-detection of CO⁺ in our ~2.5″ beam may be explained if CO⁺ is destroyed rapidly with H₂ or electrons (Stäuber & Bruderer 2009, and the references therein) in the high density (≥10⁹ cm^-3) material near SMM1 and SMM2 identified by our data.

Fig. 12 Hanning-smoothed spectrum taken from the compact-configuration dataset toward zone A (offset position: − 0.42″, − 0.43″). The black lines are the observed data, while the red lines are the model spectrum of CO+N = 3–2, F = 3/2–1/2 (353.0589 GHz), 5/2–3/2 (353.7413 GMHz) and 7/2–5/2 (354.0142 GHz) under LTE assumption. The LTE intensity ratio of the hyperfine transitions is ~1:14:20 (increasing in frequency). The F = 7/2–5/2 transition at 354.0142 GHz is not detected in our SMA data.

6. Conclusions

Our conclusions can be summarized as follows.

• 1.

  The 1″ resolution continuum image at 353.6 GHz shows 5 compact submillimeter sources SMM1–5 (Fig. 1). SMM1 and SMM2 have embedded massive stars (~20 M_⊙), judging from the detection of compact cm-wave emission. The other sources do not contain massive stars currently, and may be forming low-mass stars or simply be starless.

• 1.

  The gas densities of SMM1–5 are high as 10⁷ cm^-3, but their core masses are surprizingly low, 0.2 − 0.6 M_⊙.

• 2.

  If the massive stars in W3 IRS5 accrete locally from the surrounding core, the derived small core masses of the SMM sources imply that the main accretion phase is almost concluded.

• 3.

  From Jeans analysis, the cloud structures from 1′ scales down to 1″ scales are likely determined by gravitational fragmentation in turbulent environment.

• 4.

  The free-fall accretion rates toward SMM1 and SMM2 are about 5 × 10^-4 M_⊙ yr^-1, while SMM3–5 have accretion rates of less than 1–2 × 10^-4M_⊙ yr^-1 if their stellar masses are 8 M_⊙ or 2–4 × 10^-5 M_⊙ yr^-1 if they are starless.

• 5.

  If star formation in W3 IRS5 follows the competitive accretion model, global collapse of the core and subsequent accretion of the material from the collapse onto the stars can sustain infall rates in excess of 10^-4M_⊙ yr^-1 for SMM1/2 and are sufficient to sustain high-mass star formation.

• 6.

  From the molecular line images, we identified four molecular zones (Fig. 3). Zone A and B are the major places where sulfur chemistry and hot core chemistry are taking place. Zone C and D, however, seem kinematically unrelated to zone A and B. There is a “hot core” traced by CH₃OH and CH₃CN toward the line of sight of zone B.

• 7.

  The large abundances of SO and SO₂ with respect to H₂ (few 10^-7 to 10^-6) derived in zone A and B are indicative of active sulfur chemistry.

• 8.

  The low abundance of CH₃OH (5 × 10^-8) in the hot core toward zone B implies that it may be recently heated either internally or externally by feedback from zone A. Zone B seems chemically younger than zone A.

• 9.

  The non-detection of CO⁺J = 3–2 on small scales supports the idea that CO⁺ is formed in the FUV-irradiated outflow cavity walls, and implies that X-rays either do not contribute significantly to the CO⁺ production or that this molecule is rapidly destroyed again in the dense material around SMM1 and SMM2.

Acknowledgments

We thank the anonymous referee and the editor Malcolm Walmsley for reviewing our paper. We also thank the SMA staff for conducting our observations. The research of K.-S.W. at Leiden Observatory is supported through a Ph.D. grant from the Nederlandse Onderzoekschool voor Astronomie (NOVA). This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

All Tables

All Figures

Fig. 1 a) Continuum emission of W3 IRS5 region at 353.6 GHz (black contours) overplotted with the NICMOS 2.22 μm emission (color scale). The synthesized beam size is , PA − 14°. The contour levels are –3, 3, 5, 7, ..., 19, 23, 27, ... σ with 1σ of 10 mJy beam-1. 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.

In the text

Fig. 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 , (PA − 12°) 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 (PA − 17°). The five crosses represent the positions of SMM1 to SMM5 (cf. 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-1 km s-1 for panel a)–k) and m)–r), and Jy beam-1 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.

In the text

	Fig. 3 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.
In the text

Fig. 4 Velocity channel maps of CS J = 7–6, SO NJ = 88–77 and HCN J = 4–3 toward W3 IRS5. Robust weighting (robust = 0) is used for the imaging, resulting resolutions of , PA − 15°, , PA − 13°, and , PA − 14°, for CS, HCN and SO, respectively. The x- and y-axes are RA offset and Dec offset relative to the phase center α(2000) = 02h25m4078, 50 in arcseconds, respectively. The upper three rows shows the velocity range between − 49 km s-1 to − 30 km s-1. The lower three rows show a narrower velocity range from − 41 km s-1 to − 35 km s-1. 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-1). The absolute contour levels for HCN in the upper panel are 3, 8, 13,... σ (1σ = 0.12 Jy beam-1). The absolute contour levels for SO in the upper panel are 3, 13, 23, 43, 63... σ (1σ = 0.11 Jy beam-1). The absolute contour levels for CS in the lower panel are 3, 8, 13,... σ (1σ = 0.18 Jy beam-1). The absolute contour levels for HCN in the lower panel are 3, 7, 11,... σ (1σ = 0.16 Jy beam-1). The absolute contour levels for SO in the lower panel are 3, 13, 23, 43, 63,... σ (1σ = 0.14 Jy beam-1).

In the text

Fig. 5 Hanning-smoothed molecular spectra toward molecular zone A (thick line; offset position (− 0.42″, − 0.43″)) and B (thin line; offset position (− 1.89″, − 1.78″)). The compact-configuration dataset and natural weighting is used for all molecules except CS 7–6, SO NJ = 88–77 and HCN 4–3 which are taken from the combined dataset using natural weighting but convolved to the beam size of the compact dataset (, PA − 12°). The apparent offset in velocity seen in 33SO and 33SO2 is due to the blending of hyperfine transitions for which the reference velocity is set to the transition with lowest frequency (see Table 2).

In the text

	Fig. 6 a) Rotation diagram of SO2 (filled squares) and 34SO2 (filled triangles) toward zone A. The straight lines are the model fit. The derived rotational temperatures are indicated. b) Rotation diagram of CH3CN J = 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).
In the text

	Fig. 7 RADEX excitation analysis toward zone A traced by 34SO2 (left four panels) and zone B traced by CH3CN (right four panels). The χ2 surfaces near the best-fit solution on each parameter axis are displayed.
In the text

Fig. 8 RADEX analysis of the 34SO2 line ratios toward zone A and B assuming a FWHM line width 5.0 km s-1. The intensity ratios 191,19 − 180,18/104,6 − 103,7 are plotted in color scale. The observed integrated intensities of the 191,19 − 180,18 transition are plotted in curves (zone A: ~63 Jy beam-1 km s-1, zone B: ~15 Jy beam-1 km s-1). Toward both zones, the line ratio is about 1.8 (Table 3). Our results suggest that a high H2 volume density (> 109 cm-3) is needed to reproduce the observed line ratios.

In the text

Fig. 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.

In the text

	Fig. 10 Channel-by-channel locations of the emission peaks of the SO2 343,31 − 342,32 and HCN 4–3 v2 = 1. The emission peaks are color-coded with respect to their VLSR in km s-1 shown in the color wedge. Two distinct, orthogonal velocity gradients are seen.
In the text

Fig. 11 Integrated intensity maps at different velocity ranges (blue: −50 to −43 km s-1; green: −43 to −35 km s-1; red: −35 to −28 km s-1). 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 ≥3σ noise level. The markers are identical to the ones in Fig. 1.

In the text

Fig. 12 Hanning-smoothed spectrum taken from the compact-configuration dataset toward zone A (offset position: − 0.42″, − 0.43″). The black lines are the observed data, while the red lines are the model spectrum of CO+N = 3–2, F = 3/2–1/2 (353.0589 GHz), 5/2–3/2 (353.7413 GMHz) and 7/2–5/2 (354.0142 GHz) under LTE assumption. The LTE intensity ratio of the hyperfine transitions is ~1:14:20 (increasing in frequency). The F = 7/2–5/2 transition at 354.0142 GHz is not detected in our SMA data.

In the text