Magnetic Fields in Massive Star-forming Regions (MagMaR). II. Tomography through Dust and Molecular Line Polarization in NGC 6334I(N)

Figure 1 shows the total intensity (Stokes I) thermal dust continuum emission map from NGC 6334I(N) along with the main sources identified by Hunter et al. (2014) from their Sub Millimeter Array (SMA) data. The total intensity dust emission map shows an elongated filament with two cores (1b and 1c) dominating the emission and what appears to be a cavity in the dust emission toward the southern part of the filament. In this work we focus on the magnetic field, leaving the mass of the cores and density statistical analysis for further work (Cortes et al. 2021, in preparation). Figure 2 shows the magnetic field morphology on the plane of the sky as derived from polarized dust emission. The magnetic field morphology is derived by assuming grain alignment by magnetic fields, where the polarization position angles are rotated by 90° to obtain the field direction. The field pattern covers most of the NGC 6334I(N) filament, showing a clear indication of an "hourglass" shape over the 1b and 1c cores (inside the purple oval in Figure 2). Additionally, to the south of the main two cores we see that the field is also pinched over the third brightest core in the region (1a), with the field smoothly connecting to the aforementioned "hourglass" component. A pinched field morphology in NGC 6334I(N) has been suggested from cloud to core scales by Li et al. (2015). In fact, Li et al. traced an "hourglass" morphology with the SMA at core scales, which we reproduced here in Figure 3 by using data from Zhang et al. (2014). In this work we are further tracing the magnetic field morphology with ALMA from core to envelope scales. Note, we are referring to a pinched morphology instead of "hourglass" for the whole set of scales; we will discuss this in Section 4.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Recently, observations with the James Clerk Maxwell Telescope (JCMT) of polarized dust emission at 850 μm revealed a detailed field morphology of NGC 6334I(N) at clump scales (14'' resolution; Arzoumanian et al. 2021). To compare with ours, we show a zoomed map of the JCMT data (see Figure 3, left panel), where the pinched morphology is seen over a broader region encompassing both NGC 6334I(N) and NGC 6334I sources (NGC 6334I is not covered by the ALMA data presented here) and outlined by translucent red lines. Over the NGC 6334I(N) filament, the JCMT data show a mostly uniform pattern covering the region mapped by ALMA (see the yellow circle in Figure 3). Furthermore, Arzoumanian et al. compared the JCMT data to Planck data, where pinching of the field is only seen at the northeast edge of the cloud at the scales traced by Planck. Although these new data show a slightly different scenario from the one proposed by Li et al. (2015), the magnetic field appears to evolve coherently from clump to core-envelope scales (see Section 4.2 for a discussion).

To the SW of the filament, a cavity is seen in the dust emission traced by ALMA (see Figure 2). This cavity is well encircled by the field, which covers most of its perimeter. Magnetic fields along cavity walls produced by outflows have been seen in a number of low-mass star-forming regions (Hull et al. 2017; Maury et al. 2018; Le Gouellec et al. 2019; Hull et al. 2020a). Coincidentally, the blue lobe of the CS outflow, previously discovered by McCutcheon et al. (2000) and also reported here, appears to be cospatial with the cavity.

To ascertain the importance of the magnetic field in NGC 6334I(N), we estimate its strength onto the plane of the sky component, B_pos, using a number of variants of the Davis, Chandrasekhar, and Fermi method (or DCF; Davis 1951; Chandrasekhar & Fermi 1953; Heitsch et al. 2001; Falceta-Gonçalves et al. 2008) by using the angle dispersion function method (or ADF; Hildebrand et al. 2009; Houde et al. 2009, 2013a, 2016) and by using a recently derived approach for DCF that considers magnetosonic perturbations instead of Alfvén waves (Skalidis & Tassis 2021). We will discuss the applicability of such methods to regions such as NGC 6334I(N) in Section 4.3. The computations were executed by following Cortes et al. (2019) for the DCF ²⁰ variants and following Houde et al. (2009, 2016) for the dispersion function analysis. We obtained field strength estimates by considering only the emission within the purple ellipse shown in Figure 2 (see Table 3 for the results). This is justified because it is within this region that we have obtained sufficient overlap between the polarized dust and CS emission tracing the hourglass shape of the magnetic field (see Section 3.2). The B_pos estimates range between 1.4 and 23.6 mG, with an average of $\left\langle {B}_{\mathrm{pos}}\right\rangle =16$ mG. In contrast with previous works, here we estimate the field strength in a self-consistent manner by using parameter values derived directly from our data. For instance, we derive the velocity dispersion from our C³³S spectrum, which, by being optically thin (see Section 3.2), traces the turbulent motions inside the region. The column and volume densities are also derived directly from the Stokes I dust and CS emissions; we compute all values within the same region used to derive the polarization position angle dispersion (δ ϕ). To derive column density from dust emission, we followed the standard approach (Hildebrand 1983) assuming a dust opacity of k_1.3mm = 0.01 cm² g⁻¹ (Ossenkopf & Henning 1994), which assumes a gas to dust mass ratio of 100:1 and an average dust temperature of T_dust = 50 K (Sadaghiani et al. 2020). We estimate the volume density by assuming a cylindrical ellipsoid as the geometrical shape of the purple ellipse in Figure 2, where the height of the cylinder is taken as the mean of the major and minor axes of the ellipse, which equals 36, or 22 mpc, at the distance of NGC 6334I(N). We also use a mean molecular weight of μ = 2.8, which assumes that the gas has a 70% H₂ content and it is not stratified (Kirk et al. 2013). As a result, the uncertainties in the magnetic field strength estimation result primarily from the assumptions behind the validity of the DCF method and the geometrical assumption used to compute the density. The dispersion angle analysis is discussed in Section 4.3.

