Magnetic Field Alignment Relative to Multiple Tracers in the High-mass Star-forming Region RCW 36
Akanksha Bij Laura M. Fissel Lars Bonne Nicola Schneider Marc Berthoud Dennis Lee Giles A. Novak Sarah I. Sadavoy Thushara G. S. Pillai Maria Cunningham Paul Jones Robert Simon


0000-0001-7505-5223]Akanksha Bij
Department of Physics, Engineering Physics and Astronomy, Queen’s University, Kingston, ON K7L 3N6, Canada 

0000-0002-4666-609X]Laura M. Fissel
Department of Physics, Engineering Physics and Astronomy, Queen’s University, Kingston, ON K7L 3N6, Canada 

0000-0002-0915-4853]Lars Bonne
SOFIA Science Center, NASA Ames Research Center, Moffett Field, CA 94045, USA 

0000-0003-3485-6678]Nicola Schneider
I. Physik. Institut, University of Cologne, Zülpicher Str. 77, D-50937 Cologne, Germany 


Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), 1800 Sherman Avenue, Evanston, IL 60201, USA
Engineering + Technical Support Group, University of Chicago, Chicago, IL 60637, USA

0000-0002-3455-1826]Dennis Lee
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), 1800 Sherman Avenue, Evanston, IL 60201, USA
Department of Physics & Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA 

0000-0003-1288-2656]Giles A. Novak
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), 1800 Sherman Avenue, Evanston, IL 60201, USA
Department of Physics & Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA

0000-0001-7474-6874]Sarah I. Sadavoy
Department of Physics, Engineering Physics and Astronomy, Queen’s University, Kingston, ON K7L 3N6, Canada 

0000-0003-2133-4862]Thushara G. S. Pillai
MIT Haystack Observatory, 99 Millstone Road , Westford, MA, 01827 

0000-0001-7020-6176]Maria Cunningham
School of Physics, University of New South Wales, Sydney NSW 2052, Australia 

0000-0001-9429-9135]Paul Jones
School of Physics, University of New South Wales, Sydney NSW 2052, Australia 

0000-0003-2555-4408]Robert Simon
I. Physik. Institut, University of Cologne, Zülpicher Str. 77, D-50937 Cologne, Germany 

Akanksha Bij
a.bij@queensu.ca



§ ABSTRACT 
We use polarization data from SOFIA HAWC+ to investigate the interplay between magnetic fields and stellar feedback in altering gas dynamics within the high-mass star-forming region RCW 36, located in Vela C. This region is of particular interest as it has a bipolar  region powered by a massive star cluster which may be impacting the surrounding magnetic field. To determine if this is the case, we apply the Histogram of Relative Orientations (HRO) method to quantify the relative alignment between the inferred magnetic field and elongated structures observed in several datasets such as dust emission, column density, temperature, and spectral line intensity maps. The HRO results indicate a bimodal alignment trend, where structures observed with dense gas tracers show a statistically significant preference for perpendicular alignment relative to the magnetic field, while structures probed by photo-dissociation region (PDR) tracers tend to align preferentially parallel relative to the magnetic field. Moreover, the dense gas and PDR associated structures are found to be kinematically distinct such that a bimodal alignment trend is also observed as a function of line-of-sight velocity. This suggests that the magnetic field may have been dynamically important and set a preferred direction of gas flow at the time that RCW 36 formed, resulting in a dense ridge developing perpendicular to the magnetic field. However on filament-scales near the PDR region, feedback may be energetically dominating the magnetic field, warping its geometry and the associated flux-frozen gas structures, causing the observed the preference for parallel relative alignment.






§ INTRODUCTION


Observations and simulations suggest that star formation occurs when density fluctuations undergo gravitational collapse in molecular clouds <cit.>. The interstellar magnetic field is thought to influence the structure and evolution of these molecular clouds through regulating the rate and efficiency at which gas is converted into pre-stellar structures by providing support against collapse and/or directing gas flow <cit.>. In the vicinity of massive stars, stellar feedback in the form of winds, outflows, and radiation pressure can further alter the dynamical and chemical evolution of the molecular cloud <cit.>. 

Stellar feedback has also been observed to reshape magnetic field geometries around expanding ionized bubbles <cit.>, which is consistent with predictions from magnetohydrodynamic (MHD) simulations <cit.>. However the combined impact of both magnetic fields and stellar feedback in high-mass star forming regions remains poorly understood due to various constraints. For instance, simulating the effect of stellar feedback on the parent molecular cloud requires complex sub-grid physics and  demanding computational resources <cit.>. Furthermore measuring the magnetic field strength through observations is challenging. 
While numerous observational techniques such as Zeeman splitting <cit.>, Faraday rotation <cit.> and the Davis-Chandrasekhar-Fermi (DCF) method <cit.> applied to polarized light have been used in the past <cit.>, each technique has limitations and/or only provides partial information about the magnetic field structure. 

Of all methods to study magnetic fields, polarized dust emission is the most commonly used observational tracer in dense molecular clouds. The plane-of-sky magnetic field orientation can be inferred from the linearly polarized emission of non-spherical dust grains, which are thought to align their long axes perpendicular to the local magnetic field lines on average <cit.>. Dust polarization angle maps can therefore be used as a proxy for the magnetic field orientation weighted by density, dust grain efficiency, and dust opacity. Various comparisons between the orientation of the magnetic field lines to the orientation of molecular cloud structures have been studied to gain insight into the role of the magnetic field in the star-forming process <cit.>.

A numerical method known as the Histogram of Relative Orientations (HRO) was developed by <cit.> to statistically characterize this comparison by measuring the relative alignment between the magnetic field orientation, and the orientation of iso-contours of elongated structures measured from a gradient field. Several HRO studies have found that the relative alignment between interstellar structures and the plane-of-sky magnetic field orientation is dependent on density and column density <cit.>. Most notably,  <cit.> implemented the HRO method for 10 nearby (< 400 pc) molecular clouds using polarimetry data with a resolution of 10 at 353 GHz from the Planck satellite. They found that the overall alignment of elongated structures transitioned from either random or preferentially parallel relative to the magnetic field at lower column densities, to preferentially perpendicular at higher column densities, with the switch occurring at different critical column densities for each cloud. In simulations this signature transition to perpendicular alignment for an increasing column density has been seen for strong magnetic fields that are significant relative to turbulence and able to influence the gas dynamics (i.e. dynamically important)
 <cit.>. 

The HRO method has since been applied to younger and more distant giant molecular clouds, such as Vela C (distance of ∼ 900 pc), whose magnetic field morphology was inferred by the BLASTPol (Balloon-borne Large-Aperture Submillimetre Telescope for Polarimetry) instrument at 250, 350 and 500 μm <cit.>. The BLASTPol-inferred magnetic field was compared to column densities derived from Herschel observations by <cit.>, and to the integrated intensities of molecular gas tracers by <cit.>. Both studies found a similar tendency for elongated structures to align preferentially parallel relative to the magnetic field for low column densities or low-density gas tracers, which then switched to preferentially perpendicular for higher column densities or high-density gas tracers. However with a full width at half maximum (FWHM) resolution of ∼3, the BLASTPol observations are only able to probe the Vela C magnetic field geometry on cloud-scales (> 1 pc).

In this work, we extend these HRO studies to filament-scales (∼0.1–1 pc) by using higher resolution polarimetry data from the Stratospheric Observatory for Infrared Astronomy (SOFIA) High-resolution Airborne Wide-band Camera (HAWC+) instrument at 89 μm (Band C) and 214 μm (Band E), with angular resolution of 7.8 (0.03 pc) and 18.2 (0.08 pc), respectively <cit.>. Moreover, we focus on the role of magnetic fields in high-mass star formation by targeting the densest region within Vela C, known as RCW 36 (visual extinction of  > 100 mag) which is within a parsec of a ionizing young (∼ 1 Myr) OB cluster responsible for powering a bipolar  nebula within the region <cit.>.  

Previous HRO studies have mostly compared the relative orientation of the magnetic field to molecular gas structures. Since this study aims to understand the role of stellar feedback, we wish to additionally apply the HRO analysis to structures associated with the photodissociation region (PDR). The PDR is the interfacing boundary between the ionized  region and surrounding molecular cloud where far-ultraviolet (FUV) photons with energies in the range of 6–13.6 eV dominate and dissociate H_2 and CO molecules <cit.>. Several recent  region studies have used observations of  since it traces the PDR <cit.>. For RCW 36, the bipolar  region was investigated by <cit.> who examined the kinematics of  158 μm and  63 μm data from the SOFIA legacy project FEEDBACK <cit.>. 

