Magnetic Field Alignment Relative to Multiple Tracers in the High-mass Star-forming Region RCW 36

The linearly polarized intensity P can be found from the linear polarization Stokes parameters Q and U using

$$\begin{eqnarray}&&P=\sqrt{{Q}^{2}+{U}^{2}},\end{eqnarray} \tag{ 1 }$$

while the polarization fraction is given by p = P/I, where I is the total intensity. The polarized intensity and polarization fraction are both constrained to be positive quantities, which results in a positive bias at low total intensities (S. O. Rice 1945; K. SERKOWSKI 1962). This can be corrected for with "debiasing" (J. F. C. Wardle & P. P. Kronberg 1974), using

$$\begin{eqnarray}&&{P}_{{db}}=\sqrt{{Q}^{2}+{U}^{2}-\displaystyle \frac{1}{2}(\delta {Q}^{2}+\delta {U}^{2})},\end{eqnarray} \tag{ 2 }$$

where δ Q and δ U are the measurement uncertainties in Q and U, respectively. The polarization angle $\hat{E}$ can be calculated using $\hat{E}=0.5\,\arctan (U,Q)$. The orientation of the plane-of-sky magnetic field ${\hat{B}}_{\perp }$ is then inferred to be orthogonal to $\hat{E}$ such that

$$\begin{eqnarray}&&{\hat{B}}_{\perp }=\displaystyle \frac{1}{2}\arctan (U,Q)+\displaystyle \frac{\pi }{2},\end{eqnarray} \tag{ 3 }$$

where an angle of 0° points toward Equatorial north and increases eastward.

2.1.1. SOFIA HAWC+

In this paper, we publish for the first time, observations of RCW 36 using publicly available archival data from HAWC+ (D. A. Harper et al. 2018), the far-IR polarimeter on board SOFIA. The RCW 36 region was observed by SOFIA/HAWC+ on 2018 June 6 and 14 as part of the Guaranteed Time Observing program (AOR: 70_0609_12), in both Band C (89 μm) and Band E (214 μm), at nominal angular resolutions at FWHM of 78 and 182, respectively. The observations were done using the matched-chop-nod method (described in R. H. Hildebrand et al. 2000) with a chopping frequency of 10.2 Hz, chop angle of 1124, nod angle of −675, and chop throw of 240''. Each observing block consisted of four dithered positions, displaced by 12'' for Band C and 27'' for Band E. The total observation times for Band C and Band E were 2845 and 947 s, respectively. The total intensity maps for each band are shown in the top row of Figure 1.

To reduce this data, we used the HAWC+ Data Reduction Pipeline, which is described in detail in F. P. Santos et al. (2019) and D. Lee et al. (2021) and summarized here as follows. The pipeline begins by demodulating the data and discarding any data points affected by erroneous telescope movements or other data acquisition errors. This demodulated data are then flat-fielded to calibrate for gain fluctuations between pixels and combined into four sky images per independent pointing. The final Stokes I, Q, and U maps are generated from these four maps after performing flux calibrations accounting for the atmospheric opacity and pointing offsets. Next, the polarization is debiased (see Equation (2)), and the polarization percentage and the polarization angle are calculated. A χ² statistic is then computed by comparing the consistency between repeated measurements to estimate additional sources of uncertainties such as noise correlated across pixels, which can lead to an underestimation of errors in the I, Q, U maps (J. A. Davidson et al. 2011; N. L. Chapman et al. 2013). This underestimation can be corrected for using an excess noise factor (ENF) given by $\mathrm{ENF}=\sqrt{{\chi }^{2}/{\chi }_{\mathrm{theo}}}$ where χ_theo is theoretically expected, and χ² is measured. The ENF is estimated in the HAWC+ Reduction Pipeline by fitting two parameters I₀ and C₀ using

$$\begin{eqnarray}&&\mathrm{ENF}={C}_{0}\sqrt{1+{\left(\displaystyle \frac{I}{5{I}_{0}}\right)}^{2}},\end{eqnarray} \tag{ 4 }$$

where I is the total intensity (Stokes I) of a pixel in units of Jy pixel⁻¹. The errors of the I, Q, U maps are then multiplied by the ENF. The values of I₀, C₀ used for I, Q, and U maps in Band C and E are summarized in Table 2. We note that the Stokes I errors for Band E were a special case where the pipeline ENF fitting routine failed (giving a value of I₀ = 0, which resulted in a nonphysical ENF), likely due to large intensities from bright emission. In this case, we forced the ENF to be about 1 (by manually setting I₀ to be a sufficiently large number such as 100 and C₀ = 1) such that the errors were neither underestimated nor overestimated.

Table 2. Excess Noise Factor (ENF) Fitting Parameters I₀ and C₀ Found by the HAWC+ Data Reduction Pipeline for Stokes I, Q, U Maps in Band C and E (See Equation (4))

Band	Stokes	I0	C0
C 89 μm	I	0.332	25.929
	Q	0.716	1.723
	U	0.661	1.769
E 214 μm	I a	100	1
	Q	1.067	1.544
	U	0.696	1.739

Note.

^a Fitting routine failed, and values were manually set to get an ENF ∼1 (see text for details).

Download table as:  ASCIITypeset image

After correcting the errors, the pipeline rejects any measurements falling below the 3σ cutoff in the degree of polarization p to the associated uncertainty σ_p (p < 3σ_p ), which roughly corresponds to a 10° uncertainty in the polarization angle.

After running the pipeline, we also applied a 3σ signal-to-noise threshold on the total intensity flux and polarized flux to further remove noisy polarization vectors. As a final diagnostic, we checked for potential contamination from the reference beam position due to dithering the data in Chop-Nod polarization observations (following the method described in the Appendix section of G. Novak et al. 1997). For this test, we used Herschel far-IR intensity maps (described in Section2.2.2) since they cover both the RCW 36 region and the surrounding Vela C cloud-scale region, which include the HAWC+ chop reference beam positions that are outside the HAWC+ maps. We used Herschel maps of comparable wavelengths to the HAWC+ data (PACS 70 μm to compare with HAWC+ 89 μm and PACS 160 and SPIRE 250 μm to compare with HAWC+ 214 μm) and found the ratio of the total intensity of the HAWC+ region compared to the chop reference beam region in the Herschel data. To estimate the ratio of polarization flux, we conservatively assume that there is a 10% polarization at the reference beam positions and remove points where the estimated polarized flux in the reference beam is more than 1/3 times the polarized flux in the HAWC+ map.

2.2.1. ALMA

2.2.2. Herschel SPIRE and PACS

2.2.3. Spitzer IRAC

The gas structure of the region is also of significant interest in understanding the dynamic importance of the magnetic field and how it has been affected by stellar feedback. To this end, we use a myriad of archival spectroscopic line data to probe different chemical, thermal, and density conditions within the RCW 36 region. Table 3 summarizes the lines of interest, including their transitions, rest frequencies, velocity resolution, and the channels used to make the integrated intensity maps. For our analysis of the spectroscopic data, we use both an integrated intensity map for a wide velocity range (column (5) of Table 3) as well as channel maps over narrower velocity ranges (column (6) of Table 3). Our spectral data cube analysis is described in Sections 4.2.1 and 4.2.2.

2.3.1. SOFIA upGREAT Feedback

2.3.2. APEX LAsMA

2.3.3. Mopra

On cloud scales of >1 pc, the RCW 36 region is located within the Vela C giant molecular cloud. Vela C consists of a network of filaments, ridges, and nests, which were identified by T. Hill et al. (2011) using Herschel data. The densest and most prominent of the ridges is the Centre-Ridge, with column densities of A_V > 30, and a length of roughly ∼10 pc (T. Hill et al. 2011). The Centre-Ridge contains the RCW 36 region. It has a bipolar nebula morphology (V. Minier et al. 2013) with two fairly symmetric lobes oriented in the east–west direction that are traced well by the green PACS 70 μm emission in the lower left panel of Figure 3 (see also the dotted green ellipses). This bipolar nebula is roughly centered around a young (1.1 ± 0.6 My) massive cluster with two late-type O-type stars and ∼350 members (L. E. Ellerbroek et al. 2013). The position of the most massive star (spectral type O9 V) is indicated by a white star-shaped marker in Figure 3. The ionizing radiation from this cluster is powering an expanding H ii gas shell traced by Hα (shown in blue, Figure 3, left panel; V. Minier et al. 2013).

