Magnetized converging flows towards the hot core in the intermediate/high-mass star-forming region NGC 6334 V

TL;DR

This study uses SMA observations to reveal converging magnetic field patterns and filamentary flows towards a hot core in NGC 6334 V, illustrating gravity-driven gas infall in high-mass star formation.

Contribution

It provides the first detailed observational evidence of magnetic field convergence and filamentary inflows in a high-mass star-forming region, supported by MHD simulation comparisons.

Findings

Magnetic field shows bimodal converging pattern towards hot core.

Two filamentary structures at different velocities converge to the hot core.

Gravity appears to dominate over magnetic field strength in gas infall.

Abstract

We present Submillimeter Array (SMA) observations at 345 GHz towards the intermediate/high-mass cluster-forming region NGC 6334 V. From the dust emission we spatially resolve three dense condensations, the brightest one presenting the typical chemistry of a hot core. The magnetic field (derived from the dust polarized emission) shows a bimodal converging pattern towards the hot core. The molecular emission traces two filamentary structures at two different velocities, separated by 2 km/s, converging to the hot core and following the magnetic field distribution. We compare the velocity field and the magnetic field derived from the SMA observations with MHD simulations of star-forming regions dominated by gravity. This comparison allows us to show how the gas falls in from the larger-scale extended dense core (~0.1 pc) of NGC 6334 V towards the higher-density hot core region (~0.02 pc) through two distinctive converging flows dragging the magnetic field, whose strength seems to have been overcome by gravity.

Tables

Table 1. Table 1 SMA observations summary

Array configuration	subcompact	compact	extended
Observing date	2011 July 2	2011 June 18	2011 July 21
${\tau }_{225\text{text}GHz}$	0.09	0.06–0.1	0.05
Number of antennas	7	7	8
Flux calibrator	Callisto	Uranus	Callisto
Passband calibrator	3C 279/3C 454.3	3C 279/3C 454.3	3C 279/3C 454.3
Gain calibrator	1733-130	1733-130	1733-130
Polarization calibrator	3C 279/3C 454.3	3C 279/3C 454.3	3C 279/3C 454.3
uv range ($k\lambda$)	6–42	8–80	24–206

Table 2. Table 2 Parameters of the sources detected with the SMA at 870  μ 𝜇 \mu m dust continuum emission.

	Positiona		Deconv.ang.sizea	P.A.a	$I_{\nu }^{\mathrm{peak}}$ a	$S_{\nu }$ a,b	Massc
Source	$\alpha \left(\mathrm{J2000}\right)$	$\delta \left(\mathrm{J2000}\right)$	${(}^{\prime \prime }{\times }^{\prime \prime })$	(∘)	(Jy Beam-1)	(Jy)	(M⊙)
HC	17:19:57.364	$-35$:57:53.47	$2.9\pm 0.6\times 1.66\pm 0.6$	$-89\pm 41$	$1.27\pm 0.16$	$3.6\pm 0.7$	4-9
N core	17:19:57.408	$-35$:57:48.55	$6.2\pm 2.3\times 2.62\pm 1.2$	$85\pm 18$	$0.50\pm 0.02$	$3.3\pm 0.7$	16
SE core	17:19:57.792	$-35$:57:53.50	$2.4\pm 0.4\times 1.78\pm 0.3$	$45\pm 20$	$0.80\pm 0.04$	$1.9\pm 0.4$	9

Table 3. Table 3 Molecular lines

Molecular line	Transition	Frequency (GHz)	E${}_{\text{text}U}$ (K)
Hot core tracers
SO2	21(2,20)$-$21(1,21)	332.09143	220
SO2	16(4,12)$-$16(3,13)	346.52388	164
SO2	19(1,19)$-$18(0,18)	346.65217	168
34SO2	16(4,12)$-$16(3,13)	332.83622	163
34SO2	19(1,19)$-$18(0,18)	344.58104	168
34SO2	11(4,8)$-$11(3,9)	344.99816	99
34SO2	8(4,4)$-$8(3,5)	345.16866	71
34SO	7(8)$-$6(7)	333.90098	80
CH3OH	18(2,16)$-$17(3,14)	344.10913	419
CH3OH	19(1,19)$-$18(2,16)++	344.44390	451
CH3OCHO	31(0,31)$-$30(1,30)	333.44902	260
CH3OCH3	19(1,19)$-$18(0,18)AA	344.35806	167
CH3OCH3	11(3,9)$-$10(2,8)EE	344.51538	73
HC3N	J=38$-$37	345.60901	323
Dense core tracers
SO	8(8)$-$7(7)	344.31061	87
SO2	4(3,1)$-$3(2,2)	332.50524	31
SO2	13(2,12)$-$12(1,11)	345.33854	93
CH3OH	7(1,7)$-$6(1,6)	335.58200	79
HC15N	(4$-$3)	344.20011	41
NS	J=15/2$-$13/2	346.22116	71
H13CO+	(4$-$3)	346.99835	41
Outflow tracers
CO	(3$-$2)	345.7959	33
SiO	(8$-$7)	347.33058	75

Table 4. Table 4 Angular dispersion function fit parameters

Parameters	SMA observations	Simulations
${\delta }^a$ (mpc)	17, 25	17, 25
$f_{\mathrm{NC}}^b$	0.29, 0.52	0.13, 0.28
${\Delta }^{\prime }$ (pc)	0.07	0.07
$n$ (cm-3)	$6.1\times {10}^5$	$6.1\times {10}^5$
${\sigma }_v$ (km s-1)	0.6	0.6
$⟨B_t^2⟩/{⟨B_0^2⟩}^b$	$\simeq 0.67,0.68$	$\simeq 0.31,0.37$
${⟨B⟩}_{pos}$ (mG)	$\simeq 0.7$	$\simeq 1.1$

Equations

$$M=R\frac{d^2{S}_{\nu }}{B_{\nu }(T_d){\kappa }_{\nu }},$$

$$\frac{\delta B}{B}\simeq \frac{{\sigma }_v}{V_A},$$

$$1-⟨\cos [\Delta \Phi (l)]⟩\simeq \frac{⟨B_t^2⟩}{⟨B_0^2⟩}\frac{1}{N}\left[1-{\mathrm{e}}^{-l^2/2({\delta }^2+2W^2)}\right]+\sum _{j=1}^{\infty }{a}_{2j}^{\prime }{l}^{2j},$$