We build upon these studies by applying the HRO method to multiple complementary observations tracing column density, temperature, molecular gas, as well as the PDR, to construct a more complete picture of how the magnetic field is affecting star formation within RCW 36. The paper is organized as follows: Section <ref> describes the observations. Section <ref> details the physical structure of the RCW 36 and its magnetic field morphology, including noteworthy regions. Section <ref> discusses the HRO method, with the results presented in Section <ref>. In Section <ref>, we interpret the results and compare this work to other studies. Finally in Section <ref>, we summarize our main conclusions.



§ DATA





In this study we use several types of data for the HRO analysis, which can be separated into three categories: polarized emission from dust (Section <ref>), thermal continuum emission from dust (Section <ref>), and  spectroscopic lines (Section <ref>), each of which are summarized in Table <ref>. The observations from SOFIA HAWC+ and ALMA are presented for the first time. Figure <ref> shows the new HAWC+ and ALMA data and Figure <ref> shows archival data. A detailed overview of the each dataset is provided in the following subsections.



 §.§ Dust Polarization Maps





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


    P = √(Q^2 + U^2),


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 <cit.>. This can be corrected for with `debiasing' <cit.>, using:

    P_db = √(Q^2 + U^2 - 1/2(δ Q^2 + δ U^2)),


where δ Q and δ U are the measurement uncertainties in Q and U, respectively. The polarization angle Ê can be calculated using Ê=0.5 arctan(U,Q). The orientation of the plane-of-sky magnetic field  is then inferred to be orthogonal to  Ê such that:


    B̂_⊥=1/2arctan(U,Q) + π/2,


where an angle of 0^∘ points towards Equatorial north and increases eastward. 





  §.§.§ SOFIA HAWC+


In this paper we publish, for the first time, observations of RCW 36 using publicly available archival data from HAWC+ <cit.>, the far-IR polarimeter onboard SOFIA. The RCW 36 region was observed by SOFIA/HAWC+ on 6 and 14 June 2018 as part of the Guaranteed Time Observing (GTO) program (AOR: 70_0609_12), in both Band C (89 μm) and Band E (214 μm), at nominal angular resolutions at FWHM of 7.8 and 18.2, respectively. The observations were done using the matched-chop-nod method <cit.> with a chopping of frequency of 10.2 Hz, chop angle of 112.4^∘, nod angle of -67.5^∘ 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 time for Band C and Band E was 2845s and 947s, respectively. The total intensity maps for each band are shown in the top row of Figure <ref>.

To reduce this data, we used the HAWC+ Data Reduction Pipeline which is described in detail in <cit.> and <cit.> 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 is 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 <ref>) and the polarization percentage and the polarization angle are calculated.  
A χ^2 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 IQU maps <cit.>. This underestimation can be corrected for using an Excess Noise Factor (ENF) given by ENF = √(χ^2/χ_theo) where χ_theo is theoretically expected and χ^2 is measured. The ENF is estimated in the HAWC+ Reduction Pipeline by fitting two parameters I_0 and C_0 using


    ENF = C_0√(1 + (I/5I_0)^2),


where I is the total intensity (Stokes I) of a pixel in units of Jy pixel^-1. The errors of the IQU maps are then multiplied by the ENF. The values of I_0, C_0 used for I, Q, and U maps in Band C and E are summarized in Table <ref>. 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=0 which resulted in a non-physical ENF), likely due to large intensities from bright emission. In this case we forced the ENF to be about 1 (by manually setting I_0 to be a sufficiently large number such as 100 and C_0=1) such that the errors were neither under nor overestimated. 

After correcting the errors, the pipeline rejects any measurements falling below the 3-σ cutoff in the degree of polarization p to 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 <cit.>. For this test we used Herschel far-IR intensity maps (described in Section <ref>) 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.





 §.§ Dust Emission Maps


 


  §.§.§ ALMA

In this work, we present two new interferometric data sets from the Atacama Large Millimeter/submillimeter Array (ALMA) and Atacama Compact Array (ACA). 

The first dataset includes observations of dust continuum and line emission in Band 6 (1.1–1.4 mm) using the ACA with 7-m dishes (ID 2018.1.01003.S, Cycle 6, PI: Fissel, Laura). The observations took place from 22 April to 14 July of 2019 with a continuum sensitivity of 2.15 mJy beam^-1. The configuration resulted in a minimum baseline of 9 m and a maximum baseline of 48 m. The angular resolution is approximately 4.9 and the maximum recoverable scale is 28, which corresponds to spatial scales ranging from ∼0.02 to 0.12 pc (4410 to 25200 au) using a distance estimate of ∼ 900 pc for Vela C from <cit.>. The imaged area was 108 × 324, with 87 mosaic pointings and a mosaic spacing of 21.8. The average integration time was 20 minutes per mosaic pointing.  

The second ALMA program used the 12-m array in the C43-1 configuration to observe both polarization mosaics of some of the dense cores identified in the ACA data, as well as larger spectroscopic and continuum observations with the same correlator configuration as the ACA observations (ID 2021.1.01365.S, Cycle 8, PI: Bij, Akanksha). The spectroscopic and continuum observations took place on 19–24 March 2022 with a continuum sensitivity of 0.2 mJy beam^-1. The configuration resulted in a minimum baseline of 14.6 m and a maximum baseline of 284 m. The angular resolution is approximately 1.2 and the maximum recoverable scale is 11.2, which corresponds to ∼0.0052 to 0.05 pc (1080 to 10080 au) resolution at the distance to Vela C. The imaged area was 42 × 85 in size, with 26 mosaic pointings and a mosaic spacing of 12.3. The average integration time was 9.5 minutes per mosaic pointing. 

In this work, we present only the total intensity dust continuum maps from both datasets, as delivered by the QA2 (Quality Assurance 2) process, with no further data reduction performed. These maps are shown in the bottom row of Fig. <ref>. We do not analyze the spectral data and polarization mosaics from these observations in this work, as further data reduction is required and left for future studies. 



  §.§.§ Herschel SPIRE and PACS


To study the cloud structures probed by thermal dust emission, we use publicly available archival maps from the Herschel Space Observatory, which observed Vela C on 18 May 2010 <cit.> as part of the Herschel OB Young Stars (HOBYS) key programme <cit.>. The observations were conducted using the SPIRE instrument at 500, 350, and 250 μm <cit.>, and the PACS instrument at 70 and 160 μm <cit.>. Additionally, we include a Herschel-derived temperature map at an angular resolution of 36 and a 18 column density map. The column density map is derived from a spectral energy distribution fit to the 160, 250, 350 and 500 μm flux maps, following the procedure described in detail in Appendix A of . We show the 250 μm, 70 μm, and temperature maps in panels b, c, and d of Figure <ref>, respectively. 



  §.§.§ Spitzer IRAC

To trace warmer dust grains, we use archival mid-IR maps from the Spitzer Space Telescope, obtained from the publicly available ISRA NASA/IPAC Infrared Science Archive[Available at https://irsa.ipac.caltech.edu/about.html]. Spitzer observed RCW 36 in May 2006, employing the four-channel camera IRAC to capture simultaneous broadband images at channels 1–4, covering bands centered at 3.6, 4.5, 5.8, and 8.0 μm, respectively <cit.>. IRAC uses two 5.2 × 5.2 fields of view, where one field simultaneously images at 3.6 and 5.8 μm and the other at 4.5 and 8.0 μm. All four detector arrays are 256 × 256 pixels, with 1.2 square pixels. This dataset was published in <cit.>. In this work, we only use data from 3.6 μm (channel 1) and 4.5 μm (channel 2) for our HRO analysis, as channels 3 and 4 have artifacts and saturation. Channels 1 and 2 have resolutions of 1.66 and 1.78, respectively. We show the Spitzer 3.6 μm map in panel j of Figure <ref>. 



 §.§ Atomic and Molecular Line Maps





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 <ref> 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 <ref>) as well as channel maps over narrower velocity ranges (Column 6 of Table <ref>). Our spectral data cube analysis is described in Sections <ref> and <ref>.



  §.§.§ SOFIA upGREAT FEEDBACK

To trace the PDR, we use a 158 μm map of the  ^3P_3/2→ ^3P_1/2 transition (at native resolution of 14.1) and a 63 μm map of the  ^3P_1→ ^3P_2 transition (at native resolution of 6.3) from the SOFIA FEEDBACK C^+ legacy survey <cit.>. The survey was conducted by the upGREAT (upgraded German REceiver for Astronomy at Terahertz frequencies) heterodyne spectrometer <cit.> onboard the SOFIA aircraft <cit.> on 6 June 2019 from New Zealand. The upGREAT receivers use a low frequency array (LFA) to cover the 1.9-2.5 THz band with 14 pixels and a high frequency array (HFA) covering the 4.7 THz line with 7 pixels. The on-the-fly observing strategy and data reduction for the survey is discussed in <cit.>. In this work, we use the data reduced by <cit.> who smoothed the  map to 20 and  map to 30 to reduce noise. They also applied principal component analysis (PCA) to identify and remove systematic components of baseline variation in the spectra. Both maps are 14.4 × 7.2 in size, with a spectral binning of 0.2 km/s, for which the typical rms noise is 0.8-1.0 K for  and ∼0.8-1.5 K for  <cit.>. The integrated intensity maps for  and  and have been integrated from -20 to +20 km/s are shown in panels k and l of Fig. <ref>, respectively.



  §.§.§ APEX LAsMA

To trace the molecular gas regions in RCW 36, we use observations of ^12CO (3–2) and ^13CO (3–2) obtained on 27 September 2019 and 21 July 2021 with the heterodyne spectrometer LAsMA (Large APEX Sub-Millimetre Array), which is a 7 pixel receiver on the APEX telescope <cit.>. The maps were scanned in total power on-the-fly mode and are sized 20 × 15, with a beamsize of 18.2 at 345.8 GHz. <cit.> reduced the data to produce the baseline-subtracted spectra presented here, which have a spectral resolution of 0.2 km/s, pixel size of 9.1 and rms noise of 0.45 K. The integrated intensity maps have been integrated from -20 to +20 km/s for ^12CO and ^13CO are shown in panels e and f of Fig. <ref>, respectively.



  §.§.§ Mopra

In our analysis, we also utilize complementary molecular line surveys from the 22-m Mopra Telescope, which observed the large-scale dense gas over the entire Vela C molecular cloud from 2009 to 2013. In this work, we use only the (1–0) transitions of C^18O as well as HNC and N_2H^+ which were originally presented in <cit.>. The C^18O observations were performed by scanning long rectangular strips of 6 height in the galactic latitude and longitude directions using Mopra's fast-scanning mode. The HNC and N_2H^+ observations used overlapping 5 square raster maps. The data reduction procedure, performed by <cit.>, includes bandpass correction using off-source spectra, bandpass polynomial fitting and Hanning smoothing in velocity. The resulting FWHM angular resolution and velocity resolution is 33 and 0.18 km/s for C^18O, and 36 and ∼0.2 km/s for both HNC and N_2H^+. The integrated intensity maps for HNC, N_2H^+ and C^18O have been integrated over the velocity range given in Column 5 of Table <ref> and are shown in panels g–i of Fig. <ref>, respectively.




§ RCW 36 STRUCTURE


    


In this section, we give an overview of the morphological structure and magnetic field geometry of RCW 36 on varying spatial scales based on previous studies of the region, as well as inferences based on our observations. Figure <ref> showcases various continuum and polarization observations for RCW 36 that will be used to describe its general structure in the following subsections.  



 §.§ Cloud-Scale

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 <cit.> using Herschel data. The densest and most prominent of the ridges is the Centre-Ridge, with column densities of  > 30, and a length of roughly ∼10 pc <cit.>. The Center-Ridge contains the RCW 36 region. It has a bipolar nebula morphology <cit.> 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 <ref> (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 <cit.>. The position of the most massive star (spectral type O9V) is indicated by a white star-shaped marker in Fig. <ref>. The ionizing radiation from this cluster is powering an expanding  gas shell traced by Hα (shown in blue, Fig. <ref>, left panel) <cit.>. 

Bipolar  regions are of great interest because, though they have been observed in other high-mass star forming regions such as S106 <cit.> and G319.88+00.79 <cit.>, they seem to be more rare than single  bubbles <cit.>.

Within the bipolar cavities, <cit.> identified blue-shifted  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, <cit.> 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 <cit.>. 

The magnetic field geometry on > 1 pc cloud scales is traced by the BLASTPol 500 μm polarization map <cit.> which has a FWHM 2.5 resolution, corresponding to 0.65 pc at the distance of Vela C. The BLASTPol magnetic field orientation is shown by vectors in top-left panel of Fig. <ref>, 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 towards the center, following the bipolar shape of the structure. 



 §.§ Filament-Scale

The middle panels of Figure <ref>  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, <cit.> identify a ring-like 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 Fig. <ref>). The majority of the dense material, as traced by the column density contours, is contained within this ring. <cit.> also model the kinematics of the ring and find that the north-eastern (NE) half is mainly blue-shifted while the south-western (SW) half is red-shifted, consistent with an expanding cloud with speeds of 1–2 km/s.
To trace the ionized gas, i.e.   we use archival  data from the SuperCOSMOS H-alpha Survey <cit.>. From the SuperCOSMOS map, we note that the eastern side of the ring is seen in  absorption, signifying that it is in front of the ionizing gas and associated massive star cluster. Whereas,  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 SW 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, Fig <ref>). <cit.> estimate that the mass per unit length of the Main-Fil region is 400±85 M_⊙ pc^-1. 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 Fig <ref>.

Several diffuse filamentary structures traced by ^12CO and  (shown in red and blue, respectively in middle panel of Fig. <ref>) can also be observed in the ambient cloud surrounding the ring. <cit.> have reasoned that due to the curved shape of these filaments, they are not part of the larger expanding  shells but may have been low density filaments originally converging towards the centre dense ridge (similar to the converging filaments seen in Musca, B211/3, and DR 21; ) 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 E-W geometry seen on cloud scales, with some interesting exceptions. The most striking deviation is the region located just east of the ionising stars, hereafter referred to as `Flipped-Fil' (labeled in Fig <ref>, middle and right panels). Here the magnetic field morphology, as traced by SOFIA/HAWC+, deviates from the E-W trend and abruptly flips almost by 90^∘ to follow a more N-S configuration. This geometry appears to follow the elongation of a lower density filament traced by ^12CO (red, middle panel),  (blue, middle panel), and ^13CO emission  (blue, right panel).



 §.§ Core-Scale