We detect CS (J = 5 → 4) emission, in total intensity, along the NGC 6334I(N) filament across a velocity range from –40 to 20 km s⁻¹. Here, we focus primarily on the polarization properties of the emission and leave a detailed analysis of the gas kinematics to future work (C. Cortes et al. 2021 in preparation). We analyze the CS (J = 5 → 4) and C³³S (J =5 → 4) spectra from the same region used to estimate the magnetic field strength from polarized dust emission (see Figure 4 for the spectra). Because C³³S is an isotopologue of the CS molecule, we can assume that both species are cospatially located and thus we can use the C³³S emission, likely optically thin, to estimate the properties of the CS gas. To estimate the column density and optical depth of the lines, we used the MADCUBA software package to model the CS and C³³S line profiles (Martín et al. 2019); the model is shown in Figure 4, right panel. Under local thermodynamic equilibrium (LTE) conditions, the three spectral features of the classic asymmetric top methyl formate (CH₃OCHO) detected close to the C³³S line (see Figure 4) allow us to put a good constraint on the excitation temperature, which was found to be T_ex = 220 ± 80 K. Note, this temperature is likely probing a gas kinetic temperature, which is significantly higher than the assumed dust temperature, where the CS emission probably arises from a cooler layer than CH₃OCHO. Although this temperature seems high, it is not uncommon to find such excitation temperatures when a hot molecular core (HMC) has developed, which is the case in NGC 6334I(N) 1b (Hunter et al.2014). There, complex organic molecules such as methyl formate act as excellent thermometers for the gas temperature. For instance, in W43-Main MM1, Sridharan et al. (2014) found an excitation temperature close to 400 K when considering spectral features of methyl cyanide in their data. Assuming that the C³³S emission is thermalized to this temperature and a source size of 1'' as derived from the fit to the integrated emission, we obtain a C³³S total column density ${{\rm{N}}}_{{C}^{33}S}\,=1.0\pm 0.1\times {10}^{15}$ cm⁻², with a peak optical depth of ${\tau }_{{{\rm{C}}}^{33}{\rm{S}}}\,=0.18\pm 0.02$. Scaling the C³³S emission using the sulfur 32/33 relative abundance ratio reported by Chin et al. (1996) yields a CS total column density of N_CS = 2.9 ± 1.1 × 10¹⁷ cm⁻² with a peak optical depth of τ_CS = 32 ± 12, and therefore strongly optically thick. For completeness, we modeled the CS emission with the parameters above, together with a foreground component under T = 50 K that absorbs both the background line and continuum emission. We note the good agreement between the velocity fit of CS, C³³S, and CH₃OCHO of −2.46 ± 0.15 km s⁻¹, −2.3 ± 0.16 km s⁻¹, and 2.0 ± 0.3 km s⁻¹, with the absorption layer slightly blueshifted to −3.4 ± 0.2 km s⁻¹. The errors are obtained from the model fit to the lines.

Zoom In Zoom Out Reset image size