$$N=\left[\frac{({\delta }^2+2W^2){\Delta }^{\prime }}{\sqrt{2\pi }{\delta }^3}\right]$$

$$\frac{⟨B_t^2⟩}{⟨B_0^2⟩}=N{f}_{\mathrm{NC}}(0).$$

Full text

1. Introduction

Previous observations have led to the consensus that the magnetic (B) field strength is non-negligible during the formation of molecular gas clouds and dense gas filaments, for low- and high-mass star- and cluster-forming regions (e.g., Matthews et al. 2009; Crutcher et al. 2010; Koch et al. 2014; Fissel et al. 2016; Soler et al. 2016; see Li et al. 2014 for a review). However, the role of the magnetic field in the formation and evolution of dense molecular cores (which have sizes of $\sim 0.1$ pc) is still a matter of debate (e.g., Vázquez-Semadeni et al., 2011; Bertram et al., 2012; Crutcher, 2012; Li et al., 2014; Körtgen & Banerjee, 2015). Its role can be studied observationally by combining high angular resolution ($\sim$few arcseconds) observations of molecular lines with polarimetric observations of dust continuum (e.g., Girart et al., 2013; Hull et al., 2014; Zhang et al., 2014). However, this type of studies has been done only towards a limited sample of star-forming cores. Within this sample, there are cases where the magnetic field is consistent with an hourglass morphology, which is an indication of strong magnetic fields in infalling cores (Fiedler & Mouschovias, 1993; Galli & Shu, 1993a, b; Nakamura & Li, 2005; Frau et al., 2011). Examples include the low- and high-mass star-forming regions NGC1333 IRAS4A (Girart et al., 2006; Frau et al., 2011; Liu et al., 2016), G31.41+0.31 (Girart et al., 2009), W51e2 (Tang et al., 2009), W51 North (Tang et al., 2013), G35.2-0.74N (Qiu et al., 2013), G240.31+0.07 (Qiu et al., 2014), W43-MM1 (Cortes et al., 2016).

NGC 6334 is one of the nearest high-mass star-forming regions, located at a distance of 1.3 kpc (Chibueze et al., 2014; Reid et al., 2014). The dense molecular gas is dominated by a $\sim$20 pc long filamentary structure oriented along northeast-southwest direction. Along the filament there are several parsec-scale sites of active star formation (e.g., Li et al., 2015, and references therein). Previous optical and submillimeter polarimetric observations found a coherent magnetic field from cloud to core scales with its main direction approximately perpendicular to the long axis of the dense filament. This led Li et al. (2015) to suggest that NGC 6334 undergoes self-similar fragmentation regulated by the magnetic field.

NGC 6334 V is the southernmost active star-forming site along the large-scale northeast-southwest dense gas filament of NGC 6334. This region was not covered by the multi-scale studies of Li et al. (2015). NGC 6334 V is located at the intersection between the main filament and a smaller filament that extends along the northwest-southeast direction (e.g., André et al., 2016). This characteristic makes the NGC 6334 V star-forming core different from other cores along the main filament. Therefore, the study of this core, at scales similar to Li et al. (2015), can be used to test whether the physical properties of NGC 6334 V are affected by the convergence of these two filaments.

Signposts of recent and ongoing star formation in NGC 6334 V include its high far-infrared luminosity (L${}_{\text{bol}}\sim 10^{4}$ L*⊙*), a CO molecular outflow (Fischer et al., 1982; Kraemer & Jackson, 1995; Zhang et al., 2014) and H2O and OH masers (Forster & Caswell, 1989; Raimond & Eliasson, 1969; Brooks & Whiteoak, 2001). It also harbors several infrared sources (Simon et al., 1985; Harvey & Gatley, 1983; Harvey & Wilking, 1984; Kraemer et al., 1999) and three faint associated radio sources (Rengarajan & Ho, 1996).

Previous studies (Hashimoto et al., 2007; Simpson et al., 2009) suggest the existence of two independent outflows in NGC 6334 V, one oriented east-west along the line of sight and the other with a north-south orientation almost in the plane of the sky. These two outflows could be powered by the infrared sources KDJ3 and KDJ4 (Kraemer et al., 1999).

In this paper, we present Submillimeter Array (SMA)111The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics (Ho et al., 2004). observations towards NGC 6334 V at 345 GHz. These observations covered a frequency range that allows us to detect the dust emission with a high signal-to-noise ratio together with linear polarization from the dust emission. Within the frequency range there are various molecular lines (see section below for more details on the specific lines): dense gas tracers (i. e., sensitive to densities of $\gtrsim 10^{6}$ cm*-3* and temperature $\gtrsim 100$ K) and tracers of outflows or shocks. This allows us to characterize the physical properties of the core, and in particular to shed light on the role of the magnetic field.

The details of our observations are given in Section 2. The observational results are shown in Section 3. In Section 4 we present magneto-hydrodynamic (MHD) simulations. Radiative transfer analysis and convolution with the SMA response are performed to properly compare with the observations. Discussion and conclusions are given in Section 5 and 6, respectively.

2. SMA observations and data reduction

NGC 6334 V was observed with the SMA in the subcompact, compact and extended configurations. These observations were part of the polarization legacy project carried out between 2008 and 2012 in the 345 GHz band to observe the polarization emission of a sample of 21 massive star-forming regions (for an overview of this survey see Zhang et al., 2014). Table 1 shows the basic observational parameters and lists the calibrators used during the three tracks. The combined three configurations sample visibilities between $\sim$6 k$\lambda$ to $\sim$200 k$\lambda$. This makes the observations sensitive to scales between $1^{\prime\prime}$ and $15^{\prime\prime}$. The phase and pointing center of the observations was R.A.=17:19:57.40, DEC=$-$35:57:46.0 (J2000). The frequency was approximately centered on the CO (3$-$2) line (345.79599 GHz) in the upper sideband (USB). The correlator consisted of 48 chunks with a bandwidth of 104 MHz each. This provided a total bandwidth (including the USB and the lower side band - LSB) of $\sim 8$ GHz. We set the correlator to yield a uniform spectral resolution of 128 channels per chunk, providing a channel width of 0.8 MHz or 0.7 km s*-1* at the observing frequency.