The right panels of Figure <ref> show emission on sub-clump 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 <cit.> (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 Fig. <ref>. 

<cit.> recognized that the bright rims appear near the end of the pillar-like structures. The bright-rims are traced in the right panel of Fig.<ref> by Spitzer/IRAC 3.6 μm emission, which mainly traces hot dust, found at the edges of the PDR <cit.>. 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 Fig. <ref>. The curved morphology of the Bent-Fils is noted by <cit.> 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 E-W cloud scale magnetic field. 



§ METHODS





 §.§ Histogram of Relative Orientations


In this section we discuss the procedure of the Histogram of Relative Orientation (HRO) method <cit.> 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  at each pixel. The steps for this procedure are outlined in the following subsections.



 §.§ Preparing the Structure Map


In this work, we apply two different methods to obtain a two-dimensional structure map M(x,y). In the first approach, henceforth referred to as Single Map HRO (described further in <ref>), we compare the orientation of local structure at every location using one map M(x,y), for each of the datasets listed in Table <ref>, to the magnetic field orientation measured by dust polarization as was done by previous HRO studies <cit.>. In the second approach, applied only to the spectroscopic data cubes listed in Table <ref>, 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 <ref>).



  §.§.§ Single Map HRO


The dust emission, column density, and temperature maps are already in the 2-D 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 ). To calculate the integrated line intensity, a velocity range v_0–v_1 is selected for each molecular line over which the the radiation temperature T_R in a given velocity channel v is integrated, using


    M(x,y) = ∫_v_0^v_1 T_R  dv.


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

 


  §.§.§ 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_0–v_1(specified in Column 6 of Table <ref>) 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+0.5,   v_0+1,  …,   v_1-1,   v_1-0.5}. For every velocity v_i in the set, we generate an integrated intensity map M(x,y)_i using:


    M(x,y)_i = ∫_v_i-0.5^v_i+0.5 T_R  dv.


The 1 km/s width of the slabs are 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 <ref>) 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. 



 §.§ Projection and Masking