We detect polarized emission from the CS (J = 5 → 4) molecular line toward NGC 6334I(N). The bottom-left panel of Figure 4 shows the polarized intensity spectrum and Figure 5 shows the channel maps. In the channel maps, we display the polarized CS emission as orange pseudo-vectors superposed on a coarsely plotted magnetic field morphology as inferred from polarized dust emission (shown as blue pseudo-vectors). Note, we are not applying a 90° rotation to the CS pseudo-vectors as we do to the polarized dust emission. The channel maps reveal good agreement between the hourglass magnetic field morphology seen in the dust and the polarized CS emission observed as a function of velocity. For a more quantitative comparison between the field morphology and the polarized CS emission, we compute histograms of the differences in polarization position angles between the CS and the dust emission for different velocity channels (see Figure 6). The differences are calculated only at locations where the polarized CS emission overlaps with the polarized dust emission. Furthermore, we fit Gaussian profiles to histograms to derive probability density functions. All of the histograms are well centered around zero, within ±10° for the ±4 km s⁻¹ range, with mostly symmetric Gaussian distributions suggesting that the CS polarized emission is correctly tracing the "hourglass" field morphology derived from dust (see Table 2 for the statistics). The best match is found at V = −2 km s⁻¹ where $\left\langle {\rm{\Delta }}\theta \right\rangle =0^\circ$ and σ_Δθ = 6° from the Gaussian fit. The V_lsr of NGC 6334I(N) is about ∼−3 km s⁻¹ and thus the best match is close to the systemic velocity of the source. Although there are deviations between the magnetic field and the CS position angles, these are observed at the line wings of the line, suggesting a departure from the "hourglass" where the emission from the outflow becomes dominant. The polarized CS emission also seems to trace the magnetic field morphology along the dust cavity: see channel maps v = −10 to −2 km s⁻¹, which also show the blueshifted lobe of the CS outflow. We note that the CS outflow appears to be orthogonal to the symmetry axis of the "hourglass" magnetic field. We will explore this finding in upcoming work (Cortes et al. 2021, in preparation). Additionally, we present channel maps of polarized emission from C³³S(J = 5 → 4) in Figure 7. Although the number of independent detections of polarization in the C³³S maps is smaller than for CS, the agreement with the CS polarization position angle appears consistent across velocity space.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Polarization Angle Statistics

Velocity	$\left\langle {\rm{\Delta }}\phi \right\rangle$ a	std(Δϕ) b	Δϕ0 c	σΔϕ d
(km s−1)	(°)	(°)	(°)	(°)
−12.0	8	25.9	21	−11.5
−10.0	9	20.1	16	17.4
−8.0	10	18.0	11	−17.4
−6.0	4	13.8	4	-13.1
−4.0	2	9.0	1	−8.4
−2.0	1	7.5	0	6.4
0.0	3	12.3	2	11.5
2.0	0	14.4	−2	7.6
4.0	−2	14.8	−5	−10.1
6.0	−2	15.6	−5	−7.4
8.0	−8	19.0	−6	6.2
10.0	−8	17.4	−8	6.9
12.0	−2	13.5	−4	8.7
14.0	−2	12.2	−7	8.0

Notes.

^a Here $\left\langle {\rm{\Delta }}\phi \right\rangle$ represents the mean of the difference values between dust and CS polarization angles. ^b Here std(Δϕ) corresponds to the standard deviation of the difference values between the dust and the CS polarization angles. ^c Here Δϕ₀ corresponds to the center of a Gaussian fit to the distribution of polarization angle differences. ^d Here σ_Δϕ corresponds to the width of a Gaussian fit to the distribution of polarization angle differences.

Download table as:  ASCIITypeset image

The number of independent polarization detections in the CS channel maps is large enough that we can also estimate the field strength using the DCF technique in a number of those channels (between −6 and −2 km s⁻¹; for a description of what can be considered to be a sufficient number of channels, see Appendix A in Cortes et al. 2019). Over this channel range, we obtain an average field strength estimate of ∼2 mG (see Table 3 for the channel-by-channel estimates). Polarized emission from molecular lines has a 90° ambiguity with respect to the ambient magnetic field direction (Goldreich & Kylafis 1981). However, this does not affect the estimation method because δϕ will not change if the data are rotated by 90°. We will address the relation between the CS polarization pattern and the magnetic field in Section 4.1.1.

Table 3. Magnetic Field Strength Estimations

Tracer	Velocity	Na	n	$\left\langle \phi \right\rangle$ b	δ ϕb	Δ V	B1 c	B2 d	B3 e	B4 f	B5 g	${\left\langle {B}_{t}\right\rangle }^{2}/{\left\langle B\right\rangle }^{2}$
	(km s−1)	(1024 cm−2)	(105 cm−3)	(°)	(°)	(km s−1)	(mG)	(mG)	(mG)	(mG)	(mG)
Dust	⋯	2.0	423.9	29.3	29.4	5.3	21.9	19.9	1.4	23.6	11.1	0.23 ± 0.01
CS	−6.0	2.0	2.0	179.3	6.7	5.3	3.24	3.18	0.060	2.75	1.58	0.08 ± 0.001
CS	−4.0	2.0	2.0	179.2	6.5	5.3	2.49	2.40	0.116	2.75	1.39	0.08 ± 0.001
CS	−2.0	2.0	2.0	179.3	6.2	5.3	2.39	2.31	0.165	2.75	1.36	0.08 ± 0.001

