© ESO, 2014

1. Introduction

NGC 7538 is a well-studied optically visible H II region that is located at a distance of 2.65 kpc (Moscadelli et al. 2009), and surrounded by very massive star–forming complexes (Minn & Greenberg 1975; Werner et al. 1979; Read 1980; Rots et al. 1981; Campbell & Thompson 1984; Reid & Wilson 2005). NGC 7538 IRS 1–3 form a cluster of embedded massive protostars with luminosities of ~10⁴L_⊙ located 2 arcmin (1.5 pc) southeast of the center of the optical H II region (Wynn-Williams et al. 1974). IRS 1, 2, and 3 have luminosities of 8, 5, and 0.6 × 10⁴L_⊙, respectively, and all of them have ultra compact (UC) H ii regions (Campbell & Thompson 1984; Campbell & Persson 1988).

IRS 1, whose luminosity is equivalent to a O7.5 main sequence star, is surrounded by a UC H ii region with a double-lobed structure along the N−S direction with a size of 0′′.2 (≃500 AU) (Campbell 1984; Gaume et al. 1995). On these scales, the radio emission is dominated by a collimated, ionized wind that exhibits time variability (Franco-Hernández & Rodríguez 2004; Sandell et al. 2009). This is supported by the extremely broad line widths of radio recombination lines (Gaume et al. 1995; Keto et al. 2008). A bipolar molecular outflow extending in the NW−SE direction is detected in CO (Kameya et al. 1989; Davis et al. 1998; Qiu et al. 2011, hereafter QZM11). The outflow has created a cavity that is detected in the near and mid-IR on scales of a few arcseconds N of the protostar (De Buizer & Minier 2005; Kraus et al. 2006). Despite the existence of the UC H ii region, the protostar is probably still actively accreting gas with a high infall rate (mm inverse P−Cygni profiles yield infall rates of ~10^-3M_⊙ yr^-1), and it probably has a massive circumstellar disk (QZM11; Klaassen et al. 1992; Sandell et al. 2009; Beuther et al. 2012). This source is associated with a rich variety of masers, most of them arising from the interaction zone of the dense molecular gas with the ionized gas and from the outflow (Rots et al. 1981; Dickel et al. 1982; Johnston et al. 1989; Schilke et al. 1991; Minier et al. 1998; Hoffman et al. 2003; Hutawarakorn & Cohen 2003; Pestalozzi et al. 2006; Galván-Madrid et al. 2010; Surcis et al. 2011; Hoffman & Seojin Kim 2011).

IRS 2 and IRS 3 are not associated with a dusty envelope and are probably in a more evolved phase. IRS 2 is associated with an O5 star, and it is probably the most evolved source. Bloomer et al. (1998) propose that its stellar winds shock the surrounding material, generating a “stellar wind bow shock” visible as a shell in H II and Fe ii. IRS 3 is associated with an O6−O9 star that might power one or more CO outflows (QZM11; Ojha et al. 2004).

The dense molecular gas around IRS 1−3 appears to have a filamentary morphology, with an arc-like shape SE of IRS 1 (QZM11; Pratap et al. 1990). Kawabe et al. (1992) interpret this dense molecular gas filament as part of an expanding ring-like structure with a radius of 0.25 pc and a mass of 230 M_⊙, created or piled up by the strong protostellar winds. QZM11 present a study of the NGC 7538 IRS 1−3 region using the Submillimeter Array (SMA) in the 1.3 mm wave band. They found star-forming cores deeply embedded within the filamentary, dense molecular gas cloud of IRS 1−3 and multiple molecular outflows.

Previous submm polarization observations by Momose et al. (2001) were done with single-dish bolometers toward IRS 1−3 and IRS 11 (a younger source, located S of IRS 1). The two sources show striking differences in the polarization properties. Thus, while IRS 11 exhibits an extremely well-ordered magnetic field and high degrees of polarization, in IRS 1 the field appears locally disturbed, and the degrees of polarization are much lower than those of IRS 11. They interpret this as an evolutionary effect (more ordered fields are observed in younger sources), which has also been observed in other high-mass star-forming regions (Girart et al. 2009, 2013; Tang et al. 2009a,b).

Polarimetric observations allow us to study the magnetic fields on the relevant scales (100−10⁴ AU) where the star formation takes place (Girart et al. 1999, 2006; Rao et al. 2009, 2014; Hull et al. 2013; Liu et al. 2013; Qiu et al. 2013). Recently, Zhang et al. (2014) have presented a large sample of massive, clustered, star-forming clumps and showed that magnetic fields play an important role during the formation of dense cores on spatial scales of 0.01−0.1 pc. In this paper, we present SMA observations carried out at 345 GHz toward the massive cluster NGC 7538 IRS 1. Section 2 briefly describes observations and data reduction. Section 3 presents the results of the observations. Sections 4 and 5 contain the analysis and discussion. Finally, in Sect. 6 we draw our main conclusions.

Fig. 1 Contour map of the SMA dust continuum emission map towards NGC 7538 IRS 1 at 878 μm overlapped with the gray image of the polarized intensity. The blue and red segments show the magnetic field direction obtained from the polarization data with cutoff at rms of 2 and 3σ, respectively (see Sect. 3.2). Solid contours go from 3σ to 33σ in steps of 3σ, where σ = 0.017 Jy beam-1 is the rms noise of dust continuum. The dashed contours around IRS 1 (MM1) go from 50σ to 300σ in steps of 50σ. The scale of the polarized intensity image is shown in the righthand side of the figure (the units are Jy beam-1). The black and white crosses show the peak positions of the dust continuum sources (Table 1). The name of the dust continuum sources are also shown. Red circles mark the positions of IRS 1−3. The solid circle shows the FWHM of the SMA primary beam at the observed frequency. The physical scale of the map and synthesized beam are shown in the bottom left corner of the panel.

2. Observations and data reduction

The SMA observations were undertaken on 2005 July 7 with seven antennas in the compact configuration. The weather was good during the observations, with system temperatures in the range of 200−300 K. The phase center used for NGC 7538 IRS 1 was α(J2000.0) $\mathrm{=}{\mathrm{23}}^{\mathrm{h}}{\mathrm{13}}^{\mathrm{m}}\mathrm{43}\stackrel{\mathrm{s}}{.}\mathrm{359}$, δ(J2000.0) (Davis et al. 1998). A single receiver was used and tuned to a local oscillator frequency of 341.6 GHz (878 μm), with a total bandwidth of 2 GHz per sideband, covering a frequency range of 335.51−337.49 GHz in the lower sideband and about 345.51−347.89 GHz in the upper sideband. At these frequencies, the full width at half maximum (FWHM) of the antennae’s primary beam is . The 2 × 2 GHz correlator was configured to sample the aforementioned frequency ranges with a uniform spectral resolution of 0.81 MHz (≃0.7 km s^-1). The standard reduction of the data was done by using the IDL MIR package and selecting the QSO 3C454.3 for the bandpass calibration and BLLAC for the gain calibration. The polarization calibration was performed by observing 3C454.3 over a wide parallactic angle range (−70° to +70°). This allowed us to correct for the instrumental polarization using the MIRIAD software package at an accuracy of 0.1% (Marrone & Rao 2008). We used the MIRIAD task UVLIN to separate the continuum and the line data in the u,v domain. The line-free channels were selected by inspecting MM1 because it is the core with the most molecular lines. These channels were used for the final maps, minimizing the possible contamination from molecular line emission. Self-calibration on NGC 7538 IRS 1 was performed using the Stokes I continuum data for each baseline independently. The derived gain solutions were applied to the molecular line data.