To directly compare the structure map M(x,y) and the plane-of-sky magnetic field map   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 on to the grid of the map with the finer pixel scale (i.e., if the M(x,y) map has a lower resolution than the  map then M(x,y) is projected on the  map). The pixel sizes of each dataset are given in Table <ref> for comparison. When the  data has lower resolution of the two, rather than directly projecting the orientation of the magnetic field lines inferred from Eq. <ref>, we instead project the Stokes Q and U intensity maps separately and then recalculate the inferred . 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  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. 



 §.§ Calculating the Relative Orientation Angle


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


    ψ=arctan( δ M / δ x/δ M / δ y),


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


    δ M/δ x = M(x,y) ∗δ/δ x G(x,y),


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 <ref>. The same smoothing lengths listed in Table <ref> for the molecular line data are applied for the Velocity Dependent HRO analysis.

The relative angle ϕ(x,y) between the iso-contour direction ψ(x,y) and the plane-of-sky magnetic field  can then be computed with:


    ϕ≡arctan( |ψ×B̂_⊥|/ψ·B̂_⊥).


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 summarise the overall trend within the map. We place the ϕ(x,y) measurements into 20 bins over [0^∘, 90^∘], where each bin is of 4.5^∘ in size.



 §.§ Projected Rayleigh Statistic



While the histogram is a useful tool for checking if there is a preference towards 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 ). 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:


    Z = (Σ_i^n cosθ_i)^2 + (Σ_i^n sinθ_i)^2/n,


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 non-uniform. 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, <cit.> modify this statistic by calculating only the horizontal displacement  in the hypothetical random walk:


    Z'_x = Σ_i^n cosθ_i/√(n/2).


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


    Z'_x = Z'_x/√(2n),


where a value of Z'_x=±1 would correspond to perfectly parallel or perpendicular alignment. 

If the n data points are independent, then  should have an uncertainty of 1. However, most of the maps used in this HRO study are over-sampled and adjacent pixels are not entirely independent. Since the magnitude of  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  magnitude. To determine the statistical significance of   we follow the methodology in <cit.> 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 <ref> followed by the same mask which was applied to the real data M(x,y). We calculate the corresponding PRS, Z'_WN for each white noise map and determine the mean ⟨ Z'_WN⟩ and standard deviation σ_Z'_WN of the PRS in the 1000 runs. The value of σ_Z'_WN estimates the amount of oversampling in the map. We can then correct the PRS for oversampling using:


    Z_x = Z^'_x/σ_Z^'_WN.


After correction, the error on the corrected PRS  is now =1, and a magnitude of | Z_x| > 3 is considered statistically significant. The number of independent samples can be estimated as


    n_ind = n/(σ_Z'_WN)^2.


For the Velocity Dependent HRO analysis, a corrected PRS 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.



§ RESULTS





 §.§ HRO between Band C and Band E



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 <ref> shows the relative angles between the Band C (B̂_⊥ C) and Band E (B̂_⊥ E) magnetic field vectors at each pixel, calculated using Equation <ref>. 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 ⟨ϕ⟩ = 6.6^∘ as well as a a large positive PRS of Z_x=69.8 and normalized statistic of Z'_x=0.95, indicating a strong preference for parallel relative alignment. A discussion of this result is given in Section <ref>. 

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 <ref>. The Band C Single Map and Velocity Dependent HRO results are presented in Sections <ref> and <ref>, respectively.



 

 §.§ Single Map HRO Results



Table <ref> summarizes the Single Map HRO results, where the SOFIA/HAWC+ Band C data has been used to infer the magnetic field orientation. Most tracers have a negative PRS (), 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.



The single-dish dust emission maps show a distinct variation in the magnitude of the PRS values with wavelength. The top panel of Figure <ref> shows the oversampling-corrected  values, which indicate the statistical significance of the PRS, while the bottom panel shows the normalized uncorrected  values, which indicate the degree of alignment. Both panels show that  and  are negative for the sub-mm 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  (i.e. | Z_x| > 3, where  = 1). A notable trend is seen for the normalized statistic where the  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  decreases from 500–250 μm and peaks in magnitude at 214 μm. This is because the oversampling correction factor () is proportional to the number of independent samples in the masked area and lower  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 <ref>). 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 in an uncertainty of ±0.2 for  and ±0.002 for . These tests and a discussion of our error propagation methods are described in Appendix <ref>.

Figure <ref> summarizes the oversampling-corrected and normalized PRS for the total integrated intensity spectral line maps. Unlike the  values for the dust map  shown in Fig. <ref>, 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 <ref>, 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  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 , signifying a preference for perpendicular alignment, with the exception of  which has a positive .

In Figure <ref> we identify which structures are aligned with the magnetic field for a select number of datasets. Similar plots for the remaining datasets in Table <ref> are shown and discussed in Appendix <ref>, Figures <ref>–<ref>. 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.







From Fig. <ref>, we first  compare the relative angle maps for dust emission at wavelengths of 500 μm 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 <ref>. 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 μm and 89 μm wavelengths also show near-parallel relative angle measurements (ϕ∼ 0^∘) within the roughly N-S oriented Flipped-Fil structure, oriented parallel to the local N-S magnetic field. The main difference between the 500 μm and 89 μm, however appears to be the south Bent-Fil structure (labeled in Figure <ref>). This structure is traced at the shorter 89 μm wavelength, but not at 500 μm. Since the Bent-Fil structures are elongated E-W, parallel to the local magnetic field, this results in the 89 μm having an overall lower  magnitude, indicating less of a statistical preference for perpendicular alignment relative to the magnetic field. This general trend is noted for all other sub-mm and far-IR the dust maps from 350–70 μm as well (see Appendix <ref>, Fig. <ref>), where the emergence of the southern Bent-Fil structure at wavelengths < 160 μm results in less negative  values as observed in Fig. <ref>. 

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

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 <ref> and Figure <ref>). 

Examining the molecular gas maps, we find that ^13CO, 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  value comparable to the 89 μm dust map. 

The  relative angle map shows that the E-W Bent-Fil structures contribute about the same number of parallel aligned relative orientations (ϕ∼ 0^∘) as the perpendicular ϕ∼ 90^∘ relative orientations near the N-S 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  also does not peak at near-parallel or near-perpendicular angles, but rather close to ϕ∼ 40^∘. Though it is not statistically significant, the positive  result of  is interesting as it is in contrast with the negative  results for all other of all other spectral line tracers, which predominantly trace the molecular dense ring (Fig. <ref>). Furthermore, we see that the emission in the integrated intensity  map correlates with the Spitzer emission, which probes warm dust likely found near PDRs. It is therefore noteworthy that the  results for both the  and Spitzer maps are positive, indicating that structures associated with PDRs have a preference towards parallel alignment relative to the magnetic field.  

In summary, we find a fairly bimodal trend in the Single Map HRO results. Maps which predominately 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 which trace more diffuse structures near or within the PDR (such as mid-infrared dust maps and ) show more of a tendency towards parallel alignment relative to the magnetic field. Maps which 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 of our HRO analysis in Appendix <ref>. 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.
 



 



 §.§ Velocity Dependent HRO







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 Section <ref>. Figure <ref> shows the corrected PRS results for different spectral lines as a function of velocity. We note that while the magnitude of the  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^18O, HNC, N_2H^+ typically have statistically significant negative  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 <cit.>. In contrast, the  PRS results switch from negative  values around 5–6 km/s to statistically significant positive  values 9–11 km/s, indicating a preference for parallel alignment at higher line-of-sight velocities. 

Figure <ref> also demonstrates the limitations of using a single integrated intensity map for a spectroscopic cube in the HRO analysis, as was done in Section <ref>. 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 <ref>–<ref> show channel emission maps from 4–10 km/s. 

In the ^13CO channel maps (shown in Fig.  <ref>), 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 N-S, approximately perpendicular to the E-W magnetic field, the overall  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 larger southwestern component of the ring, resulting in a less negative . 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  channel maps (Fig. <ref>) 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  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  in these velocity channels. Since the Bent-Fil features are the more prominent in  than other gas tracers, it has the largest positive  magnitude. 

In contrast, HNC (Fig. <ref>) 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 . The channel maps for the rest of the molecular lines are shown in Appendix <ref>, Figures <ref>–<ref>, all of which show emission from the dense ring from 4–8 km/s, and the Flipped-Fil and Bent-Fil features from 8–10 km/s.

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 <ref> 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 line-of-sight velocity. In the next section, we interpret these results and suggest potential physical mechanisms that may be causing the observed trends.