Polarization observations were carried out using quarter-wave plates which convert linear polarization signals to circular ones (see Marrone et al., 2006; Marrone & Rao, 2008, for a more detailed description of the SMA polarimeter and the calibration process). The polarization leakages were measured to an accuracy of 0.1% (Marrone & Rao, 2008). The calibration for absolute flux, passband, and gains were carried out using the mir idl software package (Scoville et al., 1993). The images were created using the Miriad software package (Sault et al., 1995). We averaged the line-free channels in the LSB and USB to generate the 870 $\mu$m continuum. The 870 $\mu$m maps were obtained combining the three configurations. We set the Robust visibility weighting parameter to 0.5 (Briggs, 1995), which provides a good compromise between angular resolution and sensitivity. The rms noise level of the resulting continuum map (Fig. 1) is 5.5 mJy beam*-1* and the synthesized beam is $2\farcs 15\times 1\farcs 49$, P.A. $=-2.92^{\circ}$. To generate the polarization maps we used the impol task in Miriad which computes the linearly bias-corrected polarized intensity and position angle from Stokes Q and U (QU) images. We used 3$\sigma$ as the cutoff for Stokes QU (Vaillancourt, 2006) and I images. The rms level for Stokes QU is 3.1 mJy beam*-1*. The different molecular line maps were also generated using a Robust parameter of 0.5, except for the SiO (8$-$7) line. This line is strong and compact. Thus, we used a Robust parameter of $-2$ (uniform weighting), which provides the best possible angular resolution.

3. Results

3.1. Morphology

3.1.1 Dust continuum

The upper panel of Figure 1 shows the 870 $\mu$m dust continuum emission image. Overplotted are associated infrared and radio sources, as well as OH and H2O masers. This continuum image reveals asymmetric dense gas structures. They form a well resolved $\sim 0.07$ pc ($\sim$14,000 AU) scale arc-like structure along the east-west direction (measured from the $8\sigma$ contour from the 870 $\mu$m dust continuum emission), with three distinguishable peaks surrounding a lower density cavity, labeled as N (North), SE (South-East) and HC (Hot Core; see Section 3.1.3). Positions, deconvolved sizes, peak intensities, flux densities and mass estimates of the three main peaks are listed in Table LABEL:tab:dustpar. Most of the known infrared and radio sources are located either around the 870 $\mu$m peak (hot core; e.g., IRS V3 and radio source R E-3), or are inside the cavity.

Based on the VLA observations of the centimeter free-free continuum emission presented in Rengarajan & Ho (1996), we constrain the 870 $\mu$m free-free emission to be < 4.7 mJy towards R-E1, R-E2 and R-E3, which is negligible as compared to the detected emission level in Figure 1 (see table LABEL:tab:dustpar). Assuming that the 870 $\mu$m dust emission is optically thin, we estimated the gas+dust masses according to

[TABLE]

where $d$ is the distance, $S_{\nu}$ is the flux density, $B_{\nu}(T_{d})$ is the Planck function at the isothermal dust temperature $T_{d}$, $\kappa_{\nu}$ is the dust opacity and the factor R is the gas-to-dust ratio which we assumed to be 100 (Hildebrand, 1983; Lis et al., 1998). We adopted $\kappa_{\nu}=2.03$ cm2g*-1* for thin ice mantles after $10^{5}$ years of coagulation at a gas number density of $10^{6}$ cm*-3* and for the frequency of our observations (Ossenkopf & Henning, 1994) and we assumed $T_{d}=50-100$ K (hot core temperatures $\gtrsim 100$ K; e.g., Palau et al., 2011) for the hot core and 30 K for the north (N) and southeast (SE) cores (the average temperature in protostellar cores, e.g., Faúndez et al., 2004; Sánchez-Monge et al., 2013, see Section 3.2.2). The resulting masses are given in Table LABEL:tab:dustpar. The mass derived from the total 870 $\mu$m flux density detected by the SMA is $\sim 50$ M*⊙* for a $T_{d}=30$ K. This value is <10% of the mass estimated by Muñoz et al. (2007) based on 1.2 mm dust continuum emission obtained with the SEST SIMBA bolometer array. There are two reasons for this difference. We might be missing some flux due to the interferometric filtering of extended emission, and Muñoz et al. (2007) use an integration area for the flux calculation approximately ten times larger than our work (using an effective radius of 0.86 pc).

3.1.2 Polarization data

It is generally accepted that non-spherical interstellar dust grains are locally aligned with the B-field, with their longest axes being perpendicular to the B-field directions (see Hildebrand et al., 2000, for a review). Therefore, the projected orientation of the B-field can be inferred by rotating the dust linear polarization (i.e., electric field) by 90 degrees (hereafter B-segment for the 90*∘* rotated linear polarization).

The lower panel of Figure 1 shows the polarization map from the 870 $\mu$m continuum observations. The projected magnetic field shows an east-west orientation following the dust morphology, from the south-east (SE) continuum peak towards the hot core region. There are also several detections at the northern part of the dust continuum emission showing a near north-south orientation. This northern component of the magnetic field seems to connect to the hot core position through several detections presenting a northeast-southwest orientation.

The histogram of the B-segment position angles (see upper panel of Figure 2) shows a disperse distribution due to the changes in the morphology of the magnetic field. We identified two major orientations, at $\sim 90^{\circ}$ (east-west component) and $\sim 50^{\circ}$ (northeast-southwest component) which seem to be converging at the hot core region. The northern component at $\sim 160^{\circ}$ (north-south) is also apparent as a separate peak.

3.1.3 Molecular lines

We have detected 24 molecular lines in the 8 GHz band. Table 3 lists the frequency and the energy of the upper level (EU) for each transition. We found that the emission from several molecular transitions with EU between $\sim 100-400$ K spatially coincides with the strongest dust continuum peak (see Fig. 3 and Fig. 5), including transitions from SO2, 34SO2, HC3N, CH3OH, CH3OCH3 and CH3OCHO. Thus, these molecular transitions likely trace a hot molecular core.