Bipolar H ii regions are of great interest because, although they have been observed in other high-mass star-forming regions such as S106 (N. Schneider et al. 2018) and G319.88+00.79 (L. Deharveng et al. 2015), they seem to be more rare than single H ii bubbles (e.g., E. Churchwell et al. 2006; L. Deharveng et al. 2010; L. D. Anderson et al. 2011; S. Kendrew et al. 2012; M. R. Samal et al. 2018).

Within the bipolar cavities, L. Bonne et al. (2022) identified blueshifted [C ii] shells with a velocity of 5.2 ± 0.5 km s⁻¹, likely driven by stellar winds from the massive cluster. Additionally, they find diffuse X-ray emission (observed with the Chandra X-ray Observatory) in and around the RCW 36 region, which is tracing hot plasma created by the winds. However, L. Bonne et al. (2022) estimate that the energy of the hot plasma is 50%–97% lower than the energy injected by stellar winds and reason that the missing energy may be due to plasma leakage, as has been previously suggested for RCW 49 (M. Tiwari 2021).

The magnetic field geometry on >1 pc cloud scales is traced by the BLASTPol 500 μm polarization map (L. M. Fissel et al. 2016), which has an FWHM 25 resolution, corresponding to 0.65 pc at the distance of Vela C. The BLASTPol magnetic field orientation is shown by vectors in the top left panel of Figure 3, which follow a fairly uniform east–west morphology that is mostly perpendicular to the orientation of the dense ridge. However, around the north and south "bends" of the bipolar structure, the magnetic field lines also appear to bend inward toward the center, following the bipolar shape of the structure.

The middle panels of Figure 3 highlight structures on filament scales of ∼0.1–1 pc, and the cyan contours in all panels represent Herschel-derived column densities to show the filament. At the waist of the bipolar nebula, V. Minier et al. (2013) identify a ringlike structure that extends ∼1 pc in radius and is oriented perpendicular to the bipolar nebula lobes, in the north–south direction (labeled in yellow in both the left and middle panels of Figure 3). The majority of the dense material, as traced by the column density contours, is contained within this ring. V. Minier et al. (2013) also model the kinematics of the ring and find that the northeastern (NE) half is mainly blueshifted while the southwestern (SW) half is redshifted, consistent with an expanding cloud with speeds of 1–2 km s⁻¹. To trace the ionized gas, i.e., H ii, we use archival Hα data from the SuperCOSMOS H-alpha Survey (N. C. Hambly et al. 2001; Q. A. Parker et al. 2005). From the SuperCOSMOS map, we note that the eastern side of the ring is seen in Hα absorption, signifying that it is in front of the ionizing gas and associated massive star cluster. Whereas, Hα emission is seen across the western region of the ring, and therefore, this part of the ring is likely behind the cluster.

The highest column density contours are observed within the southwestern half of the ring, where most of the next-generation star formation appears to be taking place. We henceforth refer to this region as the "Main-Fil" (labeled in white, middle panel, Figure 3). T. Hill et al. (2012) estimate that the mass per unit length of the Main-Fil region is 400 ± 85 M_⊙ pc⁻¹. The Main-Fil is seen to host multiple star-forming cores and/or clumps, which are shown by the black ALMA Band 6 continuum contours in the middle panel of Figure 3.

Several diffuse filamentary structures traced by ¹²CO and [C ii] (shown in red and blue, respectively in the middle panel of Figure 3) can also be observed in the ambient cloud surrounding the ring. L. Bonne et al. (2022) have reasoned that, due to the curved shape of these filaments, they are not part of the larger expanding [C ii] shells but may have been low-density filaments originally converging toward the center dense ridge (similar to the converging filaments seen in Musca, B211/3, and DR 21; P. F. Goldsmith et al. 2008; N. Schneider et al. 2010; P. Palmeirim et al. 2013; N. L. J. Cox et al. 2016; L. Bonne et al. 2023) that have instead been swept away at velocities >3 km s⁻¹ due to stellar feedback.

The magnetic field geometry on filament scales, traced by SOFIA/HAWC+ 214 μm, is fairly consistent with the east–west geometry seen on cloud scales, with some interesting exceptions. The most striking deviation is the region located just east of the ionizing stars, hereafter referred to as "Flipped-Fil" (labeled in Figure 3, middle and right panels). Here, the magnetic field morphology, as traced by SOFIA/HAWC+, deviates from the east–west trend and abruptly flips almost by 90° to follow a more north–south configuration. This geometry appears to follow the elongation of a lower-density filament traced by ¹²CO (red, middle panel), [C ii] (blue, middle panel), and ¹³CO emission (blue, right panel).

The right panels of Figure 3 show emission on subclump and core scales (<0.1 pc) in RCW 36. These data reveal complex substructure within the Main-Fil region. The Main-Fil is clumpy with several bright-rims, voids, and pillar-like structures identified by V. Minier et al. (2013; see their Figure 3). This matches our ALMA band 6 continuum observations as shown by the five near-round clumps and associated elongated pillars, as seen in the right panel of Figure 3.

V. Minier et al. (2013) recognized that the bright-rims appear near the end of the pillar-like structures. The bright-rims are traced in the right panel of Figure 3 by Spitzer/IRAC 3.6 μm emission, which mainly traces hot dust, found at the edges of the PDR (V. Minier et al. 2013). These rims appear to wrap around the cold ALMA clumps, without covering them completely, in a manner resembling bow shocks (though actual bow shocks are unlikely in this region). These bright-rims are of great interest in this work and are therefore collectively labeled as "Bent-Fils" as they will be referred to in later sections. There is a prominent northern Bent-Fil and southern Bent-Fil shown by the two yellow-dotted ovals in the lower right panel of Figure 3. The curved morphology of the Bent-Fils is noted by V. Minier et al. (2013) to be likely due to tracing the inner border of the dense ring, which is being progressively photoionzed by the star cluster. Interestingly, the HAWC+ magnetic field morphology seems to follow the Bent-Fil features, deviating once again from the general east–west cloud-scale magnetic field.

In this section, we discuss the procedure of the HRO method (see J. D. Soler et al. 2013; J. D. Soler & P. Hennebelle 2017, for a more detailed description), which computes the relative angle ϕ(x, y) between the gradient vector field of a structure map M(x, y) and the plane-of-sky magnetic field orientation ${\hat{B}}_{\perp }$ at each pixel. The steps for this procedure are outlined in the following subsections.

In this work, we apply two different methods to obtain a 2D structure map M(x, y). In the first approach, henceforth referred to as single map HRO (described further in 4.2.1), we compare the orientation of local structure at every location using one map M(x, y), for each of the data sets listed in Table 1, to the magnetic field orientation measured by dust polarization as was done by previous HRO studies (Planck Collaboration et al. 2016; J. D. Soler et al. 2017; L. M. Fissel et al. 2019; D. Lee et al. 2021, e.g.). In the second approach, applied only to the spectroscopic data cubes listed in Table 3, we slice the spectral line cube into multiple velocity slabs v_i and compare the orientation of structures in the integrated intensity of each slab M(x, y)_i to the inferred magnetic field. This quantifies the relative alignment as a function of line-of-sight velocity, and will therefore be referred to as the velocity dependent HRO (see 4.2.2).

4.2.1. Single Map HRO

The dust emission, column density, and temperature maps are already in the 2D spatial format and are thus used directly as the M(x, y) structure map in the single map HRO analysis. For the spectroscopic cubes, we generate a single integrated line intensity map as M(x, y) (following the procedure of L. M. Fissel et al. 2019). To calculate the integrated line intensity, a velocity range v₀–v₁ is selected for each molecular line over which the radiation temperature T_R in a given velocity channel v is integrated, using