The final maps were obtained using natural weighting, which yielded a synthesized beam of with a position angle of 34°. The continuum map sensitivity is σ_I = 0.017 Jy beam^-1 for Stokes I and σ_pol = 0.010 Jy beam^-1 for Stokes Q and U. For the molecular line data, we made channel maps with a spectral resolution of 1.4 km s^-1 that resulted in a sensitivity of σ = 0.25 Jy beam^-1 per channel.

Fig. 2 Moment maps of molecular transitions with extended emission. Rows: moments 0, 1, and 2, in descending order, labeled on the lefthand side of the figure. Columns: from left to right, C17O 3−2, H13CO+ 4−3, SO 8−7, and C34S 7−6, respectively, labeled on the top of the figure. Color map: moment maps. Common scale is shown on the righthand side of the figure. Contours: 3σ, 9σ, and 27σ continuum emission levels, where σ = 0.017 Jy beam-1.

3. Results

3.1. Continuum emission

Figure 1 presents the 878 μm continuum SMA map of NGC 7538 IRS 1. A total of fourteen sources are detected, including IRS 1 (MM1 in this figure). MM1 is the brightest source at 878 μm, with an intensity of 5.24 Jy beam^-1. The rest of the continuum sources have intensities in the 0.11−0.56 Jy beam^-1 range. Our mass sensitivity at a 3σ level for a typical T_dust = 40 K is ~0.9 M_⊙ beam^-1, which is twice better than for QMZ11. (Our second contour in Fig. 1 roughly compares to their first contour in Fig. 1.) As a result, only nine sources were previously detected by QMZ11. MM8, detected by QZM11, is not detected at 878 μm because it falls outside the FWHM of the primary beam at this wavelength. The 878 μm position of MM9 is offset by with respect to the value given at 1.3 mm by QZM11. This is probably due to the relative weakness of this source at the two wavelengths and it being close to the primary beam edge. MM3 and MM7 show a secondary peak at 878 μm, MM3b and MM7b, which were not detected previously. Both new sources are located westward of the main component and have peak fluxes of ~55% that of the main component. In addition, four new sources are detected at 878 μm, named MM10−MM13 following the convention of QZM11. All the dust continuum sources except MM11 and MM13 appear to be located in a diffuse arc-like filamentary structure, resembling a “spiral arm” that encloses IRS 1. (Hereafter we refer to this filament as “spiral-arm”.) This spiral arm was first reported by Kawabe et al. (1992) from CS 2−1 observations. At the center of the complex, MM1 and MM2 sources seem to be closer to each other and embedded in a dense environment. The other sources, located mainly E and S of IRS 1, are embedded in more diffuse medium and form a C-shaped structure that seems to arise from MM1-MM2.

The measured and derived physical parameters of the sources are listed in Table 1. Temperatures are assumed to be those used by QZM11: 245 K for MM1 and 40−58 K for the rest of the sample. For the new sources (MM10−13), we adopted a temperature of 40 K. (These sources are weaker and relatively far from IRS 1.) The diameter of the sources (FWHM of the dust emission) is relatively homogeneous with an average value of 3′′.0 ± 0′′.5 (≃8 × 10³ AU), slightly larger than the beam, and thus are poorly resolved. The densest object is MM1 (2.8 × 10⁷ cm^-3), but the most massive one is MM2 (≃37 M_⊙). The total mass is ≃160 M_⊙, where 50% is contained in the MM1-MM2-MM4 region (≃80 M_⊙) and 40% in the C-shaped filament SE to MM1 (≃60 M_⊙). In addition, the central sources appear to be denser than those in the C-shaped filament, with the exception of MM3.

Fig. 3 Spectrum toward the MM1 peak. Top panel: lower side band. Bottom panel: upper side band. A common temperature range from −2 to 14 K is displayed for better visualization. Additional temperature ranges are shown for the upper (−2 to −20 K) and lower (14 to 42 K) sidebands to display the entire spectrum. A dashed lines marks the 0 K level. For both sidebands, the frequencies of the main molecular transitions are noted by a dotted vertical line and the transition is specified.

3.2. Linearly polarized continuum emission

The polarized emission is broadly detected along the dust filamentary structure, with a polarized intensity between ≃0.017 and 0.059 Jy beam^-1. Two sets of polarization segments were computed. Firstly, the high−significance set (Fig. 1) is computed using 3.0σ_pol cutoff in the Stokes Q and U maps and 6σ_I cutoff in the Stokes I maps. Secondly, the low−significance set (Fig. 1) is computed using 2σ_pol and 3σ_I cutoffs, respectively. The agreement in the magnetic field direction is remarkable between both sets, thus the low cutoff values deliver realistic information (Sect. 4.2). The overall morphology of the magnetic field segment directions¹ is not uniform across the region, unlike the ordered directions detected in other regions (Girart et al. 2006, 2009), and seem to roughly follow the arc−like filament direction as traced by the continuum emission.

3.3. Molecular line emission

In this section we present the molecular line data toward NGC 7538 IRS 1 at 878 μm that can be summarized by two different behaviors: (i) four molecular transitions that have extended emission arising from the diffuse material, and (ii) many spatially unresolved hot-core lines only present toward the chemically differentiated MM1 source.

The four molecular transitions that exhibit extended emission are shown in Fig. 2, which presents the zero−, first−, and second-order moment maps (integrated intensity, velocity, and velocity dispersion maps, respectively) of the emission. The different transitions are ordered by increasing critical density²: ~10⁵ cm^-3 for C¹⁷O 3−2, a few 10⁶ cm^-3 for H¹³CO⁺ 4−3 and SO 8−7, and ~10⁷ cm^-3 for C³⁴S 7−6. These moment maps cover the velocity range between −66 and −60 km s^-1, where the emission is detected. The two most extended molecular transitions, C¹⁷O 3−2 and H¹³CO⁺ 4−3 (first and second columns, respectively) trace the dust emission with high fidelity. For both transitions, the first−order moment maps show an almost equivalent velocity pattern, which strongly suggests that the complex velocity structure is real. MM1−MM2 sources have a velocity around −64.5 km s^-1, while the rest of the sources seem to be closer to − 61 km s^-1. These features are not observed as well in the C¹⁷O because it does not trace MM1 and MM4 very well. The “spiral arm” possibly starting in the MM1−MM2 region and ending up in MM9, is traced very well by the C¹⁷O and H¹³CO⁺ emission. The projected velocity pattern shows an increase in velocity from MM1 to MM6, and then a slight decrease along the filament down to MM7. An interesting feature of this filament is that it appears to be a marginal velocity gradient across the filament, with the inner edge of the arc−like structure blueshifted with respect to the outer edge.

The SO 8−7 and C³⁴S 7−6 emission (third and fourth columns, respectively, of Fig. 2) appears to be less extended than that of the C¹⁷O and H¹³CO⁺ lines, only partially tracing the dusty arc−like filament. This can be due to their higher critical density, i.e., they trace the densest parts of the filament. Both are present around MM1 (including MM2, MM3, MM3b, and MM5) and MM7. The SO line is also detected toward MM6 and MM11. MM1 is particularly bright in these two molecular transition maps and shows a large velocity dispersion of 2−3 km s^-1, in contrast to the rest of the filament that hardly achieves a dispersion of 1.5 km s^-1. The velocity pattern of the SO and C³⁴S is compatible with the more extended, C¹⁷O and H¹³CO⁺ lines.