In Fig. 3 we present the moment 0 (integrated intensity) maps of H13CO*+* (4–3), CH3OH 7(1,7)$-$6(1,6), and CH3OCHO 31(0,31)$-$30(1,30), which have upper level energies of 41 K, 79 K, and 260 K, respectively. The very different excitation conditions for these molecular line transitions trace different regions in the molecular gas, which show a variety of morphologies. H13CO+ (4$-$3) (upper panel) is a good tracer of the dense and relatively cold gas and traces the extended dense core of NGC 6334 V. It shows two main peaks forming a ring-like structure with lack of emission at the center of the ring, coinciding with the dust continuum cavity. CH3OH 7(1,7)$-$6(1,6) (middle panel) traces warmer gas, and it is generally a good hot core tracer. It presents strong emission at the hot core position and also traces a more compact region of the dense core, forming a ring-like structure, as the H13CO+ molecular line.

Many of the infrared and radio sources found in the region are located where we find a deficit of molecular emission in both H13CO+ and CH3OH (center of the ring). In addition, the cavity is filled with 4.5 $\mu$m emission (see Fig. 4). The origin of this extended emission is unclear. Possible emission mechanisms could be scattered continuum emission in outflow cavities (e.g., Qiu et al., 2008; Takami et al., 2012), H2 emission and CO emission from shocks (e.g., Noriega-Crespo et al., 2004; Smith et al., 2006; Davis et al., 2007; Takami et al., 2012). In any case, the 4.5 $\mu$m emission is likely to be a signpost of a cavity generated by a still active outflow.

Finally, in the lower panel of Fig. 3, we present a high-excitation transition (EU=260 K) of CH3OCHO which shows more compact emission elongated along an east-west direction, centered only at the hot core region.

3.2. Kinematics

A Gaussian fit to the averaged spectrum of the entire emission of the dense core tracer H13CO+ (4$-$3) suggests that the systemic velocity of NGC 6334 V is $-5.7$ km s*-1*. Therefore, we use this value for the following kinematic analysis.

3.2.1 Velocity gradients and mass estimates

A comparison between the different gas tracers can provide clues about the kinematics from larger (extended dense core) to smaller (hot core) scales. Figure 5 shows the intensity-weighted averaged velocity maps (i.e., moment 1 maps) of the same molecular lines shown in Fig. 3. The extended dense core tracers H13CO+ and CH3OH (upper and middle panels, respectively) show a velocity field dominated by a clear east-west gradient from $-8$ to $-5$ km s*-1*. However, the overall velocity field has a quadrupolar pattern centered at the hot core region. This excludes a simple interpretation for the velocity field, such as a rotating core.

To know if the observed gas belongs to a gravitationally bound system, we estimated the needed minimum mass associated to the observed velocity gradient. The maximum velocity for an object accelerated in a gravitationally bound system is given by $v_{\mathrm{max}}=(2GM/R)^{1/2}$ where $G$ is the gravitational constant, $R$ is the distance to the gravitational center and $M$ is its mass. Adopting the gravitational center as the center of the velocity gradient located at the hot core position and using a distance of $6^{\prime\prime}$ and a velocity of 2 km s*-1* from the observed velocity gradient, we derive a mass of $\sim 18$ M*⊙*. Note that the radial velocity is a lower limit of the total velocity. This mass is lower than the mass estimated from the dust continuum emission, which is a lower limit considering that we are filtering extended emission with the SMA (see Section 3.1.1), and thus, we consider the system to be gravitationally bound.

The hot core tracer CH3OCHO also presents a clear velocity gradient following the orientation of the larger-scale velocity field. The intermediate velocity $-6.5$ km s*-1* is shifted from the dust peak position $0.5^{\prime\prime}$ or 650 AU towards the east (several transitions of CH3OH, CH3OCH3, 34SO2 and SO2 also present the same velocity field). The velocity gradient is also compatible with a gravitationally bound system, with an enclosed mass of $\sim 6$ M*⊙* (using $R=2^{\prime\prime}$ and $v=2$ km s*-1*). The lower limit of the hot core mass obtained from the dust continuum emission is $\sim$4–9 M*⊙*.

3.2.2 High-velocity gas

Figure 6 shows the red- and blue-shifted CO (3$-$2) emission at velocities of $\pm 18$, $\pm 24$, $\pm 30$ and $\pm 36$ km s*-1* from the systemic velocity, $-5.7$ km s*-1*. At these high velocities, the CO (3$-$2) is tracing molecular outflows. The first panel ($\pm 18$ km s*-1*) shows a very complicated structure with red- and blue-shifted emission spread all over the dust continuum. In the second panel ($\pm 24$ km s*-1*), we can distinguish collimated blue-shifted emission along the dust continuum cavity. This emission suggests that the high-velocity gas may have produced the observed cavity. Furthermore, in this panel, we observe blue- and red-shifted emission extended towards the north and south of the hot core position, respectively, in agreement with the larger-scale previously reported CO outflow (Kraemer & Jackson, 1995; Zhang et al., 2014). This emission is oriented almost perpendicular to the velocity gradient seen at the hot core region (see bottom panel of Fig. 5). Note that outflow powering source candidates, KDJ3 and KDJ4 (Kraemer et al., 1999; Hashimoto et al., 2007; Simpson et al., 2009), are located close to the hot core position (see Fig. 1). KDJ3 does not have dust continuum emission associated as it is located near the dust cavity. On the other hand, KDJ4 seems to be a younger object than KDJ3 and a better outflow powering source candidate as it coincides well with the hot core dust continuum peak.

The highest outflow velocities, at $\pm 30$ and $\pm 36$ km s*-1* in the last two panels of Fig. 6, show red-shifted emission near the southeast continuum peak (SE) and blue-shifted emission towards the hot core region, the northern continuum peak and the dust cavity. Especially in the third panel (at $\pm 30$ km s*-1*), it seems that the southeast (SE) and north (N) cores are associated with high-velocity gas. Thus, even though the CO high-velocity emission presents a very complicated structure, it suggests the presence of at least two or three independent outflows. Multiple molecular outflows associated with a massive dense core have been previously reported in several high-mass star-forming regions (e.g. Naranjo-Romero et al., 2012; Girart et al., 2013; Fernández-López et al., 2013).

4. Analysis

4.1. Converging flows

As we have seen in Section 3.2.1, the molecular dense core traced by H13CO+ and CH3OH presents a velocity gradient with a certain symmetry centered at the hot core region (see upper and middle panels in Fig. 5). Figure 7a shows the blue- and red-shifted integrated emission of H13CO+ (4–3) showing two filamentary structures at the eastern and northern parts of the hot core. Both velocity structures seem to converge to the hot core. At the hot core region, the velocity gradient follows the orientation of the larger-scale velocity field (see Fig. 5). As the center of the velocity gradient is shifted from the dust peak, the observed velocity gradient is not likely due only to rotational motions but additional motions may be also affecting it (e.g., infalling motions).