$$\begin{eqnarray}&&M(x,y)={\int }_{{v}_{0}}^{{v}_{1}}{T}_{R}\,{dv}.\end{eqnarray} \tag{ 5 }$$

The velocity ranges used to calculate the integrated intensity map in the single map HRO analysis for each spectral cube are specified in column (5) of Table 3.

4.2.2. Velocity Dependent HRO

Additionally, for the molecular line data, we also perform a velocity dependent HRO analysis. Here, we slice a selected velocity range v₀–v₁ (specified in column (6) of Table 3) into narrower velocity slabs with a width of 1 km s⁻¹. We increment the center velocities v_i of each slab by 0.5 km s⁻¹, such that v_i = {v₀ + 0.5, v₀ + 1, ..., v₁ − 1, v₁ − 0.5}. For every velocity v_i in the set, we generate an integrated intensity map M(x, y)_i using

$$\begin{eqnarray}&&M{(x,y)}_{i}={\int }_{{v}_{i}-0.5}^{{v}_{i}+0.5}{T}_{R}\,{dv}.\end{eqnarray} \tag{ 6 }$$

The 1 km s⁻¹ width of the slabs is chosen to be roughly a factor of 5 larger than the ∼0.2 km s⁻¹ velocity resolution of the data cubes (listed in column (4) of Table 3) such that enough velocity channels are included in the integrated intensity map. This ensures that there is sufficient signal-to-noise in each slab, and small local fluctuations are averaged over. The HRO analysis is then repeated for each M(x, y)_i in the set.

To directly compare the structure map M(x, y) and the plane-of-sky magnetic field map ${\hat{B}}_{\perp }$, the next step is to ensure that both maps share the same spatial coordinate grid such that there is one-to-one mapping between the pixels. To do this, the map with the coarser pixel scale is projected onto the grid of the map with the finer pixel-scale (i.e., if the M(x, y) map has a lower resolution than the ${\hat{B}}_{\perp }$ map, then M(x, y) is projected on the ${\hat{B}}_{\perp }$ map). The pixel sizes of each data set are given in Table 1 for comparison. When the ${\hat{B}}_{\perp }$ data have the lower resolution of the two, rather than directly projecting the orientation of the magnetic field lines inferred from Equation (3), we instead project the Stokes Q and U intensity maps separately and then recalculate the inferred ${\hat{B}}_{\perp }$. This avoids an incorrect assignment of vector orientation angles to the newly sized pixels.

Next, a 3σ signal-to-noise cut is applied to the data points in the structure map M(x, y). For the single map HRO analysis, all of the M(x, y) maps of the RCW 36 region are above this threshold for every point that overlaps with the ${\hat{B}}_{\perp }$ data, with the exception of the ALMA continuum maps for which the signal is concentrated near the dense clumps. For the velocity dependent HRO analysis, the M(x, y)_i integrated intensity map of each velocity slab is masked individually.

To determine the orientation of elongated structures in M(x, y), we calculate the direction of the isocontours ψ(x, y) (which is by definition perpendicular to the gradient vector field ∇M), given by

$$\begin{eqnarray}&&\psi =\arctan \left(\displaystyle \frac{\delta M/\delta x}{\delta M/\delta y}\right),\end{eqnarray} \tag{ 7 }$$

where ψ is calculated at each pixel (x, y). The partial derivatives are calculated by convolving M(x, y) with Gaussian derivative kernels G, using

$$\begin{eqnarray}&&\displaystyle \frac{\delta M}{\delta x}=M(x,y)* \displaystyle \frac{\delta }{\delta x}G(x,y),\end{eqnarray} \tag{ 8 }$$

and similarly δM/δ y = M(x, y) ⋆ δ G(x, y)/δ y. This reduces noise and avoids erroneous relative angle measurements due to map pixelization. The size of the Gaussian kernels in angular units θ_G is chosen to be one-third of the FWHM angular resolution θ_beam of the M(x, y) map, using θ_G = θ_beam/3. If this kernel size θ_G is less than 3 pixels, then a minimum kernel size of 3 pixels is used instead. A summary of all the kernel sizes in angular units θ_G and pixel units l_G is provided in columns (5) and (6), respectively, of Table 1. The same smoothing lengths listed in Table 1 for the molecular line data are applied for the velocity dependent HRO analysis.

The relative angle ϕ(x, y) between the isocontour direction ψ(x, y) and the plane-of-sky magnetic field ${\hat{B}}_{\perp }$ can then be computed with

$$\begin{eqnarray}&&\phi \equiv \arctan \left(\displaystyle \frac{| \psi \times {\hat{B}}_{\perp }| }{\psi \cdot {\hat{B}}_{\perp }}\right).\end{eqnarray} \tag{ 9 }$$

The resulting ϕ falls within the range [0°, 180°], but since ϕ measures only the orientation and not direction, the angles ϕ and 180−ϕ are redundant. The range can therefore be wrapped on [0°, 90°] as we are only concerned with angular distance, such that ϕ = 0° (and equivalently ϕ = 180° before wrapping) corresponds to the local structures being aligned parallel relative to the magnetic field orientation, while ϕ = 90° corresponds to perpendicular relative alignment. A histogram can then be used to combine the relative angle measurements across all pixels in order to summarize the overall trend within the map. We place the ϕ(x, y) measurements into 20 bins over [0°, 90°], where each bin is of 45 in size.

While the histogram is a useful tool for checking if there is a preference toward a particular relative angle, we can go a step further and quantify the statistical significance of such a preference by calculating the projected Rayleigh statistic (PRS; as described in D. L. Jow et al. 2018). The PRS is a modified version of a classic Rayleigh statistic, which tests for a uniform distribution of angles using a random walk. The classic Rayleigh statistic characterizes the distance Z from the origin if one were to take unit-sized steps in the direction determined by each angle. Given a set θ_i of n angles within the range [0°, 360°], this distance Z can be calculated as follows:

$$\begin{eqnarray}&&Z=\displaystyle \frac{{\left({{\rm{\Sigma }}}_{i}^{n}{\rm{\cos }}{\theta }_{i}\right)}^{2}+{\left({{\rm{\Sigma }}}_{i}^{n}{\rm{\sin }}{\theta }_{i}\right)}^{2}}{n},\end{eqnarray} \tag{ 10 }$$

where n is the number of data samples. To use the set of relative angles ϕ(x, y) in the range [0°, 90°] determined from HROs, we can map each angle θ = 2ϕ. The range of possible Z then is [0, n], where Z = 0 is expected if the angles θ_i are distributed randomly, and Z = n is expected if all angles are the same. Any significant deviation from the origin would signify that the angles θ_i have a directional preference and are nonuniform. While this statistic is useful for testing for uniformity, it cannot differentiate between the preference for parallel versus perpendicular alignment, which is what we would like to measure in the context of HROs. To achieve this, D. L. Jow et al. (2018) modify this statistic by calculating only the horizontal displacement ${Z}_{x}^{{\prime} }$ in the hypothetical random walk:

$$\begin{eqnarray}&&{Z}_{x}^{{\prime} }=\displaystyle \frac{{{\rm{\Sigma }}}_{i}^{n}\cos {\theta }_{i}}{\sqrt{n/2}}.\end{eqnarray} \tag{ 11 }$$

Now, a parallel relative angle ϕ_i = 0 will map to θ = 2ϕ = 0 and give a positive $\cos (0)=1$ contribution to ${Z}_{x}^{{\prime} }$, while a perpendicular relative angle ϕ_i = π/2 will map to θ = 2ϕ = π and give a negative $\cos (\pi )=-1$ contribution to ${Z}_{x}^{{\prime} }$. Therefore, a statistic of ${Z}_{x}^{{\prime} }$ ≫ 0 indicates strong parallel alignment, while ${Z}_{x}^{{\prime} }$ ≪ 0 indicates strong perpendicular alignment. Since the value of ${Z}_{x}^{{\prime} }$ is within the range $[-\sqrt{2n},+\sqrt{2n}]$, the statistic can be normalized by $\sqrt{2n}$ to give a measure of the degree of alignment:

$$\begin{eqnarray}&&{\tilde{Z}}_{x}^{{\prime} }=\displaystyle \frac{{Z}_{x}^{^{\prime} }}{\sqrt{2n}},\end{eqnarray} \tag{ 12 }$$

where a value of ${\tilde{Z}}_{x}^{{\prime} }=\pm 1$ would correspond to perfectly parallel or perpendicular alignment.

If the n data points are independent, then ${Z}_{x}^{{\prime} }$ should have an uncertainty of 1. However, most of the maps used in this HRO study are oversampled, and adjacent pixels are not entirely independent. Since the magnitude of ${Z}_{x}^{{\prime} }$ is proportional to the number of data points n^1/2 and the relative alignment ϕ is measured at each pixel (x, y), oversampling within the map can result in a misleadingly large ${Z}_{x}^{{\prime} }$ magnitude. To determine the statistical significance of ${Z}_{x}^{{\prime} }$, we follow the methodology in L. M. Fissel et al. (2019) and correct for oversampling by repeating the HRO analysis on 1000 independent white noise maps M_WN, smoothed to the same resolution as M(x, y), and compared to the magnetic field orientation.

The white noise maps M_WN are generated to be the same size as the real data (M(x, y)), using Gaussian noise with a mean and standard deviation of M(x, y). In these maps, the orientation of the gradient will be a uniformly random distribution. The white noise M_WN map then follows the same projection procedure that was applied to M(x, y) in Section 4.3 followed by the same mask, which was applied to the real data M(x, y). We calculate the corresponding PRS, ${Z}_{\mathrm{WN}}^{{\prime} }$, for each white noise map and determine the mean $\langle {Z}_{\mathrm{WN}}^{{\prime} }\rangle$ and standard deviation ${\sigma }_{{Z}_{\mathrm{WN}}^{{\prime} }}$ of the PRS in the 1000 runs. The value of ${\sigma }_{{Z}_{\mathrm{WN}}^{{\prime} }}$ estimates the amount of oversampling in the map. We can then correct the PRS for oversampling using

$$\begin{eqnarray}&&{Z}_{x}=\displaystyle \frac{{Z}_{x}^{{\prime} }}{{\sigma }_{{Z}_{\mathrm{WN}}^{{\prime} }}}.\end{eqnarray} \tag{ 13 }$$

After correction, the error on the corrected PRS Z_x is now ${\sigma }_{{Z}_{x}}$ = 1, and a magnitude of ∣Z_x ∣ > 3 is considered statistically significant. The number of independent samples can be estimated as

$$\begin{eqnarray}&&{n}_{\mathrm{ind}}=\displaystyle \frac{n}{{\left({\sigma }_{Z{{\prime} }_{\mathrm{WN}}}\right)}^{2}}.\end{eqnarray} \tag{ 14 }$$

For the velocity dependent HRO analysis, a corrected PRS ${{\rm{Z}}}_{{x}_{i}}$ is measured for each integrated intensity map M(x, y)_i found at the velocity slab centered at the velocity v_i . This generates a PRS as a function of velocity.

In this section, we use the HRO method to compare the magnetic field orientations inferred by SOFIA at 89 μm (Band C) to 214 μm (Band E). Figure 4 shows the relative angles between the Band C (${\hat{B}}_{\perp C}$) and Band E (${\hat{B}}_{\perp E}$) magnetic field vectors at each pixel, calculated using Equation (9). We find that the relative angles are near parallel ϕ ∼ 0° at most locations in RCW 36, signifying that the magnetic field orientations in the two different bands are highly consistent. This gives a mean relative angle of 〈ϕ〉 = 66 as well as a large positive PRS of Z_x = 69.8 and normalized statistic of ${\tilde{Z}}_{x}^{{\prime} }=0.95$, indicating a strong preference for parallel relative alignment. A discussion of this result is given in Section 6.1.2.

Zoom In Zoom Out Reset image size

Since the magnetic field orientations in the two bands are very similar, we present only the Band C HRO results in the main text and defer the Band E results to Appendix A. The Band C single map and velocity dependent HRO results are presented in Sections 5.2 and 5.3, respectively.

Table 4 summarizes the single map HRO results, where the SOFIA/HAWC+ Band C data have been used to infer the magnetic field orientation. Most tracers have a negative PRS (Z_x ), indicating a statistical preference for perpendicular alignment. There are also some notable exceptions that have a positive PRS, indicating a preference for parallel alignment. We discuss the results from the various tracers below.

Table 4. PRS Results for the Single Map HRO Analysis Using Magnetic Field Orientation Inferred from HAWC+ 89 μm Data

Instrument	Map	Zx a	${\sigma }_{{Z}_{\mathrm{WN}}^{{\prime} }}$ b	${Z}_{x}^{{\prime} }$ c	${\tilde{Z}}_{x}^{{\prime} }$ d	ne
SPIRE	500 μm	−8.8	7.9	−69.5	−0.48	10,447
SPIRE	350 μm	−11.2	5.7	−63.9	−0.44	10,447
SPIRE	250 μm	−12.4	4.2	−52.2	−0.36	10,447
HAWC+	214 μm	−13.3	4.1	−54.3	−0.39	9798
PACS	160 μm	−9.0	4.1	−37.0	−0.26	10,447
HAWC+	89 μm	−9.3	3.7	−34.7	−0.27	8564
PACS	70 μm	−5.7	3.7	−21.3	−0.15	10,447
IRAC	4.5 μm	+5.9	4.1	+24.3	+0.06	85,309
IRAC	3.6 μm	+7.7	4.1	+31.1	+0.08	85,010
ALMA ACA	1.4 mm	−0.6	3.5	−2.2	−0.03	2638
ALMA 12 m	1.4 mm	−1.0	3.8	−3.8	−0.01	32,739
Herschel	N(H2)	−11.0	3.9	−42.7	−0.30	10,447
Herschel	Temp	−2.8	7.0	−19.7	−0.14	10,447
LAsMA	12CO	−2.2	4.7	−10.3	−0.07	10,447
LAsMA	13CO	−6.9	4.6	−32.0	−0.22	10,447
upGREAT	[C ii]	+0.5	4.4	+2.2	+0.02	9071
upGREAT	[O i]	−3.0	7.2	−21.2	−0.15	9634
Mopra	HNC	−2.2	7.6	−16.4	−0.11	10,447
Mopra	C18O	−2.7	7.2	−19.3	−0.13	10,447
Mopra	N2H+	−3.5	7.7	−27.0	−0.19	10,447

Notes. The structure map used for each molecular line is a single intensity map integrated over the velocity range specified in Table 1 using Equation (5). A negative Z_x corresponds to an overall preference for perpendicular alignment, while a positive value corresponds to a preference for parallel alignment. The larger the magnitude of Z_x , the stronger the statistical significance of the preferred alignment.

^a PRS corrected for oversampling using Equation (13). ^b The standard deviation for 1000 white noise runs used to correct oversampling in PRS from Equation (13). ^c The uncorrected PRS from Equation (11). ^d The normalized uncorrected PRS from Equation (12). ^e The number of relative angle points used to calculate the uncorrected PRS.

Download table as:  ASCIITypeset image