The set of molecules detected only toward MM1 continuum emission peak are all spatially unresolved. Figure 3 presents the spectrum toward the peak of MM1 with most of the lines labeled. This spectra is typical of an evolved hot core (see Sect. 5): it contains a few transitions of SO, SO₂, ³⁴SO₂, NS, HC₃N, H₂CO, and NH₂CHO, as well as many transitions of CH₃OH and CH₃OCHO. The upper energies of these transitions span an order of magnitude, ranging from ~50 K to ~500 K. Several unidentified lines were detected as well.

Fig. 4 Velocity structure of NGC 7538 IRS 1. Contours: observed dust continuum map, contours are 3σ to 21σ in steps of 6σ, where σ = 0.025 Jy beam-1. First column: order 1 moment maps, i.e. velocity structure. The common scale is shown at the top of the column. Second to fifth columns: channel maps with the velocity labeled at the top. Rows a and b: observed H13CO+ 4−3 a) and C17O 3−2 b) maps. The grayscale for the channel maps is shown in the righthand side of the figure. Rows c and d: synthetic maps generated with RATPACKS (Sect. 4.1.1) for a logarithmic spiral with radial expanding motions c) and with rotational motions d).

4. Analysis

4.1. The kinematics of the “spiral arm” around NGC 7538 IRS 1

The velocity pattern found toward the “spiral arm” of NGC 7538 IRS 1 is consistent among the tracers and shows very smooth variations (Fig. 2). This pattern is present towards the S and E of the UC H II region and includes MM3, MM3b, MM5, MM6, MM7, MM7b, MM10, and MM12. The total amount of mass in the “spiral arm” from dust continuum is ≃60 M_⊙.

4.1.1. Kinematic models

We used a set of simplified models with different geometries and different velocity structures. For the geometrical structure, we used two-dimensional Archimedean spirals r = a + bθ and logarithmic spirals r = ae^θ/b, both expressed in polar coordinates. A certain thickness was applied to the spiral to form a tubular three-dimensional structure.

For the kinematics, we took radial and rotational motions into account, both with speed following a potential law of the form v(r)_rad/rot = v_{rad / rot}(r_out)[r/r_out]^α. At any point in space, the radial velocity vector is tangential to the radial direction with respect to the source center, while the rotation velocity vector is simultaneously perpendicular to the rotation axis and the radial direction. Positive v(r_out) means expansion and counter-clockwise rotation for the cases of radial and rotational motions, respectively. However, changing the signs of the angle with respect to the LOS and of the velocity would produce the same map, and there is therefore an uncertainty on the direction of the gas flow. The parameter α can be used in radial motions to accelerate or decelerate the gas as a function of the radius. For rotating motions, α can be used to simulate rigid body rotation, constant speed rotation, and Keplerian rotation with α = 1, α = 0, and α = −1/2, respectively.

Finally, we developed a simple RAdiative Transfer Package for Adaptable Construction of Kinematical and Structural models (RATPACKS) to generate synthetic velocity cubes using any combination of geometric and kinematic input as a synthetic source. The synthetic source can be rotated around any axis allowing any orientation in three-dimensional space. The velocities are projected on the plane-of-the-sky (POS) according to the velocity pattern chosen and to the three-dimensional orientation given to the synthetic source. Then, a simple radiative transfer routine assuming optically thin emission is used to derive noiseless synthetic channel maps, and order 0, 1, and 2 moment maps.

4.1.2. Application to NGC 7538 IRS 1

We explored the parameter space to best reproduce the H¹³CO⁺ and C¹⁷O maps of the entire “spiral arm”. We adopted a depth of the spiral arm equal to the average observed width: 5′′ (64 mpc). Since we are only interested in the general kinematics, we assumed uniform gas density, constant molecular abundance, and optically thin emission. In Fig. 4 we present the models with the best-fit parameters, together with the H¹³CO⁺ and C¹⁷O data for comparison. Logarithmic spirals fit the data better than do Archimedean spirals and are the ones presented in Fig. 4. For the models presented, we used a =1.4 × 10⁴ AU and b = 42 rad. The parameters used for the kinematics are listed in Table 2. For both velocity cases, we set v_LSR = −63 km s^-1 and r_out = 1.5 × 10⁴ AU (~5′′.7), and ${\mathit{\theta }}_{\mathrm{LOS}}^{\mathrm{spir}}$ represents the angle between the plane containing the spiral and the LOS, where 0° means edge-on and 90° means face-on.

Row c in Fig. 4 shows a synthetic source with radial expanding motions, while row d shows counter-clockwise rotating motions. The radial expanding motions seem to reproduce the channel maps best in both molecules. This velocity pattern produces channel maps with extended emission. The best fit is achieved using α = 0 with a constant radial velocity of 9 km s^-1 and ${\mathit{\theta }}_{\mathrm{LOS}}^{\mathrm{spir}}\mathrm{=}\mathrm{70}$°, close to face-on. The rotational pattern produces channel maps with concentrated emission, not seen in the data. The best fit is achieved with a rigid-body rotation of 2 km s^-1 at r_out and ${\mathit{\theta }}_{\mathrm{LOS}}^{\mathrm{spir}}\mathrm{=}\mathrm{45}$°. The values of the velocities depend strongly on the orientation angles of the synthetic source, and therefore velocities should be taken as upper limits.

4.2. Statistical derivation of the magnetic field strength

Fig. 5 Angular dispersion function of the magnetic field segments detected toward the spiral arm (panels a) and b) and the central region (panels c) and d)). Top subpanels (a) and c)): dots represent the data with uncertainty bars, dashed line marks the zero value, dotted vertical line notes the beam size, dotted horizontal line shows the expected value for a randomic magnetic field, red line shows the best fit to the large-scale magnetic field (summation in Eq. (2)), and blue line shows the best fit to the data (Eq. (2)). Bottom subpanels (b) and d)): dots represent the correlated component of the best fit to the data, dashed line marks the zero value, dotted vertical line notes the beam size, red line shows the correlation due to the beam, and blue line shows the correlation due to the beam and the turbulent component of the magnetic field.

4.2.1. Formalism

Based on observational data, a widely used method of estimating the magnetic field strength of the POS component of the large-scale magnetic field is the Chandrasekhar-Fermi (hereafter CF; Chandrasekhar & Fermi 1953) equation $\frac{\mathit{\delta B}}{{\mathit{B}}_{\mathrm{0}}}\mathrm{\simeq }\frac{{\mathit{\sigma }}_{\mathit{v}}}{{\mathit{V}}_{\mathrm{A}}}\mathit{,}$(1)

where B₀ = |B₀| is the large-scale component of the magnetic field, δB the variation about B₀, ${\mathit{V}}_{\mathrm{A}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \rho }}$ is the Alfvén speed at density ρ, and σ_v the velocity dispersion of an appropriate spectral line. Recently, different statistical methods have been developed to avoid some of the CF method caveats (Hildebrand et al. 2009; Houde et al. 2009, 2011; Koch et al. 2010). These methods rely on studying the extended polarized emission in observational maps.