Notes.

^a The column density corresponds to the region used to extract δ ϕ. ^b $\left\langle \phi \right\rangle$ is the average polarization position angle (PA) and δ ϕ is the PA dispersion (calculated using circular statistics). ^c Estimations of the magnetic field, in the plane of the sky, done with the original DCF method. ^d Estimations of the magnetic field in the plane of the sky, done using the corrections implemented in Equation (9) in Falceta-Gonçalves et al. (2008). ^e Estimations of the magnetic field in the plane of the sky, done using the corrections implemented in Equation (12) in Heitsch et al. (2001). ^f Magnetic field strength estimated by using the ratio of turbulent to total magnetic energy. ^g Magnetic field strength estimated using the DCF modification proposed by Skalidis & Tassis (2021).

Download table as:  ASCIITypeset image

Linearly polarized emission from molecular lines was first detected in CS emission by Glenn et al. (1997) toward the IRC +10216 evolved star. Since then, it has been detected toward a number of sources, particularly in HMSFRs (Girart et al. 1999; Lai et al. 2003; Cortes et al. 2005; Beuther et al. 2010; Hirota et al. 2020), and toward evolved stars (Vlemmings et al. 2012; Girart et al. 2013). Linearly polarized emission from molecular lines is expected when the magnetic sublevels are unevenly populated because of anisotropies in the medium. Initial models assumed that large velocity gradients resulting from the kinematics were the dominant source of anisotropy (Goldreich & Kylafis 1982; Deguchi & Watson 1984). However, anisotropies in the radiation field can also affect the level populations and may be due to a geometrical distribution of the gas that produces optical depths that are not the same in all directions, or from embedded sources such as a protostellar core (in the case of star-forming regions; Cortes et al. 2005) or a star (in the case of a circumstellar shell around evolved sources; Vlemmings et al. 2012). Although it is expected that the amount of polarized emission will decrease with increasing optical depth as a result of photon trapping, we nevertheless find significant amounts of polarized emission across the CS spectrum in NGC 6334I(N) (between 2% and ∼10% within −20 to 20 km s⁻¹; see the bottom-left panel of Figure 4), which our calculations show to be optically thick. Deguchi & Watson (1984) performed multilevel calculations for the CS (J = 1 → 0) and (J = 2 → 1) transitions and found linearly polarized emission at a level of ∼1% for τ ∼ 10 with a significant decrease in fractional polarization for increasing τ. However, this computation is only one-dimensional and represents the best-case scenario where the optical depth is taken along the velocity gradient (see Figure 3 in Deguchi & Watson 1984). Cortes et al. (2005) also performed multilevel radiative transfer calculations for both CO (J = 1 → 0) and (J = 2 → 1) transitions, but this time adding a blackbody to the computation in order to introduce anisotropies in the radiation field. They found an increasing amount of polarization as the line becomes optically thin (see Figure 7 in Cortes et al. 2005). Recently, Lankhaar & Vlemmings (2020) presented a comprehensive quantum mechanical treatment of the alignment of the molecule's angular momentum by assuming only an anisotropic radiation field. Their treatment is fully three dimensional, allowing for simulations of both radiative transfer and gas dynamics. Although they did not model CS emission in particular, their simple example of a collapsing spherical cloud produced polarization fractions for HCO⁺ (J = 3 → 2) and (J = 2 → 1) at levels over 1% for radial distances beyond 600 au (∼3 mpc) for the (J = 3 → 2) and 900 au (∼4.5 mpc) for the (J = 2 → 1) transition. However, all of these calculations find that polarized line emission significantly decreases with higher optical depths. A possible explanation for this may be related to interferometric filtering and missing flux. Fractional polarization with interferometers has to be analyzed with care because of spatial filtering effects (Le Gouellec et al. 2020). The maximum recoverable angular scale for the configuration used to acquire these data is ∼5'', which suggests that a significant fraction of the extended total intensity CS emission seen from the single dish (McCutcheon et al. 2000) might be filtered out by ALMA, and thus the fractional polarization values presented here are overestimated. This is because linearly polarized emission tends to be more compact than the total intensity emission. Estimating how far we are from the true fractional polarization will require sampling larger angular scales in full polarization mode, which the ALMA compact array (ACA) can do.