The single-dish dust emission maps show a distinct variation in the magnitude of the PRS values with wavelength. The top panel of Figure 5 shows the oversampling-corrected Z_x values, which indicate the statistical significance of the PRS, while the bottom panel shows the normalized uncorrected ${\tilde{Z}}_{x}^{{\prime} }$ values, which indicate the degree of alignment. Both panels show that Z_x and ${\tilde{Z}}_{x}^{{\prime} }$ are negative for the submillimeter and far-IR Herschel and SOFIA data indicating a preference for perpendicular alignment, and positive for the mid-IR Spitzer data, indicating a preference for parallel alignment. All the maps have statistically significant values Z_x (i.e., ∣Z_x ∣ > 3${\sigma }_{{Z}_{x}}$, where ${\sigma }_{{Z}_{x}}$ = 1). A notable trend is seen for the normalized statistic where the ${\tilde{Z}}_{x}^{{\prime} }$ values roughly increase (i.e., becomes less negative) as the wavelength decreases, suggesting that successively more structures within the maps align parallel to the magnetic field at shorter wavelengths. In comparison, the trend for the oversampling-corrected statistic Z_x decreases from 500 to 250 μm and peaks in magnitude at 214 μm. This is because the oversampling correction factor (${\sigma }_{{Z}_{\mathrm{WN}}^{{\prime} }}$) is proportional to the number of independent samples in the masked area, and lower Z_x values are expected for the same degree of alignment at longer wavelengths with larger beams, where there are fewer pixels over the same area (see Equation (14)). We also ran Monte Carlo simulations to test whether measurement uncertainties in the magnetic field orientation could affect our measured PRS values. We find that the uncertainty in the relative angles ϕ(x, y) has a negligible effect on the PRS, resulting in an uncertainty of ±0.2 for Z_x and ±0.002 for ${\tilde{Z}}_{x}^{{\prime} }$. These tests and a discussion of our error propagation methods are described in Appendix B.

Zoom In Zoom Out Reset image size

Figure 6 summarizes the oversampling-corrected and normalized PRS for the total integrated intensity spectral line maps. Unlike the Z_x values for the dust map shown in Figure 5, many of the spectral line intensity maps show no overall preference for alignment, or only show a statistically insignificant alignment trend. These maps show different alignment preferences relative to magnetic field in different regions, or overlapping along the line of sight. In Section 5.3, we will show that some of these structures that overlap in the integrated intensity map can be decomposed into different line-of-sight velocity channels. In some cases, particularly for the Mopra observations, the low Z_x values are also in part due to lower resolution and higher noise levels of the spectroscopic data. Overall, we note that all gas tracers have a negative ${\tilde{Z}}_{x}^{{\prime} }$, signifying a preference for perpendicular alignment, with the exception of [C ii], which has a positive ${\tilde{Z}}_{x}^{{\prime} }$.

Zoom In Zoom Out Reset image size

In Figure 7, we identify which structures are aligned with the magnetic field for a select number of data sets. Similar plots for the remaining data sets in Table 4 are shown and discussed in Appendix C, Figures 14–16. The right column shows the structure map M(x, y) overlaid with the magnetic field orientation. The middle column shows the relative angle ϕ(x, y) calculated at each location in the region, where purple (ϕ ∼ 0°) is associated with local parallel alignment, and orange (ϕ ∼ 90°) is associated with local perpendicular alignment relative to the magnetic field. The left column summarizes the alignment trend using a histogram of the relative angles.

Zoom In Zoom Out Reset image size

From Figure 7, we first compare the relative angle maps for dust emission at wavelengths of 500 and 89 μm. We note that, for both maps, the majority of the relative angles are near-perpendicular and are concentrated at the left and right sides of the dust map, which correspond to the east and west halves of the dense molecular ring labeled in Figure 3. This is consistent with the visual observation that the ring is elongated approximately along the north–south direction, which is oriented roughly perpendicular to the mostly east–west magnetic field morphology from HAWC+ Band C observations. Both 500 and 89 μm wavelengths also show near-parallel relative angle measurements (ϕ ∼ 0°) within the roughly north–south oriented Flipped-Fil structure, oriented parallel to the local north–south magnetic field. The main difference between the 500 and 89 μm, however, appears to be the south Bent-Fil structure (labeled in Figure 3). This structure is traced at the shorter 89 μm wavelength, but not at 500 μm. Since the Bent-Fil structures are elongated east–west, parallel to the local magnetic field, this results in the 89 μm having an overall lower ${\tilde{Z}}_{x}^{{\prime} }$ magnitude, indicating less of a statistical preference for perpendicular alignment relative to the magnetic field. This general trend is noted for all other submillimeter and far-IR the dust maps from 350 to 70 μm as well (see Appendix C.0.1, Figure 14), where the emergence of the southern Bent-Fil structure at wavelengths <160 μm results in less negative Z_x values as observed in Figure 5.

However, unlike the far-IR and submillimeter dust maps, the 3.6 μm Spitzer data are less sensitive to the high column density ring structure and instead predominantly trace emission near the north and south Bent-Fils, which are oriented parallel to the east–west magnetic field. The lack of perpendicular relative angles from the dense ring results in an overall positive Z_x .

The ALMA continuum maps show no significant preference for either parallel or perpendicular alignment. This is likely because ALMA interferometer is resolving out many of the large-scale dense ring, Main-Fil, Flipped-Fil, and Bent-Fils structures (see full discussion in Appendix C.0.1 and Figure 16).

Examining the molecular gas maps, we find that ¹³CO, which is sensitive to intermediate-density gas, is able to trace both the dense molecular ring and the south Bent-Fil structure, resulting in a Z_x value comparable to the 89 μm dust map.

The [C ii] relative angle map shows that the east–west Bent-Fil structures contribute about the same number of parallel aligned relative orientations (ϕ ∼ 0°) as the perpendicular ϕ ∼ 90° relative orientations near the north–south ring. This results in a PRS (Z_x = +0.5) that is close to 0 and has neither a statistical preference for perpendicular nor parallel alignment relative to the magnetic field. Distinctively, the histogram of [C ii] also does not peak at near-parallel or near-perpendicular angles, but rather close to ϕ ∼ 40°. Although it is not statistically significant, the positive Z_x result of [C ii] is interesting as it is in contrast with the negative Z_x results for all other spectral line tracers, which predominantly trace the molecular dense ring (Figure 17). Furthermore, we see that the emission in the integrated intensity [C ii] map correlates with the Spitzer emission, which probes warm dust likely found near PDRs. It is therefore noteworthy that the Z_x results for both the [C ii] and Spitzer maps are positive, indicating that structures associated with PDRs have a preference toward parallel alignment relative to the magnetic field.

In summary, we find a fairly bimodal trend in the single map HRO results. Maps that predominantly trace the high column density ridge and ring structure (such as longer wavelength dust maps and high-density gas tracers) show an overall preference for perpendicular alignment relative to the magnetic field. Whereas maps that trace more diffuse structures near or within the PDR (such as mid-infrared dust maps and [C ii]) show more of a tendency toward parallel alignment relative to the magnetic field. Maps that show a combination of the two types of structures (such as 160–70 μm dust maps and low-to-intermediate-density gas tracers) show both regions of parallel and perpendicular alignment relative to the local magnetic field, resulting in a final PRS of lower magnitude. We discuss some caveats and considerations from our HRO analysis in Appendix C.0.3. In the next section, we discuss the HRO analysis for different line-of-sight velocity ranges in the spectral line cubes. We use this velocity dependent HRO approach to examine the relationship between the orientation of the different line-of-sight gas structures and the magnetic field.

In this section, we present the results of the velocity dependent HRO analysis, which measures the PRS of a spectral cube as a function of line-of-sight velocity, using the method described in Section 4.2.2. Figure 8 shows the corrected PRS results for different spectral lines as a function of velocity. We note that, while the magnitude of the Z_x is not always statistically significant (<3σ) over the velocity range for all tracers, the overall trend is interesting and consistent with the single map HRO results. The intermediate and dense gas tracers, C¹⁸O, HNC, N₂H⁺, typically have statistically significant negative Z_x values, especially around 5–8 km s⁻¹, implying a preference for perpendicular alignment. This velocity range matches the mean line-of-sight velocity of the cloud of around 7 km s⁻¹ (V. Minier et al. 2013; L. Bonne et al. 2022). In contrast, the [C ii] PRS results switch from negative Z_x values around 5–6 km s⁻¹ to statistically significant positive Z_x values 9–11 km s⁻¹, indicating a preference for parallel alignment at higher line-of-sight velocities.

Zoom In Zoom Out Reset image size