Houde et al. (2009) assume two statistically independent components of B, the large-scale magnetic field B₀(x), and the turbulent magnetic field B_t(x). Then, they derive the turbulent to large-scale magnetic field strength ratio from the angular dispersion function that accounts for the polarization angle differences as a function of the distance between the measured positions. The analysis is based on an analytical derivation for a turbulent cloud (see their Eq. (4)), including the effect of beam and LOS averaging. They further assume a stationary, homogeneous, and isotropic magnetic field strength, an isotropic but negligible turbulent polarized emission, and a magnetic field turbulent correlation length δ that is much smaller than the thickness of the cloud Δ′ (δ ≪ Δ′). Applying all these simplifications, the angular dispersion function can be written as $\begin{array}{ccc}\mathrm{1}\mathrm{-}\mathrm{⟨}\cos \left[\mathrm{\Delta \Phi }\mathrm{(}\hspace{0.17em}\mathit{l}{\hspace{0.17em}}^{\mathrm{)}}\right]\mathrm{⟩}\mathrm{\simeq }\frac{\mathrm{⟨}{\mathit{B}}_{\mathrm{t}}^{\mathrm{2}}\mathrm{⟩}}{\mathrm{⟨}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{⟩}}\hspace{0.17em}\frac{\mathrm{1}}{\mathit{N}}\hspace{0.17em}\left[\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{l}}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{(}{\mathit{\delta }}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{W}}^{\mathrm{2}}\mathrm{)}}\right]\mathrm{+}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}{\mathit{a}}_{\mathrm{2}\mathit{j}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{l}}^{\hspace{0.17em}\mathrm{2}\mathit{j}}\mathit{,}& & \end{array}$(2)where $\mathit{N}\mathrm{=}\left[\frac{\mathrm{(}{\mathit{\delta }}^{\mathrm{2}}\mathrm{+}\mathrm{2}\hspace{0.17em}{\mathit{W}}^{\mathrm{2}}\mathrm{)}{\mathrm{\Delta }}^{\mathrm{\prime }}}{\sqrt{\mathrm{2}\hspace{0.17em}\mathit{\pi }}\hspace{0.17em}{\mathit{\delta }}^{\mathrm{3}}}\right]$(3)

is the number of independent turbulent cells along the LOS, W the standard deviation (σ = FWHM/$\sqrt{\mathrm{8}\ln \mathrm{2}}$) of the Gaussian beam, and the summation is a Taylor expansion representing the large-scale magnetic field component that does not involve turbulence. The first term in the square brackets contains the integrated turbulent magnetic field contribution, while the exponential term represents the correlation by the combined effect of the beam (W) and the turbulent magnetic field (δ). The intercept of the fit to the data of the uncorrelated part at l = 0, f_NC(0), and the assumption of a cloud thickness Δ′ allow us to estimate the turbulent to large-scale magnetic field strength ratio as $\frac{\mathrm{⟨}{\mathit{B}}_{\mathrm{t}}^{\mathrm{2}}\mathrm{⟩}}{\mathrm{⟨}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{⟩}}\mathrm{=}\mathit{N}\hspace{0.17em}{\mathit{f}}_{\mathrm{NC}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}$(4)

Finally, identifying $\mathrm{⟨}{\mathit{B}}_{\mathrm{t}}^{\mathrm{2}}\mathrm{⟩}\mathrm{\equiv }\mathit{\delta }{\mathit{B}}^{\mathrm{2}}$, one can apply the CF equation (Eq. 1) to derive the large-scale component of the magnetic field as