We also see significant amounts of fractional polarization at higher velocities (polarization levels ≳10% at v < −20 and v > 20 km s⁻¹). It is likely that in these velocity ranges where the CS emission is tracing the outflow, the CS emission might be optically thin and thus the values estimated might be closer to the "true" fractional polarization values. The large velocity gradient and the radiation field from the embedded protostars are, most likely, the sources of anisotropy necessary to produce the large polarization fractions that we see there (see Cortes et al. 2005, Figure 7, for a CO computation). We also note that strong polarization from CS is expected because of its large dipole moment (μ= 1.96 D) compared with, for example, that of CO (μ = 0.12 D), because of the μ² dependence of the radiative rates. Thus, it seems that the combination of anisotropies in both the radiation field and the gas kinematics is causing the strong CS polarized emission that we see in NGC 6334I(N), where the filtering effect might explain the large fractional polarization values seen when comparing to radiative transfer calculations. Nonetheless, three-dimensional modeling of polarized CS emission will be needed to further understand the emission detected here, which is beyond the scope of this work.

Finally, if a foreground screen of molecular gas is present between the source and the telescope, some amount of linear polarization might be converted into circular polarization through anisotropic resonant scattering, or ARS (Houde et al.2013b). The ARS will systematically corrupt the linear polarization position angle from the CS emission and change its relation with respect to the ambient magnetic field. The ARS will manifest itself through the presence of a statistically significant signal in Stokes V, which we do not detect. The peak emission is Stokes V is between about −2.1 and 3.3 mJy beam⁻¹ with a fractional level between 0.6% and 1.0%, which is below the statistical uncertainty of ∼2% that ALMA can measure (Cortes et al. 2021). The morphology of the Stokes V velocity channel maps is consistent with noise where the peaks alternate between positive and negative between channels, which is inconsistent with a coherent conversion of linear to circular polarization by ARS (see Table 1 and Figure 9 in Appendix A for the Stokes V channel maps). Thus, we conclude that polarized emission from CS is purely linear and its relation to the field is subject to the known 90° ambiguity (see the next section for a discussion).

4.1.1. The 90° Ambiguity in the Position Angle of the Polarized CS Emission

Hourglass magnetic field morphologies have long been predicted by magnetically regulated star formation models (e.g., Mouschovias 1976; Mouschovias & Morton 1985) and may be the result of ambipolar diffusion (Mouschovias 1991). This field morphology has been seen in a number of low- and high-mass star-forming regions and at different length and mass scales (Schleuning 1998; Girart et al. 2006, 2009; Maury et al. 2018; Beltrán et al. 2019). In the magnetically regulated star formation scenario, a molecular cloud will initially be supported against gravitational collapse by the magnetic field, which is initially thought to be uniform. Because the neutrals are only weakly coupled to the charge carriers, which are held in place by the field, they will slowly diffuse past the magnetic field via ambipolar diffusion. Because the field is frozen into the charge carriers, the field will be pinched, attaining the hourglass shape. Alternatively, a self-gravitating core can also pull the field at the core center, which will create an "hourglass" shape, but in this case ambipolar diffusion might be a by-product of the process.

Although the "hourglass" shape seen in NGC 6334I(N) is clear, the exquisite sensitivity and resolution of the ALMA data allow us to see local deviations from which the "hourglass" is a simple model. The density regime traced by ALMA is likely showing the result of complex physics that is difficult to explain without detailed numerical modeling and outside the scope of this paper. Furthermore, a clear "hourglass" field morphology in NGC 6334I(N) is not seen at all length scales. From the envelope scales traced by ALMA to the core scales imaged by the SMA, and to the clump scales traced by the JCMT, the hourglass appears to be the shape of the magnetic field when we also consider NGC 6334I (see Figures 2 and 3). However, at the cloud scales traced by Planck (see Figure 5 in Arzoumanian et al. 2021), the field is only pinched at the northeast side of the cloud. Whether this is a projection effect or not is not clear from the figure alone, where a detailed analysis of the Planck data would be required. Even though we cannot state with certainty that we see the "hourglass" field morphology as the ubiquitous shape across all length scales in NGC 6334, the pinching in the field seen from the cloud to the envelope scales is quite remarkable, suggesting that the magnetic field is strong in this region. Furthermore, the field pattern also seems preserved in velocity space from –10 to 4 km s⁻¹ at the ALMA scales, which is a 14 km s⁻¹ range, almost three times the 5.3 km s⁻¹ FWHM line width of C³³S (our proxy for the turbulent motions).