§ DISCUSSION



The goal of this section is to better understand the physical processes behind the observed magnetic field geometry and morphology of the star-forming structures within the RCW 36 region. We are particularly interested in understanding the energetic impact of the magnetic field and stellar feedback in shaping the gas dynamics. To do this, we examine the magnetic field observations inferred from HAWC+ in Section <ref>, followed by an interpretation of the HRO results and discussion of the origins of the Flipped-Fil structure in Sections <ref> and <ref>, and finally comment of the energetic balance of the region in Section <ref>.



 §.§ The Magnetic Field Structure of RCW 36



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 <ref> and compare the polarization data and magnetic field morphologies at the different HAWC+ wavelengths in Section <ref>.



  §.§.§ Optical Depth of Dust Emission


To better understand the location of the dust grains contributing to SOFIA/HAWC+ polarized emission maps, we estimate the optical depth τ_ν to check whether the magnetic field is being inferred from the average column of material along the entire line of sight or if it is tracing only the outer surface layer of an opaque dust cloud. The full method is discussed in Appendix <ref> and summarized here. We use τ_ν = N_H_2μ m_H R_dgκ_ν, where N_H_2 is the column density, μ is the mean molecular weight, m_H is the mass of Hydrogen, R_dg is the dust-to-gas ratio, and κ_ν is the dust opacity. We adopt the same dust opacity law (given in Appendix <ref>) as in previous HOBYS and Herschel Gould Belt Survey <cit.> studies <cit.>. The opacity law is independent of temperature and assumes a dust-to-gas fraction of 1%. We use a Herschel column density map derived by <cit.> which has angular resolution of 36, and is different than the 18 column density map listed in Table <ref> used for the HRO analysis. We choose to use the 36 column density map since the resolution matches the temperature map. The assumed dust opacity law from <cit.> and spectral index is also consistent with <cit.> and <cit.>. 

We estimate that the optical depth for 214 μm (Band E) is optically thin τ_ν≪ 1 everywhere within RCW 36 (see Fig. <ref> in Appendix <ref>). At 89 μm, we also find that the emission is fairly optically thin τ_ν < 1 at 89 μm (Band C) for most regions, except for certain locations within the Main-Fil where the τ_ν can reach values of ∼1.4. It should be noted that these optical depth estimates are uncertain due to the difference in resolution and the possibility of emission from very small dust grains (see Appendix <ref> for details). As a first approximation however, we find that for most regions in RCW 36 the dust emission should be optically thin in HAWC+ Band C and E, and we should therefore be able to probe the full dust column. 



  §.§.§ Magnetic Field Comparison





In this section, we discuss the wavelength dependence of the polarization data from HAWC+ since the 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 <ref>, 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 is 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  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   and anti-correlate with the ALMA ACA continuum which traces cold dense cores (as shown in Appendix <ref>). Therefore the HAWC+ magnetic field is likely weighted towards 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 datasets, however, was found by comparing the total polarized intensity given by Eq. <ref> in Band C  (P_C; smoothed to the resolution of Band E) to Band E (P_E). Figure <ref> 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 E-W 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 <cit.> and <cit.> 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 <cit.>. In this case, we may expect to see the polarization efficiency increase towards 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. 



 §.§ Interpretation of HRO Results




  §.§.§ Preferentially Perpendicular Alignment for Dense Tracers


In Section <ref>, we found that the structure maps M(x,y) which predominantly trace dense structures such as the ring and Main-Fil showed a statistically significant preference for perpendicular alignment relative to the filament-scale (7.8 FWHM) magnetic field probed by HAWC+ Band C. This result is consistent with previous large-scale HRO studies which compared the alignment of structures within the Central-Ridge of Vela C relative to the cloud-scale magnetic field probed by BLASTPol at 250, 350 and 500 μm (2.5–3 FWHM) <cit.>. These studies found that the relative alignment between large-scale structures in the Vela C cloud and the magnetic field is column density and density dependent.  

<cit.> showed that for both the entire Vela C cloud and the Centre-Ridge region (containing RCW 36), the relative alignment trend transitions from preferentially parallel or no preferential alignment at low column densities, to preferentially perpendicular at high column densities. They find that the transition occurs at a threshold column density of N_H∼ 10^21.8. Additionally, <cit.> compared the orientation of the magnetic field inferred from BLASTPol 500 μm to integrated line intensity maps of different molecular lines tracing low, intermediate and high density structures averaged over the entire Vela C Cloud. They found that the low density gas tracers were more likely to align parallel to the magnetic field while intermediate to high density tracers were more likely to align perpendicular, with the transition occurring at a density of n_H_2∼ 10^3 cm^-3. This signature transition to preferentially perpendicular alignment at a critical column density has been observed for several other molecular clouds as well <cit.>.

In this work, we do not see this transition as a function of column density and only observe a preference for perpendicular relative alignment for different column density bins. This could be because RCW 36 region is the densest region within Vela C (with N_H≳ 10^22.4) and most of its structures are above the critical column density.  

We also compare our work to magnetic field studies done on comparable small-scales. <cit.> find that the Musca filament is oriented roughly perpendicular to the surrounding magnetic field morphology, as traced by SOFIA/HAWC+ 214 μm observations. Moreover, <cit.> applied the HRO method to dense cores in the Ophiuchus molecular cloud and found similar results of preferential perpendicular alignment between high column density, elongated filament and core-scale structures in ρ Oph A and ρ Oph E relative to the magnetic field traced by SOFIA/HAWC+ 154 μm observations. The prevalence of this perpendicular relative alignment trend across different star-forming regions in varying molecular cloud environments suggests that shared physical processes may be underlying the observations.

Possible interpretations of such processes have been explored by comparing observations to simulations. For instance, <cit.> propose that the degree of magnetization of a cloud impacts the trend in relative alignment, where the high-magnetization case specifically reproduces the transition from preferentially parallel to preferentially perpendicular at a critical density. Other studies such as <cit.> reason that the preferentially parallel relative alignment occurs in magnetically dominated (sub-Alfvénic) gas, while preferentially perpendicular relative alignment occurs in turbulence dominated (super-Alfvénic) gas with the transition occurring when the kinetic energy due to gravitational contraction becomes larger than the magnetic energy. This connection to the energy balance is consistent with <cit.>, who demonstrate that a transition from parallel to perpendicular relative alignment can occur as a result of convergent velocity flows, which could be due to gravitational collapse. They also find that the transition in alignment can occur when the large-scale magnetic field is strong enough to impose an anisotropic velocity field and set an energetically preferred direction in the gas flow. However, simulations also caution that projection effects are an important consideration in the interpretation of HRO results as <cit.> showed that the relative orientation trends also strongly dependent on viewing angle.

Based on these studies, we propose that the large-scale magnetic field surrounding RCW 36 may have been dynamically important during it's formation, allowing gas to flow preferentially parallel to the E-W magnetic lines. This may have resulted in the formation of an elongated molecular gas sheet or filament (currently the Centre-Ridge) since material could have been inhibited in collapsing perpendicular to the magnetic field lines in the N-S direction. As the region went on to form stars, <cit.> suggest that ionizing radiation from the massive star cluster would have then reshaped the surrounding gas into a bipolar nebula, forming a ring of dense material at the center as an  region expanded into the elongated structure. Both the BLASTPol and HAWC+ maps show that the magnetic field lines pinch near the waist of the bipolar nebula which could be evidence that the ram pressure may be overpowering the local magnetic pressure in that region, as the magnetic field lines are being warped along with the flux-frozen gas. So while the magnetic field may have set a preferred direction of gas flow during the formation of Centre-Ridge and Main-Fil, it may no longer be energetically significant across all of the RCW 36 region since the birth of the massive stars.



  §.§.§ Parallel Alignment for PDR Tracers


Section <ref> also showed that some regions and tracers had a preferential parallel alignment between elongated structures and the inferred magnetic field. The decrease in the statistic magnitude || from dust map wavelengths shorter than 214 μm, was found to be largely due to the gradual emergence of the north and south Bent-Fil features (labeled in Fig. <ref>), which are elongated along the orientation of the HAWC+ Band C magnetic field lines. The emergence of Bent-Fils towards shorter wavelength  (70–214 μm) Herschel and SOFIA dust maps, implies that these features are likely tracing cloud structures with warmer dust populations, near the PDR. The north and south Bent-Fils are also traced by the Spitzer mid-infrared 3.6–4.5 μm maps, which are sensitive to emission from hot dust found near the PDR. 