$\mathrm{⟨}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\sqrt{\mathrm{4}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}\mathit{\rho }}\hspace{0.17em}{\mathit{\sigma }}_{\mathit{v}}\hspace{0.17em}{\left[\frac{\mathrm{⟨}{\mathit{B}}_{\mathrm{t}}^{\mathrm{2}}\mathrm{⟩}}{\mathrm{⟨}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{⟩}}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}$(5)We note that the magnetic field component labeled as “turbulent” describes, more generally, any contribution to the total magnetic field other than the uniform large-scale one. Therefore, when we refer in the next sections to the “turbulent” magnetic field, we are discussing the nonuniform magnetic field contribution.

4.2.2. Application to NGC 7538 IRS 1

The magnetic field segments in this complex region do not follow a defined homogeneous pattern along the observed field (see Fig. 1) such as the hour-glass shape reported and modeled in simpler sources (Girart et al. 2006, 2009; Frau et al. 2011; Padovani et al. 2012). No analytical models are available that can be compared to this complex source. Thus, to extract physical information, a statistical approach seems to be the best strategy. The “spiral arm” and the central sources seem to have different kinematics and different directional patterns of the segments, probably related to the YSO embedded in the central region (Sect. 4.1 and Fig. 1). Consequently, we analyzed the magnetic field for each of the regions independently. Figure 5 shows the angular dispersion function for both structures. Bins are equally spaced by 1′′. Data points represent the mean within the bin, with uncertainties that are smaller than the point size. We use the nonlinear least-squares Marquardt-Levenberg algorithm to fit Eq. (2) to the data. The best fit is shown in Fig. 5 and the parameters are listed in Table 3.

Panels (a) and (b) show the results for the spiral arm. The uncorrelated large-scale component is fitted with a j = 1 polynomial following Eq. (2). The correlated component dominates at small distances (~5′′−6′′ or ~64−77 mpc at 2.65 kpc). The turbulent magnetic field effect in correlating the segments is significantly more important than the beam effect. The correlation length is (33 ± 4 mpc at 2.65 kpc), almost three times larger than the beam correlation distance.

Panels (c) and (d) show the results for the central sources. The correlated component is only important in a distance half that of the spiral arm, and is mostly due to the effect of the beam. For larger distances, the data flattens to a value compatible with a random magnetic field ((1 − cos ( ≃52°) ≃0.384: Poidevin et al. 2010, see also Girart et al. 2013). Therefore, the summation in Eq. (2) (large-scale magnetic field) was dropped in our analysis, and only the correlated component was used (see Table 3). The best fit leads to a turbulent magnetic field correlation length of δ = 1′′.0 ± 0′′.6 (13 ± 8 mpc at 2.65 kpc), roughly equal to (the telescope beam correlation length, Eq. (2)).

To estimate the magnetic field strength, one has to assume either a certain cloud thickness Δ′ or a certain number of turbulent cells N (see Eq. (3)). Since the magnetic field that we are tracing mostly comes from a filamentary structure, an educated guess is to assume that the thickness is that of the filament width, 5′′ (64 mpc). Under this assumption, we find that the spiral arm contains one turbulent cell along the LOS, while the result for the central sources is poorly constrained to 6 ± 8 cells. In both cases, the local turbulent field is more important than the large-scale, ordered field. However, the central sources have a magnetic turbulence two times more important than the spiral arm. The strength of the field is ~2.5 mG, similar for both regions.

Fig. 6 Comparison of the dusty filament orientation to the magnetic field segments orientation toward NGC 7538 IRS 1. Pixels have been resized to the beam size to ensure statistical independence (see Fig. 1 for Nyquist sampling). Top panel: contours are 3, 9, and 27 times 0.02 Jy beam-1 continuum emission levels. Blue segments are derived as for Fig. 1. The red thick line is the axis of the filament (see Sect. 4.3). Red segments show the orientation of the filament corresponding to each pixel with polarization detection. Bottom panel: histogram of the angle difference between the polarization segment and the filament axis for each pixel. The correlated and uncorrelated distributions used for the χ2 test in Sect. 4.3 are shown as dotted and dashed lines, respectively.

4.3. Comparison of dust and magnetic field structures

The magnetic field segments seem to roughly follow the direction of the “spiral arm”. To examine this, we defined the dusty “spiral arm” axis as the line connecting the dusty sources (nodes) within the structure (shown in Fig. 6). Then, we obtained new maps with pixels of the size of the beam to ensure statistical independence. Each independent pixel with polarized emission was assigned a segment representing the local direction of the filament. This direction was defined as the line connecting the two nearest nodes (red vectors in Fig. 6). Finally, the difference between the magnetic field segment direction and the filament local direction was computed and binned. As shown in the histogram in Fig. 6, nine of the fourteen segments (~64%) have differences of less than 20°, and none have an angular difference larger than 50°. The number of independent measurement is relatively small, so we performed a χ² test to assess the statistical significance of the results. On the one hand, we compared the data to a flat distribution representing uncorrelated orientations. We found a <1% probability of obtaining higher χ² values with random data, so the null hypothesis of uncorrelated orientations was rejected. On the other hand, we compared the data to a simple distribution representing correlated orientations. We simplified this distribution to a linearly decreasing function that evolves from total correlation to none in half of the angular range covered (see Fig. 6). We found a 98% probability of obtaining a higher χ² value with random data, well above the standard 5% rejection threshold of the χ² test, so the null hypothesis of correlated orientations cannot be rejected. Based on this analysis, we conclude that the orientations of filament and polarization segments are correlated.

4.4. Energy state of the individual sources: the “mass balance”

We analyzed the main causes that are in interplay: gravity, magnetic field, thermal pressure, and internal dynamics. On the one hand, gravity has the effect of bringing mass together. On the other hand, the rest of the causes exert the opposite effect by either stopping the mass from accreting or dispersing it. We computed a series of meaningful parameters that relate these physical quantities. We also compared the relative strength of these causes in terms of the mass supported against gravitational collapse.

4.4.1. Formulae

We used the viral theorem to check whether the different cores are gravitationally bounded, and to estimate the maximum mass supported by the thermal and nonthermal motions. These take the different internal pressure components into account. In Sect. 5, we discuss the effect of the external pressure. The virial mass for the thermal component, M_T, is ${\mathit{M}}_{\mathrm{T}}\mathrm{=}\frac{\mathit{k}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{R}}{\mathit{G}}\mathit{,}$(6)where c_s, R and G are the sound speed, the core radius and the gravitational constant, respectively. The parameter k takes the specific density distribution of the core into account. We use k = 1, which is the value for a density profile ρ ∝ r^-2 (MacLaren et al. 1988). Similarly, for the nonthermal component the viral mass term, M_NT, is ${\mathit{M}}_{\mathrm{NT}}\mathrm{=}\frac{\mathit{k}{\mathit{\sigma }}_{\mathrm{NT}}^{\mathrm{2}}\mathit{R}}{\mathit{G}}\mathit{,}$(7)where σ_NT is the full three-dimensional velocity dispersion of the gas due to nonthermal motions.

The support of the magnetic field can be included as an additional component in the virial mass. Thus, the mass for a critical mass-to-flux ratio is given by Nakano & Nakamura (1978): ${\mathit{M}}_{\mathrm{mag}}\mathrm{=}\frac{\mathit{\pi }{\mathit{R}}^{\mathrm{2}}\mathit{B}}{\sqrt{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}\mathit{G}}}\mathrm{·}$(8)where B is the field strength.

The nonthermal kinetic energy can be compared to thermal kinetic energy by the square of the turbulent Mach number ${\mathrm{ℳ}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}{\left(\frac{{\mathit{\sigma }}_{\mathrm{NT}}}{{\mathit{c}}_{\mathrm{s}}}\right)}^{\mathrm{2}}\mathit{,}$(9)where ℳ_s> 1 means that nonthermal motions are supersonic, hence more dynamically important than thermal motions. The thermal energy is compared to magnetic field energy by the plasma β_T${\mathit{\beta }}_{\mathrm{T}}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{therm}}}{{\mathit{P}}_{\mathrm{mag}}}\mathrm{=}\mathrm{2}{\left(\frac{{\mathit{c}}_{\mathrm{s}}}{{\mathit{v}}_{\mathrm{A}}}\right)}^{\mathrm{2}}\mathit{,}$(10)where ${\mathit{v}}_{\mathrm{A}}\mathrm{=}{\mathit{B}}_{\mathrm{3}\mathrm{D}}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi }\mathit{\rho }}$ is the Alfvén speed. Similarly, magnetic fields compare to nonthermal motions by ${\mathit{\beta }}_{\mathrm{NT}}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{NT}}}{{\mathit{P}}_{\mathrm{mag}}}\mathrm{=}\mathrm{2}{\left(\frac{{\mathit{\sigma }}_{\mathrm{NT}}^{\mathrm{mol}}}{{\mathit{v}}_{\mathrm{A}}}\right)}^{\mathrm{2}}\mathit{,}$(11)where β_T < 1 or β_NT < 1 implies that magnetic pressure overcomes thermal or kinetic pressure, respectively.

4.4.2. Energy ratios in NGC 7538 IRS 1 sources

The parameters described in the previous section that measure the energy balance among forces are listed in Table 4. The sound speed ranges from 0.49 km s^-1 to 0.59 km s^-1 except for MM1 that hosts an O7.5 star and is significantly warmer. Nonthermal velocity dispersion σ_NT range from 0.23 km s^-1 to 1.27 kms. The mean value is 0.8 ± 0.3 km s^-1, ~60% greater that the typical sound speed at 40 K. As a result, 70% of the sources show supersonic gas motions ℳ_s > 1.

We use for each source the magnetic field strength derived in Sect. 4.2.2: 2.64 mG and 2.3 mG for the spiral arm and the central sources, respectively. For the isolated cores we used the average of 2.5 mG. This assumption implies that, within each region, the variation in the derived magnetic quantities depends on clump properties: ${\mathit{v}}_{\mathrm{A}}\mathrm{\propto }{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}^{-1}$ and M_mag ∝ r², where r is the clump radius derived as half the diameter from Table 1 (see discussion in Sect. 5.1). The Alfvén speed ranges from 0.96 km s^-1 to 5.3 km s^-1, with mean value of 3.0 ± 1.1 km s^-1. All sources but MM1 have v_A>c_s, so that magnetic pressure dominates locally over thermal pressure (β < 1). Nonthermal kinetic energy is comparable to the magnetic energy in four sources: MM1, MM2, MM4, and MM7. For the rest of the sample, magnetic pressure locally dominates nonthermal pressure.

4.4.3. “Mass balance” in NGC 7538 IRS 1 sources

A similar analysis can be performed in terms of the maximum mass supported by each force, listed in Table 4. The mass balance accounts for all the available information at once. MM1 hosts an embedded O7.5 star whose mass is taken into account (30 M_⊙: Pestalozzi et al. 2004). In Fig. 7, we compare collapse forces versus support forces to derive the individual “mass balance”. This analysis shows clearly two groups of sources in terms of stability, well correlated with their location in either (i) the central structure or (ii) the spiral arm and isolated sources.

The sources in the central structure seem to have more mass than what can be supported: MM1 and MM2 have masses that are higher than the combined maximum supported mass, while in MM4 masses are comparable. In contrast, all sources located either at the filament or in isolation have virial masses higher than measured masses. In general, the main agent against gravity for sources at the filament is the internal dynamics, while it is the magnetic field for the isolated sources. See Sect. 5.4 for a discussion on the implications.

5. Discussion

NGC 7538 IRS 1−3 as a whole presents a complex and rich structure, velocity field, and magnetic field. In addition, a number of individual cores can be identified. In this section, we first discuss the global properties, and then proceed to analyze the state of the individual cores.

Fig. 7 “Mass balance” analysis. Comparison of the measured mass to the maximum supported mass by different forces. Cores are ordered according to their location in the central massive structure, in the filament, or isolated. Black: measured mass from continuum maps. Gray: mass of the star embedded in the MM1 clump. Red: mass supported by magnetic fields assuming a uniform value across the source. Light and dark blue: mass supported by virialized gas motions due to internal dynamics and thermal dispersion, respectively.

