We use polarization data from SOFIA HAWC to investigate the interplay between magnetic fields and stellar feedback in altering gas dynamics within the highmass starforming region RCW 36 located in Vela C. This region is of particular interest as it has a bipolar H II 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 data sets 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 the photodissociation 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 lineofsight 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 fluxfrozen gas structures causing the observed preference for parallel relative alignment.
['Star formation', 'Star forming regions', 'Interstellar medium', 'Interstellar magnetic fields', 'H II regions', 'Photodissociation regions', 'Molecular clouds', '1569', '1565', '847', '845', '694', '1223', '1072', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Solar and Stellar Astrophysics']
Magnetic Field Alignment Relative to Multiple Tracers in the Highmass Starforming Region RCW 36
{'Magnetic Field Alignment Relative to Multiple Tracers in the High-mass Star-forming Region RCW 36': "<!-- image -->  \nROBERT SIMON  \n1 Department of Physics, Engineering Physics and Astronomy, Queen's University, Kingston, ON K7L 3N6, Canada 2 SOFIA Science Center, NASA Ames Research Center, Moffett Field, CA 94045, USA 3 I. Physik. Institut, University of Cologne, Z¨ulpicher Str. 77, D-50937 Cologne, Germany Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), 1800 Sherman Avenue, Evanston, IL 60201, USA 5 Engineering + Technical Support Group, University of Chicago, Chicago, IL 60637, USA 6 Department of Physics & Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA 7  \n4  \nMIT Haystack Observatory, 99 Millstone Road , Westford, MA, 01827  \n8 School of Physics, University of New South Wales, Sydney NSW 2052, Australia", 'ABSTRACT': 'Weuse 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 H II 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.', '1. INTRODUCTION': "Observations and simulations suggest that star formation occurs when density fluctuations undergo gravitational collapse in molecular clouds (McKee & Ostriker 2007). 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 prestellar structures by providing support against collapse and/or directing gas flow (Pattle et al. 2023). 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 (e.g., Whitworth 1979; Gritschneder et al. 2009; Lopez et al. 2011; Chuss et al. 2019).  \nCorresponding author: Akanksha Bij  \na.bij@queensu.ca  \nStellar feedback has also been observed to reshape magnetic field geometries around expanding ionized bubbles (e.g., Heiles 1989; Soler et al. 2018; Tahani et al. 2019), which is consistent with predictions from magnetohydrodynamic (MHD) simulations (e.g. Krumholz et al. 2007). 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 (e.g., Dale et al. 2014; Geen et al. 2020). Furthermore measuring the magnetic field strength through observations is challenging. While numerous observational techniques such as Zeeman splitting (Crutcher 2012; Crutcher & Kemball 2019), Faraday rotation (Tahani et al. 2018) and the Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951; Chandrasekhar & Fermi 1953) applied to polarized light have been used in the past (e.g., Girart et al. 2006; Pillai et al. 2015; Planck Collaboration et al. 2016),  \neach technique has limitations and/or only provides partial information about the magnetic field structure.  \nOf 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 (Davis & Greenstein 1951; Andersson et al. 2015). 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 starforming process (e.g., Goldsmith et al. 2008; Tassis et al. 2009; Li et al. 2013; Planck Collaboration et al. 2016; Pillai et al. 2020).  \nA numerical method known as the Histogram of Relative Orientations (HRO) was developed by Soler et al. (2013) 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 (e.g., Soler & Hennebelle 2017; Pillai et al. 2020; Seifried et al. 2020). Most notably, Planck Collaboration et al. (2016) 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) (e.g., Soler & Hennebelle 2017; Kortgen & Soler 2020).  \nThe 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 (Fissel et al. 2016). The BLASTPol-inferred magnetic field was compared to column densities derived from Herschel observations by Soler et al. (2017), and to the integrated intensities of molecular gas tracers by Fissel et al. (2019). 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  \nof ∼ 3 ' , the BLASTPol observations are only able to probe the Vela C magnetic field geometry on cloud-scales ( > 1 pc).  \nIn this work, we extend these HRO studies to filamentscales ( ∼ 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 (Harper et al. 2018). 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 A V > 100 mag) which is within a parsec of a ionizing young ( ∼ 1 Myr) OB cluster responsible for powering a bipolar H II nebula within the region (Minier et al. 2013; Ellerbroek et al. 2013).  \nPrevious 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 H II 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 (Hollenbach & Tielens 1997). Several recent H II region studies have used observations of [C II] since it traces the PDR (e.g., Schneider et al. 2018; Anderson et al. 2019; Pabst et al. 2019; Luisi et al. 2021; Beuther et al. 2022; Tram et al. 2023). For RCW 36, the bipolar H II region was investigated by Bonne et al. (2022) who examined the kinematics of [C II] 158 µ mand [O I] 63 µ mdata from the SOFIA legacy project FEEDBACK (Schneider et al. 2020).  \nWe 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 2 describes the observations. Section 3 details the physical structure of the RCW 36 and its magnetic field morphology, including noteworthy regions. Section 4 discusses the HRO method, with the results presented in Section 5. In Section 6, we interpret the results and compare this work to other studies. Finally in Section 7, we summarize our main conclusions.", '2. 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 2.1), thermal continuum emission from dust (Section 2.2), and spectroscopic lines (Section 2.3), each of which are summarized in Table 1. The observations from SOFIA HAWC+ and ALMA are presented for the first time. Figure 1 shows the new HAWC+ and ALMA data and Figure 2 shows archival data. A detailed overview of the each dataset is provided in the following subsections.  \n2.1. Dust Polarization Maps  \nFigure 1. Far-infrared and millimeter intensity maps of the RCW 36 region presented for the first time in this work. Top Left: SOFIA HAWC+ 89 µ m(Band C) total intensity in Jy pixel -1 . Top Right: SOFIA HAWC+ 214 µ m(Band E) total intensity in Jy pixel -1 . Bottom Left: ALMA ACA 1.1-1.4 mm (Band 6) continuum in Jy beam -1 Bottom Right: ALMA 12m 1.1-1.4 mm (Band 6) continuum in Jy beam -1 . All color maps are in linear scale. A 0.25 pc scale bar is shown on the bottom right. The FWHM beam size is shown on bottom left (beam and pixel sizes are given in Table 1). The gray and cyan contours show Herschel -derived column densities at values of 1.5 × 10 22 and 4.7 × 10 22 cm -2 , respectively. These contours are shown as a position reference on all maps of RCW 36 presented in this work.  \n<!-- image -->  \nTable 1. Summary of all datasets used for Single Map HRO analysis described in section 4.2.1  \nThe linearly polarized intensity P can be found from the linear polarization Stokes parameters Q and U using:  \nP = √ Q 2 + U 2 , (1)  \nwhile 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 (Rice 1945; Serkowski 1962). This can be corrected for with 'debiasing' (Wardle & Kronberg 1974), using:  \nP db = √ Q 2 + U 2 -1 2 ( δQ 2 + δU 2 ) , (2)  \nwhere δQ and δU are the measurement uncertainties in Q and U , respectively. The polarization angle ˆ E can be calculated using ˆ E = 0 . 5 arctan( U, Q ) . The orientation of the plane-of-sky magnetic field ˆ B ⊥ is then inferred to be orthogonal to ˆ E such that:  \nˆ B ⊥ = 1 2 arctan( U, Q ) + π 2 , (3)  \nwhere an angle of 0 · points towards Equatorial north and increases eastward.", '2.1.1. SOFIA HAWC+': "In this paper we publish, for the first time, observations of RCW36using publicly available archival data from HAWC+ (Harper et al. 2018), 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 (described in Hildebrand et al. 2000) 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 1.  \nTo reduce this data, we used the HAWC+ Data Reduction Pipeline which is described in detail in Santos et al. (2019) and Lee et al. (2021) and summarized here as follows. The pipeline begins by demodulating the data and discarding any data points affected by erroneous telescope movements or other data acquisition errors. This demodulated data 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 2) 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  \nof uncertainties such as noise correlated across pixels which can lead to an underestimation of errors in the IQU maps (Davidson et al. 2011; Chapman et al. 2013). 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  \nENF = C 0 √ 1 + ( I 5 I 0 ) 2 , (4)  \nwhere 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 2. We note that the Stokes I errors for Band E were a special case where the Pipeline ENF fitting routine failed (giving a value of I 0 = 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.  \nAfter 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.  \nAfter running the pipeline, we also applied a 3σ signal-tonoise threshold on the total intensity flux and polarized flux to further remove noisy polarization vectors. As a final diagnostic we checked for potential contamination from the reference beam position due to dithering the data in Chop-Nod polarization observations (following the method described in the Appendix section of Novak et al. 1997). For this test we used Herschel far-IR intensity maps (described in Section 2.2.2) since they cover both the RCW 36 region and the surrounding Vela C cloud-scale region, which include the HAWC+ chop reference beam positions that are outside the HAWC+ maps. We used Herschel maps of comparable wavelengths to the HAWC+ data (PACS 70 µ m to compare with HAWC+ 89 µ m and PACS 160 and SPIRE 250 µ m to compare with HAWC+ 214 µ m) and found the ratio of the total intensity of the HAWC+ region compared to the chop reference beam region in the Herschel data. To estimate the ratio of polarization flux, we conservatively assume that there is a 10% polarization at the reference beam positions and remove points where the estimated polarized flux in the reference beam is more than 1/3 times the polarized flux in the HAWC+ map.", '2.2. Dust Emission Maps': "2.2.1. ALMA  \nIn this work, we present two new interferometric data sets from the Atacama Large Millimeter/submillimeter Array (ALMA) and Atacama Compact Array (ACA).  \nFigure 2. Maps of RCW 36 derived from different infrared and spectral line tracers. Panels a -i share the same RA and Dec axes, for which a scale bar is given in the bottom right of panel i . Panels j -l share another set of RA and Dec axes, for which a scale bar is given in panel l . The FWHMbeam size is shown on the bottom left (beam and pixel sizes are given in Table 1). For BLASTPol, the beam size is 2.5 ' and pixel size is 4.6 '' . The contours show the same Herschel -derived column densities as Fig. 1. The Herschel far-IR flux maps in panels b and c are plotted with a log scale color map, the rest of the color maps are in linear scale. Panels e -i and k -l are integrated over the velocity ranges specified in Column 5 of Table 3.  \n<!-- image -->  \nTable 2. Excess Noise Factor (ENF) fitting parameters I 0 and C 0 found by the HAWC+ Data Reduction Pipeline for Stokes I , Q , U Maps in Band C and E (see Eq. 4).  \nThe 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 Zucker et al. (2020). 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.  \nThe 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.  \nIn 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. 1. 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.", '2.2.2. Herschel SPIRE and PACS': "To study the cloud structures probed by thermal dust emission, we use publicly available archival maps from the Her-  \nschel Space Observatory, which observed Vela C on 18 May 2010 (Hill et al. 2011) as part of the Herschel OB Young Stars (HOBYS) key programme (Motte et al. 2010). The observations were conducted using the SPIRE instrument at 500, 350, and 250 µ m (with FWHM angular resolutions of 36 '' , 25 '' , and 18 '' , respectively; Griffin et al. 2010), and the PACS instrument at 70 and 160 µ m (with resolutions of 8 '' and 13 '' , respectively; Poglitsch et al. 2010; Giannini et al. 2012). 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 Palmeirim et al. 2013. We show the 250 µ m, 70 µ m, and temperature maps in panels b, c, and d of Figure 2, respectively.", '2.2.3. 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 1 . Spitzer observed RCW 36 in May 2006, employing the fourchannel 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 (Fazio et al. 2004). IRAC uses two 5.2 ' × 5.2 ' fields of view, where one field simultaneously images at 3.6 and 5.8 µ mand 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 (Ellerbroek et al. 2013). 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 µ mmap in panel j of Figure 2.", '2.3. 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 3 summarizes the lines of interest, including their transitions, rest frequencies, velocity resolution, and the channels used to make the integrated intensity maps. For our analysis of the spectroscopic data, we use both an integrated intensity map for a wide velocity range (Column 5 of Table 3) as well as channel maps over narrower velocity ranges (Column 6 of Table 3). Our spectral data cube analysis is described in Sections 4.2.1 and 4.2.2.', '2.3.1. SOFIA upGREAT FEEDBACK': "To trace the PDR, we use a 158 µ m map of the [C II] 3 P 3 / 2 → 3 P 1 / 2 transition (at native resolution of 14.1 '' ) and  \nTable 3. Summary of Spectroscopic Data Cubes  \nThe systemic velocity for the Centre-Ridge is ∼ 7 km/s (Minier et al. 2013).  \na 63 µ m map of the [O I] 3 P 1 → 3 P 2 transition (at native resolution of 6.3 '' ) from the SOFIA FEEDBACK C + legacy survey (Schneider et al. 2020). The survey was conducted by the upGREAT (upgraded German REceiver for Astronomy at Terahertz frequencies) heterodyne spectrometer (Risacher et al. 2018) onboard the SOFIA aircraft (Young et al. 2012) 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 Schneider et al. (2020). In this work, we use the data reduced by Bonne et al. (2022) who smoothed the [C II] map to 20 '' and [O I] 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 [C II] and ∼ 0.8 -1.5 K for [O I] (Bonne et al. 2022). The integrated intensity maps for [C II] and [O I] and have been integrated from -20 to +20 km/s are shown in panels k and l of Fig. 2, respectively.", '2.3.2. APEX LAsMA': "To trace the molecular gas regions in RCW 36, we use observations of 12 CO (3-2) and 13 CO (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 (Gusten et al. 2006). The maps were scanned in total power on-thefly mode and are sized 20 ' × 15 ' , with a beamsize of 18.2 '' at 345.8 GHz. Bonne et al. (2022) 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 12 CO and 13 CO are shown in panels e and f of Fig. 2, respectively.  \nIn 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 18 O as well as HNC and N 2 H + which were originally presented in Fissel et al. (2019). The C 18 O 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 2 H + observations used overlapping 5 ' square raster maps. The data reduction procedure, performed by Fissel et al. (2019), 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 18 O, and 36 '' and ∼ 0.2 km/s for both HNC and N 2 H + . The integrated intensity maps for HNC, N 2 H + and C 18 O have been integrated over the velocity range given in Column 5 of Table 3 and are shown in panels g-i of Fig. 2, respectively.", '3. 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 3 showcases various continuum and polarization observations for RCW 36 that will be used to describe its general structure in the following subsections.', '3.1. 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 Hill et al. (2011) using Herschel data. The densest and most prominent of the ridges is the Centre-Ridge, with column densities of A V > 30 , and a length of roughly ∼ 10 pc (Hill et al. 2011). The Center-Ridge contains the RCW 36 region. It has a bipolar nebula morphology (Minier et al. 2013) with two fairly symmetric lobes oriented in the east-  \nFigure 3. Three color-composite images of the RCW 36 region on different scales. All images use either linear intensity or log intensity to highlight emission structures. The white star marker in all panels shows the location of the brightest O9 V star in RCW 36 (Object ID number 1 from Minier et al. 2013, and number 462 from Bik et al. 2005, 2006). The cyan contours represent the same Herschel -derived column densities as Figure 1. Left: Cloud-scales . The red-blue-green colours are the same for both top and bottom panels with Herschel /PACS 160 µ m (red), Herschel /PACS 70 µ m (green), and SuperCOSMOS H α (blue). The top panel vectors show the magnetic field orientation inferred from BLASTPol 500 µ m polarized data. The bottom panel has labels for the bipolar nebula (green dotted), ring (yellow dashed) and ionized gas shell (blue dashed). Middle: Filament-scales The red-blue-green colours are the same for both top and bottom panels with APEX 12 CO integrated intensity(red), APEX 13 CO integrated intensity (green), and SOFIA [C II] integrated intensity (blue). Note that the [C II] emission does not cover the northern part of the image. The top panel vectors show the magnetic field orientation inferred from HAWC+ 214 µ m(Band E) polarized data. The bottom panel labels the Main-Fil (white dashed) and Flipped-Fil (pink dotted). The Main-Fil contains five dense starforming clumps outlined by black contours showing ALMA ACA continuum at an intensity of 0.018 Jy beam -1 . The Flipped-Fil corresponds to the region where the magnetic field orientation appears to abruptly change by almost 90 · . Right: Clump-scales The red-blue-green colours are the same for both top and bottom panels with ALMA ACA 1.1-1.4 mm continuum (red), Spitzer 3.6 µ mintensity (green), and APEX 13 CO integrated intensity (blue). Top panel vectors show the magnetic field orientation inferred from HAWC+ 89 µ m (Band C) polarized data. The Bottom panel shows the North and South Bent-Fils (yellow-dotted). Each of the lower panels also shows a scale bar.  \n<!-- image -->  \nst direction that are traced well by the green PACS 70 µ m emission in the lower-left panel of Figure 3 (see also the dotted green ellipses). This bipolar nebula is roughly centered around a young (1.1 ± 0.6 My) massive cluster with two late-type O-type stars and ∼ 350 members (Ellerbroek et al. 2013). The position of the most massive star (spectral type O9V) is indicated by a white star-shaped marker in Fig. 3. The ionizing radiation from this cluster is powering an expanding H II gas shell traced by H α (shown in blue, Fig. 3, left panel) (Minier et al. 2013).  \nBipolar H II regions are of great interest because, though they have been observed in other high-mass star forming regions such as S106 (Schneider et al. 2018) and G319.88+00.79 (Deharveng et al. 2015), they seem to be  \nmore rare than single H II bubbles (e.g., Churchwell et al. 2006; Deharveng et al. 2010; Anderson et al. 2011; Kendrew et al. 2012; Samal et al. 2018).  \nWithin the bipolar cavities, Bonne et al. (2022) identified blue-shifted [C II] shells with a velocity of 5.2 ± 0.5 km/s, likely driven by stellar winds from the massive cluster. Additionally, they find diffuse X-ray emission (observed with the Chandra X-ray Observatory) in and around the RCW 36 region which is tracing hot plasma created by the winds. However, Bonne et al. (2022) estimate that the energy of the hot plasma is 50-97% lower than the energy injected by stellar winds and reason that the missing energy may be due to plasma leakage, as has been previously suggested for RCW 49 (Tiwari 2021).  \nThe magnetic field geometry on > 1 pc cloud scales is traced by the BLASTPol 500 µ m polarization map (Fissel et al. 2016) 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. 3, which follow a fairly uniform east-west morphology that is mostly perpendicular to the orientation of the dense ridge. However, around the north and south 'bends' of the bipolar structure, the magnetic field lines also appear to bend inward towards the center, following the bipolar shape of the structure.", '3.2. Filament-Scale': "The middle panels of Figure 3 highlight structures on filament scales of ∼ 0 . 1-1 pc and the cyan contours in all panels represent Herschel -derived column densities to show the filament. At the waist of the bipolar nebula, Minier et al. (2013) 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. 3). The majority of the dense material, as traced by the column density contours, is contained within this ring. Minier et al. (2013) 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. HII we use archival H α data from the SuperCOSMOS H-alpha Survey (SHS; Hambly et al. 2001; Parker et al. 2005). From the SuperCOSMOS map, we note that the eastern side of the ring is seen in H α absorption, signifying that it is in front of the ionizing gas and associated massive star cluster. Whereas, H α emission is seen across the western region of the ring and therefore this part of the ring is likely behind the cluster.  \nThe 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 3). Hill et al. (2012) 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 3.  \nSeveral diffuse filamentary structures traced by 12 CO and [C II] (shown in red and blue, respectively in middle panel of Fig. 3) can also be observed in the ambient cloud surrounding the ring. Bonne et al. (2022) have reasoned that due to the curved shape of these filaments, they are not part of the larger expanding [C II] shells but may have been low density filaments originally converging towards the centre dense ridge (similar to the converging filaments seen in Musca, B211/3, and DR 21; Goldsmith et al. 2008; Schneider et al. 2010; Palmeirim et al. 2013; Cox et al. 2016; Bonne et al. 2023) that have instead been swept away at velocities > 3 km/s due to stellar feedback.  \nThe magnetic field geometry on filament scales, traced by SOFIA/HAWC+ 214 µ m is fairly consistent with the E-W  \ngeometry 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 3, 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 12 CO (red, middle panel), [C II] (blue, middle panel), and 13 CO emission (blue, right panel).", '3.3. Core-Scale': "The right panels of Figure 3 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 MainFil is clumpy with several bright rims, voids, and pillar-like structures identified by Minier et al. (2013) (see their Figure 3). This matches our ALMA band 6 continuum observations as shown by the five near-round clumps and associated elongated pillars, as seen in the right panel Fig. 3.  \nMinier et al. (2013) recognized that the bright rims appear near the end of the pillar-like structures. The bright-rims are traced in the right panel of Fig.3 by Spitzer /IRAC 3.6 µ m emission, which mainly traces hot dust, found at the edges of the PDR (Minier et al. 2013). These rims appear to wrap around the cold ALMA clumps, without covering them completely, in a manner resembling bow-shocks (though actual bow-shocks are unlikely in this region). These bright rims are of great interest in this work and are therefore collectively labeled as 'Bent-Fils' as they will be referred to in later sections. There is a prominent northern Bent-Fil and southern Bent-Fil shown by the two yellow-dotted ovals in the lower right panel of Fig. 3. The curved morphology of the Bent-Fils is noted by Minier et al. (2013) to be likely due to tracing the inner border of the dense ring which is being progressively photoionzed by the star cluster. Interestingly, the HAWC+ magnetic field morphology seems to follow the Bent-Fil features, deviating once again from the general E-W cloud scale magnetic field.", '4.1. Histogram of Relative Orientations': 'In this section we discuss the procedure of the Histogram of Relative Orientation (HRO) method (see Soler et al. 2013; Soler & Hennebelle 2017, for a more detailed description) which computes the relative angle ϕ ( x, y ) between the gradient vector field of a structure map M ( x, y ) and the planeof-sky magnetic field orientation ˆ B ⊥ at each pixel. The steps for this procedure are outlined in the following subsections.', '4.2. 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 4.2.1), we compare the orientation of local structure at every location using one map M ( x, y ) , for each of  \nthe datasets listed in Table 1, to the magnetic field orientation measured by dust polarization as was done by previous HRO studies (Planck Collaboration et al. 2016; Soler et al. 2017; Fissel et al. 2019; Lee et al. 2021, e.g.). In the second approach, applied only to the spectroscopic data cubes listed in Table 3, we slice the spectral line cube into multiple velocity slabs v i and compare the orientation of structures in the integrated intensity of each slab M ( x, y ) i to the inferred magnetic field. This quantifies the relative alignment as a function of line-of-sight velocity, and will therefore be referred to as the Velocity Dependent HRO (see 4.2.2).', '4.2.1. Single Map HRO': 'The dust emission, column density, and temperature maps are already in the 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 Fissel et al. 2019). 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  \nM ( x, y ) = ∫ v 1 v 0 T R dv. (5)  \nThe velocity ranges used to calculate the integrated intensity map in the Single Map HRO analysis for each spectral cube are specified in Column 5 of Table 3.', '4.2.2. Velocity Dependent HRO': 'Additionally for the molecular line data, we also perform a Velocity Dependent HRO analysis. Here we slice a selected velocity range v 0 -v 1 (specified in Column 6 of Table 3) into narrower velocity slabs with a width of 1 km/s. We increment the center velocities v i of each slab by 0.5 km/s, such that v i = { v 0 +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:  \nM ( x, y ) i = ∫ v i +0 . 5 v i -0 . 5 T R dv. (6)  \nThe 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 3) such that enough velocity channels are included in the integrated intensity map. This ensures that there is sufficient signal-to-noise in each slab and small local fluctuations are averaged over. The HRO analysis is then repeated for each M ( x, y ) i in the set.', '4.3. Projection and Masking': 'To directly compare the structure map M ( x, y ) and the plane-of-sky magnetic field map ˆ B ⊥ 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  \ndo 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 ˆ B ⊥ map then M ( x, y ) is projected on the ˆ B ⊥ map). The pixel sizes of each dataset are given in Table 1 for comparison. When the ˆ B ⊥ data has lower resolution of the two, rather than directly projecting the orientation of the magnetic field lines inferred from Eq. 3, we instead project the Stokes Q and U intensity maps separately and then recalculate the inferred ˆ B ⊥ . This avoids an incorrect assignment of vector orientation angles to the newly sized pixels.  \nNext, 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 ˆ B ⊥ 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.', '4.4. 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:  \nψ = arctan ( δM/δx δM/δy ) , (7)  \nwhere ψ is calculated at each pixel ( x, y ) . The partial derivatives are calculated by convolving M ( x, y ) with Gaussian derivative kernels G, using  \nδM δx = M ( x, y ) ∗ δ δx G ( x, y ) , (8)  \nand similarly δ M /δy = M( x, y ) ⋆ δG ( x, y ) /δy . This reduces noise and avoids erroneous relative angle measurements due to map pixelization. The size of the Gaussian kernels in angular units θ G is chosen to be one third of the FWHM angular resolution θ beam of the M ( x, y ) map, using θ G = θ beam / 3 . If this kernel size θ G is less than 3 pixels, then a minimum kernel size of 3 pixels is used instead. A summary of all the kernel sizes in angular units θ G and pixel units l G is provided in Columns 5 and 6, respectively, of Table 1. The same smoothing lengths listed in Table 1 for the molecular line data are applied for the Velocity Dependent HRO analysis.  \nThe relative angle ϕ ( x, y ) between the iso-contour direction ψ ( x, y ) and the plane-of-sky magnetic field ˆ B ⊥ can then be computed with:  \nϕ ≡ arctan ( | ψ × ˆ B ⊥ | ψ · ˆ B ⊥ ) . (9)  \nThe 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  \nwrapped 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.', '4.5. 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 Jow et al. 2018). The PRS is a modified version of a classic Rayleigh statistic which tests for a uniform distribution of angles using a random walk. The classic Rayleigh statistic characterizes the distance Z from the origin if one were to take unit-sized steps in the direction determined by each angle. Given a set θ i of n angles within the range [0 · , 360 · ] , this distance Z can be calculated as follows:  \nZ = (Σ n i cosθ i ) 2 +(Σ n i sinθ i ) 2 n , (10)  \nwhere 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, Jow et al. (2018) modify this statistic by calculating only the horizontal displacement Z ' x in the hypothetical random walk:  \nZ ' x = Σ n i cosθ i √ n/ 2 . (11)  \nNow a parallel relative angle ϕ i = 0 will map to θ = 2 ϕ = 0 and give a positive cos (0) = 1 contribution to Z ' x , while a perpendicular relative angle ϕ i = π/ 2 will map to θ = 2 ϕ = π and give a negative cos ( π ) = -1 contribution to Z ' x . Therefore, a statistic of Z ' x ≫ 0 indicates strong parallel alignment, while Z ' x ≪ 0 indicates strong perpendicular alignment. Since the value of Z ' x is within the range [ -√ 2 n, + √ 2 n ] , the statistic can be normalized by √ 2 n to give a measure of the degree of alignment:  \n˜ Z ' x = Z ' x √ 2 n , (12)  \nwhere a value of ˜ Z ' x = ± 1 would correspond to perfectly parallel or perpendicular alignment.  \nIf the n data points are independent, then Z ' x 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 Z ' x is proportional to the number of data points n 1 / 2 and the relative alignment ϕ is measured at each pixel ( x, y ) , oversampling within the map can result in a misleadingly large Z ' x magnitude. To determine the statistical significance of Z ' x we follow the methodology in Fissel et al. (2019) and correct for oversampling by repeating the HRO analysis on 1000 independent white noise maps M WN , smoothed to the same resolution as M ( x, y ) , and compared to the magnetic field orientation.  \nThe white noise maps M WN are generated to be the same size as the real data ( M ( x, y ) ), using Gaussian noise with a mean and standard deviation of M ( x, y ) . In these maps the orientation of the gradient will be a uniformly random distribution. The white noise M WN map then follows the same projection procedure that was applied to M ( x, y ) in Section 4.3 followed by the same mask which was applied to the real data M ( x, y ) . We calculate the corresponding PRS, Z ' 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:  \nZ x = Z ' x σ Z ' WN . (13)  \nAfter correction, the error on the corrected PRS Z x is now σ Z x = 1 , and a magnitude of | Z x | > 3 is considered statistically significant. The number of independent samples can be estimated as  \nn ind = n ( σ Z ' WN ) 2 . (14)  \nFor 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.", '5.1. 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 4 shows the relative angles between the Band C ( ˆ B ⊥ C ) and Band E ( ˆ B ⊥ E ) magnetic field vectors at each pixel, calculated using Equation 9. We find that the relative angles are near parallel ϕ ∼ 0 · at most locations in RCW 36, signifying that the magnetic field orientations in the two different bands are highly consistent. This gives a mean relative angle of ⟨ ϕ ⟩ = 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 rela-  \nFigure 4. Map of the relative angles ϕ ( x, y ) measured between the magnetic field orientation inferred from HAWC+ Band C (89 µ m) and Band E (214 µ m). The line segments show the plane-of-sky magnetic field orientation inferred from SOFIA HAWC+ Band C polarization data.  \n<!-- image -->  \nssion of this result is given in Section 6.1.2.  \nSince the magnetic field orientations in the two bands are very similar, we present only the Band C HRO results in the main text and defer the Band E results to Appendix A. The Band C Single Map and Velocity Dependent HRO results are presented in Sections 5.2 and 5.3, respectively.", '5.2. Single Map HRO Results': "Table 4 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 ( Z x ), indicating a statistical preference for perpendicular alignment. There are also some notable exceptions that have a positive PRS, indicating a preference for parallel alignment. We discuss the results from the various tracers below.  \nThe single-dish dust emission maps show a distinct variation in the magnitude of the PRS values with wavelength. The top panel of Figure 5 shows the oversampling-corrected Z x values, which indicate the statistical significance of the PRS, while the bottom panel shows the normalized uncorrected ˜ Z ' x values, which indicate the degree of alignment. Both panels show that Z x and ˜ Z ' x are negative for the submmand far-IR Herschel and SOFIA data indicating a preference for perpendicular alignment, and positive for the mid-IR Spitzer data, indicating a preference for parallel alignment. All the maps have statistically significant values Z x (i.e. | Z x | > 3 σ Z x , where σ Z x = 1 ). A notable trend is seen for the normalized statistic where the ˜ Z ' x values roughly increase  \nTable 4. PRS Results for the Single map HRO Analysis using magnetic field orientation inferred from HAWC+ 89 µ mdata.  \nNotes: The structure map used for each molecular line is a single intensity map integrated over the velocity range specified in Table 1 using Equation 5. A negative Z x corresponds to an overall preference for perpendicular alignment, while a positive value corresponds to a preference for parallel alignment. The larger the magnitude of Z x , the stronger the statistical significance of the preferred alignment.  \n- a PRS corrected for oversampling using Eq. 13\n- b The standard deviation for 1000 white noise runs used to correct oversampling in PRS from Eq. 13\n- c The uncorrected PRS from Eq. 11\n- d The normalized uncorrected PRS from Eq. 12\n- e The number of relative angle points used to calculate the uncorrected PRS.  \n(i.e becomes less negative) as the wavelength decreases, suggesting that successively more structures within the maps align parallel to the magnetic field at shorter wavelengths. In comparison, the trend for the oversampling-corrected statistic Z x decreases from 500-250 µ mandpeaks in magnitude at 214 µ m. This is because the oversampling correction factor ( σ Z ' WN ) is proportional to the number of independent samples in the masked area and lower Z x values are expected for the same degree of alignment at longer wavelengths with larger beams, where there are fewer pixels over the same area (see Equation 14). We also ran Monte Carlo simulations to  \ntest 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 Z x and ± 0.002 for ˜ Z ' x . These tests and a discussion of our error propagation methods are described in Appendix B.  \nFigure 6 summarizes the oversampling-corrected and normalized PRS for the total integrated intensity spectral line maps. Unlike the Z x values for the dust map shown in Fig. 5, many of the spectral line intensity maps show no overall preference for alignment, or only show a statistically insignificant alignment trend. These maps show different alignment preferences relative to magnetic field in different regions, or overlapping along the line of sight. In Section 5.3, we will show that some of these structures that overlap in the integrated intensity map can be decomposed into different line of sight velocity channels. In some cases, particularly for the Mopra observations, the low Z x values are also in part due to lower resolution and higher noise levels of the spectroscopic data. Overall, we note that all gas tracers have a negative ˜ Z ' x , signifying a preference for perpendicular alignment, with the exception of [C II] which has a positive ˜ Z ' x .  \nIn Figure 7 we identify which structures are aligned with the magnetic field for a select number of datasets. Similar plots for the remaining datasets in Table 4 are shown and discussed in Appendix C, Figures 14-16. The right column shows the structure map M ( x, y ) overlaid with the magnetic field orientation. The middle column shows the relative angle ϕ ( x, y ) calculated at each location in the region, where purple ( ϕ ∼ 0 · ) is associated with local parallel alignment and orange ( ϕ ∼ 90 · ) is associated with local perpendicular alignment relative to the magnetic field. The left column summarizes the alignment trend using a histogram of the relative angles.  \nFrom Fig. 7, 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 3. This is consistent with the visual observation that the ring is elongated approximately along the north-south direction which is oriented roughly perpendicular to the mostly eastwest 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 3). This structure is traced at the shorter 89 µ mwavelength, 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 ˜ Z ' x 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 Ap-  \nFigure 5. The PRS results from the Single Map HRO analysis which compares the orientation of elongated structures in dust maps of varying observing wavelengths relative to the magnetic field orientation inferred from SOFIA/HAWC+ 89 µ m. The top panel shows the Z x values which have been corrected for oversampling (using Equation 13), and indicate the statistical significance of a preference for parallel ( Z x > 0 ) or perpendicular alignment ( Z x < 0 ) of map structures relative to the magnetic field. The Z x values have errors of 1, as shown by the error bars. The bottom panel shows the normalized uncorrected ˜ Z ' x values which indicate the degree of alignment. A maximum value of ˜ Z ' x = +1 corresponds to complete parallel alignment while ˜ Z ' x = -1 corresponds to complete perpendicular alignment.  \n<!-- image -->  \nFigure 6. Same as Figure 5 but now showing the total integrated intensity spectral line maps.  \n<!-- image -->  \nFigure 7. HRO results for selected maps. Right: Intensity maps (labeled in the left column text) that are projected onto the HAWC+ 89 µ m grid. The colorbar is units of MJy sr -1 for the Herschel and Spitzer maps, Jy pixel -1 for 89 µ m map, and K km/s for the integrated intensity [CII] and 13 CO maps. The Herschel 500 µ m and SOFIA 89 µ m maps are in log-scale, while the rest are linear-scale. The vectors show the orientation of the magnetic field inferred from HAWC+ 89 µ m polarized emission, where detections below 3 σ have been masked. Middle: Spatial distribution of relative angles ϕ (x,y), sharing the same RA/Dec axes as the right column. Only pixels where the inferred magnetic field orientation is not masked have ϕ (x,y) values. Contours for right and middle panel show column densities of 1.5 × 10 22 and 4.7 × 10 22 cm -2 (same as Fig 1). Left: Histogram density for relative angles in 5 · wide bins. The color map for angles is the same as the middle panel colorbar, where purple signifies a parallel relative angle of ϕ ( x, y ) = 0 · and orange signifies a perpendicular relative angle of ϕ ( x, y ) = 90 ·  \n<!-- image -->  \npendix C.0.1, Fig. 14), where the emergence of the southern Bent-Fil structure at wavelengths < 160 µ m results in less negative Z x values as observed in Fig. 5.  \nHowever 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 Z x .  \nThe 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 BentFils structures (see full discussion in Appendix C.0.1 and Figure 16).  \nExamining the molecular gas maps, we find that 13 CO, which is sensitive to intermediate density gas, is able to trace both the dense molecular ring and the south Bent-Fil structure, resulting in a Z x value comparable to the 89 µ m dust map.  \nThe [C II] relative angle map shows that the E-W BentFil 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 [C II] also does not peak at near-parallel or nearperpendicular angles, but rather close to ϕ ∼ 40 · . Though it is not statistically significant, the positive Z x result of [C II] is interesting as it is in contrast with the negative Z x results for all other of all other spectral line tracers, which predominantly trace the molecular dense ring (Fig. 17). Furthermore, we see that the emission in the integrated intensity [C II] map correlates with the Spitzer emission, which probes warm dust likely found near PDRs. It is therefore noteworthy that the Z x results for both the [C II] and Spitzer maps are positive, indicating that structures associated with PDRs have a preference towards parallel alignment relative to the magnetic field.  \nIn 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 [C II]) 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 C.0.3. In the next section, we discuss the HRO analysis for different line-ofsight velocity ranges in the spectral line cubes. We use this Velocity Dependent HRO approach to examine the relation-  \nhip between the orientation of the different line-of-sight gas structures and the magnetic field.  \nFigure 8. Corrected PRS Z x as a function of velocity for different molecular lines.  \n<!-- image -->", '5.3. 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 4.2.2. Figure 8 shows the corrected PRS results for different spectral lines as a function of velocity.  \n13  \nFigure 9. The integrated intensities of 13 CO (3-2) for 1 km/s wide velocity slabs (labeled on the top-left of each panel) from 4-10 km/s. The maps have been projected onto the HAWC+ Band C grid. The color-scale indicates integrated intensity in units of K km s -1 . The contours are the same as Figure 1. The vectors show the magnetic field orientation inferred from HAWC+ 89 µ mdata  \n<!-- image -->  \nWe note that while the magnitude of the Z x is not always statistically significant ( < 3σ ) over the velocity range for all tracers, the overall trend is interesting and consistent with the Single Map HRO results. The intermediate and dense gas tracers, C 18 O, HNC, N 2 H + typically have statistically significant negative Z x values, especially around 5-8 km/s, implying a preference for perpendicular alignment. This velocity range matches the mean line-of-sight velocity of the cloud of around 7 km/s (Minier et al. 2013; Bonne et al. 2022). In contrast, the [C II] PRS results switch from negative Z x values around 5-6 km/s to statistically significant positive Z x values 9-11 km/s, indicating a preference for parallel alignment at higher line-of-sight velocities.  \nFigure 8 also demonstrates the limitations of using a single integrated intensity map for a spectroscopic cube in the HROanalysis, as was done in Section 5.2. Since there can be multiple overlapping elongated structures at different line-ofsight 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  \nvelocities, Figures 9-11 show channel emission maps from 4-10 km/s.  \nIn the 13 COchannel maps (shown in Fig. 9), we notice that emission from the northeastern half of the ring structure is most prominent at ∼ 4-6 km/s, while the southwestern section of the ring is most prominent at ∼ 6-8 km/s. Since the ring including the Main-Fil is oriented N-S, approximately perpendicular to the E-W magnetic field, the overall Z x at these velocities is preferentially perpendicular and therefore negative. However, the area of the northeastern region of the ring is smaller and contains fewer HAWC+ Band C polarization detections, leading to less ϕ ( x, y ) ∼ 90 · pixels at 4-5 km/s compared to larger southwestern component of the ring, resulting in a less negative Z x . At line-of-sight velocities > 8 km/s, the gas traces the Flipped-Fil and north Bent-Fil, causing a weak preference for parallel alignment around 10 km/s.  \nSimilarly, in the [C II] channel maps (Fig. 10) we notice that the eastern half of the ring can be seen in the 4-5 and 5-6 km/s maps, followed by the western half of the ring at 68 km/s along with the southern Bent-Fil, all of which result in an overall negative and therefore preferentially perpendicular, negative Z x for these velocity channels. At higher ve-  \nFigure 10. Same as Fig. 9 but for SOFIA [C II] data.  \n<!-- image -->  \nlocities we see emission from the Flipped-Fil and northern Bent-Fil at 8-10 km/s, resulting in an overall positive Z x in these velocity channels. Since the Bent-Fil features are the more prominent in [C II] than other gas tracers, it has the largest positive Z x magnitude.  \nIn contrast, HNC (Fig. 11) which is tracing denser gas, mostly shows emission tracing the dense ring (eastern half at 4-6 km/s and western half at 6-8 km/s) has a mostly negative Z x . The channel maps for the rest of the molecular lines are shown in Appendix D, Figures 19-21, all of which show emission from the dense ring from 4-8 km/s, and the FlippedFil and Bent-Fil features from 8-10 km/s.  \nThis switch from negative to positive PRS as a function of line-of-sight velocity is consistent with our Single Map HRO findings. In Section 5.2 we noted a bimodal trend in the PRS, where high column density gas tracers showed a preference of perpendicular relative alignment while tracers associated with the PDR and warmer dust showed more parallel relative orientations. From the Velocity Dependent HRO analysis, we now learn that the dense gas and PDR structures are also kinematically distinct, such that the same PRS bimodality is also observed as a function line-of-sight velocity. In the next section, we interpret these results and suggest potential physical mechanisms that may be causing the observed trends.', '6. 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 6.1, followed by an interpretation of the HRO results and discussion of the origins of the Flipped-Fil structure in Sections 6.2 and 6.3, and finally comment of the energetic balance of the region in Section 6.4.', '6.1. 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 6.1.1 and compare the polarization data and magnetic field morphologies at the different HAWC+ wavelengths in Section 6.1.2.  \nMopra HNC  \nFigure 11. Same as Fig. 9 but for Mopra HNC data.  \n<!-- image -->  \n6.1.1. Optical Depth of Dust Emission  \nTo 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 E 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 E) as in previous HOBYS and Herschel Gould Belt Survey (HGBS: Andr'e et al. 2010) studies (e.g. Motte et al. 2010; Roy et al. 2014). 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 Hill et al. (2011) which has angular resolution of 36 '' , and is different than the 18 '' column density map listed in Table 1 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 Hildebrand (1983) and spectral index is also consistent with Hill et al. (2011) and Roy et al. (2014).  \nWe estimate that the optical depth for 214 µ m (Band E) is optically thin τ ν ≪ 1 everywhere within RCW 36 (see Fig. 23 in Appendix E). 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 E 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.", '6.1.2. 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 5.1, we used the HRO method to statistically show that the magnetic field  \nFigure 12. The total polarized intensity measured for HAWC+ Band C (89 µ m) divided by the total polarized intensity measured for HAWC+ Band E (214 µ m). The ratio of intensity P C /P E is shown by the colorbar, where a ratio ∼ 1 corresponds to where the intensities are roughly equal. The Band E polarized intensity was projected on to the Band C grid, and the Band C data was smoothed to the Band E resolution. The Band C polarized emission detections below a 3 σ signal-to-noise cutoff have been masked.  \n<!-- image -->  \nmorphologies inferred from Band C and Band E are almost identical.  \nThis 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 starforming region like RCW 36, are likely being heated from the radiation of the massive stellar cluster. This warm dust is therefore probably located near the H II expanding gas shell and associated PDR. This is also supported by the observation that dust polarized intensity of Band C appears to correlate with the PDR-tracing [C II] and anti-correlate with the ALMA ACA continuum which traces cold dense cores (as shown in Appendix F). Therefore the HAWC+ magnetic field is likely weighted towards the surface of the cloud near the PDR, rather than the colder denser dust structures.  \nAside 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  \nconsistent magnetic field morphology can be expected at the different wavelengths if the dust emission is optically thin, as previously suggested.  \nOne noteworthy difference between the Band C and Band E datasets, however, was found by comparing the total polarized intensity given by Eq. 1 in Band C ( P C ; smoothed to the resolution of Band E) to Band E ( P E ). Figure 12 shows that the ratio of P C /P E is close to unity for the majority of RCW 36 except for certain regions. These regions have higher polarized intensities in the Band C map than they do in Band E by a factor of ∼ 2-4 . Interestingly, the regions also overlap with where the HAWC+ magnetic field is seen to deviate from the general 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 Vaillancourt et al. (2008) and Santos et al. (2019) who also compared the polarization ratio in different bands. Radiative alignment torques, which are thought to be responsible for the net alignment of the dust grains with their short axes parallel to the magnetic field, require anisotropic radiation fields from photons of wavelengths comparable or less than the grain size (Andersson et al. 2015). 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.', '6.2.1. Preferentially Perpendicular Alignment for Dense Tracers': "In Section 5.2, 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) (Soler et al. 2017; Fissel et al. 2019). 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.  \nSoler et al. (2017) showed that for both the entire Vela C cloud and the Centre-Ridge region (containing RCW 36), the  \nrelative 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, Fissel et al. (2019) compared the orientation of the magnetic field inferred from BLASTPol 500 µ mto 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 (e.g., Planck Collaboration et al. 2016, at 10 ' FWHM).  \nIn 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.  \nWe also compare our work to magnetic field studies done on comparable small-scales. Kaminsky et al. (2023) find that the Musca filament is oriented roughly perpendicular to the surrounding magnetic field morphology, as traced by SOFIA/HAWC+ 214 µ m observations. Moreover, Lee et al. (2021) 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.  \nPossible interpretations of such processes have been explored by comparing observations to simulations. For instance, Soler et al. (2013) 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 Chen et al. (2016) reason that the preferentially parallel relative alignment occurs in magnetically dominated (sub-Alfv'enic) gas, while preferentially perpendicular relative alignment occurs in turbulence dominated (super-Alfv'enic) 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 Soler &Hennebelle (2017), 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  \nstrong 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 Seifried et al. (2020) showed that the relative orientation trends also strongly dependent on viewing angle.  \nBased 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, Minier et al. (2013) 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 H II 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.", '6.2.2. Parallel Alignment for PDR Tracers': 'Section 5.2 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 | Z x | 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. 3), 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 µ mmaps, which are sensitive to emission from hot dust found near the PDR.  \nThe 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 fluxfrozen 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  \nBent-Fil regions traced by hot dust and [C II] and [O I] PDR tracers.  \nFurthermore, the Velocity Dependent HRO results (see Section 5.3) 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. Bonne et al. (2022) estimate an expansion velocity of 1-2 km/s for the dense molecular ring and and [C II] 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.  \nWhile 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. 7 and Appendix C.0.1, Fig. 14). Furthermore, the Flipped-Fil is not traced by the Spitzer maps (see Fig. 7 and Appendix C.0.1, Fig. 15), 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.', '6.3. Origins of the Flipped-Fil': "One region of particular interest throughout this study has been the Flipped-Fil (labeled in Fig. 3) 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.  \nObservational 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 neardiscontinuous 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 (Le Gouellec et al. 2023).  \nThe change in magnetic field orientation within the Flipped-Fil can also be explained through physical origins. One plausible formation scenario was presented by Pattle et al. (2018) studying the magnetic field morphology of the Pillars of Creation seen in M16, which resembles the morphology of the Flipped-Fil. The scenario (detailed schematically in Figure 5 of Pattle et al. 2018) is summarized here. An ionization front fueled by photon flux from a massive radiating star or cluster is envisioned to approach molecular gas, which may have regions of varying density. The gas being dragged by the ionization front may bow-around the an overdensity to form an elongated pillar. The flux-frozen magnetic field lines within the pillar would then follow the gas motion and end up perpendicular relative to the background magnetic field orientation. Such a structure could remain stable as the compressed magnetic field lines would provide support against radial collapse since gas flow perpendicular to the field lines would be inhibited. The pillar may gradually erode in the length-wise direction however, as gas flow parallel to the field lines would still be allowed.  \nWhile 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 H II 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 2Dprojection 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 FlippedFil is directly associated with the PDR as it is not seen in the Spitzer maps which trace hot dust, but is seen in the [C II] and [O I] integrated intensity maps, which traces irradiated dense gas. The lack of mid-infrared emission towards the FlippedFil is likely not due to absorption from foreground structures as the region is associated with H α emission (see Fig. 3). Additionally, at a column density of A V ∼ 13 , the FlippedFil is traced by low and intermediate gas tracers such as 12 CO and 13 CO indicating that it is has molecular gas component, but is not quite dense or cold enough to be traced as clearly in N 2 H + and HNC.  \nThe Flipped-Fil thus shows clear differences from the Bent-Fils which are likely associated with the warm dust structures and traced by shorter wavelength maps (3.5-160 µ m) as well as PDR tracers such as [C II] and [O I]. The Pillars of Creation formation scenario suggested for the Flipped-Fil can also be applicable to the Bent-Fils. In this picture, the star-forming clumps seen in ALMA ACA continuum data (see Fig. 3) could be the over-dense structures envisioned Figure 5 of Pattle et al. (2018), around which the bright rimmed Bent-Fil structures are being bowed around. The orientation of the bow-shapes then may suggest the di-  \nection of these ionization fronts. This model is similar to Minier et al. (2013) who suggest from comparisons of numerical simulations by Tremblin et al. (2012a,b) 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.  \nWhile 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.", '6.4. Energetic Balance': "In this section we examine the energetic balance of the RCW36 region in light of the new HAWC+ polarization observations presented in this work.  \nThe 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. Bonne et al. (2022) estimate the ram pressure energy density within the ring (labeled in Fig. 3) 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:  \nu ram ≥ u B = B 2 8 π (cgs) , (15)  \nand 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 Kusune et al. (2016) for the Centre-Ridge using the Davis-Chandrasekhar-Fermi method. Furthermore, our estimate is also lower in comparison to similar high-mass star-forming regions. For instance, DR 21 is measured to have a magnetic field strength of 130 µ G(Koley et al. 2021) and RCW 120 is estimated to have 100 µ G (Chen et al. 2022). Since our upper-limit is crude and based on the assumption that the feedback is ram pressure dominated, we may be underestimating the magnetic field strength.  \nBonne et al. (2022) 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'enic conditions (where u B ≥ u turb ) (Stone et al. 1998). If the turbulent energy was dominant, we would expect the magnetic field orientation to have more random vari-  \nns, which would decrease any alignment trend and result in Z x values with smaller magnitudes than our current measurements. Alternatively, the effects of turbulence on the magnetic field morphology may not be visible on the spatial scales probed by SOFIA/HAWC+ if the size of the turbulent eddies are smaller than the size of the beam, such that the polarization component from turbulent motion cancels out along the line-of-sight to give the appearance of low dispersion, as demonstrated in Jones (1989). On filament-scales, the ordered magnetic field observations from HAWC+ suggest a near-equipartition between the magnetic field energy and turbulent energy, with the ram pressure from stellar feedback dominating in certain regions.  \nThis interpretation is different from the large-scale HRO results using BLASTPol, which suggested that the magnetic field may have been dominating the energetic balance, setting a preferred direction of gas motion on cloud scales for the dense Centre-Ridge to form preferentially perpendicular relative to the magnetic field (see Section 6.2.1). This indicates that the dynamic importance of the magnetic field may be scale dependent in this region or that the energetic balance has changed since the formation of the original generation of stars which is currently driving the feedback within RCW 36. It should also be noted that, since the filament-scale magnetic field traced by HAWC+ Band C is weighted 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.", '7. 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 mmdata at clump scales.  \nWe 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 µ mto 1.4 mm, for a multiscale, 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  \nimportance of the magnetic field in RCW 36. The main conclusions of this analysis are:  \n- 1. 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.\n- 2. 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.\n- 3. Repeating the HRO analysis for different of line-ofsight 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 811 km/s. This technique allows us to disentangle otherwise overlapping structures in the single integrated intensity map.\n- 4. 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 (Soler et al. 2017; Fissel et al. 2019). Comparing this result to studies which applied the HRO method to synthetic observations of MHD simulations (e.g., Chen et al. 2016; Kortgen & Soler 2020) 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.\n- 5. 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 ra-  \ndiation 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 Pattle et al. (2018) 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.  \nIn 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 cloudscale 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 fur 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  \nvia 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/062662020, and funding through the Queen's University Research Initiation Grant. G.N. is grateful for financial support for  \nthis 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.  \nFacilities: SOFIA, ALMA, BLAST, APEX, Herschel, Spitzer, Mopra  \nSoftware: Astropy (Astropy Collaboration et al. 2013, 2018, 2022), NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020), APLpy (Robitaille & Bressert 2012), Spectral-Cube (Ginsburg et al. 2019)", 'A. HRO RESULTS USING HAWC+ BAND E': 'Table 5. Same as Table 4 but using the magnetic field inferred from SOFIA/HAWC+ Band E (214 µ m).  \nIn this section, we present the HRO Single Map results (method described in Section 4.2.1) using the HAWC+ Band E data to infer the magnetic field orientation. Table 5 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  \nthe results found for Band C (see Table 4). The consistency of the results is due to the similarity of the magnetic field morphologies traced in Band C and Band E (see Figure 4).  \nFigure 13. Same as Fig. 5 but showing HRO results using HAWC+ Band E (214 µ m) polarization data to infer the magnetic field  \n<!-- image -->  \nFigure 13 shows the PRS ( Z x ) for the different single-dish dust map wavelengths. In both Band C and Band E HRO analyses, the resulting sign (positive or negative) of Z x 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 Z x 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  \ntracers show a preference for perpendicular alignment, with the exception of [C II], 12 CO, and [O I] which have an overall statistically insignificant Z x . The only difference between Band C and Band E results is the magnitude of the Z x values. This is mostly due to the difference in resolution of the Band C and Band E data, since the magnitude of Z x depends on the number of independent samples (as discussed in Section 5.2 and 4.5).  \nThus 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 6.2 are therefore also applicable to Band E.', 'B. UNCERTAINTY ESTIMATION FOR THE PRS': "In Section 4.5, we discussed how the oversamplingcorrected Projected Rayleigh statistic Z x is expected to have an uncertainty of 1 (Jow et al. 2018). These uncertainties do not, however, account for the measurement uncertainties in the magnetic field orientation ˆ B ⊥ 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 Z x , Z ' x , and ˜ Z ' x calculations.  \nTo estimate the impact of the uncertainties on the PRS due to ˆ B ⊥ measurement uncertainties, we repeat the HRO analysis for 1000 magnetic field map iterations ( ˆ B ⊥ + σ ˆ B ⊥ ), where a magnetic field orientation error term σ ˆ B ⊥ is added to the measured HAWC+ 89 µ m(Band C) polarization angles. The error σ ˆ B ⊥ 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 2.1.1). The uncertainty of the uncorrected statistic is then determined from the distribution of Z ' x values. We perform this test for two selected maps, the HAWC+ Band C intensity and the [C II] 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 [C II] has no clear statistical preference for parallel or perpendicular alignment ( Z x = 0 . 5 ). For both maps, the Monte Carlo analysis with ˆ B ⊥ uncertainties result in Z ' x and ˜ Z ' x 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 Z ' x and ˜ Z ' x values, which did not include measurement uncertainties (given in Table 4).  \nTo 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  \n4.3 with the same kernel sizes listed in Table 1) 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 Z x results.  \nThe Monte Carlo tests that have been presented in this section thus far have assumed that each pixel samples the probability distribution function of ˆ B ⊥ 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 Z x , 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 ˆ B ⊥ measurement uncertainties ( σ ˆ B ⊥ ) gives a standard deviation of σ Z x = 0 . 2 for both the Band C and [C II] intensity maps, which is the same as the distribution in Z ' x 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, Z x Z ' x , and ˜ Z ' x 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.", 'C. SINGLE MAP HRO': "C.0.1. Dust Emission  \nIn this section, we show the remainder of the relative angle maps and histograms from the Band C (89 µ m) Single Map HRO analysis that were not presented in Figure 7.  \nFigure 14 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 ˜ Z ' x values for the longer wavelength dust maps in Table 4.  \nOne similarity for all dust maps at 350-70 µ m in Fig. 14 is the alignment measured for the Flipped-Fil (labeled in Fig. 3). All relative angle maps find ϕ ∼ 0 · within the Flipped-Fil structure, even though the region is not always particularly noticeable in the intensity maps. This result of a mostly parallel relative alignment within the Flipped-Fil is consistent with the visual observation of the filament being elongated approximately in the direction of the N-S local magnetic field orientation, as seen in the middle and right panels of Fig. 3. This is in contrast to the otherwise E-W orientation of the  \nFigure 14. Same as Fig. 7. The right column is showing dust intensity with a log-scale colorbar in units of MJy sr -1 for the 350 µ m and 250 µ mmaps, and Jy pixel -1 for the 214, 160 and 70 µ mmaps.  \n<!-- image -->  \noverall magnetic field morphology in the surrounding RCW 36 region.  \nThe main difference between the different wavelengths is the emergence of the bright-rimmed Bent-Fil structures (labeled in Fig. 3), 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 Z x to be less negative than at 350-214 µ m. This trend can be also be seen in the left column HROs, which show a decreasing histogram density near ϕ ∼ 90 · for the shorter wavelengths.  \nWe contrast these results to the Spitzer data at 4.5 µ m, shown in the two row of Figure 15. Similar to the 3.6 µ m map shown in Fig. 7, the 4.5 µ mmap 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 Z x . In Section 5.2, we made note of the apparent correlation between the emission of the Spitzer 3.6 µ m map and the [C II] total integrated intensity. Here, we once again note the similarities between the Spitzer emission and the [O I] total integrated intensity map (shown in the second row of Fig. 15), which is also a PDR tracer like [C II]. The HRO results for the [O I] data are further discussed in Section C.0.2.  \nWe now analyze the HRO results of the Herschel -derived column density and temperature maps, shown in Figure 15. Similar to longer wavelength dust maps, the column density map is also mostly sensitive to the dense molecular ring elongated in the roughly 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 Z x . In contrast, the temperature map shows only a slight preference for perpendicular relative alignment as it appears to be mostly tracing the bipolar morphology of the region. The HAWC+ magnetic field follows the curvature of the bipolar nebula and thus results in more parallel relative angles between the magnetic field and structures in the temperature map and a lower magnitude Z x .  \nFigure 16 shows the results for the ALMA data. The HRO analysis for both the ALMA 12-m and ACA maps finds a statistically insignificant PRS of Z x ∼ 0 , implying that the structures traced by ALMA do not have a preferred direction of orientation relative to the HAWC+ Band C inferred magnetic field. ALMA being an interferometer, resolves out much of the large scale structure such as the dense ring, Main-Fil, Flipped-Fil and Bent-Fils, which are the main features observed by the other dust maps. The HAWC+ Band C data is also likely not probing the magnetic field within the dense cores detected by ALMA as is further discussed in Appendix F. This may explain why there is no correlation between the structures traced by the ALMA data and the  \nmagnetic 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 ALMApolarization mosaics. This is outside the scope of this paper.", 'C.0.2. Spectral lines': 'Figure 17 shows the Single Map HRO for the different molecular gas tracers. We compare these results with the atomic gas data of [O I] in Figure 15 and [C II] in Figure 7. For the atomic gas, the [C II] and [O I] data from SOFIA probe the transition from molecular to ionized gas in the PDR. The regions of parallel (purple) and perpendicular (orange) relative orientation angle in the [C II] relative angle map is similar to [O I]. The main difference between the two is that [C II] (with a critical density for collisional deexcitation of ∼ 2 . 6 × 10 3 cm -3 ; Rollig et al. 2006) traces more of the E-W oriented Bent-Fils in the west half of the ring as compared to the [O I]. This may be because the BentFils features in the western half are more diffuse, while [O I] tends to trace denser PDR regions (critical density of 5 × 10 5 cm -3 ; Rollig et al. 2006) and hotter gas (typically ∼ 200 K). As such, [O I] emission appears to better trace the N-S MainFil, resulting in a slight overall preference for perpendicular alignment as compared to [C II].  \nNext, we examine the Single Map HRO results for the molecular line data shown in Figure 17. The high-density gas tracers such as HNC and N 2 H + (critical densities > 10 4 cm -3 ; Shirley 2015; Fissel et al. 2019) 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 Z x value. The C 18 O data traces intermediate densities ( ∼ 10 3 cm -3 ; Fissel et al. 2019) similar to 13 CO, but the C 18 OMopra data has a lower signalto-noise ratio and lower resolution, resulting in a lower Z x than the APEX 13 CO data. Interestingly we find that the south Bent-fil traced by 13 CO and HAWC+ 89 µ m intensity (shown in Fig. 7), and can be somewhat seen in the maps of C 18 Oand HNC, but is not seen in the N 2 H + map which tends to trace only high density and colder molecular gas. These observations are contrasted with the 12 CO integrated intensity map which traces even lower density molecular gas and shows very different structure compared to the other spectral lines. While it also maps parts of the ring, 12 CO shows bright emission around the Flipped-Fil and the north Bent-Fil which are elongated along the direction of the magnetic field lines resulting in an overall weak preference for perpendicular alignment relative to the magnetic field ( Z x = -2 . 2 ). Though 12 CO is detected throughout the entire RCW 36 region, it is optically thick at the densest regions. Spectra of 12 CO, 13 CO, [C II] and [O I] for the Flipped-Fil are given in  \nFigure 15. Same as Fig. 7. The right column colobar is in units of Jy pixel -1 for the Spitzer map, K km/s for the SOFIA [O I] map, cm -2 for the column density map and K for the temperature map. All colorbars are linear-scale.  \n<!-- image -->  \nFigure 16. Same as Fig. 7 but for ALMA data. The right column is showing a log-scale colorbar in Jy beam -1 units for both continuum maps. Only ALMA data points above a 3σ signal-to-noise threshold are used for the HRO analysis.  \n<!-- image -->  \nAppendix D for reference (for spectra in other regions, see Fig. 2 of Bonne et al. 2022).', 'C.0.3. 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. 7), the Gaussian kernel smoothing removes some of the relative angle measurements near the south edge of the map boundary which is not the case for the Herschel 70 µ mwhich covers a larger area on the sky (see Figure 14). Another caveat to note for all our HRO analysis plots, is that some of the ϕ ∼ 0 · relative angle points are due to the gradient amplitude approaching zero as the gradient changes direction at the peak of the iso-contours. For example in Fig. 7 at the center of the highest column density contour within the Main-Fil, a thin row of purple ϕ ∼ 0 · pixels can be seen along the 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 Z x is insignificant. Furthermore it should also be mentioned that in addition to the hot dust detected by Spitzer , the instrument is also clearly detecting starlight from the massive stellar cluster. Since this emission is not extended but rather from point sources, it is  \ngenerally not elongated in a particular direction and therefore also does not significantly affect Z x .', 'D. VELOCITY DEPENDENT HRO': 'In this section, we show the velocity channel maps for the gas tracers not included in Figures 9-11. For reference, Figure 18 shows the spectra for 12 CO, 13 CO, [C II] and [O I] at the Flipped-Fil region. Figures 19-22 show the integrated intensities for 1 km/s wide velocity slabs from 4-10 km/s for 12 CO, [O I], C 18 O, and N 2 H + , 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.', 'E. 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 6.1.1 from the main text.  \nFigure 17. Same as Fig. 7. The right column is showing a integrated intensity in units of K km/s using a linear-scale colorbar for all maps.  \n<!-- image -->  \nThe 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 Hildebrand (1983) that was applied by Hill et al. (2011) to originally derive the Vela C column density map, which is:  \nκ λ = 0 . 1 × ( λ 300 µm ) -β (E1)  \nwhere we have used β = 2 , to match Hill et al. (2011). The resulting estimates for optical depth τ are shown in Figure 23. 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  \nFigure 18. Spectra showing the antenna temperature of 12 CO, 13 CO, [C II] [O I], for the Flipped-Fil region.  \n<!-- image -->  \nwould 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.  \nHowever, 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. 23) 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 (Roy et al. 2014). 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 H II 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) (Greenberg 1968), which do not emit at sub-millimeter wavelengths but can at 70-100 µ m(e.g., Draine & Anderson 1985; Walterbos & Schwering 1987). 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.", 'F. 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 1. 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.  \nWe do this by individually overlaying contours of the different tracers on the Band C total polarization intensity map and visually comparing the emission. Figure 24 shows that the 89 µ m polarized emission correlates well with the [C II] 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. 24. 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.  \nFigure 19. Same as Fig. 9 but for APEX 12 CO data.  \n<!-- image -->  \nFigure 20. Same as Fig. 9 but for SOFIA [O I] data.  \n<!-- image -->  \nFigure 21. Same as Fig. 9 but for Mopra C 18 O data.  \n<!-- image -->  \n+  \nFigure 22. Same as Fig. 9 but for Mopra N 2 H + data.  \n<!-- image -->  \nFigure 23. The optical depth τ (shown by the colorbar) of dust emission at 89 µ m ( left ) and 214 µ m ( right ), where τ < 1 corresponds to optically thin emission and τ > 1 corresponds to optically thick emission. The optical depth is estimated from Herschel -derived column density and temperature maps at 36 '' .  \n<!-- image -->  \nFigure 24. The colormap shows the polarized flux for HAWC+ Band C (89 µ m) in both panels. Left: shows the Band C polarized flux overlaid with SOFIA [C II] integrated intensity contours at the levels of 300, 400, 450, 500 K km s -1 . Right: shows the Band C polarized flux overlaid with ALMA ACA Band 6 1.1-1.4 mm continuum contours at the level of 0.018 Jy beam -1 .  \n<!-- image -->", 'REFERENCES': 'Tassis, K., Dowell, C. D., Hildebrand, R. H., Kirby, L., &  \nVaillancourt, J. E. 2009, MNRAS, 399, 1681,  \ndoi: 10.1111/j.1365-2966.2009.15420.x  \n- Tiwari, M. 2021, in American Astronomical Society Meeting Abstracts, Vol. 53, American Astronomical Society Meeting Abstracts, 137.03\n- Tram, L. N., Bonne, L., Hu, Y., et al. 2023, ApJ, 946, 8, doi: 10.3847/1538-4357/acaab0  \nTremblin, P., Audit, E., Minier, V., Schmidt, W., & Schneider, N.  \n2012a, A&A, 546, A33, doi: 10.1051/0004-6361/201219224  \n- Tremblin, P., Audit, E., Minier, V., & Schneider, N. 2012b, A&A, 538, A31, doi: 10.1051/0004-6361/201118031'}