The observation of a preferential parallel relative alignment between the direction of elongation of the Bent-Fils and the local Band C magnetic field orientation can be explained by the coupling of the gas and the magnetic field. We propose that the stellar feedback in the form of ionizing radiation from the high-mass cluster may be warping the flux-frozen gas, thereby dragging the magnetic field lines along with it. The higher resolution and/or shorter wavelengths of the SOFIA/HAWC+ observations are able to trace the regions where the magnetic field orientation appears to be altered from the otherwise uniform east-west geometry traced by 500 μm BLASTPol observations. The altered field lines follow the warped morphology seen for the bright-rimmed Bent-Fil regions traced by hot dust and  and  PDR tracers. 

Furthermore, the Velocity Dependent HRO results (see Section <ref>) show that the Flipped-Fil and Bent-Fil structures had a line-of-sight velocity of ∼8–10 km/s, while the ring and Main-Fil structures were seen at velocities of ∼5–7km/s. If the Bent-Fil features are in fact being warped by expansion pressures from the ionization front, then it may be expected that these features have different velocities as compared to the dense structures which may be more shielded. <cit.> estimate an expansion velocity of 1–2 km/s for the dense molecular ring and and  expanding shells in the bipolar cavities with velocities of ∼5 km/s. However, expansion is only one explanation and there are others plausible reasons as to why dense ring and PDR regions have different line-of-sight velocities such as rotations, tidal forces, etc. If the magnetic field lines are indeed being altered by the radiation from the massive stars, then this may suggest that the magnetic field pressure may not be sufficient to the support the cloud structures against the kinetic energy injected by stellar feedback. 

While the Flipped-Fil also shows a strong preference for parallel alignment relative to the local magnetic field and similar line-of-sight velocities as the Bent-Fils, it is not as clear at this stage whether the Flipped-Fil is an irradiated structure associated with warped gas near the PDR. Unlike the Bent-Fils, the Flipped-Fil is not preferentially observed at shorter wavelength dust maps but rather, appears faintly in dust emission across the wavelengths 500–70 μm (see Fig. <ref> and Appendix <ref>, Fig. <ref>). Furthermore, the Flipped-Fil is not traced by the Spitzer maps (see Fig. <ref> and Appendix <ref>, Fig. <ref>), which may be expected if the structure was associated with warmer dust grains. A full discussion of the origins of the Flipped-Fil is presented in the next section. 



 §.§ Origins of the Flipped-Fil


One region of particular interest throughout this study has been the Flipped-Fil (labeled in Fig. <ref>) due to the N-S orientation of the magnetic field lines locally within the filament, which is in stark contrast with the general E-W orientation of the surrounding HAWC+ Band C magnetic field morphology. While the magnetic field lines appear to deviate slightly from the E-W trend in several regions such as that 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 2-D orientation, which may not be as drastic in 3-D. A change in the grain alignment mechanism of the dust grains could also cause the near-discontinuous behaviour 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 <cit.>. 

The change in magnetic field orientation within the Flipped-Fil can also be explained through physical origins. One plausible formation scenario was presented by <cit.> studying the magnetic field morphology of the Pillars of Creation seen in M16, which resembles the morphology of the Flipped-Fil. The scenario <cit.> 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 an over-density 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 length-wise 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 spectral line data of the Flipped-Fil is observed at line-of-sight velocities of ∼8–10 km/s, which is red-shifted compared to the northern and southern halves of the dense ring. This arrangement could have occurred if the expansion of the  region swept up the Flipped-Fil and pushed it behind the massive cluster such that it is currently at a further 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 2-D 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  and  integrated intensity maps, which traces irradiated dense gas. The lack of mid-infrared emission towards the Flipped-Fil is likely not due to absorption from foreground structures as the region is associated with  emission (see Fig. <ref>). Additionally, at a column density of  ∼ 13, the Flipped-Fil is traced by low and intermediate gas tracers such as ^12CO and ^13CO indicating that it is has molecular gas component, but is not quite dense or cold enough to be traced as clearly in N_2H^+ 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  and . 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 Fig. <ref>) could be the over-dense structures envisioned Figure 5 of <cit.>, 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 <cit.> who suggest from comparisons of numerical simulations by <cit.> 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 favoured explanation. Higher resolution infrared observations may help distinguish these structures better, giving more insight to their morphology and origin. 



 §.§ Energetic Balance


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 E-W 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. <cit.> estimate the ram pressure energy density within the ring (labeled in Fig. <ref>) to be u_ram = 0.41–3.67× 10^-10 erg cm^-3. Assuming equipartition, we set the ram pressure energy density u_ram equal to the magnetic energy density u_B so that:

    u_ram≥ u_B = B^2/8π(cgs),


and the 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 <cit.> 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 <cit.> and RCW 120 is estimated to have 100 μG <cit.>. 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. 

<cit.> also calculate that the turbulent energy density within the ring is u_turb= 4.1–5.1× 10^-10 erg cm^-3 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) <cit.>. 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  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 <cit.>. 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 <ref>). 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 towards 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. 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 dataset for future work. 



§ CONCLUSION


The motivation of this work was to better understand the combined influence of stellar feedback and magnetic fields on high-mass star formation. To do this, we targeted the extensively studied region RCW 36 in the Vela C giant molecular cloud, which has been previously observed using many different tracers. Adding to this suite of complimentary data, we presented new, higher resolution observations of the magnetic field morphology inferred from SOFIA/HAWC+ linearly polarized dust emission maps at 89 and 214 μm at filament-scales as well as ALMA Band 6 continuum 1.1-1.4 mm data at clump scales.

We then employed the Histogram of Relative Orientations (HRO) method to compare the orientation of the HAWC+ magnetic field to the orientation of physical structures in RCW 36 as traced by 7 spectral lines and dust emission and continuum maps ranging from 3.6 μm to 1.4 mm, for a multi-scale, multi-wavelength study. Comparing our HRO results to previous larger cloud-scale studies and simulations, we discussed the implications of our findings on the energetic importance of the magnetic field in RCW 36. The main conclusions of this analysis are:


    
  * We find that the inferred filament-scale magnetic field from HAWC+ generally matches the east-west morphology of the cloud-scale magnetic field inferred from BLASTPol, except for a few notable regions of interest. One exception we call the Flipped-Fil region, where the field switches to a roughly north-south orientation and the other exception we call the Bent-Fils region, where the field follows a bent shape around star forming clumps. We also find that the magnetic field morphologies inferred by Band C (89 μm) and Band E (214 μm) are highly similar, indicating that they may be tracing similar dust grain populations, scales and densities.
    
    
  * The HRO analysis between the inferred magnetic field and single intensity maps show differences in orientation between dense gas tracers and PDR tracers. Structures observed in dense gas tracers show a preference for perpendicular alignment relative to the magnetic field, whereas the tracers of warm dust and the PDR show a preference for parallel relative alignment. The aforementioned Flipped-Fil region, however, tends to be preferentially parallel in most tracers for which it is well detected, indicating that this is a special case.
    
    
  * Repeating the HRO analysis for different of line-of-sight velocities in the spectroscopic data cubes shows that the relative alignment of structures also varies with velocity. Structures associated with dense gas show a preference for perpendicular alignment relative to the magnetic field at line-of-sight velocities of 4–7 km/s, while structures associated with the PDR show a preference for parallel alignment at velocities of 8–11 km/s. This technique allows us to disentangle otherwise overlapping structures in the single integrated intensity map.  
    
    
  * The finding the dense ridge of RCW 36 is oriented perpendicular to the magnetic field is consistent with previous cloud-scale HRO studies of the Centre-Ridge within Vela C <cit.>. Comparing this result to studies which applied the HRO method to synthetic observations of MHD simulations <cit.> suggests that the magnetic field may have been dynamically important on cloud scales when the dense ridge of RCW 36 region first formed, however this may no longer be the case after the formation of the massive stars.     
    
    
  * The HRO results from the warm dust and PDR tracers suggest that the magnetic field lines are perhaps being altered near the ionization front such that they align parallel relative to gas warped by stellar feedback. This could indicate that ram pressure and radiation from the nearby massive cluster may be dominating the energetic balance on filament-scales. This is potentially causing the flux-frozen magnetic fields to be bent in directions which follow the elongation of the bright-rimmed Bent-fil structures. The parallel relative alignment observed for the Flipped-Fil may have resulted from a formation scenario similar to what has been suggested by <cit.> for the Pillars of Creation where gas bows around an over-density creating a pillar-like structure in which the local magnetic field is rotated orthogonally in comparison to the background magnetic field orientation. 