In Figure 7b and c we show two position-velocity diagrams, along the east-west and north-south direction centered at the hot core region, to compare the kinematics from the extended dense core (traced by H13CO+) and the hot core region (traced by CH3OCHO). Both dense core and hot core tracers show the same two distinctive velocity components, one at $-5.0$ km s*-1* and another one at $-7.4$ km s*-1*. In addition, we did the same position-velocity cuts towards the shock tracer SiO (8–7). The SiO emission is a signpost of strong shocks. It has, however, also been detected at velocities close to the systemic velocity without a clear association to protostellar outflow activity (Jiménez-Serra et al., 2010; Nguyen-Lu’o’ng et al., 2013). In these cases, SiO can be produced in low-velocity shocks (e.g. Duarte-Cabral et al., 2014; Girart et al., 2016). As shown in grey scale in Fig. 7b and c, SiO presents a velocity component just at an intermediate velocity of $\sim-6.0$ km s*-1* and at the dust continuum peak position, suggesting that there could be interaction between the two filamentary structures precisely at the hot core region. Such a converging flow toward a hot core was also reported in W33A by Galván-Madrid et al. (2010).

Furthermore, the magnetic field distribution described in Section 3.1.2 also shows a “bimodal converging” pattern towards the hot core. In the upper panel of Figure 8 we overlaid the moment 1 map of H13CO+ with the magnetic field distribution. We used the polarization detections from only compact and subcompact configurations as the polarized emission is stronger than with all three configurations combined. The morphology of the magnetic field is very well correlated with the two velocity components shown by the H13CO+ velocity field, suggesting that the magnetic field could be being dragged by the gas dynamics. Previous observations of the B-type binary forming region G192.16-3.84 (Liu et al., 2013) have shown that the magnetic field is also dragged by the gravitationally dominant infall and rotational motions, where the magnetic field orientation is approximately parallel to the velocity gradient of the gas.

The combination of these results shows how the gas seems to be infalling from the larger-scale dense core of NGC 6334 V towards the higher-density hot core region through two distinctive converging flows dragging along the magnetic field whose strength seems to have been overcome by gravity.

4.2. Synthetic observations

Our observational analysis so far points out that the gas is converging towards the hot core through two main infalling flow components, suggesting that gravity could play a major role in the dynamics of intermediate/massive star-forming regions. However, our interpretation can be biased due to several observational effects. First, certain information has been lost when projecting a 3-dimensional density distribution onto a 2-dimensional plane. In addition, the vector fields (e.g., B-field, mean velocity, etc) were smeared due to the integration along the line-of-sight. Moreover, the interferometric observations of molecular lines are biased by the excitation conditions and the incomplete sampling in the Fourier domain. To facilitate the test of whether or not the observational results are indeed consistent with our interpretation, we have carried out 3-D numerical magnetohydrodynamic (MHD) simulations for massive star-forming regions which are dominated by gravity. We produced synthetic observational results based on these simulations, to make a consistent comparison with the observational results of SMA.

4.2.1 Numerical simulation

We carry out a 3-D numerical simulation of a massive, compact molecular cloud, similar to that presented in Ballesteros-Paredes et al. (2015), although with a different numerical scheme: rather than the Lagrangian smoothed particle hydrodynamical code used before, we use the Eulerian adaptive mesh refinement (AMR) FLASH (version 2.5) code (Fryxell et al., 2000) in a magneto-hydrodynamical regime. As in Ballesteros-Paredes et al. (2015), we include self-gravity, and the formation of sink particles, in an initially turbulent velocity field which is let to evolve freely, i.e., no turbulent forcing is imposed. The simulation is isothermal in order to isolate the thermal effects on the cloud evolution. We also did not include radiative feedback or stellar heating from the sink particles in order to understand the very early stages of the molecular cloud collapse.

The initial conditions are as follows. We consider a numerical box of 1 pc per side containing 1000 M⊙ of molecular gas. The mass in the box is distributed homogeneously, with a number density of $n_{0}=1.7\times 10^{4}\,\mathrm{cm}^{-3}$. The free-fall time for this configuration is 256,000 yr. Following the prescription by, e.g., Stone et al. (1998), we included a pure rotational velocity spectrum with random phases and amplitudes that peak at wave numbers $k=4\pi/L_{0}$, where $L_{0}$ is the linear size of the box. The resulting initial velocity field was, thus, an incompressible supersonic turbulent fluid, with an rms Mach number of $\mathcal{M}_{\rm rms}=8$. No forcing at later times is imposed.

The magnetic field is initially uniform along the $x$-direction with a strength of $50\,\mu\mathrm{G}$, which is consistent with magnetic field intensities at densities of the order of $n_{0}$ (see e.g., Crutcher, 2012). Our box is, thus, magnetically supercritical, with a mass-to-flux ratio, $\mu$, 6.7 times the critical value, $\mu_{\mathrm{crit}}=(4\pi^{2}G)^{-1/2}$ (see e.g., Nakano & Nakamura, 1978), and prone to gravitational collapse as soon as the initial turbulent field is dissipated.

For the dynamical mesh refinement, we use a Jeans criterion in which we resolve the local Jeans length with at least 8 grid cells in order to prevent spurious fragmentation (Truelove et al., 1997). Thus, the numerical grid is refined to reach a maximum resolution of $\Delta x\sim 9.8\times 10^{-4}\,\rm{pc}$ ($\sim 201\,{\rm AU}$). Once the maximum refinement level is reached in a given cell, no further refinement is performed and a sink particle can be formed when the density in this cell exceeds a threshold density, $n_{\rm thr}=1.7\times 10^{7}\,{\rm cm}^{-3}$. The sink particles can then accrete mass from their surroundings (within an accretion radius of $\sim 2.5\Delta x$), increasing their mass. The numerical scheme for creation of sinks and further mass accretion is that of Federrath et al. (2010).