Figure 8 also demonstrates the limitations of using a single integrated intensity map for a spectroscopic cube in the HRO analysis, as was done in Section 5.2. Since there can be multiple overlapping elongated structures at different line-of-sight velocities, measuring the PRS from only one integrated intensity map may result in a loss of information on the alignment preferences of kinematically distinct structures. To differentiate which structures are being observed at different velocities, Figures 9–11 show channel emission maps from 4 to 10 km s⁻¹.

Zoom In Zoom Out Reset image size

In the ¹³CO channel maps (shown in Figure 9), we notice that emission from the northeastern half of the ring structure is most prominent at ∼4–6 km s⁻¹, while the southwestern section of the ring is most prominent at ∼6–8 km s⁻¹. Since the ring including the Main-Fil is oriented north–south, approximately perpendicular to the east–west magnetic field, the overall Z_x at these velocities is preferentially perpendicular and therefore negative. However, the area of the northeastern region of the ring is smaller and contains fewer HAWC+ Band C polarization detections, leading to less ϕ(x, y) ∼ 90° pixels at 4–5 km s⁻¹ compared to the larger southwestern component of the ring, resulting in a less negative Z_x . At line-of-sight velocities >8 km s⁻¹, the gas traces the Flipped-Fil and north Bent-Fil, causing a weak preference for parallel alignment around 10 km s⁻¹.

Similarly, in the [C ii] channel maps (Figure 10), we notice that the eastern half of the ring can be seen in the 4–5 and 5–6 km s⁻¹ maps, followed by the western half of the ring at 6–8 km s⁻¹ along with the southern Bent-Fil, all of which result in an overall negative and therefore preferentially perpendicular, negative Z_x for these velocity channels. At higher velocities, we see emission from the Flipped-Fil and northern Bent-Fil at 8–10 km s⁻¹, resulting in an overall positive Z_x in these velocity channels. Since the Bent-Fil features are more prominent in [C ii] than other gas tracers, it has the largest positive Z_x magnitude.

Zoom In Zoom Out Reset image size

In contrast, HNC (Figure 11), which is tracing denser gas, mostly shows emission tracing the dense ring (eastern half at 4–6 km s⁻¹ and western half at 6–8 km s⁻¹) has a mostly negative Z_x . The channel maps for the rest of the molecular lines are shown in Appendix D, Figures 19–21, all of which show emission from the dense ring from 4 to 8 km s⁻¹, and the Flipped-Fil and Bent-Fil features from 8 to 10 km s⁻¹.

Zoom In Zoom Out Reset image size

This switch from negative to positive PRS as a function of line-of-sight velocity is consistent with our single map HRO findings. In Section 5.2, we noted a bimodal trend in the PRS, where high column density gas tracers showed a preference of perpendicular relative alignment while tracers associated with the PDR and warmer dust showed more parallel relative orientations. From the velocity dependent HRO analysis, we now learn that the dense gas and PDR structures are also kinematically distinct, such that the same PRS bimodality is also observed as a function of line-of-sight velocity. In the next section, we interpret these results and suggest potential physical mechanisms that may be causing the observed trends.

In this section, we discuss the polarization data from SOFIA/HAWC+ in more detail to try and infer the density scales for which the RCW 36 magnetic field is being traced, i.e., the population of the dust grains, which contribute to the majority of the polarized emission. To do this, we first estimate the optical depth of the dust emission in Section 6.1.1 and compare the polarization data and magnetic field morphologies at the different HAWC+ wavelengths in Section 6.1.2.

6.1.1. Optical Depth of Dust Emission

6.1.2. Magnetic Field Comparison

In this section, we discuss the wavelength dependence of the polarization data from HAWC+ since Band C (89 μm) and Band E (214 μm) may be sensitive to different dust grain populations. Emission at 89 μm is typically more sensitive to warm (T ≥ 25 K) dust grains and less sensitive to cold (T ≤ 15 K) dust grains. This is in contrast to 214 μm, which can also probe the magnetic field orientation in colder, more shielded dust columns. However, in Section 5.1, we used the HRO method to statistically show that the magnetic field morphologies inferred from Band C and Band E are almost identical.

This high degree of similarity could suggest that the Band E observations may be measuring polarized radiation mostly emitted by warm dust grains, similar to Band C. Alternatively, the magnetic field morphology in regions where the dust grains are warmer (T > 25 K) may be similar to the field morphology over a wider range of dust grain temperatures. In either case, the Band C polarization data are tracing dust grains with higher temperatures, which, in a high-mass star-forming region like RCW 36, are likely being heated from the radiation of the massive stellar cluster. This warm dust is therefore probably located near the H ii expanding gas shell and associated PDR. This is also supported by the observation that dust polarized intensity of Band C appears to correlate with the PDR-tracing [C ii] and anticorrelate with the ALMA ACA continuum, which traces cold dense cores (as shown in Appendix F). Therefore, the HAWC+ magnetic field is likely weighted toward the surface of the cloud near the PDR, rather than the colder denser dust structures.

Aside from dust temperature, the similarity between Band C and Band E magnetic field orientations may further indicate that the magnetic fields are likely being traced at comparable scales and densities in the two bands. Moreover, a consistent magnetic field morphology can be expected at the different wavelengths if the dust emission is optically thin, as previously suggested.

One noteworthy difference between the Band C and Band E data sets, however, was found by comparing the total polarized intensity given by Equation (1) in Band C (P_C ; smoothed to the resolution of Band E) to Band E (P_E ). Figure 12 shows that the ratio of P_C /P_E is close to unity for the majority of RCW 36 except for certain regions. These regions have higher polarized intensities in the Band C map than they do in Band E by a factor of ∼2–4. Interestingly, the regions also overlap with where the HAWC+ magnetic field is seen to deviate from the general east–west trend of the BLASTPol magnetic field, i.e., the Flipped-Fil and north Bent-Fil. This may be because the dust grains traced by the Band C map produce radiation with a higher degree of linear polarization, due to higher grain alignment efficiency based on a change in temperature and/or emissivity. A similar analysis has been performed by J. E. Vaillancourt et al. (2008) and F. P. Santos et al. (2019) who also compared the polarization ratio in different bands. Radiative alignment torques, which are thought to be responsible for the net alignment of the dust grains with their short axes parallel to the magnetic field, require anisotropic radiation fields from photons of wavelengths comparable or less than the grain size (B. G. Andersson et al. 2015). In this case, we may expect to see the polarization efficiency increase toward regions where the dust has been heated by the young star cluster, such as the PDR as was noted for the Bent-Fils. Another possibility is that the magnetic field lines are more ordered in the gas traced by the warm dust, which is being preferentially traced in Band C. More ordered fields could mean less cancellation of the polarized emission and therefore a higher polarized intensity in comparison to a sight line with more tangled fields. The geometry of the region is also a consideration. The warm dust structures could be inclined at a different angle compared to the cooler layers, as the dust polarized emission is only sensitive to the plane-of-sky magnetic field component.

Zoom In Zoom Out Reset image size

6.2.1. Preferentially Perpendicular Alignment for Dense Tracers

6.2.2. Parallel Alignment for PDR Tracers

One region of particular interest throughout this study has been the Flipped-Fil (labeled in Figure 3) due to the north–south orientation of the magnetic field lines locally within the filament, which is in stark contrast with the general east–west orientation of the surrounding HAWC+ Band C magnetic field morphology. While the magnetic field lines appear to deviate slightly from the east–west trend in several regions such as that in Bent-Fils, the Flipped-Fil region is the most striking feature as the magnetic field lines appear to change direction more abruptly and are almost orthogonal to the magnetic field of the surroundings.

Observational effects like projection may be contributing to the abrupt 90° change in 2D orientation, which may not be as drastic in 3D. A change in the grain alignment mechanism of the dust grains could also cause the near-discontinuous behavior of the Flipped-Fil if the reference direction for grain alignment changed from the magnetic field to the radiation field, as has been theorized for other high-mass star-forming regions such as the Orion Bar (V. J. M. Le Gouellec et al. 2023).