In conclusion, the SOFIA/HAWC+ polarization data provided new insights into the RCW 36 region, particularly regarding how the magnetic field may have been altered near the PDR region due to ionization from the massive stellar cluster. The filament-scale HRO analysis highlighted structures showing parallel alignment relative to the local magnetic field which were not observed in previous HRO cloud-scale studies. This altered magnetic field near the PDR may impact the formation of next-generation stars by influencing gas dynamics. Thus comparing the magnetic field from higher resolution, shorter wavelength polarization data to PDR tracers may offer useful insight when studying the impact of feedback on the magnetic field in other high mass star-forming regions as well.




§ ACKNOWLEDGMENTS

We thank the referee for their discerning feedback which has greatly improved the presentation of this work. This research has made use of data from the Herschel Gould Belt survey (HGBS) project (http://gouldbelt-herschel.cea.fr). The HGBS is a Herschel Key Programme jointly carried out by SPIRE Specialist Astronomy Group 3 (SAG 3), scientists of several institutes in the PACS Consortium (CEA Saclay, INAF-IFSI Rome and INAF-Arcetri, KU Leuven, MPIA Heidelberg), and scientists of the Herschel Science Center (HSC). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01003.S, ADS/JAO.ALMA#2021.1.01365.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This study was partly based on observations made with the NASA/DLR  SOFIA. SOFIA is jointly operated by the Universities Space Research Association Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI), under DLR contract 50 OK 0901 to the University of Stuttgart. upGREAT is a development by the MPIfR and the KOSMA/University Cologne, in cooperation with the DLR Institut für Optische Sensorsysteme. Financial support for FEEDBACK at the University of Maryland was provided by NASA through award SOF070077 issued by USRA. The FEEDBACK project is supported by the BMWI via DLR, project number 50 OR 2217 (FEEDBACK-plus).  The BLASTPol telescope was supported by through grant numbers NNX13AE50G, 80NSSC18K0481,
NAG5-12785, NAG5-13301, NNGO-6GI11G, NNX0-9AB98G, the Illinois Space Grant Consortium, the Canadian Space Agency, the Leverhulme Trust through the Research Project Grant F/00 407/BN, the Natural Sciences and Engineering Research Council of Canada, the Canada Foundation for Innovation, the Ontario Innovation Trust, and the US National Science Foundation Office of Polar Programs. L.M.F acknowledges support from the National Science and Engineering Research Council (NSERC) through Discovery Grant RGPIN/06266-2020, and funding through the Queen’s University Research Initiation Grant. G.N. is grateful for financial support for this work from NASA via award SOF06-0183 issued by USRA to Northwestern University. N.S. and R.S. acknowledge support by the SFB 1601, sub-project B2, funded by the DFG. T.G.S.P gratefully acknowledges support by the National Science Foundation under grant No. AST-2009842 and AST-2108989 and by NASA award #09-0215 issued by USRA.

SOFIA, ALMA, BLAST, APEX, Herschel, Spitzer, Mopra
Astropy <cit.>, NumPy <cit.>, Matplotlib <cit.>, Scipy <cit.>, APLpy <cit.>, Spectral-Cube <cit.> 





§ HRO RESULTS USING HAWC+ BAND E





In this section, we present the HRO Single Map results (method described in Section <ref>) using the HAWC+ Band E data to infer the magnetic field orientation. Table <ref> gives the corrected and uncorrected PRS values for the different tracers. We find that the general trend of the Single Map results from the HAWC+ Band E data are fairly similar to the results found for Band C (see Table <ref>). The consistency of the results is due to the similarity of the magnetic field morphologies traced in Band C and Band E (see Figure <ref>).



Figure <ref> shows the PRS () for the different single-dish dust map wavelengths. In both Band C and Band E HRO analyses, the resulting sign (positive or negative) of  as a function of dust wavelength is the same. The longer wavelength (500–70 μm) dust maps show a statistically significant (| Z_x| > 3) preference for perpendicular alignment, while the Spitzer maps show a statistically significant preference for parallel alignment. Similar to Band C, the Band E HRO results also show insignificant  values for the ALMA continuum data. Furthermore, the Band E Single Map HRO results for the column density, temperature, and atomic and molecular lines are all consistent with the Band C results. Most gas tracers show a preference for perpendicular alignment, with the exception of , ^12CO, and  which have an overall statistically insignificant . The only difference between Band C and Band E results is the magnitude of the  values. This is mostly due to the difference in resolution of the Band C and Band E data, since the magnitude of  depends on the number of independent samples (as discussed in Section <ref> and <ref>).

Thus the Band E HRO results similarly find that tracers which are mostly sensitive to the N-S dense molecular ring and Main-Fil features show a preference for perpendicular alignment, while the tracers mostly sensitive to the E-W Bent-Fil features show a preference for parallel alignment. The interpretations of the results for Band C presented in Section <ref> are therefore also applicable to Band E.



§ UNCERTAINTY ESTIMATION FOR THE PRS



In Section <ref>, we discussed how the oversampling-corrected Projected Rayleigh statistic  is expected to have an uncertainty of 1 <cit.>. These uncertainties do not, however, account for the measurement uncertainties in the magnetic field orientation  and the structure map M(x,y), which may contribute additional sources of error in the HRO analysis. In this section, we perform Monte Carlo tests to propagate measurement uncertainties to the , , and  calculations.

To estimate the impact of the uncertainties on the PRS due to  measurement uncertainties, we repeat the HRO analysis for 1000 magnetic field map iterations (+ ), where a magnetic field orientation error term  is added to the measured HAWC+ 89 μm (Band C) polarization angles. The error  is drawn from a normal distribution centered at 0 with a standard deviation equal to the polarization angle error, which is estimated from the HAWC+ Data Reduction Pipeline (discussed in Section <ref>). The uncertainty of the uncorrected statistic is then determined from the distribution of  values. We perform this test for two selected maps, the HAWC+ Band C intensity and the  integrated intensity map. We choose these maps since they have different alignment trends. The Band C intensity has a strong preference for perpendicular alignment (Z_x=-9.3), while  has  no clear statistical preference for parallel or perpendicular alignment (Z_x=0.5). For both maps, the Monte Carlo analysis with  uncertainties result in  and  distributions with standard deviations of σ_Z'_x=0.2 and σ_Z'_x=0.002. The mean of the distributions is the same values as the Single Map HRO  and  values, which did not include measurement uncertainties (given in Table <ref>).

To test whether measurement errors σ_M in the intensity maps M(x,y) could increase the uncertainty of our HRO analysis, we ran a Monte Carlo simulation which generates 1000 structure map iterations (M(x,y) + σ_M). We select the Band C intensity map for M(x,y) since the uncertainties in total intensity have been estimated by the HAWC+ Data Reduction pipeline. We then smooth the maps by the Gaussian gradient kernel (using the method described in Section <ref> with the same kernel sizes listed in Table <ref>) and calculate the relative angles for each iteration. We find that the standard deviation of the PRS values for this test to be to be σ_Z'_x=0.8 and σ_Z'_x=0.006. We note that the uncertainties in the statistics are slightly larger in comparison to the polarization angle error propagation. However, neither the magnetic field orientation uncertainties nor the polarization angle uncertainties have a large impact on the final results. 

The Monte Carlo tests that have been presented in this section thus far have assumed that each pixel samples the probability distribution function of  independently from neighbouring pixels. However, since the FWHM beam area spans many pixels, the measurement errors are correlated between adjacent pixels. To estimate the uncertainties on the oversampling corrected , we re-calculate the PRS using only independent relative angle pixels (i.e one pixel per FWHM beam area). Using this approach, our Monte Carlo test with  measurement uncertainties () gives a standard deviation of σ_Z_x=0.2 for both the Band C and  intensity maps, which is the same as the distribution in  found in the over-sampled case. The Monte Carlo test with M(x,y) measurement uncertainties (σ_M) gives a standard deviation of σ_Z_x=0.3 for Band C intensity.  In all cases,   , and  have standard deviations less than 1. Therefore, the uncertainties in the PRS statistics are primarily due to the distribution in the relative orientation angles sampled at different locations in the map, rather than by measurement uncertainties in the maps or inferred magnetic field orientation. 



§ SINGLE MAP HRO





  §.§.§ Dust Emission



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

Figure <ref> shows the results for the 350–70 μm dust maps. We see that all maps trace the east and west halves of the N-S oriented dense ring which contribute most of the 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  values for the longer wavelength dust maps in Table <ref>.



One similarity for all dust maps at 350–70 μm in Fig. <ref> is the alignment measured for the Flipped-Fil (labeled in Fig. <ref>). 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 N-S local magnetic field orientation, as seen in the middle and right panels of Fig. <ref>. This is in contrast to the otherwise E-W 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 Fig. <ref>), particularly the south Bent-Fil which begins to appear over the south-western 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 E-W elongation, resulting in increasingly parallel relative orientations at 160 μm and 70 μm, which decreases the overall magnitude of the  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 row of Figure <ref>. Similar to the 3.6 μm map shown in Fig. <ref>, the 4.5 μm map is also highly sensitive to the E-W Bent Fils structure, but does not trace the dense molecular ring. Since the Bent-Fils are oriented roughly parallel to the E-W magnetic field, this results in a statistically significant positive . In Section <ref>, we made note of the apparent correlation between the emission of the Spitzer 3.6 μm map and the  total integrated intensity. Here, we once again note the similarities between the Spitzer emission and the  total integrated intensity map (shown in the second row of Fig. <ref>), which is also a PDR tracer like . The HRO results for the  data are further discussed in Section <ref>. 