Besides tracing a pinched field shape through multiple orders of magnitude in spatial scales, we are also tracing the field at different densities as well (from 10 to 10⁷ cm⁻³ when considering density estimates from Li et al. 2015; Arzoumanian et al. 2021). The critical density for the CS (J = 5 → 4) transition is 9 × 10⁶ cm⁻³, which is obtained by assuming that the emission is optically thin. Because of the high optical depth estimated for the CS (J = 5 → 4) transition, that number density is most likely smaller as the line may be subthermally excited. Shirley (2015) accounted for optical depth effects by computing the critical density considering photon trapping. Thus, the number density can be approximated as ${n}_{\mathrm{crit}}^{\mathrm{thick}}={n}_{\mathrm{crit}}^{\mathrm{thin}}/{\tau }_{{\nu }_{\mathrm{jk}}}$ when the optical depth is much larger than one, where ${n}_{{\rm{crit}}}^{{\rm{thick}}}$ is the optically thin critical density and τ_νjk is the line optical depth. From this assumption, we obtain ${n}_{\mathrm{CS},\mathrm{crit}}^{\mathrm{thick}}=2.8\times {10}^{5}$ cm⁻³. Note, here the ${n}_{\mathrm{CS},\mathrm{crit}}^{\mathrm{thick}}$ refers to the collisional partners involved in the excitation of the CS molecule, which corresponds to H₂. Although an abuse of notation, we use CS in the subscript to indicate that the density is derived from the CS emission. Furthermore, Shirley (2015) tabulated effective excitation densities for the CS molecule by considering a number of rotational transitions. The effective excitation density is an empirical quantity defined by considering a 1 K km s⁻¹ integrated emission. This quantity also takes into account optical depth effects such as radiation trapping (for a review, see Evans 1999). For the CS (J = 5 → 4) transition, Shirley (2015) obtains n_eff = 7.6 × 10⁴ cm⁻³, assuming a kinetic temperature of 50 K, the dust temperature that is a good estimation of the kinetic temperature in the region of interest (see the purple oval in Figures 2 and 5). Although these two estimates of the number density are close in value, detailed numerical radiative transfer modeling is required to obtained more accurate population numbers, which we leave for future work. Nonetheless, the two criteria are reasonable approximations to the number density of the molecular hydrogen CS collisional partners under high optical depth conditions, and thus we use the average as the estimate of the "true" number density, or n_CS ∼ 2 × 10⁵ cm⁻³.

A value of n_CS = 2 × 10⁵ cm⁻³ is ∼2 orders of magnitude less than the volume density estimated from the dust emission, which we calculated to be 4.2 × 10⁷ cm⁻³. Thus, the polarized CS emission appears to be tracing the field at lower densities than the dust and at a level comparable with the JCMT observations. Therefore, within a single ALMA data set, we are not only tracing the "hourglass" shape as a function of velocity but also as a function of density: i.e., we are effectively performing magnetic field tomography. When we consider the evolution of the field morphology along these three axes (length scale, velocity, and density), the striking coherence seen in the field structure strongly suggests that the magnetic field remains dynamically important from the diffuse to the high-density regime in NGC 6334I(N).

The DCF method does not consider the effects of finite angular resolution or integration along the line of sight, which at lower resolutions smooth out the polarized emission, reducing its dispersion, and therefore it might overestimate the field strength. Moreover, the effect of self-gravity in bending the field lines, which will affect the dispersion in the polarization position angle, is also not considered by DCF and its variants. Thus and because of the exquisite uv-coverage, resolution, and sensitivity of ALMA, local deviations in the position angle from the main field model, due to gravity, make the applicability of the DCF method more challenging to the data. Here, we discuss how the different DCF variants used in this work attempt to correct the polarization position angle dispersion. To account for poor resolution in polarization maps, Heitsch et al. (2001) used the geometric mean between two modified DCF equations. The first modification attempts to address the small-angle approximation by replacing the polarization angle dispersion by the dispersion of the tangent values of the position angle. The second modification attempts to deal with the case where the dispersion in the field lines is larger than the mean field. They do this by considering the three-dimensional expansion of the field where all of the random components are assumed to be the same. This yields an equation that is also dependent on the dispersion of the tangent of the position angle values (see Equation (11) in Heitsch et al. 2001). Although consistent with their own simulations, their variant appears to underestimate the field strength when compared to other numerical results (Falceta-Gonçalves et al.2008). In our data, this method yielded an estimate that is a factor of ∼10 smaller than the other estimates. It is likely that this is because local deviations from the "hourglass" main field morphology produced larger values when using the tangent function, which yields a larger dispersion and therefore a smaller estimate. In contrast, the modification proposed by Falceta-Gonçalves et al. (2008) was implemented by assuming that δ B/B is a global relation and thus they modified DCF by taking the tangent of the dispersion angle instead of δ ϕ (see Equation (9) in Falceta-Gonçalves et al. 2008). This was proposed to address larger angle dispersions due to an increasing turbulent component in the field. However, this will rapidly decrease the field strength as the tangent function quickly diverges for δ ϕ > 60°. Recently, the DCF method was revisited and another correction that considers the compressible modes from small amplitude MHD waves (magnetosonic) was proposed, instead of purely Alfvén transverse waves (Skalidis & Tassis 2021). This method implements a substitution in the DCF equation in the form of $\delta \phi \to \sqrt{\delta \phi }$, which seems to improve the recovery of the field strength when applied to their numerical simulations. Although it is not clear if the dominant mode perturbing the main field component is transverse (Alfvénic) or compressible (small amplitude magnetosonic waves), this variant will yield estimates that are smaller in value than the original DCF. The three previously described DCF variants used simulations to test the proposed modifications. All of the simulations assumed ideal MHD, which might not be representative of the physical conditions in NGC 6334I(N). In fact, ideal MHD may produce artificially tangled magnetic field morphologies that will affect the polarization position angle dispersion, particularly in the line of sight. Furthermore and because of the low ionization rates in dense molecular clumps (∼10⁻⁷), nonideal MHD effects, such as ambipolar diffusion, are unavoidable (Hennebelle & Inutsuka2019).