We run the simulation until 0.9 times the free-fall time ($t_{\rm ff}$) and we compare the numerical model with the observations at this point, because gravity dominates the gas dynamics. In general terms, the evolution can be considered quite similar to that shown in Ballesteros-Paredes et al. (2015): the initial supersonic velocity fluctuations shock and rapidly dissipate their kinetic energy, the gas then collapses in a global but hierarchical and chaotic way forming a complex network of cores and filaments, carrying with them the magnetic field lines. The time evolution of the simulation is illustrated in Fig. 9, which shows column density maps in the $x-y$ plane at $t=$ 0.1, 0.3, 0.6 and 0.9 $t_{\rm ff}$. As it can be seen, the first time steps are dominated by density fluctuations that grow in mass and merge as collapse proceeds. These density fluctuations are elongated mostly in the $y$-direction because the initial magnetic field is aligned along constant values of $x$, allowing the mass to flow in the $y$-direction more freely than in the perpendicular directions (in contrast to the initial fluctuations of the non-magnetic simulations presented earlier in Ballesteros-Paredes et al., 2015). Nevertheless, as the collapse proceeds, no preferred elongation is seen in the column density maps, and the final structure is quite similar to the non-magnetic case (see Fig. 1 in Ballesteros-Paredes et al., 2015), suggesting that, indeed, gravity is driving the dynamics of the dense collapsing core.

The onset of star formation (sink formation) occurs at $0.75\,t_{\rm ff}$, once the cloud is globally collapsing. We compare the stellar content from both the simulation and observations in Sec. 5.1.

Finally, to compare the numerical model with the simulations, we extract, at $0.9\,t_{\rm ff}$, the central cubic subregion of 0.3 pc per side in order to match the angular size with the observations (see Fig. 11).

It is important to notice that even though the magnetic-to-flux ratio should be conserved in MHD simulations, this is valid only in an average sense over the whole box. However, as a consequence of the collapse the magnetic field lines increase their values locally, as in the portion of the simulation that we have analyzed (Fig. 11), although it starts magnetically supercritical, it ends up with equipartition values. The fact that cores and or clumps are observed in equipartition (virial) between gravity and their internal energies (magnetic, kinetic) has been interpreted as a signature of collapse: these energies do not have to know what are their detailed values to keep the cloud in equipartition, and the more reasonable physical reason for this is that gravity dominates bulk, chaotic and hierarchical motions, dragging the magnetic field as collapse proceeds (e.g., Ballesteros-Paredes, 2006; Ballesteros-Paredes et al., 2011).

4.2.2 ARTIST

We utilize ARTIST222URL: http://youngstars.nbi.dk/artist/Welcome.html (Adaptable Radiative Transfer Innovations for Submillimetre Telescopes; Brinch & Hogerheijde, 2010) to generate a synthetic intensity cube (RA, DEC, $v_{\mathrm{LSR}}$) of the (4–3) transition of H13CO+ from the 3-D structure generated in the numerical simulation. We take a fractional abundance of $3\times 10^{-11}$ with respect to the H2 density (see e.g., Girart et al., 2000), and we assume a dust temperature of 30 K, a source distance of 1.3 kpc, and no inclination/rotation, i.e., face-on view. The angular resolution is 0.05 arcsec.

Additionally, we utilize the DustPol module of ARTIST (Padovani et al., 2012), to generate a synthetic map of the dust polarized emission at $870\,\mu\mathrm{m}$ on the plane of the sky from the density distribution of the magnetic field generated in the numerical simulation.

4.2.3 SMA response

The numerical simulations and the H13CO+ (4$-$3) ARTIST spectral cube were convolved with the SMA interferometric response for a realistic comparison with the observations. This was done using the uvmodel task in Miriad by converting the modeled maps to visibilities using the same distribution of the visibilities in the (u, v)-plane as in the NGC 6334 V SMA observations. Thus, we filtered the same large-angular scale structures as in the real SMA observations. Later, the same parameters were used to create the final maps (e.g., pixel size, Robust parameter, …). For the simulated maps, the procedure was followed independently for Stokes I, Q, and U. Once the synthetic SMA-like Stokes I, Q, U maps were obtained, the polarization intensity and position angles were obtained in the same way as for the SMA data.

4.3. Observations vs. simulation: polarization angular dispersion function comparison

From dust polarization observations we are able to estimate the magnetic field strength in the plane of the sky using the Chandrasekhar-Fermi (CF) equation (Chandrasekhar & Fermi, 1953),

[TABLE]

where $B$ is the strength of the magnetic field, $\delta B$ the variation of $B$, $V_{A}=B/\sqrt{4\pi\rho}$ is the Alfvén speed at mass density $\rho$, and $\sigma_{v}$ is the velocity dispersion along the line of sight. Different statistical methods have been developed in order to avoid some of the CF method caveats (Hildebrand et al., 2009; Houde et al., 2009, 2011; Franco et al., 2010; Koch et al., 2010).

Assuming two statistically independent components of $B$ (the large-scale magnetic field $B_{0}$ and the turbulent magnetic field $B_{t}$), Houde et al. (2009) suggest that the ratio of $B_{t}$ to $B_{0}$ can be evaluated from the angular dispersion function that accounts for the polarization angle differences ($\Delta\Phi$) as a function of the distance ($l$) between the measured positions.

The angular dispersion function can be written as

[TABLE]

where

[TABLE]

is the number of turbulent cells along the line of sight, $\delta$ is the magnetic field turbulent correlation length (assumed to be much smaller than the thickness of the cloud $\Delta^{\prime}$), $W$ is the synthesized beam (i.e., FWHM/$\sqrt{8~{}ln(2)}$), and the summation is a Taylor expansion representing the large-scale magnetic field component which does not involve turbulence. For displacements $l$ less than a few times $W$ we keep only the first $l^{2}$ term in the Taylor expansion. We do not use Houde et al. (2016) approximation which takes into account the interferometer filtering effect, as we want to simplify the variables to perform a qualitative comparison between the observational data and the simulation. In addition, the simulation has been convolved by the SMA response (see Section 4.2.3), thus, the effect should be the same as in the observational data.

The first term on the right hand side in equation (3) contains the integrated turbulent magnetic field contribution, while the second (exponential) term represents the correlation by the combined effect of the beam ($W$) and the turbulent magnetic field ($\delta$). The value of the correlated component at the origin $l=0$, $f_{\mathrm{NC}}(0)$, allows us to estimate the turbulent to large-scale magnetic field strength ratio as