We now analyze the HRO results of the Herschel-derived column density and temperature maps, shown in Figure <ref>. Similar to longer wavelength dust maps, the column density map is also mostly sensitive to the dense molecular ring elongated in the roughly N-S oriented, particularly the western half containing the Main-Fil. This results in a majority of locally perpendicular ϕ(x,y) angles relative to the E-W magnetic field morphology, giving an overall large negative . 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 .

Figure <ref> 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 is also likely not probing the magnetic field within the dense cores detected by ALMA as is further discussed in Appendix <ref>. 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 which 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. 





  §.§.§ Spectral lines



Figure <ref> shows the Single Map HRO for the different molecular gas tracers. We compare these results with the atomic gas data of  in Figure <ref> and  in Figure <ref>. For the atomic gas, the  and  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  relative angle map is similar to . The main difference between the two is that  <cit.> traces more of the E-W oriented Bent-Fils in the west half of the ring as compared to the . This may be because the Bent-Fils features in the western half are more diffuse, while  tends to trace denser PDR regions <cit.> and hotter gas (typically ∼200 K). As such,  emission appears to better trace the N-S Main-Fil, resulting in a slight overall preference for perpendicular alignment as compared to .



Next, we examine the Single Map HRO results for the molecular line data shown in Figure <ref>. The high-density gas tracers such as HNC and N_2H^+ <cit.> clearly trace the densest N-S 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  value. The C^18O data  traces intermediate densities <cit.> similar to ^13CO, but the C^18O Mopra data has a lower signal-to-noise ratio and lower resolution, resulting in a lower  than the APEX ^13CO data. Interestingly we find that the south Bent-fil traced by ^13CO and HAWC+ 89 μm intensity (shown in Fig. <ref>), and can be somewhat seen in the maps of C^18O and HNC, but is not seen in the N_2H^+ map which tends to trace only high density and colder molecular gas. These observations are contrasted with the ^12CO 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, ^12CO 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). Though ^12CO is detected throughout the entire RCW 36 region, it is optically thick at the densest regions. Spectra of ^12CO, ^13CO,  and  for the Flipped-Fil are given in Appendix <ref> for reference (for spectra in other regions, see Fig. 2 of ).



  §.§.§ Caveats and Considerations


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 Fig. <ref>), 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 <ref>). 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 iso-contours. For example in Fig. <ref> 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 N-S direction where we would expect the gradient to change direction. This is most obvious for the 350 μm 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  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 .



§ VELOCITY DEPENDENT HRO



In this section, we show the velocity channel maps for the gas tracers not included in Figures <ref>–<ref>. For reference, Figure <ref> shows the spectra for ^12CO, ^13CO,  and  at the Flipped-Fil region. Figures <ref>–<ref> show the integrated intensities for 1 km/s wide velocity slabs from 4–10 km/s for ^12CO, , C^18O, and N_2H^+, respectively. Similar to the main text, we note that the dense ring is traced at line-of-sight velocities of 4–7 km/s, while the Bent-Fils and Flipped-Fil are seen at 8–10 km/s.















§ OPTICAL DEPTH OF HAWC+ DATA





Dust emission can be optically thin (τ≪ 1), such that the plane-of-sky magnetic field orientation is inferred from the average emission by all dust grains along the line of sight. Alternatively, the emission can be optically thick (τ > 1), such that only the flux from the outer surface layer of the cloud is traced. Understanding the optical depth helps identify the location of the dust grains dominating the emission, whether they are within a translucent cloud, or on the surface of an opaque cloud. To estimate the optical depth of RCW 36 at the HAWC+ wavelengths, 89 μm and 214 μm we use Equation <ref> from the main text.

The Herschel-derived column density (36 FWHM) map is used N_H_2. The 36 column density map is chosen over the 18 version as this resolution matches the temperature map. We use the dust opacity law from <cit.> that was applied by <cit.> to originally derive the Vela C column density map, which is:


    κ_λ = 0.1 ×(λ/300 μ m)^-β


where we have used β=2, to match <cit.>. The resulting estimates for optical depth τ are shown in Figure <ref>. From the colorbar, we see that at the longer HAWC+ wavelength of 214 μm (Band E), the dust emission is optically thin τ < 1 everywhere in the region. At the shorter wavelength in HAWC+ at 89 μm (Band C), we estimate the emission is roughly optically thin everywhere, except at the brightest peaks near the clumps, where τ∼ 1.4. While optically thick emission at an observing wavelength 89 μm would not be unexpected, in this case the maximum optical depth is still relatively close to the τ∼ 1 surface, indicating that the emission is only moderately thick, rather than very thick (e.g for say τ∼ 10). This means that we are missing some flux at the brightest peaks, but not too much.

However, we note that since we are using a 36 column density map to extrapolate the optical depth of the Band C emission, which has a native resolution of 8, the actual optical depth at 8 could be higher than what is estimated from 36 FWHM map (shown in Fig. <ref>) at the smallest scales closest to the brightest peaks. Therefore, an optical depth of τ∼ 1.4 at these points may be an underestimation. Another caveat is that a singular mass-weighted average dust temperature is assumed when generating the Herschel column density map <cit.>. This does not account for temperature gradients along the line of sight. If the 89 μm dust is probing a population of warmer dust grains, that only exist near the surface of the  region, then that dust will not be probing the entire column traced by the Herschel column density map. Furthermore, the emission at 89 μm may also be tracing very small dust grains (VSGs) <cit.>, which do not emit at sub-millimeter wavelengths but can at 70–100 μm <cit.>.  These grains are stochastically heated and are not in equilibrium, which makes inferring their properties difficult. The 160–500 μm emission used to generate the column density map are likely tracing emission from the larger dust grains, meaning the column density map derived from these wavelengths will not include VSGs. 



§ CORRELATION OF POLARIZED DUST EMISSION WITH OTHER TRACERS





In this section, we compare the 89 μm polarized emission to the other dust emission, column density, temperature and molecular line maps listed in Table <ref>. These tracers probe different physical properties of the gas, and a strong correlation between the polarized emission and a particular dataset may imply the magnetic field is primarily being traced in regions with similar density, temperature, chemical and excitation conditions. 

We do this by individually overlaying contours of the different tracers on the Band C total polarization intensity map and visually comparing the emission.  Figure <ref> shows that the 89 μm polarized emission correlates well with the  integrated intensity, shown in the left panel. This further reinforces our previous assertion that the polarization data is preferentially tracing the magnetic field from the warm dust located near the PDR, where [CII] is abundant. This is contrasted to the apparent anti-correlation of the Band C polarized data observed for the ALMA ACA continuum data, as can be seen in the right panel of Fig. <ref>. The ALMA clumps appear to be located in areas where there is a lack of polarized intensity. This finding is consistent with Band C emission being sensitive to warmer dust, rather than the cold dense structures traced by ALMA. Polarization measurements at longer wavelengths and higher resolution (e.g. with ALMA) would be needed to probe the magnetic field within these colder dense structures.



aasjournal