5.1. Magnetic field properties

NGC 7538 IRS 1−3 contains two different regions in terms of magnetic field properties as shown in Sect. 4.2. These differences are important in two different but related aspects: the relative importance of unordered magnetic fields with respect to ordered ones, and the relative dynamical importance of magnetic fields in the overall picture.

For the central sources, the magnetic field segments are only correlated at slightly larger distances than the correlation distance of the beam. The turbulence of the magnetic field has a mild effect, but the data does not allow us to accurately constrain the contribution. In any case, the transitions to values for the angular dispersion function compatible to a random field happens at a very short distance (~3′′, ~39 mpc). The number of turbulent cells along the LOS seems to be large, which is a clue that the field is severely distorted. This suggests that not ordered magnetic fields are more important than ordered magnetic fields. In fact, the ratio of energies $\mathrm{⟨}{\mathit{B}}_{\mathrm{t}}^{\mathrm{2}}\mathrm{⟩}\mathit{/}\mathrm{⟨}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{⟩}$ is ${{\mathrm{~}}_{}}^{\mathit{>}}$3. The size of the sources is larger than the correlation distance in all cases, and thus, this region needs to be analyzed source-by-source. The study of the energetic state of the individual sources suggests that nonthermal motions are dominant. In addition, the magnetic field tends to a random configuration on relatively short scales. These facts suggest that magnetic fields are not important in the overall picture. In summary, the relevance of the magnetic field is low in this region, and the field is mostly not ordered.

The spiral arm shows different magnetic field properties with respect to the central structure. The correlation length due to the not ordered magnetic field is accurately determined since it is about three times longer than the beam correlation length. The distance between consecutive embedded sources within the filament of 2′′ (26 mpc) in average, smaller than the magnetic correlation length of 2′′.6 (33 mpc). This implies that the field properties among cores are not independent. Consequently, the filament analysis must take the whole complex into account. According to the statistical analysis, the spiral arm contains only one turbulent cell along the LOS; as a result, it has relatively little turbulence. This agrees with the well-aligned segments observed and with the sources having dependent magnetic fields. Finally, the field shows a ∝l² trend on large scales that suggests a smoothly varying field (Hildebrand et al. 2009). Moreover, the analysis in Sect. 4.3 shows that the magnetic field direction is correlated to the dust morphology along the entire 11′′-long filament. All these suggest that the magnetic field has a strong internal coherence within the filament and that it is somehow tied to the dust structure.

Fig. 8 Left panel: overlayed contours for dust continuum (gray, this work), and blue− and redshifted 13CO 2−1 outflow (blue and red, QZM11). Crosses mark the positions for dusty cores and red dots for IR sources. Right panel: schematic 3D cartoon of the proposed scenario for the NGC 7538 IRS 1−3 complex (Sect. 5.3). Yellow sphere represents MM1, blue and red cones represent the blue and red outflow lobes, and the gray structure represents the spiral arm. Faded colors represent the structures behind MM1 in the LOS direction.

5.2. Energetics of the spiral arm

The analysis of the kinematics in Sect. 4.1 suggests that the spiral arm around NGC 7538 IRS 1 is expanding, although we cannot discard a certain contribution from rotation. This result has already been suggested based on CS observations by Kawabe et al. (1992).

For the expansion to happen, the filament must be gravitationally unbound to the total mass around IRS 1. The combined mass of MM1, MM2, MM4, plus the star embedded in MM1, is 110 M_⊙. The distance of the filament with respect to IRS 1 is ~13′′ (0.17 pc at 2.65 kpc). For these quantities, the virialized rotation velocity of the filament is 1.7 km s^-1, and the escape velocity is 2.4 km s^-1. The velocity difference of the filament with respect to the central region is in the 1.5−4.5 km s^-1 range when projected in the POS, or 9 km s^-1 according to the best-fitting kinematic model. Therefore, the spiral arm appears to be gravitationally unbound with respect to the massive MM1, MM2, and MM4 cores. We note that the measured mass is a lower limit due to the filtering effect of the SMA. However, a mass of 1.6 × 10³M_⊙ is required within the central ~10′′ to gravitationally bind the gas moving at 9 km s^-1. The single-dish measured mass of the entire 1′ clump is 3.7 × 10³M_⊙ (Momose et al. 2001, corrected for the different κ_dust and distance used). It seems unlikely that 40% of the total mass is accumulated within the central 10′′, and so it is unlikely that the spiral−arm is gravitationally bound.

When we focus on the filament, the total mass combining MM3, MM6, MM7, MM7b, MM10, and MM12 is 45.8 M_⊙. The total gravitational energy of the filament using a radius of 13′′ is thus E_grav = 2.6 × 10⁴⁵ erg. In addition, we can derive relevant dynamical parameters from the kinematic model. We assumed for the calculations that the expansion velocity is constant. Then, we considered two different scenarios: (i) a conservative one where the expansion velocity is the maximum velocity projected in the POS, 4.5 kms; and (ii) the one drawn by the kinematic model that takes projection effects into account and results in a faster velocity, 9 km s^-1. The conservative scenario delivers an age of t_dyn ≃ 3.6 × 10⁴ yr and kinetic energy of E_kin ≃ 9.2 × 10⁴⁵ erg, while the deprojected scenario delivers t_dyn ≃ 1.8 × 10⁴ yr and E_kin ≃ 3.7 × 10⁴⁶ erg. The gravitational well is a factor of 3.5−14 smaller, thus confirming the plausibility of an expanding filament.

To complete the picture, a crude estimate of the magnetic field energy can be done by multiplying the volume of mass permeated by the field times the overall magnetic pressure (${\mathit{P}}_{\mathrm{B}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathrm{8}\mathit{\pi }\mathrm{\right)}$). An approximate area of 150′′² and depth of 5′′ delivers a magnetic energy of E_mag = 1.3 × 10⁴⁵ erg, which is negligible for the overall filament dynamics when compared to kinetic energy. The correlation of the field morphology to the dust morphology can now be tentatively explained: the three to nine times more energetic expansion motions seem to push away both matter and magnetic field, shaping them in a similar morphology as shown in Fig. 6.

5.3. Formation of the spiral arm through stellar feedback

The large-scale expansion motions proposed require a powerful driving source. A plausible origin is the IRS 1 feedback, in particular, from the powerful NW-SW molecular outflow powered by IRS 1 (QZM11). We show in Fig. 8 a comparison of the ¹³CO outflow from MM1 to our continuum map. Projected on the POS, it seems that the outflow is perpendicular to the spiral arm, suggesting that it could be formed by swept material.

Quantitatively, QZM11 derive E_outflow = 4.9 × 10⁴⁶ erg, which is a factor of 1.2−3.7 larger than the combined kinetic, magnetic, and gravitational energy of the filament. In addition, energy losses are expected in the form of turbulence in this complex scenario. The average nonthermal velocity dispersion is Δv^NT = 1.8 km s^-1, which implies an extra 1.5 × 10⁴⁵ erg that the outflow can provide. Therefore, it is feasible that the molecular outflow is the energy source responsible for the expansion of the spiral arm. Furthermore, the outflow momentum for one lobe, 2.3 × 10²M_⊙ km s^-1, is comparable to the momentum of the spiral arm, ≃2 × 10²M_⊙ km s^-1, suggesting that the outflow is the unique cause of the filament being set into motion. Finally, the outflow dynamical timescale of 2 × 10⁴ yr falls within the filament age range if expanding, suggesting that the two structures have contemporary births.