The ADF method provides a way to quantify the turbulent component in the field and to better estimate the value of δ ϕ. This is done by fitting a structure function of the polarization position angle data, by incorporating a turbulence model, and also by considering the effects of interferometric filtering (for a summary of the technique see Section 2 of Houde et al. 2016). This analysis allows us to associate $\delta \phi ={\left[\left\langle {B}_{t}^{2}\right\rangle /\left\langle {B}^{2}\right\rangle \right]}^{1/2}$, where the quantity on the right-hand side is the ratio of turbulent to total magnetic energy, one of the quantities derived from the dispersion function analysis. We can then use δ ϕ to estimate the field strength in the plane of the sky via the usual DCF technique, which takes the form of

$$\begin{eqnarray}&&{B}_{\mathrm{pos}}\simeq \sqrt{4\pi \rho }\sigma (v){\left[\displaystyle \frac{\left\langle {B}_{t}^{2}\right\rangle }{\left\langle {B}^{2}\right\rangle }\right]}^{-1/2},\end{eqnarray} \tag{ 1 }$$

where $\sigma (v)={\rm{\Delta }}V/2\sqrt{2\mathrm{log}2}$, ΔV is the FWHM line width of the C³³S line, and ρ is the volume density. By applying this analysis to the polarized dust emission data, considering the same region as before (see purple ellipse in Figure 2), we obtain a turbulence correlation length of δ = 0262 ± 0008, or ∼2 mpc, at the distance to NGC 6334I(N), and $\left\langle {B}_{t}^{2}\right\rangle /\left\langle {B}^{2}\right\rangle =0.23\pm 0.01$, which we use to estimate a plane-of-sky magnetic field strength B_pos = 24 mG (see Figure 8 and Table 3 for the results). Furthermore, we apply the same analysis to the CS polarized emission data considering a wide velocity range and obtain a good fit between −6 and −2 km s⁻¹. In this interval, we find a turbulence correlation length of δ = 042 ± 0011, or ∼2.6 mpc, and $\left\langle {B}_{t}^{2}\right\rangle /\left\langle {B}^{2}\right\rangle =0.08\pm 0.0013$, which yields a plane-of-sky magnetic field strength B_pos = 2.8 mG. Note, the analysis was applied to the integrated data between −6 and −2 km s⁻¹, but given the small differences in δ ϕ within the interval, we used the same $\left\langle {B}_{t}^{2}\right\rangle /\left\langle {B}^{2}\right\rangle$ value to estimate B_pos at each velocity channel (see Table 3). The differences in the field strength estimation are largely due to the differences in density used. That is, the dispersion analyses only contribute a factor of ∼1.7 (lowering the dust estimate) while the densities increase the dust value by ∼14.

Zoom In Zoom Out Reset image size

Recent work by Liu et al. (2021) analyzed ideal MHD simulations of protocluster formation at clump scales. They applied various statistical methods to synthetic magnetic field maps to study the applicability of the DCF method and the variants used here (including ADF). Because the magnetic morphology in NGC 6334I(N) has an "hourglass" shape, it is likely that we are in the strong-field regime. Thus, the Liu et al. results suggest that the magnetic field strength estimates derived here are good to a factor of a few, again, subject to the caveat of ideal MHD simulations.

We have five estimates for the field strength onto the plane of the sky. Each of them originates from modifications to the DCF method that try to address finite resolutions, the polarization angle dispersion value due to a random component of the field, and the effect of a different perturbation mode. Because none of these methods used here consider all of the relevant physics in this region, e.g., self-gravity, the "true" value for the field strength remains unconstrained. We lack an actual measurement of the field strength such as the one provided by the Zeeman effect. Although still contested (see Jiang et al. 2020), results from Zeeman measurements show that the field strength will grow with density as a power law, or ∼ n^2/3 (Crutcher & Kemball 2019). Thus, the DCF and its variants still give us a first-order statistical approximation to the "true" magnetic field strength onto the plane of the sky. We quantify a final estimate by taking the average of all five estimates, obtaining $\left\langle {B}_{\mathrm{pos}}\right\rangle =16$ mG as the average magnetic field on the plane of the sky at densities of n = 4.2 × 10⁷ cm⁻³ and $\left\langle {B}_{\mathrm{pos}}\right\rangle =2$ mG at densities of n = 2.0 × 10⁵ cm⁻³ when considering the polarized CS emission.