[TABLE]

To fit the angular dispersion function to our results and to be able to perform a qualitative comparison between the observations and the numerical simulations, we used two fixed values of the parameter $\delta$, 17 mpc (from Girart et al., 2013) and 25 mpc which produced a better fit (see Fig. 10). The derived value of the correlated component at the origin is $f_{\mathrm{NC}}\simeq 0.29,0.52$ and $\simeq 0.13,0.28$ for both values of $\delta$ and the observations and simulations respectively (see Table 4). A reasonable approximation to the core’s effective thickness, $\Delta^{\prime}$, is to assume that it is similar to the average diameter of the dense core measured in the plane of the sky with the SMA (Koch et al., 2010), which in our case is $\simeq 0.07$ pc (see Section 3.1.1). Using eqs. (4) and (5) $\langle B_{t}^{2}\rangle/\langle B_{0}^{2}\rangle$ is $\simeq 0.7$ for the observations, i.e., it is close to equipartition between the turbulent and ordered magnetic field energies. This is similar to the results obtained towards the high-mass star-forming region DR21(OH) (Girart et al., 2013). However, for the simulations the values are lower ($\simeq 0.3,0.4$). One possibility for this difference is that the size of the core in the simulations could be larger than the observed one, as $N$ is proportional to $\Delta^{\prime}$. Also, the thermal Q and U instrumental noise produced by the SMA can induce an increase of the angular dispersion which is not taken into account in the simulations (see upper panel of Fig. 10).

To estimate the large-scale magnetic field strength in the plane of the sky we used the CF equation, $\langle B_{0}^{2}\rangle^{1/2}=\sqrt{4\pi\rho}~{}\sigma_{v}~{}[\langle B_{t}^{2}\rangle/\langle B_{0}^{2}\rangle]^{-1/2}$, (e.g., see Equation (57) by Houde et al., 2009). Table 4 shows the values used for the velocity dispersion $\sigma_{v}$ and the number density $n$ (using a mean molecular weight of 2.33). To derive the averaged volume density we used the mass derived from the total 870 $\mu$m flux density detected with the SMA ($\sim 50$ M*⊙) and the size of the core ($\sim 0.07$ pc). We estimated the velocity dispersion from the H13*CO+ (4–3) data, because it is well correlated with the 870 $\mu$m dust emission (see upper panel of Fig. 3). The derived value for the ordered large-scale magnetic field strength in the plane of the sky for both observations and simulations is $\langle B_{0}^{2}\rangle^{1/2}\simeq 1$ mG. This value is in very good agreement with the average value of the magnetic field (weighted by density) measured in the simulation at 0.9 $t_{\rm ff}$, 0.9 mG.

5. Discussion

5.1. Simulation results and comparison with observations

In the bottom panel of Figure 8 we present the results from the simulations filtered out with the SMA response. The map shows different velocity structures at scales of $\sim 0.1$ pc which seem to converge towards the strongest peak of the continuum and with a similar velocity range as the one present in the observed data. As in the case of the observations, in this figure we notice that the field has three main orientations: segments that are almost horizontal, but that have two components: at +70 and +100 degrees, (also seen as a double peak in the histogram of the lower panel of Fig. 2). The third component is not highly populated, but it is clearly traced by the almost vertical white segments shown in Fig. 8. This component corresponds to the long tail towards small angles of the 70-degree component in the histogram shown in the lower panel of Fig. 2. Furthermore, note that at least the two main components, at +70 and +100 degrees, coincide with the two different velocity structures traced by H13CO+.

The magnetic field configuration shown in the lower panel of Fig. 8 can be understood in terms of the dynamics of the core. In Fig. 11 we show the gas column density (colors) and the $x-y$ velocity field of the gas (thin black arrows) after evolving the simulation. Note that this representation has no observational counterpart, and although the simulation resolution is $\sim 200$ AU, the column density map shown in this figure has been convolved with a beam of the same size of the observations. The column density map shows filamentary structures which seem to converge at the center of the box, where the highest column densities are found. The velocity field shows how the gas is converging from the dense structures at the large scale ($\sim$0.2 pc) towards the highest density region, at $\sim$0.02 pc scales, and they appear to be almost radial. Thus, it is clear that the core is collapsing. The originally horizontal magnetic field has been advected at this time by the collapsing gas, resulting in an hourglass shape with a dominant nearly horizontal component at the center of the core, with position angles of $+90\pm 20$ degrees. This dominance is just the consequence of the B-field being compressed and dragged by the collapsing gas. We notice also that in the outer parts of the core, the B-field does have components that are far from parallel to the original magnetic field. Although it is a coincidence that the simulations exhibit three main orientations of the magnetic field, as in the observations, the important point here is that the B-field, in the innermost denser part of the core, seems to be oriented mostly horizontal, while in the outermost parts it looks more radial.

It can be asked whether the selected angle of view of the simulations affects what we observe in Fig. 8. Certainly, as the magnetic field is a vector quantity, any point of view could be considered as “particular”. And we have chosen a point of view that is somehow special: it is perpendicular to the initial magnetic field (the box is seen along the $z$-axis, while the initial magnetic field runs along the $x$-axis). Given this configuration, and since the collapse occurs almost radially, one may assume that the density and magnetic field must have, in a statistical sense, an axial symmetry with respect to a declination offset=0 (Fig. 8). This means that, in principle, if we had observed the core along any other line of sight that was perpendicular to the $x$-axis, we would have observed something similar: nearly horizontal orientations of the field in the innermost parts of the core, and more radially directed orientations in the external parts of the core. However, if we were observing the simulation along the $x$-axis, the horizontal, dominant component of the B-field would almost not be observable. The B-field should look almost radial in the outskirts of the core, with a decreasing component at the center. The more general situation will be something in between these two possibilities, i.e., radial lines in the external parts of the core, and a less-dominant magnetic field in the center of the box.

As we do not include feedback processes in the numerical model, we cannot model the hot core itself (and the associated outflows). However, in the densest part of the simulation, where we locate the “Hot Core” labeled in Figure 11, we have a group of three stars (sinks) with a total mass of $\sim 9.7$ M*⊙* in a compact space of 0.03 pc ($\sim 4.7^{\prime\prime}$, assuming a distance of 1.3 kpc), which is comparable with the synthesized beam of our observations. It is remarkable that this mass is comparable with the mass inferred from kinematics arguments in the observed hot core ($\sim 6$ M*⊙*; see Section 3.2.1), which give us confidence that we are comparing our observations with an adequate model.