The change in magnetic field orientation within the Flipped-Fil can also be explained through physical origins. One plausible formation scenario was presented by K. Pattle et al. (2018) studying the magnetic field morphology of the Pillars of Creation seen in M16, which resembles the morphology of the Flipped-Fil. The scenario (detailed schematically in Figure 5 of K. Pattle et al. 2018) is summarized here. An ionization front fueled by photon flux from a massive radiating star or cluster is envisioned to approach molecular gas, which may have regions of varying density. The gas being dragged by the ionization front may bow around the overdensity to form an elongated pillar. The flux-frozen magnetic field lines within the pillar would then follow the gas motion and end up perpendicular relative to the background magnetic field orientation. Such a structure could remain stable as the compressed magnetic field lines would provide support against radial collapse since gas flow perpendicular to the field lines would be inhibited. The pillar may gradually erode in the lengthwise direction, however, as gas flow parallel to the field lines would still be allowed.

While such a physical model may be applicable to a structure similar to the Flipped-Fil, there are obvious differences between our observations and the Pillars of Creation. In the spectral line, data of the Flipped-Fil are observed at line-of-sight velocities of ∼8–10 km s⁻¹, which is redshifted compared to the northern and southern halves of the dense ring. This arrangement could have occurred if the expansion of the H ii region swept up the Flipped-Fil and pushed it behind the massive cluster such that it is currently at a farther distance away from us and thus receding at a faster line-of-sight velocity than the main ridge. It is also difficult to distinguish from a 2D projection on the plane-of-sky if the Flipped-Fil is indeed a pillar and column-like structure or whether it is a ridge of material. Moreover, while The Pillars of Creation are photoionized columns. It is not immediately obvious if the Flipped-Fil is directly associated with the PDR as it is not seen in the Spitzer maps, which trace hot dust, but is seen in the [C ii] and [O i] integrated intensity maps, which traces irradiated dense gas. The lack of mid-infrared emission toward the Flipped-Fil is likely not due to absorption from foreground structures as the region is associated with Hα emission (see Figure 3). Additionally, at a column density of A_V ∼ 13, the Flipped-Fil is traced by low and intermediate gas tracers such as ¹²CO and ¹³CO indicating that it is has a molecular gas component, but is not quite dense or cold enough to be traced as clearly in N₂H⁺ and HNC.

The Flipped-Fil thus shows clear differences from the Bent-Fils, which are likely associated with the warm dust structures and traced by shorter wavelength maps (3.5–160 μm) as well as PDR tracers such as [C ii] and [O i]. The Pillars of Creation formation scenario suggested for the Flipped-Fil can also be applicable to the Bent-Fils. In this picture, the star-forming clumps seen in ALMA ACA continuum data (see Figure 3) could be the overdense structures envisioned in Figure 5 of K. Pattle et al. (2018), around which the bright-rimmed Bent-Fil structures are being bowed around. The orientation of the bow shapes then may suggest the direction of these ionization fronts. This model is similar to V. Minier et al. (2013) who suggest from comparisons of numerical simulations by P. Tremblin et al. (2012a, 2012b) that the bright rims or what we call "Bent-Fils" are the result of density enhancements in thin shells due to gas compression around the pillar-like structures.

While this pillar formation and similar origin scenarios for the Flipped-Fil and Bent-Fils are certainly plausible, there is insufficient evidence for it to be the favored explanation. Higher resolution infrared observations may help distinguish these structures better, giving more insight to their morphology and origin.

In this section, we examine the energetic balance of the RCW 36 region in light of the new HAWC+ polarization observations presented in this work.

The suggestion that the flux-frozen magnetic field lines are transitioning from their mostly east–west cloud-scale geometry to align parallel with feedback-associated structures implies that the magnetic field pressure must be less than the ram pressure. This change in morphology indicates that the magnetic field is being altered by feedback as it is unable to support the cloud structures from warping. Setting the ram pressure in equilibrium with the magnetic field pressure would therefore give an upper-limit on the magnetic field strength. L. Bonne et al. (2022) estimate the ram pressure energy density within the ring (labeled in Figure 3) to be u_ram = 0.41–3.67 × 10⁻¹⁰ erg cm⁻³. Assuming equipartition, we set the ram pressure energy density u_ram equal to the magnetic energy density u_B so that

$$\begin{eqnarray}&&{u}_{\mathrm{ram}}\geqslant {u}_{B}=\displaystyle \frac{{B}^{2}}{8\pi }(\mathrm{cgs}),\end{eqnarray} \tag{ 15 }$$

and the upper-limit on the magnetic field strength is estimated to be B = 33–99 μG. This is lower than the magnetic field strength of 120 μG estimated by T. Kusune et al. (2016) for the Centre-Ridge using the Davis–Chandrasekhar–Fermi method. Furthermore, our estimate is also lower in comparison to similar high-mass star-forming regions. For instance, DR 21 is measured to have a magnetic field strength of 130 μG (A. Koley et al. 2021), and RCW 120 is estimated to have 100 μG (Z. Chen et al. 2022). Since our upper-limit is crude and based on the assumption that the feedback is ram pressure dominated, we may be underestimating the magnetic field strength.

L. Bonne et al. (2022) also calculate that the turbulent energy density within the ring is u_turb = 4.1–5.1 × 10⁻¹⁰ erg cm⁻³, which is comparable to the ram pressure energy density and our estimated upper-limit for the inferred magnetic field energy density. However, the magnetic field energy is likely not much weaker than the turbulent energy since a fairly ordered (rather than tangled) magnetic field geometry is observed in RCW 36, which is a signature of sub- or trans-Alfvénic conditions (where u_B ≥ u_turb; J. M. Stone et al. 1998). If the turbulent energy was dominant, we would expect the magnetic field orientation to have more random variations, which would decrease any alignment trend and result in Z_x values with smaller magnitudes than our current measurements. Alternatively, the effects of turbulence on the magnetic field morphology may not be visible on the spatial scales probed by SOFIA/HAWC+ if the size of the turbulent eddies are smaller than the size of the beam, such that the polarization component from turbulent motion cancels out along the line of sight to give the appearance of low dispersion, as demonstrated in T. J. Jones (1989). On filament scales, the ordered magnetic field observations from HAWC+ suggest a near-equipartition between the magnetic field energy and turbulent energy, with the ram pressure from stellar feedback dominating in certain regions.

This interpretation is different from the large-scale HRO results using BLASTPol, which suggested that the magnetic field may have been dominating the energetic balance, setting a preferred direction of gas motion on cloud scales for the dense Centre-Ridge to form preferentially perpendicular relative to the magnetic field (see Section 6.2.1). This indicates that the dynamic importance of the magnetic field may be scale dependent in this region or that the energetic balance has changed since the formation of the original generation of stars, which is currently driving the feedback within RCW 36. It should also be noted that, since the filament-scale magnetic field traced by HAWC+ Band C is weighted toward warm dust, the magnetic field may only be less dynamically important near the PDR. Whether this is the case within the cold dense star-forming clumps remains unclear. Presumably, gravitational in-fall will also be a strong contributing force to the energetic budget on core scales. For a more in-depth analysis of the energetic balance within the cores and clumps, polarization data at longer wavelengths and higher resolutions, such as the polarization mosaics from our ALMA 12 m Band 6 program, is needed. We leave the analysis of that data set for future work.

In this section, we show the remainder of the relative angle maps and histograms from the Band C (89 μm) single map HRO analysis that were not presented in Figure 7.

Figure 14 shows the results for the 350–70 μm dust maps. We see that all maps trace the east and west halves of the north–south oriented dense ring, which contribute most of the perpendicular ϕ(x, y) measurements. Comparing the different wavelengths maps, it can be seen that the emission from the longer wavelengths at 350 μm (first row) and 250 μm (second row) trace the ring structure most closely, particularly the denser western half, which includes the Main-Fil, as outlined by the column density contours. This results in the higher degree of perpendicular alignment relative to the magnetic field, as signified by the higher-magnitude ${\tilde{Z}}_{x}^{{\prime} }$ values for the longer wavelength dust maps in Table 4.