We can take the inclination effects into account for better determination of the outflow effect. Based on our kinematic model, we found that the plane containing the logarithmic spiral lies almost parallel to the POS, tilted by ~20° (70° with respect to the LOS, see Table 2). This information may help for understanding the three-dimensional orientation of the system. One possibility is that outflow and spiral arm flow are parallel, and thus, the gas forming the spiral arm would be pushed by the outflow end. The opposite possibility is that the spiral arm is perpendicular to the outflow, implying that the material would be blown away from the cavity (~80° wide: Kraus et al. 2006) and driven to the equatorial plane of the system. The aligned scenario implies an outflow lying almost on the POS, with an increase on de-projected outflow velocity by a factor of ~3 and on E_outflow by a factor of ~9, while the perpendicular scenario renders an outflow inclined by 20° with respect to the LOS, and a mild increase of ~6% in outflow velocity and ~13% in E_outflow. A possibly precessing outflow axis may increase the uncertainty in our analysis (Kraus et al. 2006). Also, we might consider configurations between the parallel and perpendicular scenarios. The perpendicular case would imply that both structures have opposite velocity directions when projected on the POS. Therefore, the fact that the spiral arm is mostly redshifted toward the same direction as the large scale outflow seems to favor the aligned scenario. There are other cases in the literature of sweeping up the ambient material as a snowplow and accumulating it into a shell (Anglada et al. 1995; Girart et al. 2005).

Consequently, based on the morphology considerations in this section and on the energy considerations in Sect. 5.2, we speculate that the dusty spiral arm is created by the accumulation of matter due to the IRS 1 outflow feedback.

5.4. Gravitational collapse of the individual cores: a cluster in the making

Three active, bright IR sources indicate that star formation is ongoing in the NGC 7538 IRS 1−3 cluster. Moreover, through high-resolution IR interferometry, Kraus et al. (2006) find 18 new faint stars and a NW-oriented, fan-shaped outflow arising from MM1. Interestingly, the positions of the stars are correlated well with the outflow, which they propose is precessing and triggering star formation. In this environment of high interaction, we target the cold dust to study the possible evolution of the mass reservoir in the cluster.

The three cores in the central dusty structure seem to be gravitationally dominated against the support forces (Sect. 4.4, MM4 is only at the limit). The MM1 core has already formed a still accreting protostar, which has gathered about two thirds of the locally available total mass (star plus dust core system). If the same star-to-core mass ratio applies, MM2 will form massive stars of ~24 M_⊙ and MM4 of ~17 M_⊙.

The situation is less clear for the dust cores either located in the spiral arm or else isolated. The “mass balance” analysis shows that the total mass is insufficient to gravitationally bound the cores. Figure 9 shows an extension of the traditional virial parameter analysis with inclusion of the magnetically supported mass. This parameter is usually fit by a function of the form ${\mathit{M}}_{\mathrm{obs}}^{\mathit{a}}$ that delivers a typical a = −0.68 ± 0.06, in agreement with previously reported trends (see, e.g., Bertoldi & McKee 1992, and derived works). The data can also be fitted by a function of the form b[M_vir/M_obs]^a, where the proportionality constant b carries a physical meaning, varying from 2.06 for self-graviting clumps to 2.9 for pressure-confined clumps (Bertoldi & McKee 1992). The best fit to the cores delivers b = 2.26 ± 0.12, in good agreement with the prediction of b = 2.12 for magnetized critical cores. To examine the range of applicability of the model, we show the expected results in Fig. 9 using the masses of the clumps, the theoretical b = 2.12, and the previously derived a. The prediction for magnetized critical cores is in good agreement with the magnetically dominated low-mass cores (see also Sect. 4.4 and Fig. 7). In contrast, the model prediction is less precise for the more massive cores, where the magnetic field is less important and are less likely to be magnetically critical.

The analyses performed assume cores in isolation and ignore the effects of the highly dynamical environment. In other words, the expansion powered by the outflow may help pile up material as the filament is expanding through the ISM, and more importantly, it is creating a high external pressure along the filament. This pressure can help gravitation to overcome magnetic and turbulent energy. We can estimate wether an external perturbation can have a significant impact on core evolution by using the typical crossing-time t_cross = R_core/c_s. On average, the filament cores have R_core ≃ 3.6 × 10³ AU and c_s ≃ 0.53 km s^-1, resulting in a typical t_cross ≃ 3 × 10⁴ yr. This value is comparable to the t_dyn estimates for the filament and outflow (Sect. 5.2). Therefore, a perturbation constantly acting for t_dyn, such as stellar winds, can influence the evolution of the cores in the filament. Since all measures suggest that the outflow is behind the formation of the entire filament (Sect. 5.3), comparable timescales at a core level make it reasonable for the external pressure to also influence the cores.

We speculate that the external pressure from the winds acting for t_dyn could trigger the collapse of the cores in the filament, leading to the formation of a group of low-mass stars. This triggered star formation SE of MM1 is supported by the mirrored star formation toward the NW, in a more evolved stage of evolution. These stars are older than the MM1 outflow and could have been formed through feedback from the older, more evolved IRS 2 star. This star is associated to a well-studied H II region and powers a “stellar wind bow shock” (Bloomer et al. 1998). Consequently, triggered star formation in the NGC 7538 IRS 1−3 complex could be an episodic process following the evolution of the most massive stars. Such a scenario would generate a small cluster with two stellar groups: (i) a few central high-mass stars; surrounded by (ii) a wealth of low-mass stars formed through feedback from the former group.

Fig. 9 Ratio of supported mass to observed mass as a function of the observed mass for the NGC 7538 IRS 1−3 dense cores. MM1 is represented by a red star and the other cores by black circles. The solid line shows the best fit to the starless cores while blue crosses show the expected values for critical, magnetized clumps (see Sect. 5.4). The dotted line shows the limiting value for effective support against gravitational collapse.

6. Conclusions

We have carried out a molecular, dust, and polarimetric study of the NGC 7538 IRS 1−3 star-forming cluster. We used SMA high angular resolution observations at 880 μm with the compact configuration. Here, we summarize the main results.

• 1.

  We detect up to 14 dust cores in continuum emission, six of them newly discovered, spanning one order of magnitude in mass (from 3.5 M_⊙ to 37 M_⊙). The brightest core is MM1, associated with IRS 1. IRS 2 and IRS 3 show no continuum counterpart. The dust cores are connected by diffuse gas, and are arranged in two larger scale structures: a central bar containing MM1, MM2, and MM4; and a filamentary spiral arm containing at least six cores. The total dust mass is ≃160 M_⊙, almost equally split between the two large scale structures.

• 2.

  We detect C¹⁷O 3−2 and H¹³CO⁺ 4−3 large-scale emission sharply tracing the two main large-scale structures, unveiling a velocity gradient along the spiral arm. We developed a code to generate synthetic velocity cubes, RATPACKS, and reproduced the velocity gradient through a model of a spiral expanding at 9 km s^-1 with respect to the central MM1.

• 3.

  We broadly detect polarized emission in the compact cores and in the diffuse extended structures. Based on a statistical analysis, we derived a magnetic field strength of ≃2.5 mG. The orientation of the magnetic field segments is significantly homogeneous along the spiral arm, and it is correlated at an 80% confidence level to the direction of the dust main axis. This suggests that dust and magnetic field are tightly connected.