Finally, in addition to the good qualitative comparison of the velocity and magnetic field structures between the SMA observations and the numerical simulations, the angular dispersion function analysis also reveals a very similar behavior for both (see Section 4.3 and Fig. 10). The simulations, thus, again seem to be a good approximation to the observational results.

5.2. Outflow-generated cavity

We have seen in Section 3.1.3 that the emission from the extended dense core of NGC 6334 V (traced by H13CO+ (4–3) and CH3OH 7(1,7)$-$6(1,6)) presents a clear ring-like structure around the systemic velocity, $-5.7$ km s*-1* (see Fig. 3). The infrared and radio sources found in the region are located at the center of the ring, that also coincides with the cavity of the dust continuum emission. In addition, we detected high-velocity gas traced by the CO (3$-$2) molecular line at the position of the cavity and at the center of the ring structure (see Fig. 6). 4.5 $\mu$m emission is also detected crossing through the center of the ring with an east-west elongated shape presenting a bow shock structure towards the west, suggesting the emission could be shock- or outflow-generated (see Fig. 4). These results suggest the presence of an east-west outflow (also reported by Hashimoto et al., 2007; Simpson et al., 2009) that could be affecting the central region of NGC 6334 V, causing a cavity in both the extended dense gas as well as in the dust continuum structure.

5.3. Comparison with other NGC 6334 star-forming sites

In a previous work, Li et al. (2015) studied the magnetic field properties in several star-forming regions inside NGC 6334, from 100 down to 0.01 pc scales. Source V was not included in this study as large-scale polarization observations were not available. They concluded that the well ordered magnetic field (at all scales) in sources I, I(N) and IV plays a crucial role in the fragmentation of NGC 6334.

However, in NGC 6334 V, we see a more complex magnetic field morphology governed by the gas dynamics. It seems that NGC 6334 V, which is located at the southern end of the NGC 6334 ridge and possibly at the intersection region with another northwest-southeast gas filament, with gravity already being the dominant force over the magnetic field (e.g., André et al., 2016), may be more dynamically evolved than the sources I, I(N) and IV (note that the total mass in NGC 6334 V seems to be smaller than in source I and I(N); Hunter et al., 2006).

6. Conclusions

We analyzed SMA observations at 345 GHz towards the intermediate/high-mass star-forming region NGC 6334 V. The main results are as follows:

• •

The dust continuum emission presents an arc-like structure of $\sim$0.07 pc ($\sim$14,000 AU), with three distinguishable peaks forming a cavity structure (Fig. 1). The total mass of the arc-like molecular structure is $\sim 50$ M*⊙. The main peak (HC) has a mass of $\sim 4$–9 M⊙* and exhibits the typical chemistry of a hot core with emission from 34SO, SO2, 34SO2, CH3OH, CH3OCH3, CH3OCHO and HC3N. Dense core tracers H13CO+ (4–3) and CH3OH 7(1,7)$-$6(1,6) show extended emission forming a ring-like structure at the systemic velocity of $-5.7$ km s*-1*. The presence of an outflow at the center of the region could be affecting the surrounding gas forming the cavity traced by the dust continuum and the dense gas emission from H13CO+ and CH3OH.

• •

Kinematically, the dense core has two distinctive velocity components, one at $-5.0$ km s*-1* tracing a filamentary-like arm on the southeast part of the core and another one at $-7.4$ km s*-1* arising mainly from the northern part (Figs. 5, 7). Both velocity structures converge towards the hot core. At higher density regions in the hot core, the hot core line tracers CH3OCHO 31(0,31)$-$30(1,30) and several transitions of CH3OH, CH3OCH3, 34SO2 and SO2 show the same two velocity components. The shock tracer SiO (8–7) (Fig. 7b,c) presents a velocity component just at an intermediate velocity ($\sim-6$ km s*-1*) suggesting interaction between the two flows.

• •

The magnetic field (derived from the dust thermal polarized emission at 870 $\mu$m) shows a bimodal converging pattern towards the hot core and follows the distribution of the two velocity components.

• •

We produced synthetic observations from numerical simulations of massive star-forming regions dominated by gravity. As in the SMA observations, the numerical simulations produced two distinctive velocity components (traced by ARTIST-generated H13CO+ emission) converging towards the strongest peak of the dust thermal emission, with the magnetic field being dragged by the gas. The polarization angular dispersion function comparison between the observations and simulations also reveals a very similar behavior.

• •

NGC 6334 V may be more dynamically evolved than the sources I, I(N) and IV (see Li et al., 2015), with gravity already being the dominant force over the magnetic field.

Finally, these results show how the gas is being accreted from the larger scale extended dense core ($\sim 0.1$ pc) of NGC 6334 V towards the higher density hot core region at $\sim 0.02$ pc scales, through two distinctive converging flows dragging along the magnetic field whose strength seems to have been overcome by gravity.

We would like to thank the referee for her/his useful comments. We also thank the SMA staff for their support which makes these studies possible. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics, and is funded by the Smithsonian Institution and the Academia Sinica. We thank Mauricio Tapia for thoughtful discussions about the origin of the infrared emission in NGC 6334 V. We thank Marco Padovani for his valuable support with ARTIST, and Gilberto Zavala Pérez and Alfonso H. Ginori González for their computational support. CJ acknowledges support from MINECO (Spain) BES-2012-052481 grant. CJ and JMG acknowledge support from MICINN (Spain) AYA2014-57369-C3 grant. JMG also acknowledges the support from the MECD (Spain) PRX15/00435 travel grant. AP and MZA acknowledge financial support from UNAM-DGAPA-PAPIIT IA 102815 grant (Mexico). MZA also acknowledges CONACyT for a postdoctoral fellowship at the University of Michigan. PMK acknowledges support from the Ministry of Science and Technology in Taiwan through grant MoST 103-2119-M-001-009 and an Academia Sinica Career Development Award. J.B.P. acknowledges UNAM-PAPIIT grant number IN110816, and to UNAM’s DGAPA-PASPA Sabbatical program. He also is indebted to the Alexander von Humboldt Stiftung for its invaluable support.