By comparing our results from NGC 6334I(N) to other HMSFRs, we seek to discover if there is a pattern in the physical conditions of regions where the magnetic field has a clear and distinctive shape, such as an "hourglass" morphology. Examples of similar magnetic field morphologies to NGC 6334I(N) are cores like G240, where the "hourglass" magnetic field appears as a "textbook" case for magnetically controlled star formation with a bipolar outflow closely aligned to both the rotation and magnetic field axes (Qiu et al. 2014). Note that the G240 mass, 95 M_⊙, is substantially larger than the ∼26 M_⊙ of the combined 1b and 1c core masses in NGC 6334I(N) as derived from our data inside the purple oval region shown in Figure 2. Another example is the massive core G31.41, with a mass comparable to NGC 6334I(N) (about 26 M_⊙ from Beltrán et al. 2019), has also been shown to exhibit an "hourglass" magnetic field morphology from the core (Girart et al. 2009) to envelope scales (Beltrán et al. 2019). However, the alignment between the outflow, rotation, and magnetic axes is less clear here when compared to G240. The length scales where the "hourglass" shape is traced in these two sources are similar to what we see in NGC 6334I(N). For instance, in G31.41 the field morphology is seen preserved through a scale range that matches the lower end in the NGC 6334I(N) scales. The data obtained from G240 trace the field at the core scales also where the "hourglass" is seen in NGC 6334I(N). However, not all HMSFRs show "hourglass" magnetic field morphologies. For instance, in the W43-Main molecular complex, the W43-MM1 (Cortes et al. 2016; Arce-Tord et al. 2020) and W43-MM2 clumps (Cortes et al. 2019) exhibit magnetic field morphologies that are primarily radial over their most massive cores, which is expected when gravity dominates the dynamics. W43-Main harbors some of the most massive protostars currently known (≳100 M_⊙, Cortes et al. 2016; Motte et al. 2018; Cortes et al. 2019), whereas the mass of the central cores in NGC 6334I(N) are only about ∼26 M_⊙ in total when considering our data. However and because of the angular scales sampled by our ALMA data, we might be missing flux that might make NGC 6334I(N) appear less massive than other regions. Nonetheless, this difference in core mass is not significant when comparing W43-Main with G240 where the core masses are comparable, but the field shapes are completely different. A totally different magnetic field morphology is seen in IRAS 180089–1732 where the field was found to have a spiral morphology (Sanhueza et al. 2021). Previous mapping of this source at clump scales appears to show the same field pattern (Beuther et al. 2010) as seen by ALMA at envelope scales. In this case, the total core mass is estimated to be 75 M_⊙ from the ALMA data, which is also comparable to G240 and in the lower range from the W43-Main estimates. Thus, it is also uncertain whether the core mass is a decisive factor to explain the differences in the field shape seen across these HMSFRs.

High-mass star-forming cores are usually surrounded by H ii regions, which provide significant radiative feedback. It is possible that radiation pressure coming from H ii regions may compress the field in conjunction with the effects of gravity, which the field may resist if strong enough (e.g., see Li et al. 2006; Shariff et al. 2019, for an example in the Carina nebula). For instance, W43-Main is part of a giant molecular complex, which has at its center a large H ii region powered by a number of O7 Wolf-Rayet stars, which appear to be not only ionizing the boundaries of W43-Main but also compressing the gas (Blum et al. 1999; Motte et al. 2003), while NGC 6334 contains a group of smaller H ii regions known as the "Cat's Paw," which seem distributed along the filament (Russeil et al. 2016). This also seems to be the case for G31.41, which is surrounded by both compact and extended H ii regions (J. M. Girart 2021, private communication). In contrast, for G240 and IRAS 180089–1732, the situation seems unclear as the cores appear to be more isolated than those in NGC 6334I(N), G31.41, and W43-Main. Although we note these differences, in this simple analysis we are certainly ignoring a number of other factors such as chemical diversity, stage of evolution, possible initial conditions, among many others. Thus, acquiring sufficient statistical cases is paramount to increasing our understanding of how stars form in HMSFRs and what is the role of the magnetic field. As part of this MagMaR project, we have acquired a comprehensive sample that is sufficiently large to allow us to begin addressing these questions in future work.