• 4.

  The spiral arm is gravitationally unbound with respect to the central bar. The gravitational and magnetic field energies combined are a factor of 2.3−9.5 smaller than the kinetic energy. Therefore, it is likely that the dominant expansion is shaping dust and the magnetic field into a similar morphology.

• 5.

  The total energy, linear momentum, and dynamic age (≃4.2 × 10⁴⁶ erg, ≃4 × 10²M_⊙ km s^-1, and ≃1.8 × 10⁴ yr) of the spiral arm are compatible with the values of the MM1 outflow by QZM11 when deprojected. Both spiral arm and outflow are redshifted, hence likely to flow in parallel. Consequently, it seems plausible that the dominant kinetic energy of the spiral arm has its origin in the MM1 outflow, which may be causing its formation in a snowplow fashion, in agreement with our expansion model.

• 6.

  We developed the “mass balance” analysis that compares collapse vs. support forces, accounting for all the available information on the energetics on core scales. On the one hand, the cores in the central bar seem to be gravitationally unstable and prone to forming massive stars. On the other hand, the combined support forces seem to dominate the cores located in the spiral arm or else isolated, with non-thermal motions and magnetic fields being the main agents of support, respectively. However, the dynamically important external pressure from the outflow could trigger the gravitational collapse and lead to the formation of low-mass stars as reported toward NW to MM1 (Kraus et al. 2006).

• 7.

  We speculate that the NGC 7538 IRS 1 region is forming a small cluster with a few central high-mass stars, surrounded by a number of low-mass stars formed through protostellar feedback.

Acknowledgments

We thank all members of the SMA staff that made these observations possible. This research made use of NASA’s Astrophysics Data System Bibliographic Services (http://adsabs.harvard.edu/), the SIMBAD database, operated at the CDS, Strasbourg, France (http://simbad.u-strasbg.fr/simbad/), and the Splatalogue database for astronomical spectroscopy (http://www.splatalogue.net). We thank the anonymous referee for the useful comments. PF is supported by the Spanish CONSOLIDER project CSD2009-00038. PF and JMG are supported by the Spanish MINECO AYA2011-30228-C03-02 and Catalan AGAUR 2009SGR1172 grants.

References

All Tables

All Figures

Fig. 1 Contour map of the SMA dust continuum emission map towards NGC 7538 IRS 1 at 878 μm overlapped with the gray image of the polarized intensity. The blue and red segments show the magnetic field direction obtained from the polarization data with cutoff at rms of 2 and 3σ, respectively (see Sect. 3.2). Solid contours go from 3σ to 33σ in steps of 3σ, where σ = 0.017 Jy beam-1 is the rms noise of dust continuum. The dashed contours around IRS 1 (MM1) go from 50σ to 300σ in steps of 50σ. The scale of the polarized intensity image is shown in the righthand side of the figure (the units are Jy beam-1). The black and white crosses show the peak positions of the dust continuum sources (Table 1). The name of the dust continuum sources are also shown. Red circles mark the positions of IRS 1−3. The solid circle shows the FWHM of the SMA primary beam at the observed frequency. The physical scale of the map and synthesized beam are shown in the bottom left corner of the panel.

In the text

Fig. 2 Moment maps of molecular transitions with extended emission. Rows: moments 0, 1, and 2, in descending order, labeled on the lefthand side of the figure. Columns: from left to right, C17O 3−2, H13CO+ 4−3, SO 8−7, and C34S 7−6, respectively, labeled on the top of the figure. Color map: moment maps. Common scale is shown on the righthand side of the figure. Contours: 3σ, 9σ, and 27σ continuum emission levels, where σ = 0.017 Jy beam-1.

In the text

Fig. 3 Spectrum toward the MM1 peak. Top panel: lower side band. Bottom panel: upper side band. A common temperature range from −2 to 14 K is displayed for better visualization. Additional temperature ranges are shown for the upper (−2 to −20 K) and lower (14 to 42 K) sidebands to display the entire spectrum. A dashed lines marks the 0 K level. For both sidebands, the frequencies of the main molecular transitions are noted by a dotted vertical line and the transition is specified.

In the text

Fig. 4 Velocity structure of NGC 7538 IRS 1. Contours: observed dust continuum map, contours are 3σ to 21σ in steps of 6σ, where σ = 0.025 Jy beam-1. First column: order 1 moment maps, i.e. velocity structure. The common scale is shown at the top of the column. Second to fifth columns: channel maps with the velocity labeled at the top. Rows a and b: observed H13CO+ 4−3 a) and C17O 3−2 b) maps. The grayscale for the channel maps is shown in the righthand side of the figure. Rows c and d: synthetic maps generated with RATPACKS (Sect. 4.1.1) for a logarithmic spiral with radial expanding motions c) and with rotational motions d).

In the text

Fig. 5 Angular dispersion function of the magnetic field segments detected toward the spiral arm (panels a) and b) and the central region (panels c) and d)). Top subpanels (a) and c)): dots represent the data with uncertainty bars, dashed line marks the zero value, dotted vertical line notes the beam size, dotted horizontal line shows the expected value for a randomic magnetic field, red line shows the best fit to the large-scale magnetic field (summation in Eq. (2)), and blue line shows the best fit to the data (Eq. (2)). Bottom subpanels (b) and d)): dots represent the correlated component of the best fit to the data, dashed line marks the zero value, dotted vertical line notes the beam size, red line shows the correlation due to the beam, and blue line shows the correlation due to the beam and the turbulent component of the magnetic field.

In the text

Fig. 6 Comparison of the dusty filament orientation to the magnetic field segments orientation toward NGC 7538 IRS 1. Pixels have been resized to the beam size to ensure statistical independence (see Fig. 1 for Nyquist sampling). Top panel: contours are 3, 9, and 27 times 0.02 Jy beam-1 continuum emission levels. Blue segments are derived as for Fig. 1. The red thick line is the axis of the filament (see Sect. 4.3). Red segments show the orientation of the filament corresponding to each pixel with polarization detection. Bottom panel: histogram of the angle difference between the polarization segment and the filament axis for each pixel. The correlated and uncorrelated distributions used for the χ2 test in Sect. 4.3 are shown as dotted and dashed lines, respectively.

In the text

Fig. 7 “Mass balance” analysis. Comparison of the measured mass to the maximum supported mass by different forces. Cores are ordered according to their location in the central massive structure, in the filament, or isolated. Black: measured mass from continuum maps. Gray: mass of the star embedded in the MM1 clump. Red: mass supported by magnetic fields assuming a uniform value across the source. Light and dark blue: mass supported by virialized gas motions due to internal dynamics and thermal dispersion, respectively.

In the text

Fig. 8 Left panel: overlayed contours for dust continuum (gray, this work), and blue− and redshifted 13CO 2−1 outflow (blue and red, QZM11). Crosses mark the positions for dusty cores and red dots for IR sources. Right panel: schematic 3D cartoon of the proposed scenario for the NGC 7538 IRS 1−3 complex (Sect. 5.3). Yellow sphere represents MM1, blue and red cones represent the blue and red outflow lobes, and the gray structure represents the spiral arm. Faded colors represent the structures behind MM1 in the LOS direction.

In the text

Fig. 9 Ratio of supported mass to observed mass as a function of the observed mass for the NGC 7538 IRS 1−3 dense cores. MM1 is represented by a red star and the other cores by black circles. The solid line shows the best fit to the starless cores while blue crosses show the expected values for critical, magnetized clumps (see Sect. 5.4). The dotted line shows the limiting value for effective support against gravitational collapse.

In the text