Zoom In Zoom Out Reset image size

One similarity for all dust maps at 350–70 μm in Figure 14 is the alignment measured for the Flipped-Fil (labeled in Figure 3). All relative angle maps find ϕ ∼ 0° within the Flipped-Fil structure, even though the region is not always particularly noticeable in the intensity maps. This result of a mostly parallel relative alignment within the Flipped-Fil is consistent with the visual observation of the filament being elongated approximately in the direction of the north–south local magnetic field orientation, as seen in the middle and right panels of Figure 3. This is in contrast to the otherwise east–west orientation of the overall magnetic field morphology in the surrounding RCW 36 region.

The main difference between the different wavelengths is the emergence of the bright-rimmed Bent-Fil structures (labeled in Figure 3), particularly the south Bent-Fil, which begins to appear over the southwestern region of the dense ring at 160 μm and become increasing prominent at 70 μm. The Band C HAWC+ magnetic field lines are observed to largely follow the geometry of these Bent-Fils in the direction of their east–west elongation, resulting in increasingly parallel relative orientations at 160 and 70 μm, which decreases the overall magnitude of the Z_x to be less negative than at 350–214 μm. This trend can be also be seen in the left column HROs, which show a decreasing histogram density near ϕ ∼ 90° for the shorter wavelengths.

We contrast these results to the Spitzer data at 4.5 μm, shown in the two rows of Figure 15. Similar to the 3.6 μm map shown in Figure 7, the 4.5 μm map is also highly sensitive to the east–west Bent Fils structure, but does not trace the dense molecular ring. Since the Bent-Fils are oriented roughly parallel to the east–west magnetic field, this results in a statistically significant positive Z_x . In Section 5.2, we made note of the apparent correlation between the emission of the Spitzer 3.6 μm map and the [C ii] total integrated intensity. Here, we once again note the similarities between the Spitzer emission and the [O i] total integrated intensity map (shown in the second row of Figure 15), which is also a PDR tracer like [C ii]. The HRO results for the [O i] data are further discussed in Section C.0.2.

Zoom In Zoom Out Reset image size

We now analyze the HRO results of the Herschel-derived column density and temperature maps, shown in Figure 15. Similar to longer wavelength dust maps, the column density map is also mostly sensitive to the dense molecular ring elongated in the roughly north–south oriented, particularly the western half containing the Main-Fil. This results in a majority of locally perpendicular ϕ(x, y) angles relative to the east–west magnetic field morphology, giving an overall large negative Z_x . In contrast, the temperature map shows only a slight preference for perpendicular relative alignment as it appears to be mostly tracing the bipolar morphology of the region. The HAWC+ magnetic field follows the curvature of the bipolar nebula and thus results in more parallel relative angles between the magnetic field and structures in the temperature map and a lower magnitude Z_x .

Figure 16 shows the results for the ALMA data. The HRO analysis for both the ALMA 12 m and ACA maps finds a statistically insignificant PRS of Z_x ∼ 0, implying that the structures traced by ALMA do not have a preferred direction of orientation relative to the HAWC+ Band C inferred magnetic field. ALMA being an interferometer, resolves out much of the large-scale structure such as the dense ring, Main-Fil, Flipped-Fil, and Bent-Fils, which are the main features observed by the other dust maps. The HAWC+ Band C data are also likely not probing the magnetic field within the dense cores detected by ALMA as is further discussed in Appendix F. This may explain why there is no correlation between the structures traced by the ALMA data and the magnetic field orientation inferred by SOFIA/HAWC+. Furthermore, there were not enough ALMA data points that were above a 3σ signal-to-noise threshold that also overlapped with the HAWC+ vectors, to produce a robust PRS measurement. An improved HRO study would compare the structures traced by the ALMA continuum maps to the magnetic field on similar core scales, such as that inferred from ALMA polarization mosaics. This is outside the scope of this paper.

Zoom In Zoom Out Reset image size

Figure 17 shows the single map HRO for the different molecular gas tracers. We compare these results with the atomic gas data of [O i] in Figure 15 and [C ii] in Figure 7. For the atomic gas, the [C ii] and [O i] data from SOFIA probe the transition from molecular to ionized gas in the PDR. The regions of parallel (purple) and perpendicular (orange) relative orientation angle in the [C ii] relative angle map are similar to [O i]. The main difference between the two is that [C ii] (with a critical density for collisional deexcitation of ∼2.6 × 10³ cm⁻³; M. Röllig et al. 2006) traces more of the east–west oriented Bent-Fils in the west half of the ring as compared to the [O I]. This may be because the Bent-Fils features in the western half are more diffuse, while [O I] tends to trace denser PDR regions (critical density of 5 × 10⁵ cm⁻³; M. Röllig et al. 2006) and hotter gas (typically ∼200 K). As such, [O I] emission appears to better trace the north–south Main-Fil, resulting in a slight overall preference for perpendicular alignment as compared to [C ii].

Zoom In Zoom Out Reset image size

Next, we examine the single map HRO results for the molecular line data shown in Figure 17. The high-density gas tracers such as HNC and N₂H⁺ (critical densities >10⁴ cm⁻³; Y. L. Shirley 2015; L. M. Fissel et al. 2019) clearly trace the densest north–south region in the west half of the ring where the alignment of the gas structures with respect to the magnetic field is mostly perpendicular within the Main-Fil column density contours, resulting in a negative Z_x value. The C¹⁸O data traces intermediate densities (∼10³ cm⁻³; L. M. Fissel et al. 2019) similar to ¹³CO, but the C¹⁸O Mopra data have a lower signal-to-noise ratio and lower resolution, resulting in a lower Z_x than the APEX ¹³CO data. Interestingly, we find that the south Bent-fil is traced by ¹³CO and HAWC+ 89 μm intensity (shown in Figure 7), and can be somewhat seen in the maps of C¹⁸O and HNC, but is not seen in the N₂H⁺ map, which tends to trace only high-density and colder molecular gas. These observations are contrasted with the ¹²CO integrated intensity map, which traces even lower-density molecular gas and shows very different structure compared to the other spectral lines. While it also maps parts of the ring, ¹²CO shows bright emission around the Flipped-Fil and the north Bent-Fil, which are elongated along the direction of the magnetic field lines resulting in an overall weak preference for perpendicular alignment relative to the magnetic field (Z_x = −2.2). Although ¹²CO is detected throughout the entire RCW 36 region, it is optically thick at the densest regions. Spectra of ¹²CO, ¹³CO, [C ii], and [O i] for the Flipped-Fil are given in Appendix D for reference (for spectra in other regions, see Figure 2 of L. Bonne et al. 2022).

Finally, we make note of some important considerations in our HRO analysis. We note that smoothing can reduce the number of data points near the boundaries of the map. For instance, in the HAWC+ 89 μm map (shown in Figure 7), the Gaussian kernel smoothing removes some of the relative angle measurements near the south edge of the map boundary, which is not the case for the Herschel 70 μm, which covers a larger area on the sky (see Figure 14). Another caveat to note, for all our HRO analysis plots, is that some of the ϕ ∼ 0° relative angle points are due to the gradient amplitude approaching zero as the gradient changes direction at the peak of the isocontours. For example, in Figure 7 at the center of the highest column density contour within the Main-Fil, a thin row of purple ϕ ∼ 0° pixels can be seen along the north–south direction where we would expect the gradient to change direction. This is most obvious for the 350 and 250 μm relative angle maps. Since there is a small percentage of the total ϕ(x, y) pixels in the relative angle map, which are subjected to this effect, the impact on the resulting Z_x is insignificant. Furthermore, it should also be mentioned that, in addition to the hot dust detected by Spitzer, the instrument is also clearly detecting starlight from the massive stellar cluster. Since this emission is not extended but rather from point sources, it is generally not elongated in a particular direction and therefore also does not significantly affect Z_x .