© The Authors 2024

1 Introduction

The question of what controls the star formation efficiency in molecular clouds has long been at the center of star formation research. Early studies (Zuckerman & Evans 1974) showed that if all the gas within dense interstellar clouds were to collapse freely under self-gravity, the star formation rate in the Milky Way would be two orders of magnitude higher than the observed rate of 2 M_⊙ yr⁻¹ (Robitaille & Whitney 2010). Theories proposed to explain such a low star formation rate invoke the presence of turbulence or magnetic fields supporting interstellar clouds against gravitational collapse. In some models (Tan et al. 2006), turbulent or magnetic pressure gradients are strong enough to maintain clouds in approximate hydrostatic equilibrium on all spatial scales for tens of free-fall times. In other models, high-density regions are rapidly contracting, converting a large fraction (as high as 40%; Bonnell et al. 2011) of their mass into stars in only a few free-fall times. However, the low-density parts of the cloud, which contain up to 90% of the cloud mass (Battisti & Heyer 2014), disperse over the scales of free fall time as a result of turbulence generated by galactic shear or by energy input from supernova explosions (Dobbs et al. 2011; Walch & Naab 2015), or due to radiative and mechanical stellar feedback from high-mass stars that formed early on during the evolution of the cloud (Murray 2011; Colín et al. 2013; Pabst et al. 2020; Chevance et al. 2022; Suin et al. 2022).

In addition to preventing the gas from collapsing, the turbulence and the magnetic field can also bolster the star formation processes. For example, the coupling between the magnetic field and the neutral gas can allow parts of low-density clouds to fragment and initiate star formation (Fiedler & Mouschovias 1993). In this scenario, while neutral gas is collapsing, a part of the magnetic flux is removed, which further contributes to the gravitational instability (Mouschovias & Ciolek 1999; Lazarian et al. 2012; Priestley & Whitworth 2022; Tritsis et al. 2023).

However, it has been difficult to test these star formation models observationally. Despite all of the effort made over the past 40 yr to characterize interstellar turbulence, we still struggle to understand the different energy sources that contribute to the observed line-of-sight gas velocity dispersion. The only available way to measure turbulent motions is by observing the gas velocity along the line of sight (Larson 1981; Hennebelle & Falgarone 2012), which does not provide a complete picture. As noted by Ballesteros-Paredes et al. (2011), a collapsing cloud has energetic properties (at least in terms of the total kinetic energy) similar to those of an identical cloud supported against gravity by turbulence. In addition, measurements of magnetic fields are observationally extremely challenging, as they typically involve relatively weak signals, as is the case for Zeeman gas line splitting (Crutcher 2012) or (sub)millimeter/far-infrared (FIR) dust polarization measurements (see, Pattle et al. 2019, for a recent review).

Interstellar dust thermal emission is polarized (Hall & Mikesell 1949; Hiltner 1949) due to the presence of B fields and non-spherical dust grains in the interstellar medium (ISM). The explanation for this phenomenon was firstly proposed by Davis & Greenstein (1951) as a paramagnetic alignment with the magnetic field. Later on, this process was described by the radiative alignment torques (RAT) theory (Hoang & Lazarian 2014; Andersson et al. 2015, and references therein), in which the minor axis of dust grains is aligned parallel to the direction of the B field. Consequently, the (sub)millimeter dust continuum emission is polarized perpendicular to the direction of the component of the B field in the plane of the sky (POS).

Regardless of the observational challenges, technological advances in the past decade are now allowing for polarized dust emission to be measured over increasingly extended regions. The Planck satellite provided an all-sky map of the magnetic field geometry in the diffuse ISM, albeit at a low angular resolution of 5–15′. These observations revealed that the galactic B field is intertwined with the filamentary structure of the ISM. In particular, the POS orientation of the B field and the filamentary structures is correlated with the observed column density, with the field largely parallel to diffuse structures (striations) while perpendicular to dense filaments (Planck Collaboration Int. XXXII 2016; Planck Collaboration Int. XXXV 2016). While supported by independent studies (McClure-Griffiths et al. 2006; Goldsmith et al. 2008; Peretto et al. 2012; Clark & Hensley 2019), this picture is complicated by potential projection effects (Panopoulou et al. 2016) or the presence of stellar feedback (H II bubble shells, Chapman et al. 2011).

The Balloon-borne Large Aperture Submillimeter Telescope for Polarimetry (BLASTPol) has also observed a similar B-field morphology at a higher angular resolution (2.5′, Fissel et al. 2019), which has been interpreted as evidence of the magnetic field dominating the energy balance in the diffuse gas. In contrast, gravity is dominant in the dense, star-forming regions. This interpretation has recently been reinforced by high angular resolution (14″) ground-based observations using the POL-2 instrument on the James Clerk Maxwell Telescope (JCMT; Pattle et al. 2017; Liu et al. 2019; Wang et al. 2019), which have shown that field lines bend along dense star-forming filaments, possibly as a result of gravitationally driven flows (Goldsmith et al. 2008; Chapman et al. 2011; Hwang et al. 2023).

A transition from magnetically dominated to gravitation-ally dominated gas requires a redistribution of magnetic flux (Tritsis et al. 2022). In principle, such a transition occurs before the cores form (Ching et al. 2022), which implies that cores should be generally super-critical. Nevertheless, a change in the gravitational stability of the gas should be accompanied by a corresponding change in its kinematic properties, arguing for the necessity of combining high angular resolution dust polarization and velocity-resolved molecular emission imaging. Such studies have only recently started and cover relatively small isolated fields. A recent study (Tang et al. 2019) has indeed investigated the relation between dense gas velocity gradients (as imaged in N₂H⁺) and the magnetic field direction in massive star-forming infrared dark clouds and concluded that the two are strongly correlated as a result of gravity dragging matter toward the center of the massive ridge. Extending such studies to larger fields is of crucial importance and the High-resolution Airborne Wideband Camera-Plus (HAWC+, Harper et al. 2018) instrument on the Stratospheric Observatory for Infrared Astronomy (SOFIA) offered exceptional capabilities for mapping the magnetic field geometry, as demonstrated, for example, by the observations of the OMC-1 region at 53, 89, 154, and 214 μm (Chuss et al. 2019). The capabilities of HAWC+ for large-area magnetic field mapping have further improved with the commissioning of the on-the-fly map (OTFMAP) polarimetric mode for fast wide-field polarimetric imaging.

When estimating the strength of the POS B field, it is critical to know the level of turbulence in the gas because the small-scale turbulence causes local deviations from the mean direction of the magnetic field. The Davis-Chandrasekhar-Fermi (DCF) method (Davis & Greenstein 1951; Chandrasekhar & Fermi 1953) provides a way of calculating the POS B-field strength using information about gas density, turbulence, and changes in the direction of the magnetic field. This method assumes linear geometry and sub-Alfvénic (magnetically dominated) turbulence. Therefore, it is necessary to obtain information about the turbulence of the gas and dust polarization to use this method. Several variations of the DCF have been developed over the last few decades. The main modification lies in taking into account the number of turbulent cells present along the line of sight and those captured within the beam by including a correction factor Q, which takes a value between zero and one (see Eq. (16) in Ostriker et al. 2001). In another variation, the deviation in the mean direction of the B field is substituted by the ratio between the ordered and turbulent magnetic field component (Houde et al. 2009, and references therein). The DCF method assumes that imcompressible motions cause the dispersions of the observed polarization angles, which is not always applicable within the ISM. Skalidis & Tassis (2021) provided an alternative method for deriving the magnetic field strength from dust polarization measurements – the Skalidis-Tassis (ST) method. In addition to Alfvénic modes, this study accounts for the presence of compressible motions in the gas without discarding the anisotropic nature of the turbulence, providing a physically motivated approach to estimating the magnetic field strength (Skalidis et al. 2021). For this work, we used the ST method to investigate highly structured regions in the Flame nebula in the Orion B complex, while also presenting the results from the modified DCF method (for example, Ostriker et al. 2001; Crutcher 2004; Lyo et al. 2021) for comparison.

Although OMC-1 is the closest and best-studied high-mass star-forming region, it is not a typical cloud to study the role of the magnetic field in star formation. The reason is that OMC-1 is particularly dense (10⁵ cm⁻³) and very active in terms of star formation, affected by strong shocks (including an explosive outflow resulting from the merger of three massive protostars, see, Bally 2008). These shocks exhibit exceptionally intense UV illumination, with enhancement factors of G₀ ~10⁴, up to 10⁵ relative to the standard interstellar radiation field (ISRF) (in the vicinity of the Trapezium cluster, Goicoechea et al. 2015, 2019). The presence of a strong radiation field implies that the thermal pressure is high, and the magnetic field plays a limited role in the photon-dominated region (PDR) gas dynamics.

Fortunately, the vast quantity of molecular data already obtained in Orion B makes this region a particularly well-suited location for studying the correlation between changes in the magnetic field geometry and associated gas velocity gradients. The Orion B giant molecular cloud (GMC) hosts NGC 2024, one of the closest high-mass star-forming regions (at a distance of d = 410 pc; Cao et al. 2023), and it is overall more representative of a standard GMC in our galaxy and normal galaxies, with the far-ultraviolet (FUV) radiation field G₀ in the range of ~10³−10⁴ (Santa-Maria et al. 2023).

In a recent study, Orkisz et al. (2017) used ¹²CO(1−0) and ¹³CO(1−0) lines to characterize observationally the ratio of compressive versus solenoidal motions in the turbulent flow and to relate this to the star formation efficiency in various regions of Orion B. Orkisz et al. (2019) accurately analyzed the dynamics of the filamentary network using C¹⁸O(1−0). Most identified filaments in Orion B are low-density, thermally subcritical structures, not collapsing to form stars. Only about 1% of the Orion B cloud mass can be found in super-critical, star-forming filaments, consistent with the low overall star formation efficiency of the region (Orkisz et al. 2019).

NGC 2024 (Fig. 1) is located east of Alnitak (ζ Ori) in the Orion B complex (Meyer et al. 2008). This region contains a massive and young H II region (age 2 × 10⁵ yr, Tremblin et al. 2014, and priv. comm.), deeply embedded in dust, located in the foreground, and extending to the south (Barnes et al. 1989, see gray contrours showing the network of filaments from Orkisz et al. 2019 in the left panel and labels in the right panel in Fig. 1). The north-south filament in NGC 2024 is super-critical (Orkisz et al. 2019) and the site of ongoing star formation as recently confirmed in the southern part of this filament observed by NOEMA (Shimajiri et al. 2023).

The dust bar observed along the line of sight is visible in the ESO-Vista image (right panel in Fig. 1) as the dark extinction pattern across the image, leading to an apparent cooler dust temperature derived from Herschel observations (Lombardi et al. 2014). The central part of NGC 2024 contains warm dust and gas, heated by the H II region, as well as embedded protostellar objects (Mezger et al. 1988, 1992; Lis et al. 1991) located in the background (light blue points in the right panel in Fig. 1). The young H II region is expanding, strongly impacting its parental cloud, and creating sharp ionization fronts toward the south, which makes it a good example demonstrating how such systems can efficiently exert stellar feedback. The edges of the bubble are seen toward the west and east of the center of NGC 2024 (labeled as Eastern and Western Loop in Fig. 1).

Fig. 1 Multiwavelength image of NGC 2024. Left panel: Color-composite image of the Flame nebula (NGC 2024) showing peak intensities of 12CO (1−0) (blue) emission and isotopologues, 13CO (1−0) (green) and C18O (1−0) (red) obtained by the IRAM 30-meter telescope (image credits: Pety et al. 2017). We overlaid the SOFIA HAWC+ field of view as a white dashed rectangle. Gray contours show a network of filaments presented in Orkisz et al. (2019). We label the regions we investigate in this work: Bubble East, Bubble West, Filament, and Bubble+Filament. Right panel: ESO-VISTA image (ESO/J. Emerson/VISTA, Cambridge Astronomical Survey Unit). We sketched the environments seen across NGC 2024: the H II region (white dashed circle, see Table 2 in Gaudel et al. 2023) and the filament (orange line, Orkisz et al. 2019). Light blue points are the positions of far-infrared sources in the background (Mezger et al. 1988, 1992). In addition, we labeled the edges of the expanding H II region and the direction of the stellar feedback driven by the radiation produced by recent star formation.

2 Observations

2.1 Dust polarization measurements using SOFIA HAWC+

Our work employs the dust polarization measurements acquired using HAWC+ on SOFIA in September 2021 for program 09_0015 (PI: D. Lis). Specifically, we observed NGC 2024 at 154 μm (Band D) and 214 μm (Band E) with HAWC+ in polarization mode on flights F779 (8 September 2021), F780 (9 September 2021), F782 (11 September 2021), F783 (14 September 2021), and F784 (15 September 2021). The map obtained at each wavelength comprises ten 20′ × 7′ strips in a weave pattern: five vertical and five horizontal. These strips were designed using the On-The-Fly scan mode of the HAWC+ camera, which covers the requested area using Lissajous patterns on the sky (Harper et al. 2018). Each strip was repeated at least twice.

The level 0 data for each flight were downloaded from the NASA/IPAC Infrared Science Archive (IRSA) and reduced in July 2023 using the SOFIA Data Reduction software (SOFIA Redux version 1.3.0; Clarke & Vander Vliet 2023). The HAWC+ scan mode reduction package from SOFIA Redux was initially built in Python from the Java-based CRUSH data reduction software (Kovács 2008). All available files at a given wavelength (104 each for Band D and Band E) were loaded in SOFIA Redux to be reduced using the default parameters of the software except for two options at the Compute Scan Map step of the pipeline. Specifically, the fixjumps option was set to True to filter out flux jumps in individual detectors during observations, and the rounds option was set to 40 to improve the recovery of astronomical large scale flux from the background subtraction. Appendix A.1 shows the comparison of the Stokes I values from the resulting data to archival Herschel fluxes and the improvement relative to the Level 4 data available on IRSA at the time of writing.

The pixel scale of the HAWC+ polarization maps and the effective beam size is 3.4″ and 14.0″, respectively, in Band D, and 4.6″ and 18.7″ in Band E. The final data products for both Band D and E contain the I, Q and U Stokes parameters, the polarization fraction Ρ and angle θ, and their uncertainties (Gordon et al. 2018). Clarke & Vander Vliet (2023) gives the full calculation for each quantity, which we summarize here, assuming that the cross-terms of the error covariance matrix are negligible.

Stokes I describes the total dust thermal emission, and its polarized component I′_p is calculated using the following equation: $${I^{\prime }}_p=\sqrt{Q^2+U^2}.$$(1)

The uncertainty δI′_p on the polarized intensity I′_p is then $$\delta {I^{\prime }}_p=\frac{\sqrt{{(Q\delta Q)}^2+{(U\delta U)}^2}}{{I^{\prime }}_p},$$(2)

where δQ and δU are the uncertainties on Stokes Q and U, respectively.

The polarized intensity I′_p has to be corrected for the bias created by the quadratic addition of the noise in the Stokes Q and U maps (Wardle & Kronberg 1974; Naghizadeh-Khouei & Clarke 1993). This de-biased polarized intensity I_p (although this is not the only way to do so, Montier et al. 2015), used in this work, is calculated using $$I_p=\sqrt{{I^{\prime }}^2{}_p-\delta {I^{\prime }}^2{}_p},$$(3)

with uncertainty δΙ_p = δΙ′_p.

The polarization fraction P is then obtained from $$P=100\frac{I_p}{I},$$(4)

with the uncertainty δΡ given by $$\delta P=P\sqrt{{\left(\frac{\delta {I^{\prime }}_p}{{I^{\prime }}_p}\right)}^2+{\left(\frac{\delta I}{I}\right)}^2}.$$(5)

The polarization angle θ_p is defined using the Stokes parameters Q and U as $${\theta }_p=\frac{1}{2}\arctan \frac{U}{Q},$$(6)

with the uncertainty δθ_p given by $$\delta {\theta }_p=\frac{1}{2}\frac{\sqrt{{(Q\delta U)}^2+{(U\delta Q)}^2}}{{I^{\prime }}_p^2}.$$(7)

Since the thermal emission from interstellar dust grains is preferentially polarized perpendicular to the plane of the sky magnetic field lines (Hoang & Lazarian 2014; Andersson et al. 2015, and references therein), the direction of the magnetic field in the plane of the sky can be obtained by adding π/2 to Eq. (6).

Prior to the data analysis, we filter out low signal-to-noise ratio (S/N) data points. After masking the data, we keep pixels that show an S/N ≥ 50 in total intensity, ≤30° in polarization angle uncertainty, and ≤30% in polarization fraction. Our HAWC+ observations filter out a non-negligible fraction of low-level extended emission. We estimate the amount of missing flux in the SOFIA observations by comparing our data with PACS measurements in Appendix A.1. The data set used in this study was reduced using a larger number (40) of iterations, or rounds, than the default (15) used for the original Level 4 data available on IRSA. A larger number of iterations during data reduction typically improves the recovery of diffuse large-scale emission for scan mode maps. We compare our reduced maps of each Stokes parameter and the polarization angle with those derived from the standard pipeline setup and show them in Figs. A.2 and A.3. We find an overall good agreement between these two data sets.

2.2 IRAM 30-m observations of Orion B

In our work, we make use of information of the ¹²CN(1−0) and HCO⁺ (1−0) emission from the ongoing IRAM-30m ORION-B Large Program (PIs: M. Gerin & J. Pety, see the left panel in Fig. 1 taken from Pety et al. 2017). ORION B images a 5 square degree field (~20 pc across) in the Orion B molecular cloud, at an angular resolution of 26″ (10⁴ AU, 0.05 pc) in at least 30 molecular lines in the (71–79) and (84–116) GHz range with a spectral resolution ~0.6 km s⁻¹. These include common molecular tracers such as ¹²CO(1−0), HCO+(1−0), HCN(1−0), and CS(2−1), as well as their optically thin isotopologues, which have narrow line widths and are most sensitive to kinematic variations. The resulting wide-field hyper-spectral data cube is genuinely unique in terms of its massive information content (~820 000 pixels, ~240 000 spectral channels per pixel), enabling an unprecedented characterization of the physical structure, chemistry, and dynamics of a GMC, and their connection to its star formation activity.

The J = 1−0 lines were observed in 2013–2020 in the context of the ORION-B Large Program, using a combination of the EMIR receivers and FTS spectrometers. The data reduction is described in Pety et al. (2017). It uses the standard methods provided by the GILDAS¹/CLASS software. The data cubes were reprojected on the same astrometric grid as the SOFIA data. We use the data cubes at their original spatial resolution of 30″ and 23″ for the HCO⁺ and ¹²CN lines, respectively. The velocity spacing is 0.5 km s⁻¹. The data are calibrated on the main beam temperature scale. The achieved noise levels are 0.26 and 0.58 K (see Table 2 in Gratier et al. 2017), respectively, over the studied field of view.

We use spectroscopic data of cyanide (¹²CN) and formyl cation (HCO⁺) (J = 1−0 transition line) to trace the UV-irradiated gas (Bron et al. 2018). CN is present in UV-illuminated edges as a photodissociation product of HCN and HNC, resulting from UV-dominated chemistry, and can be collisionally excited by electrons and neutrals (Santa-Maria et al. 2023). Therefore, the CN rotational lines remain bright in dense UV-illuminated edges, which makes them excellent tracers of the UV-dominated regions, like the one seen in NGC 2024. In addition, CN shows the multiple hyperfine components (Penzias et al. 1974), spectrally resolved by our observations. We report the properties of each hyperfine component used in this work in Table 1, which allow the opacity (Milam et al. 2009) and excitation temperature of the (1−0) transition to be derived, enabling its excitation to be constrained accurately and for limits on the gas density to be obtained.

HCO⁺, similarly to CN, traces photon-dominated regions (Young Owl et al. 2000; Lis & Schilke 2003). HCO⁺ has different production pathways dominant at different physical conditions of the gas. At the edges of the dissociation region, HCO⁺ is a precursor of CO, and can be produced through the CO⁺ (Hogerheijde et al. 1995; Young Owl et al. 2000; van der Werf et al. 1996) or CH⁺ molecule (Goicoechea et al. 2016). Recent studies detected HCO⁺ emission the edge of the PDR in the Orion bar (Goicoechea et al. 2016), and the Horse-head nebula (Hernández-Vera et al. 2023). Pety et al. (2017); Gratier et al. (2017) found that CN and HCO⁺ are sensitive to the UV radiation. Moreover, by applying a clustering algorithm, Bron et al. (2018) found that CN and HCO⁺ trace the UV-radiated gas, contrary to the C¹⁸O molecule, which gets photodissociated. Therefore, using CN and HCO⁺ emission in our work, together with the FUV heated dust traced by SOFIA HAWC+ observations, we can directly assess the gas impacted by the radiative stellar feedback in NGC 2024.

3 Dust polarization

In this section, we discuss our SOFIA HAWC+ dust polarization measurements introduced in Sect. 2.1, and present the magnetic field morphology across NGC 2024.

The SOFIA HAWC+ Band D measurement is shown in Fig. 2. The map is masked, as explained above, and is presented at its native angular resolution of 13.6″, which corresponds to linear scales of ~0.027 pc. The background of this figure (all panels) shows total Stokes intensity, and the black lines indicate the orientation of the magnetic field. In our work, we mainly focus on four specific regions across NGC 2024: the edges of the bubble (middle panels in Fig. 2), the filament (bottom left panel in Fig. 2), and the intersection of filament and bubble (bottom right panel in Fig. 2). Each of these panels shows a circular-shaped, beam-sized region within which we computed the magnetic field strength. We chose a beam-sized region for data analysis rather than taking the information from a single pixel to get suitable properties of each region.

We show SOFIA HAWC+ Band E map in Fig. A.4. The overall agreement in magnetic field direction traced by dust polarization at 154 and 214 μm is observed within NGC 2024. We show the difference in measured polarization angles and their rms in Fig. A.5. While the rms is only a few degrees (~2 deg) across a large part of NGC 2024, we observe higher values (> 10 deg) toward the center of NGC 2024. Nevertheless, we consider only Band D data in this work because of an overall good correspondence (measured as a low rms) between dust polarization angles measured from Band D and E across the regions we further analyze in this work.

The dust continuum emission is strongest toward the center of NGC 2024. A significant emission is observed along the filament that is in the front of the ionization region caused by the young massive star IRS 2b (the central part of our map, Bik et al. 2003), and along the ionization fronts on the southeastern and southwestern part of NGC 2024 (labeled as Bubble West and Bubble East, and shown in the middle row in Fig. 2).

In addition, we compare our magnetic field directions to the direction inferred from the near-infrared (NIR) polarization of young stars (Kandori et al. 2007). The NIR polarization directly traces the direction of the magnetic field (that is, there is no π/2 difference as for the FIR dust polarization – Sect. 2.1). We show this comparison in the top panels of Fig. A.6. As seen from the figure, we find that magnetic field vectors derived from NIR polarization and our work are parallel to each other, indicating an agreement between these two datasets in NGC 2024 (see Fig. A.6). In addition, we compare the magnetic field direction with the 100 μm (Dotson et al. 2000) and submillimeter (850 μm) polarization measurements (Matthews et al. 2002) in the central part of NGC 2024 (labeled with the dashed black rectangle in top panels of Fig. A.6) and show zoomed-in panels in the bottom row of Fig. A.6.

The area labeled with dark blue contours in the top panels in Fig. A.6 is particularly interesting because we observe an outflowing feature in the HCO⁺ emission at velocities of ~14 km s⁻¹. The outflow originates from the FIR 5 source located in the dense molecular cloud behind the expanding H II region (Richer et al. 1992; Greaves et al. 2001; Choi et al. 2015). Dust continuum emission at 154 μm possibly traces the outflow in this region (see also the bottom middle panel in Fig. 3), as we note that magnetic field lines follow the direction of the outflow. We have not observed this behavior at 214 μm.

The magnetic flux is frozen within the molecular gas, and as a consequence, the B field will trace its morphology. Therefore, the local environmental conditions that shape the distribution of molecular gas will also impact the magnetic field. On the one hand, the magnetic field is highly ordered in some regions of NGC 2024; for instance, B field follows the dusty filament to the south of NGC 2024 (see, for instance, the bottom left panel of Fig. 2), except for the very southern part of the filament (close to the bottom left corner of the top panel in Fig. 2), where the magnetic field is perpendicular to the filament. At the northern part of NGC 2024, the magnetic field direction varies from following the filament to being perpendicular.

We observe the nearly horizontal direction of the B field at the edges of the expanding H II region (middle panels in Fig. 2). Gas affected by the stellar feedback (for example, UV radiation, stellar winds) is pushed outward from the H II region.

Consequently, magnetic field lines become parallel to the edges of such an expanding shell (for example, Tahani et al. 2023). On the other hand, the magnetic field appears chaotic in the central and densest parts of NGC 2024. This area shows great complexity, as it results from the mixture of the filament located in the front, the H II region and the dense molecular gas in the background (see Fig. 8 in Matthews et al. 2002, or Fig. 5 in Roshi et al. 2014). One possible explanation for this observed morphology of the magnetic field is that the magnetic field changes its orientation relative to the filament depending on the column density contrast between the filament and the background emission, from being parallel to perpendicular to the molecular gas distribution($\text{Δ}{N}_{{\text{H}}_2}>{10}^{20} {\text{cm}}^{-2}$, for example, Planck Collaboration Int. XXXIV 2016; Alina et al. 2019). Additionally, the magnetic field could be “pinched” due to the gravitational collapse of the gas (Basu 2000; Lai et al. 2002; Doi et al. 2020, 2021), which can also explain the magnetic field direction at the northern part of NGC 2024 and the very southern part of the filament.

Fig. 2 SOFIA HAWC+ 154 μm (Band D, all panels) dust continuum maps at 13.6″ angular resolution. The maps are masked based on the measured S/N in measured Stokes intensity and polarization angle. Black lines in all panels show the orientation of the magnetic field for every fifth pixel. We show Band D dust polarization map across NGC 2024 (top panel). Bottom rows show the zoomed-in regions we analyse in this work: egde of the bubble on the western (middle left panel), eastern (middle right panel), in the filament (bottom left panel), and in the overlap region of the filament and bubble (bottom right panel). A circle in each of these zoomed-in panels shows the region we use to compute the dispersion in the mean angle of magnetic field, as well as number density and turbulent velocity dispersion.

Fig. 3 Moment maps of 12CN(1−0) and HCO+ (1−0) emission across NGC 2024. Top: maps of the CN emission: integrated intensity (the zeroth moment, left panel), the centroid velocity (the first moment, middle panel), and the FWHM (the second moment, right). Bottom: the same as in the top row, but for HCO+. We show the beam size for both molecular lines at the bottom left corner of each panel and the 0.5 pc scalebar at the bottom right corner. All pixels shown in these maps result from the masking technique presented in Einig et al. (2023). We refer the reader to the description of the production of these moment maps in Sect. 4.1.

4 Characterizing the turbulence and density structures in NGC 2024

The DCF method (Sect. 1) requires the characterization of the turbulent velocity dispersion and the gas density in the same material traced by the dust emission. Therefore, in this section, we briefly describe our analysis of molecular data from the ORION B large program and present moment maps of ¹²CN(1−0) and HCO⁺(1−0) and then provide results on measured turbulent line widths and gas volume densities of each region shown on the left panel of Fig. 1. We derive these quantities using non-LTE radiative transfer models (RADEX, van der Tak et al. 2007).

4.1 Moment maps

Figure 3 shows moment maps for the CN(1−0) (top row) and HCO⁺(1−0) (bottom row) emission. We create these moment maps using CUBE in GILDAS. Before creating moment maps, we mask CN and HCO⁺ cubes. The mask is created using the segmentation technique in CUBE, where we identify neighboring voxels with continuous S/N; we have selected all voxels with S/N above 2 (see Einig et al. 2023, for a complete description). In this case, we integrate all molecular emission along each unmasked line of sight, without considering the possible presence of multiple velocity components. Finally, we apply to the moment maps the spatial mask where the dust polarization measurements are reliable as described in Sect. 3. Similarly, as for the intensity map, the line width map shown in this figure corresponds to the measured second moment, which does not consider possible spectral complexity (multiple velocity components or the presence of an outflow) or a correction for the opacity broadening. These corrections could bias our measurement, particularly in the central part of NGC 2024, as shown by the spectral decomposition results of CN and HCO⁺ (in Figs. B.1 and B.2 respectively).

The integrated line intensity (moment 0) is shown in the left panels of Fig. 3. Both CN and HCO⁺ show the brightest emission toward the central part of NGC 2024, at the heavily dust-obscured region, as seen in Fig. 1.

The first moment, also known as the centroid velocity map, is shown in the central panels of Fig. 3. The central velocity is at 10 km s⁻¹. However, we note the presence of multiple velocity components along the line of sight in both CN and HCO⁺, as previously identified in the ¹³CO(1−0) and C¹⁸O(1−0) emission (Gaudel et al. 2023). We discuss velocity components in Appendices B.1 and B.2.

The second moment map (the observed line width (σ_obs) map) is shown in the right panels of Fig. 3. We note that to first order the HCO⁺ line is broader than CN, which has several possible reasons. Firstly, HCO⁺ is brighter and more spatially distributed than CN and, therefore, could result in having a broader line. Second, HCO⁺ emission can be optically thick (see, Barnes & Crutcher 1990), which additionally broadens the line. Third, it is possible that HCO⁺, similarly to ¹³CO and C¹⁸O has multiple velocity layers (Gaudel et al. 2023). Nevertheless, we attempt to correct these effects in our analysis described in the following sections. The regions we analyze in this work do not contain multiple velocity features in the CN emission (see Fig. B.1). However, in the case of HCO⁺ emission, we do not spectrally resolve the multiple-component HCO⁺ emission. However, we note the possibility of HCO⁺ showing multiple velocity components (Fig. B.2) in NGC 2024. The more thorough analysis of the HCO⁺ velocity field is, therefore, beyond the scope of this paper, and will be an aim of the upcoming studies.

4.2 Measuring turbulent velocity dispersion

Here, we briefly describe steps in order to constrain the turbulent velocity dispersion. Prior to this, we note the different natures of the ¹²CN and HCO⁺ line emission profiles. For instance, the ¹²CN(1−0) has a hyperfine structure (see Fig. C.9). We do not observe any anomalous hyperfine structure emission in the CN emission, as reported, for example, for the HCN emission in Orion B (Santa-Maria et al. 2023). Therefore, all components of the multiplet have the same line width. We provide more information on the properties of the hyperfine structure of the ¹²CN(1−0) emission in Appendix B.1. The profile of the HCO⁺ can be described by a Gaussian function (see Fig. C.9) assuming its optically thin emission. In the case of the optically thick emission, the line profile will be changed and depend on the observed optical depth.

To derive the turbulent line width for both CN and HCO⁺, we correct their measured line widths (that we will infer from the RADEX modeling) as follows. First, we correct our measurements for the contribution of thermal broadening: $${\sigma }_{\text{NT},\text{mol}}=\sqrt{{\sigma }_{\text{obs},\text{mol}}^2-{\sigma }_{\text{TH},\text{mol}}^2},$$(8)

where σ_obs,mol is the measured FWHM of a molecular line, and σ_TH,mol is the line width of the thermal component. We calculate the thermal broadening using the following equation: $${\sigma }_{\text{TH},\text{mol}}=\sqrt{\frac{k{T}_{\text{k}}}{m_{\text{mol}}}},$$(9)

where k is Boltzman’s constant, T_k is the kinetic temperature obtained using the CO(1−0) measurements (Orkisz et al. 2017), and m_mol is the mass of a molecule. A final correction that we apply in our analysis is for opacity broadening of molecular lines, using the following equation: $$\text{Δ}\upsilon =\frac{{\sigma }_{\text{NT},\text{mol}}}{\beta },$$(10)

where the factor β is a function of the optical depth (τ₀) at the line center, defined as (Phillips et al. 1979; Hacar et al. 2016; Orkisz et al. 2017): $$\beta =\frac{1}{\sqrt{\ln 2}}\cdot {\left[\ln \left\{\frac{{\tau }_0}{\ln \left(\frac{2}{\exp \left(-{\tau }_0\right)+1}\right)}\right\}\right]}^{\frac{1}{2}}.$$(11)

The equations presented above are crucial at subparsec scales since these effects significantly contribute to the measured line width at these scales. For example, the optical depth effect could broaden the line by a few tens of percent, which is the case of CN and HCO⁺ . This is discussed in the following section.

The FWHM of a molecular line, corrected for the contribution mentioned above of the thermal and opacity broadening, can be derived in multiple ways. For example, the FWHM can be derived directly from measuring the second order moment map (right panel in Fig. 3, Sect. 4.1), or from the line fitting by using, for example, the spectral decomposition (see Appendices B.1 and B.2), or directly from the non-LTE modeling of the emission spectrum for a set of input parameters, that describe the physical conditions of the gas within a selected region. In this work, we use the latter method, as we aim to find a set of physical parameters that describe the edges of the bubble exposed to the FUV emission and the gas coming from the filamentary structure. Therefore, we derive the FWHM from the non-LTE modeling and RADEX analysis, as well as the gas number density described in the following section.

4.3 Radiative transfer results

To derive gas volume densities and measure the turbulence in NGC 2024, we have employed the non-LTE radiative transfer code RADEX using Python wrapper (SpectralRadex, Holdship, in prep.) because it allows the user to compute the spectrum of a line and directly compare models to observations. Using the excitation of CN and HCO⁺ and some a priori information and assumptions, we derive gas number densities, $n_{{\text{H}}_2}$ and line widths across four regions in NGC 2024: two edges of the expanding H II bubble (located to the west and the east), the filament, and the region consisting of the edge of the bubble and the filament (located in the north of NGC 2024).

In the following, we provide information about our input parameters and assumptions. We take the value of 2.73 K as the background temperature, and use the kinetic temperature (T_kin) derived from the ¹²CO(1−0) peak intensity (Orkisz et al. 2017), and shown in Table 2. We create a grid of column densities in the range of (1−4) × 10¹⁴ cm⁻² for CN and (1–4) × 10¹³ cm⁻² for HCO⁺ , following results from Bron et al. (2018). The selected grid of line widths covers the range from or 1 to 2.5 km s⁻¹. Next, we create a grid in molecular hydrogen densities from 10² to 10⁵ cm⁻³. The density of H₂ is comprised of para- and ortho-H₂, assuming the temperature dependence of the ortho-to-para ratio: $9\cdot {\text{e}}^{-170/T_{\text{kin }}}$ (Mandy & Martin 1993). Here, we assume a fixed electron fraction (electron-to-H₂ density ratios), f_e = 1.4 × 10⁻⁴, which is a maximum value found in PDRs for the case when number densities of ionized and neutral carbon (including CO) are equal (Sofia et al. 2004; Graf et al. 2012). Our choice of fe corresponds to values of the ionization fraction for translucent medium within Orion B found in Bron et al. (2021). With fixed ortho-to-para H2 ratio and electron fraction, we run a grid in three independent variables: line width, column density, and H₂ volume density.

We consider all relevant collision partners for CN and HCO⁺ in our input files. For the CN emission, we take into account its hyperfine structure (Müller et al. 2005), and include collisions with ortho- and para-H₂ (Kalugina et al. 2012), as well as electrons (Harrison et al. 2013; Santa-Maria et al. 2023). The input file with collision partners comes from the combination of two data files from the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005) and Excitation of Molecules and Atoms for Astrophysics (EMAA²; EMAA 2021) to include all relevant collisional partners and take into account the hyperfine structure of CN. In the case of HCO⁺ , we consider collisions with ortho-, para-H₂ (Denis-Alpizar et al. 2020) and electrons (Fuente et al. 2008) as well.

The output parameters from the RADEX modeling are the excitation temperature, T_ex and opacity, τ. We additionally compute the peak temperature of the spectrum generated from each model using SpectralRadex. In the case of CN, we calculate excitation temperature and the opacity of each component of the multiplet, which is scaled using information about relative intensities (Table 1).

We compute model CN(1−0) and HCO⁺(1−0) spectra assuming the optical depth to be a Gaussian function of frequency, τ(ν): $$T_{mb}(v)=\left(T_{\text{ex}}-T_{\text{bg}}\right)\cdot \left(1{\text{e}}^{-\tau (v)}\right).$$(12)

The best-fit model is found based on the modeled spectrum for each combination of ($n_{{\text{H}}_2}$, N_mol, Δυ). We perform a χ² minimization for the modeled and observed peak temperature, T_peak. Then, we calculate the turbulent line width following the prescription in Sect. 4.2. We report our turbulent line width and gas number density results in Table 2. Model spectra inferred from the RADEX analysis of CN and HCO⁺ for a selection of physical parameters that best describe the observed spectra in each region we investigate in this work are shown in Fig. C.9.

5 Magnetic field strength

5.1 Davis–Chandrasekhar–Fermi method

To derive the strength of the magnetic field, we use the Davis–Chandrasekhar–Fermi method (Davis & Greenstein 1951; Chandrasekhar & Fermi 1953), which provides a recipe for calculating the POS B-field strength as following: $$B_{\text{pos}}=\sqrt{4\pi \rho }\frac{{\sigma }_{\text{NT}}}{{\widehat{\sigma }}_{\text{c}}(\phi )},$$(13)

where ρ is the gas volume density, σ_ΝΤ is the nonthermal velocity dispersion (${\sigma }_{\text{NT}}=\text{Δ}\upsilon /\sqrt{8\ln 2}$), where Δυ is the measured FWHM of the line, corrected for the line broadening and opacity effects – Eqs. (8) and (10) – Sect. 4.2, and ${\widehat{\sigma }}_{\text{c}}(\phi )$ is the spatial dispersion of magnetic field angle.

This equation is in CGS units, and the DCF method also assumes that any perturbation in the magnetic field originates from local, small-scale turbulence. The stronger the magnetic field, the smaller will be the perturbation caused by turbulence. By calculating all constants and keeping the units of number density, line width and the angle dispersion in cm⁻³, km s⁻¹, and deg, respectively (Lyo et al. 2021), Eq. (13) can be expressed as: $$B_{\text{pos}}\approx 9.3\cdot \sqrt{n_{{\text{H}}_2}}\cdot \frac{\text{Δ}v}{{\widehat{\sigma }}_{\text{c}}(\phi )}[\text{μG}].$$(14)

We have included in the above a factor of 0.5 for overestimation of the magnetic field strength due to line of sight integration effects (see, for instance, Ostriker et al. 2001).

Fig. 4 Direction of the magnetic field (left panel), the mean angle of the magnetic field (middle panel), and the rms of the magnetic field direction (right panel). Maps in the middle and right panels have been made using a 12 × 12 pixels sliding window with weight inversely proportional to the measured uncertainty variance. All maps have been regridded to match the grid size of the CN and HCO+ data. The gray contour indicates the mask we defined for dust polarization measurements, described in Sect. 2.1.

5.2 Skalidis–Tassis method

The DCF method is widely used in the literature to derive the strength of the magnetic field (see, Pattle et al. 2019). However, it assumes that isotropic turbulence and Alfvénic waves in an incompressible medium cause the observed dispersion in polarization angle. Other mechanisms such as magneto-hydrodynamic waves (Heyvaerts & Priest 1983) and entropy modes (Lithwick & Goldreich 2001) also cause fluctuations in polarization angle. Therefore, to derive the magnetic field strength in NGC 2024, it is important to acknowledge the contribution of non-Alfvénic motions and the compressible nature of the ISM. In this work, we derive the B_POS using the prescription presented in Skalidis & Tassis (2021); Skalidis et al. (2021): $$B_{\text{pos}}=\sqrt{2\pi \rho }\frac{{\sigma }_{\text{NT}}}{\sqrt{{\widehat{\sigma }}_{\text{c}}(\phi )}},$$(15)

where ρ, σ_NT, and ${\widehat{\sigma }}_{\text{c}}(\phi )$ are the same as in Eq. (13).

Similarly as in Eq. (14), after substituting all constants, the above equation becomes: $$B_{\text{pos}}\approx 1.8\cdot \sqrt{n_{{\text{H}}_2}}\cdot \frac{\text{Δ}\upsilon }{\sqrt{{\widehat{\sigma }}_{\text{c}}(\phi )}}(\text{μG}).$$(16)

We favor using this approach rather than the classical DCF because it is physically motivated, considering the nature of the ISM in NGC 2024. The ST method takes into account the compressible nature of the gas, and Orkisz et al. (2017) found compressible, non-Alfvénic motions to dominate over solenoidal modes in NGC 2024.

5.3 Sliding window

We present the dispersion of the mean direction of the magnetic field in Fig. 4. We produced this map by using a “sliding window” to remove the impact of gradients of the large-scale magnetic field (“unsharp-masking”, Dharmawansa et al. 2009; Pattle et al. 2017). The “sliding window” technique is based on computing the mean and standard deviation of the magnetic field orientation within the window that is three times bigger than the beam and has 12 × 12 pixels. The size of the sliding window is made to ensure that we remove a large-scale magnetic field contribution without losing information on the small-scale perturbation of the magnetic field (Pattle et al. 2017). As shown in Mardia & Jupp (1999), a suitable way to estimate the mean and standard deviation from the 12 × 12 = 144 directions ϕ(l) of each window is to compute: $$z=\frac{1}{L}{\displaystyle \sum _{l=1}^L{\text{e}}^{2i\phi (l)}}\text{,}$$(17)

where L is the number of pixels within the sliding window. Defining a and b such that z = a + i · b, one computes the circular mean using $${\widehat{m}}_{\text{c}}(\phi )=\frac{1}{2}\arctan \frac{b}{a},$$(18)

and the circular standard deviation using $${\widehat{\sigma }}_{\text{c}}(\phi )=\sqrt{\frac{1}{2}(1-|z|)}.$$(19)

Moreover, to consider the uncertainty associated with the different directions φ(l), we use a weighted mean to compute z with weight inversely proportional to the variance. The error of the circular standard deviation is computed as follows: $$\text{Δ}\left({\widehat{\sigma }}_{\text{c}}(\phi )\right)=\frac{{\widehat{\sigma }}_{\text{c}}(\phi )}{\sqrt{2\cdot L-2}}.$$(20)

We show maps of the magnetic field angle, the mean angle computed using Eq. (18) and 12 × 12 window, and the circular standard deviation (Eq. (19)) in Fig. 4. As seen in the left panel of Fig. 2, the magnetic field direction differs along the edges of the expanding shell and the filament. A similar behavior we observe in the mean angle of the POS magnetic field is shown in the middle panel of Fig. 4. The standard deviation of the angle that is caused by the small-scale turbulence also shows different values within the borders of the H II region and the filament. In particular, the measured ${\widehat{\sigma }}_{\text{c}}$ is higher along the filament (more than 10 degrees), whereas it is a few degrees at the edges of the bubble. This result indicates that the magnetic field is possibly stronger at the edges of a bubble than in the filament. This result could also be a consequence of a complex geometry of NGC 2024: near the edge of a bubble, we observe a coherent limb-brightened structure, whereas, near the center of NGC 2024, we find a superposition of several components.

The ratio between the line width and the angle dispersion in Eq. (14) is shown in Fig. 5. The two panels show the measured ratios between the line width of the CN line (left) and HCO⁺ emission (right) to the dispersion of the angle of the magnetic field. We used the line width information derived from the second moment map shown in the right panels in Fig. 3. This approach gives a robust estimate of the range of the magnetic field strength assuming a specific number density of the gas because the second moment map provides information on the line width of all emission along the line of sight without considering possible impacts (see Sects. 4.1 and 4.2).

By assuming a gas number density of, for instance, 10⁴ cm⁻³, using Eq. (14), the ratio of 0.5 km s⁻¹ deg⁻¹ in Fig. 5 corresponds to 90 μG, whereas a ratio of 2 km s⁻¹ deg⁻¹ corresponds to 360 μG. We note, however, that using the second order moment for getting information about line width does not take into account the opacity of the corresponding line, or a possible spectral complexity.

Fig. 5 Ratio between measured line widths (right panels in Fig. 3) and the square root of the circular standard deviation of the magnetic field angle shown in the right panel of Fig. 4. The left panel shows the line width of CN divided by the square root of the rms of an angle, and the right panel shows the same as on the left panel, but here we use the HCO+ line width.

5.4 Magnetic field strength in NGC 2024

We calculate the strength of magnetic field across the beam-averaged area at the western (Bubble West) and eastern side (Bubble East) of the bubble, the filament (Filament), and across the region where the dust filament overlaps with the expanding shell of the H II region (Filament+Bubble) (Fig. 2, also see Sect. 4). We measure the uncertainty of magnetic field strength derived using the DCF method (Sect. 5.1) as follows: $$\text{Δ}{B}_{\text{POS}}=B_{\text{POS}}\cdot \sqrt{\frac{1}{4}\cdot \frac{{(\text{Δ}n)}^2}{n^2}+\frac{{(\text{Δ}(\text{Δ}v))}^2}{{(\text{Δ}v)}^2}+\frac{{\left(\text{Δ}\left({\widehat{\sigma }}_{\text{c}}(\phi )\right)\right)}^2}{{\widehat{\sigma }}_{\text{c}}{(\phi )}^2}},$$(21)

where Δn, Δ(Δν), and $\text{Δ}\left({\widehat{\sigma }}_{\text{c}}(\phi )\right)$ are measured uncertainties of η, Δν and ${\widehat{\sigma }}_{\text{c}}(\phi )$ (Sect. 5.3) respectively. The uncertainty of the magnetic field strength derived from Eq. (14) is similar to that given by Eq. (21), with an additional factor of $\frac{1}{4}$ in front of the last term in the equation.

We measure the uncertainties of n, σ_NT from the χ² minimization (see the bottom panel in Figs. C.1–C.8). We select the range in parameter space within which the χ² reaches 0.5 to derive corresponding n, and σ_ΝΤ. Thus, we found Δn to be 10³ cm⁻³ and Δσ_ΝΤ = 0.25 km s⁻¹ for both CN and HCO⁺.

The reported measurements of the angle dispersion and the POS magnetic field strength, including their uncertainties for both methods, are shown in Table 3. We show the comparison between the magnetic field strengths derived using the DCF and the ST method, including the use of different methods to estimate the angle dispersion in Fig. D.4. Overall, the POS magnetic field strength derived using the prescription from Skalidis & Tassis (2021); Skalidis et al. (2021) is lower than that derived using classical DCF. In the following, we comment on the magnetic field strengths derived using the Eq. (16) in Sect. 5.2. The B_POS strength is in the range of ~20–90 μG. The magnetic field strength inferred from CN is overall smaller than that derived using HCO⁺, mainly due to HCO⁺ showing slightly larger line widths than the CN in some regions, which directly impacts the strength of the magnetic field. Due to ambipolar diffusion effects (Zweibel 2004; Tritsis et al. 2023), it is expected that the emission lines of neutral molecular species (CN in our case) are wider than those of the ion molecular species (such as HCO⁺) (Li & Houde 2008; Yin et al. 2021). However, in our case, we observe that HCO⁺ is broader than the CN spectral lines, even after corrections for the opacity broadening. This result has two implications. Firstly, our observations did not spectrally resolve HCO⁺ emission, which potentially results in larger line widths caused by the blending of velocity components. Secondly, it is possible that our measurements do not probe ambipolar diffusion scale (~10⁻³ pc, Li & Houde 2008) in NGC 2024.

The strongest magnetic field, derived from CN and HCO⁺ is observed in the region at the intersection of the dusty filament and the edge of the expanding H II region and at the eastern side of the bubble. This result is driven mainly by the large densities and the broader CN and HCO⁺ lines observed in these regions (see Appendix C) found in this work. The magnetic field strength measured using both CN and HCO⁺ emission is stronger toward the eastern side of the bubble than on the western edge. We found lower line widths of CN and HCO⁺ at the western side of the bubble. The weakest magnetic field is derived toward the filament, although we note the B_POS derived from HCO⁺ in the western side of the bubble and filament to be comparable given the uncertainties.

6 Discussion

6.1 The magnetic field in NGC 2024

In this section, we discuss our results on the direction of magnetic field and its morphology, presented in Sect. 3. An overview of previous studies on the structure of magnetic field within NGC 2024 is reported in Meyer et al. (2008). Information about magnetic field in NGC 2024 were based on the thermal dust continuum emission (linear dust polarization), at 100 μm (Hildebrand et al. 1995; Dotson et al. 2000) and at 850 μm (Matthews et al. 2002), dichroic polarization of point sources (Kandori et al. 2007), and the Zeeman splitting of H I and OH lines (Crutcher et al. 1999).

All these studies found that the magnetic field shows a specific structure. In particular, the line of the sight (LOS) B-field strength in the central area of NGC 2024 dominated by the dusty filament (Crutcher et al. 1999) weakens from the northeast to the southwest. Matthews et al. (2002) investigated (sub)millimeter dust polarization and modeled magnetic field in NGC 2024 with two components. The first component is related to the dense and dust obscured part, where magnetic field lines follow locations of the FIR sources (right panel in Fig. 1). The second component of magnetic field is linked to gas affected by the stellar feedback. In our case, we observe regions where these two features are dominating: Band D dust polarization traces the filament in the central and southern part of NGC 2024, whereas it is impacted by the stellar feedback at eastern and the western parts at ionization front (Eastern and Western Loop in the right panel in Fig. 1).

We overlay dust polarization vectors at 100 μm from Hildebrand et al. (1995) and at 850 μm from Matthews et al. (2002) in the bottom row of Fig. A.6. The area where dust polarization measurements from Hildebrand et al. (1995) and Matthews et al. (2002) overlaps with our work is the central part of NGC 2024, where embedded FIR sources are located. We find the overall agreement in the direction of magnetic field inferred from dust polarization at 100, 154 and 850 μm.

POS magnetic field lines reported in our work are impacted by the stellar feedback at the surfaces of molecular cloud. Similarly, Crutcher et al. (1999) reported that the LOS magnetic field is nearly zero east of the filament in NGC 2024, which indicates a possibility that total magnetic field lines are in the POS east of the filament in NGC 2024. However, we note that observations presented in Crutcher et al. (1999) do not cover the edges of H II region, and that we require a study of LOS magnetic field component that covers larger field of view in NGC 2024.

The magnetic field lines are mainly parallel to the filament in NGC 2024, which is in the agreement with results presented in other studies (such as, Planck Collaboration Int. XXXII 2016; Santos et al. 2016; Pillai et al. 2015; Pattle et al. 2017). However, we observe a change in the magnetic field direction at the very southern part of the filament in Fig. 2, particularly in the Band D data, and also at the northern part of NGC 2024. Such changes in the direction of the POS magnetic field can indicate a few possibilites. Firstly, changes of the magnetic field often trace star formation (Pillai 2009; Ward-Thompson et al. 2017) due to gravitational collapse of the gas that cause the magnetic field lines to have “pinched” structure. Moreover, the variation of the magnetic field direction can be also a consequence of changes in the column density of the gas (Alina et al. 2019).

6.2 The magnetic field strength in the PDR and filament

In this section, we discuss our findings on the measured magnetic field strength in several distinct regions: the edge of the expanding PDR and the filament, including the region where it is not possible to clearly separate these two environments.

6.2.1 The edge of the HII region

In our work, we estimate the POS magnetic field strength at the border of the expanding H II region to be ~20–40 μG on the west side, and ~82 μG on the east side. The good agreement between the CN and HCO⁺ B_POS strengths in these regions, particularly at the eastern side, could indicate that these two molecular lines are tracing similar gas. Nevertheless, the factor of almost 2 difference in the B_POS measured in the west and the east could be due to different gas densities at these sides of the bubble. As the eastern side of the bubble is denser than the western side (see, Meyer et al. 2008), we expect that the gas here is less impacted by the incoming radiation. In addition, the opacity of the HCO⁺ emission is notably higher (τ > 4) on the western side than on the eastern side (τ = 2).

As previously presented in Sect. 5, the eastern edge of H II region has a stronger magnetic field than the western side. This result can also be due the total magnetic field changing its direction with respect to the line of sight. A previous study showed that the line of sight (LOS) magnetic field changes its strength from being zero at the eastern part of NGC 2024, to 100 μG at the western side, which is indicative of the change of a direction of the overall magnetic field (Crutcher et al. 2009). Moreover, the molecular gas content located west from the center of NGC 2024 has a lower density than the gas located at the eastern side, suggesting that stellar radiation has stronger impact on the gas on the western side of the bubble (Crutcher et al. 2009). At the eastern edge of the bubble, due to higher gas densities (Table 2), stellar feedback did not sweep up gas as far as on the other side (Barnes et al. 1989).

The values of B_POS calculated in our work are generally lower than those reported in studies of other PDRs. For example, the magnetic field strength measured from the Zeeman splitting of H I and OH in M 17 is ~750 μG and ~250 μG respectively (Brogan & Troland 2001), and ~(1000 – 1700) μG using dust polarization data (Hoang et al. 2022). The magnetic field strength of the PDR in the Horsehead nebula (SMM1, Hwang et al. 2023) is a few tens of μG, and comparable to our results in NGC 2024. That study used C¹⁸O to derive the line width and dust column density and effective radius to estimate the gas volume density. Although the edges of the bubble in our work and the Horsehead nebula have comparable densities (a few × 10³ cm⁻³), measured line widths are different. It is worth pointing out that C¹⁸O and CN and HCO⁺ are tracing different gas (Philipp et al. 2006; Bron et al. 2018). C¹⁸O traces more compact structures than the CN and HCO⁺, and it gets destructed by the stellar feedback, contrary to CN and HCO⁺, whose emission becomes enhanced in these regions (Bron et al. 2018). In addition, Hwang et al. (2023) used a modified DCF method that measures the ratio of the ordered and turbulent component of the B field (Hildebrand et al. 2009) to derive the magnetic field strength.

Other PDR regions impacted by super stellar clusters have stronger magnetic field to those reported in our work. For example, Pattle et al. (2018) reported strong magnetic fields of a few hundred μG in the Pillars of M 16. This work used the DCF method to derive the magnetic field strength. Gas density of the Pillars are somewhat higher than those in our work 5 × 10⁴ cm⁻³ (Ryutov et al. 2005). Similarly, the line width measurements used in Pattle et al. (2018) are taken from the Gaussian fitting of several molecular lines from and these are wider than those reported in our work (see Table 3 in White et al. 1999, reported line widths in the range from 1.2 to 2.2 km s⁻¹). Guerra et al. (2021) showed the POS magnetic field strength across the Orion Bar PDR in OMC-1 is of a few hundreds of μG, also using the DCF method.

The main difference in magnetic field strengths computed in our work and found within the literature comes from the methods used to derive the POS B field. The DCF method generally overestimates the magnetic field strength. On the other hand, simulations have shown that the POS B field is proportional to inverse square root of the angle dispersion. However, in this case, the prefactor is smaller by a factor of 5, resulting in generally lower B_POS. We show how B_POS varies with the method we selected in Fig. D.4.

6.2.2 Filament

The filament going across NGC 2024 is super-critical (Orkisz et al. 2019), which means it is gravitationally unstable and a potential site for star formation. We measure the lowest values of the POS magnetic field (~30–50 μG) in this region, suggesting that the magnetic field cannot dominate gas dynamics. We return to this point in Sect. 6.3. These low values in the magnetic field strength come from the largest dispersion angles and narrow spectral lines computed for this region (Tables 2 and 3). Large dispersion angles suggest that the turbulence does not impact the magnetic field, which is also supported by the narrow CN and HCO⁺ lines. Since the turbulence level in this part appears low, this could lead to star formation in the filament, previously confirmed at its southern part (Hwang et al. 2023).

Overall, our measurement of the magnetic field strength in the filament is lower than in other filamentary structures. For example, Pattle et al. (2017) measured a POS magnetic field strength of a few mG in the Orion A filament. Similarly, as discussed in the previous section, it is important to highlight that Pattle et al. (2017) used C¹⁸O emission to measure the turbulence and number density and different method to derive the B_POS. Similarly, Ching et al. (2018) investigated magnetic field toward the DRL1 filament and found B_POS = 600 μG.

6.3 Magnetic field support in NGC 2024

Next, we investigate the role of the magnetic field at the edges of the bubble, in the filament, and the overlap region (bubble and the filament) in NGC 2024. To do so, we compute several parameters that describe gas stability in relation to local environmental conditions: radiation, magnetic and turbulent field. Since we do not have information on the line of sight component of magnetic field, we cannot derive the total magnetic field strength. Therefore, the following quantities need to be treated as upper limits, as they are proportional to the magnetic field strength as B⁻ⁿ.

We make an estimate on the mass-to-flux ratio of gas in NGC 2024 $\left(\frac{M}{\text{Φ}}\right)$ and compare it to the critical mass-to-flux ratio. The critical mass-to-flux ratio ${\left(\frac{M}{\text{Φ}}\right)}_{\text{crit}}$ depends on the assumed geometry. For instance, for a uniform disk, ${\left(\frac{M}{\text{Φ}}\right)}_{\text{crit,cyl }}=\frac{1}{2\pi \sqrt{G}}$ (Nakano & Nakamura 1978; Joos et al. 2012; Hanawa et al. 2019). In addition, for a spherical geometry, the critical mass-to-flux ratio will be expressed as ${\left(\frac{M}{\text{Φ}}\right)}_{\text{crit,sph }}=\frac{1}{3\pi \sqrt{\frac{G}{5}}}$ (Mouschovias & Spitzer 1976). The μ_Φ parameter is the ratio between estimated mass-to-flux and ${\left(\frac{M}{\text{Φ}}\right)}_{\text{crit}}$ and can be computed as: $${\left({\mu }_{\text{Φ}}\right)}_{\text{cyl }}=\frac{(M/\text{Φ})}{{(M/\text{Φ})}_{\text{crit }}}=7.6\times {10}^{-21}\frac{N}{B},$$(22)

in the case of cylindrical geometry. N is the column density of molecular gas taken from Lombardi et al. (2014) in units of cm⁻², and B is measured POS magnetic field strength (Eq. (16)) in μG. In the case of spherical geometry, using the expression for ${\left(\frac{M}{\text{Φ}}\right)}_{\text{crit}}$ and the above equation: $${\left({\mu }_{\text{Φ}}\right)}_{\text{sph}}=0.67\cdot {\left({\mu }_{\text{Φ}}\right)}_{\text{cyl}}\text{. }$$(23)

Therefore, the μ_Φ will vary by a factor of 0.67 that comes from the assumed geometry. It is necessary to point out that there is not a critical mass-to-flux ratio for lateral contraction of a filament when it threaded by magnetic field parallel to its axis of symmetry (Mouschovias & Morton 1992). In our work, we use Eq. (22) and assume the cylindrical geometry, but we note that this factor can vary depending on the assumed geometry of a system. Moreover, in case of the filament, the interpretation and physical meaning of mass-to-flux ratio in not straightforward, and it should be treated with caution.

Next, we compute the Alfvénic Mach number, ℳ_A: $${ℳ}_{\text{A}}=\frac{\sqrt{3}\cdot {\sigma }_{\upsilon ,\text{NT}}}{{\upsilon }_{\text{A}}}\text{.}$$(24)

where σ_v,ΝΤ is the non-thermal velocity dispersion reported in Table 2. The $\sqrt{3}$ factor comes from the assumption of isotropic turbulence (Crutcher 2004; Stewart & Federrath 2022). υ_Α is the Alfvén velocity defined as $${\upsilon }_{\text{A}}=\frac{B}{\sqrt{4\pi \rho }},$$(25)

where B is the total magnetic field strength in units of G, and ρ is the gas mass volume density in units of g cm⁻³. Equation (24) can also be represented as a ratio between turbulent gas and magnetic energies. Therefore, the Alfvénic Mach number provides information about a dominant driver of gas flows.

Finally, we compute the plasma-beta parameter, β_B, that gives information about the ratio of thermal and magnetic pressure: $${\beta }_{\text{B}}=\frac{n{k}_{\text{B}}T}{B^2/(8\pi )},$$(26)

where n is the number density, Τ is the gas temperature (Table 2), and B magnetic field strength (Table 3).

We report our results in Table 4 for values of magnetic field strengths computed from CN and HCO⁺ respectively using approach from Skalidis & Tassis (2021); Skalidis et al. (2021). We estimate μ_Φ < 1 in each region. This result could imply the presence of magnetically supported gas (Pattle et al. 2019). In general, it is worth noting that we see a difference of one and two orders of magnitude in the μ_Φ between the edges of the bubble and the filament and the overlap region.

We find M_A greater than 1 in all regions, which suggests that the magnetic field does not govern the gas motions. Moreover, we measure a plasma-beta parameter lower than 1 in almost all environments, which agrees with the lower values of μ_Φ, particularly at the edges of the bubble. A plasma beta derived from CN measurements of magnetic field at the western side of the bubble is slightly higher than one. Low values of plasma-beta suggest that the magnetic energy dominates the thermal motions of the gas.

The results reported in Table 4 should be taken as upper limits considering we calculated them using the POS component of the magnetic field. In particular, the reported mass-to-flux ratio must be taken with caution. The complex geometry of NGC 2024, large uncertainties of reported POS magnetic field strength, and non availability of the LOS component of the B field have a big impact on values reported in Table 4. Our results imply that the edges of the bubble show different properties than the filamentary structure and the overlap region in NGC 2024. These regions are indeed different in terms of the impact of stellar feedback. However, further quantification of the gravitational stability of the gas in these regions requires systematic analysis of the magnetic field. This also includes both POS and LOS components in these regions and constraints on the geometry of NGC 2024, which is beyond the scope of this study.

The properties of the gas we study in our work vary with the environment and this gas shows different characteristics. Our results thus highlight the non-negligible role of the magnetic field in NGC 2024 in regions impacted by the stellar feedback. The edges of the bubble in NGC 2024 are magnetically dominated and gas in these regions is magnetically supported against gravitational collapse. However, the high values of the Alfvénic Mach number suggest that these regions are highly turbulent. The gravitational stability of gas located at the shell of expanding H II region is not always the case. For example, the expansion of the H II region can also trigger star formation, as seen in Galactic H II regions RCW82 (Pomarès et al. 2009) and RCW120 (Figueira et al. 2017). Since the H II region in NGC 2024 is relatively young, it is possible that the stellar feedback has not yet been strong enough to trigger star formation. Nevertheless, it is crucial to point out that our result is based on CN and HCO⁺ measurements, which trace the UV-illuminated and UV-shielded gas, which do not necessarily trace the star-forming gas.

The environment in which we observe gas coming from the edge of the bubble and the filament is also magnetically supported, which could imply that the gas in this region mainly comes from the ionization front. This result can be explained by the presence of compressive motions in NGC 2024 (Orkisz et al. 2017) that possibly originate from expansion of the H II region.

The role of the magnetic field is somewhat different in the filament. High μ_Φ > 1 implies gravitational instability that could lead to star formation, as observed at the southern part of this filament (Hwang et al. 2023). Measured Alfvénic Mach numbers are also higher than 1 in the filament. ℳ_A measured from CN and HCO⁺ is also higher than 1.

Nevertheless, we should note the following. Firstly, in this work we do not investigate the central region of NGC 2024, where the protostellar candidates and ongoing star formation is observed. Therefore, our results focus only on the gas impacted by the stellar feedback and one located in the front of the bubble and a filament. We also note that due to relatively high uncertainties (around 30%) of measured magnetic field strengths presented in Table 4, these uncertainties impact our measurements, particularly they could change our presentation of the physical conditions in the filament. Given the uncertainties, it is possible that the filament is in the transition zone between being fully magnetically supported and gravitationally unstable.

7 Conclusions

We present new SOFIA HAWC+ dust polarization measurements across the NGC 2024 H II region and associated molecular cloud in the Orion B molecular cloud. We combined these measurements with molecular data observations from the ORION-B Large Program on the IRAM 30-meter telescope. We summarize our findings obtained using these observations as follows:

• 1.

  Our results focus on a subset of environments found in NGC 2024, particularly on the shell near the edge of the H II region and the filament in the front of NGC 2024, which are not locations containing protostellar cores (Fig. 1 and sites of active star formation, located in the center of NGC 2024;

• 2.

  We investigate the magnetic field morphology traced by dust polarization from HAWC+ observations using Band D at 154 μm and BandE at 216 μm. The direction of the magnetic field derived from these two bands shows a good agreement;

• 3.

  We use HAWC+ Band D at 154 μm to characterize the geometry of the magnetic field across NGC 2024. We find that the structure of the magnetic field is ordered and follows the morphology of the expanding H II region and the direction of the filament;

• 4.

  Using the CN(1−0) and HCO⁺(1−0) molecular emission obtained using the IRAM 30-meter telescope, we characterize physical conditions (turbulence and gas density) in four specific regions in NGC 2024: the edges of the expanding H II shell located to the east and the west, the filament to the south, and the environment in the northern part of NGC 2024. In our analysis, we include collisions with electrons and the hyperfine structure of CN(1−0), which are essential aspects of our calculations. Both CN(1−0) and HCO⁺(1−0) emission lines are optically thick across NGC 2024, which broadens the observed line widths. The gas number density derived from CN(1−0) is comparable or somewhat higher than that obtained from analyzing the excitation of HCO⁺(1−0);

• 5.

  We derived the POS magnetic field strength using the ST (and classical DCF) method, with values ranging from 20 to 90 μG (50–240 μG) in the environments mentioned above. The strongest magnetic field is found in the region composed of the dusty filament and the edge of the expanding H II region, located in the northern part of NGC 2024. The high magnetic field strength derived in this area is driven by the largest gas densities and line widths;

• 6.

  Magnetic field strengths derived using CN(1–0) and HCO⁺ (1–0) are comparable within the uncertainties at the edges of the bubble, especially in the eastern side of the bubble. The magnetic field strength calculated from HCO⁺(1–0) is generally higher than that inferred from CN(1–0). We note that B_POS measured from CN(1–0) and HCO⁺(1–0) is larger on the eastern side of the bubble than in the western edge. This observed difference may result from changes in the magnetic field direction indicated in previous studies. In addition, the western side of the bubble is more impacted by stellar radiation due to its lower density;

• 7.

  By analyzing the mass-to-flux ratio, Alfvénic Mach number, and plasma-beta parameter, we find that the edges of the bubble and the filament show different properties. Gas impacted by the stellar feedback and traced by the CN(1–0) and HCO⁺ (1–0) emission seems better supported against gravitational collapse by the magnetic field than the gas in the filament, which represents a location where star formation can take place. However, we find the Alfvénic Mach number higher than one in all regions, which suggests that the magnetic field does not control the gas motions. Our results should be treated as upper limits given the use of the POS component of the B field. In addition, the estimated mass-to-flux ratio has a large uncertainty due to the specific geometry of NGC 2024. We also note that these regions of NGC 2024 could be in the transition phase regarding gravitational stability.

Our research emphasizes the importance of utilizing dust polarization measurements to characterize the structure of the magnetic field and the need to combine these measurements with molecular data to accurately infer the magnetic field’s strength. Our results demonstrate that the magnetic field plays a critical role and that its contribution and other factors cannot be ignored. We also demonstrate the significant impact that different methods used to derive the POS magnetic field can have on our results and how it affects the interpretation of the gravitational stability of the gas based on the corresponding parameters. Moreover, by removing the contribution of large-scale gradients of the magnetic field direction using techniques such as the sliding window, we can consider only the local magnetic field. In addition, we highlight the necessity of conducting a careful analysis of the line radiative transfer and using appropriate rate coefficients for inelastic collisions with H2 and electrons, as well as line opacity corrections.

Acknowledgements

We thank the anonymous referee for useful comments that have improved the manuscript. This research has made use of spectroscopic and collisional data from the EMAA database (https://emaa.osug.fr/ and https://doi.osug.fr/public/EMAA/EMAA.html). EMAA is supported by the Observatoire des Sciences de l’Univers de Grenoble (OSUG). This work was supported by the French Agence Nationale de la Recherche through the DAOISM grant ANR-21-CE31-0010. This research was carried out in part at the Jet Propulsion Laboratory, which is operated by the California Institute of Technology under a contract with the National Aeronautics and Space Administration (80NM0018D0004). MGSM and JRG thank the Spanish MICINN for funding support under grant PID2019-106110GB-I00. Based in part on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (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 2002 to the University of Stuttgart. Financial support for SC was provided by NASA through award #08_0186 issued by USRA.

Appendix A HAWC+ dust continuum and polarization

We present a comprehensive analysis of the flux that is filtered out in our dust polarization measurements using SOFIA HAWC+ in Sec. A.1. Then, we show HAWC+ Band E dust continuum data, and provide comparison with Band D in Sec. A.2. The comparison between the dust polarization shown in this work with the literature data is described in Sect. A.3.

Appendix A.1 Filtering of low-level extended flux in HAWC+ images

Dust continuum images obtained with ground-based or airborne far-infrared total power continuum cameras suffer from some degree of correlated atmospheric noise, which cannot be completely separated from the extended low-level emission of the source. The HAWC+ instrument data reduction pipeline attempts to remove correlated components from the time stream using the algorithms described in Kovács (2008). To estimate the magnitude of the resulting spatial filtering at different flux levels in our HAWC+ Stokes I images, we compared the Band D image to that obtained using the 160 μm channel of the PACS instrument on Herschel (OBSID 1342206080). Only high S/N pixels are included in the dust polarization analysis (see Sec. 2.1). Fig. A.1 shows a pixel-by-pixel comparison of the HAWC+ and PACS fluxes at the spatial resolution of the HAWC+ Band D data. Given the difference in the filter passbands and overall calibration uncertainties of the two instruments, we scaled the HAWC+ data by a factor of 1.12 to enforce a unity flux ratio at pixels with fluxes above 80,000 MJy sr⁻¹ (marked by the dotted line in the upper-right corner of the figure).

A good linear correlation between the two data sets is seen for points with fluxes above about 60,000 MJy sr⁻¹. However, at lower flux levels, the HAWC+ points consistently fall below the slope-one line marked in red. To quantify this effect, we computed average HAWC+/PACS flux ratios in two narrow (±5%) intervals centered at PACS fluxes of 30,000 and 10,000 MJy sr⁻¹, marked with vertical dotted lines in Fig. A.1. The average HAWC+/PACS flux ratios in the two intervals are 0.81 and 0.57, respectively.

The effect of the filtering of extended flux in the Stokes Q and U images is impossible to estimate, given the absence of space-based polarization data at a sufficient angular resolution. However, the filtering likely affects the polarization fraction more strongly than the polarization angle. Consequently, we do not use the polarization fraction in our analysis. Using the PACS image as a prior in the data reduction pipeline would improve the fidelity of the HAWC+ Stokes I images. However, such a future modification to the pipeline is unlikely.

Appendix A.2 SOFIA HAWC+ Band E dust continuum

We show SOFIA HAWC+ Band E dust continuum map in Fig. A.4 at its native resolution of 18.7″. This map is processed and S/N masked the same way as Band D data, as described in Sec. 2.1.

To compare dust polarization angles derived from Bands D and E observations, we first convolve Band D to match the lower spatial resolution of the Band E data. Next, we compute a difference between the two angles, θ_D – θ_Ε. Finally, we measure the rms of this difference, as described in Sec. 5.3, Eq. 19. We show a map of θ_D – θ_Ε in the left panel of Fig. A.5. The right panel in Fig. A.5 shows the corresponding rms of the map shown in the left panel.

Fig. A.1 Pixel-by-pixel comparison of the HAWC+ Band D and PACS 160 μm fluxes. The PACS data have been convolved to the spatial resolution of the HAWC+ image. Only high S/N pixels included in the polarization analysis are shown. As described in the text, the HAWC+ fluxes were scaled by a factor of 1.12.

Appendix A.3 Comparison with the literature data

In Fig. A.6, we show a comparison to our SOFIA HAWC+ data with other polarization measurements from the literature. The background in all panels is Stokes I map of Bands D (left panels) and E (right panels), whereas thin black lines show the direction of the magnetic field. We overlay NIR polarization measurements (Kandori et al. 2007) as thick black lines and the position of the outflow indicated by the dark blue contours on the top panels of Fig. A.6. The bottom panels show the zoom-in of the inner region in NGC 2024 (indicated by the black dashed rectangle in the top panels) and 100 μm (brown lines, Dotson et al. 2000) and magnetic field direction inferred from the 850 μm dust continuum observations in purple (Matthews et al. 2002).

Fig. A.2 Comparison between Stokes parameters, I, Q, and U (each row) for two setups of data reduction (the one used in this work is shown in the left panel, and the other one in the middle panel; see Sec. 2.1). The right panels show the difference between the left and middle panels (computed for each Stokes parameter).

Fig. A.3 Same as in Fig. A.2, but for the polarization angles (see Eq. 6 in Sec. 2.1).

Fig. A.4 SOFIA HAWC+ 214 μm (Band E) dust continuum map at 18.2″ angular resolution corresponding to linear scales of ~0.037 pc. The map is masked the same as the Band D map (Fig. 2, see Sec. 2.1). Black lines in both panels show the orientation of the magnetic field for every fifth pixel.

Fig. A.5 Difference in polarization angles measured from Band D and E data is shown in the left panel. The right panel shows the rms of the angle differences.

Fig. A.6 SOFIA HAWC+ dust continuum Band D and Band E measurements (left and right panels in the top row, respectively), but with overlayed positions of protostars and NIR polarization (dark purple circles and lines, Kandori et al. 2007). Dark blue contours represent the outflow observed in the HCO+ emission. Light blue points indicate the positions of FIR sources. Black dashed rectangles show the central area and the zoom-in panels on the bottom row. Here, we show magnetic field lines from the FIR dust polarization at 100 μm (Dotson et al. 2000) in dark magenta and from (sub)millimeter dust polarization at 850 μm (Matthews et al. 2002) in black.

Appendix B Spectral line fitting

Fig. B.1 Fitting parameters for the hyperfine CN emission in NGC 2024 measured for both velocity components (top and bottom row each). The first parameter, p1 is shown in the left panel, followed by the centroid velocity, FWHM and the final parameter, the opacity (p4).

We fit the CN and HCO⁺ emission lines using the CUBE software, part of GILDAS/CLASS, currently under development by IRAM in Grenoble. We describe the CN fitting procedure in Sect. B.1. Additionally, we use the Semi-automated multi-COmponent Universal Spectral-line fitting Engine (SCOUSE, Henshaw et al. 2016, 2019) for fitting the HCO⁺ emission prior to using CUBE, which is further explained in Sec. B.2.

Appendix B.1 Hyperfine structure of the CN line

The CN molecule has a hyperfine structure, and the important parameters of each component of the multiplet studied in this work are shown in Table 1 in Sec. 4. Prior to fitting the hyperfine structure, we assume the following. The first assumption is that the components of a multiplet do not overlap with each other. Second, we assume the same excitation temperature for all hyperfine components in the multiplet and the same line width.

Four parameters describe the hyperfine structure model: p₁, p₂, p₃, and p₄. The first parameter is the antenna temperature multiplied by the optical depth of the CN emission: $$p_1=T_{\text{ant }}\cdot \tau$$(B.1)

The second parameter is the centroid velocity of the main hyperfine component (the F = 5/2 − 3/2 transition): $$p_2={\upsilon }_{0,F=5/2-3/2}\text{.}$$(B.2)

The next parameter is the full width at half maximum (FWHM) of the hyperfine components: $$p_3=FWHM\text{,}$$(B.3)

and the last parameter is the opacity of all components of the studied multiplet: $$p_4=\tau .$$(B.4)

We assume that the optical depth of each component can be described using a Gaussian function of velocity: $${\tau }_{\text{i}}(\upsilon )={\tau }_{\text{i}}\cdot \exp \left[-4\ln 2{\left(\frac{\upsilon -{\upsilon }_{0,\text{i}}}{p_3}\right)}^2\right]\text{.}$$(B.5)

The opacity of the multiplet is then calculated as a sum of opacities of all hyperfine components: $$\tau ={\displaystyle \sum _{i=1}^N{\tau }_i}$$(B.6)

where τ_i is the opacity of the i-th component.

Fig. B.2 Fitting parameters for the HCO+ emission assuming a Gaussian line profile. Each row shows parameters for each fitted component. The results in the first row represent the brightest HCO+ component. The second row shows the fainter component and the outflow observed close to the center of NGC 2024.

Finally, the antenna temperature, T_ant is thus derived as the ratio between the first and the fourth fitting parameters, p₁, and ^p₄: $$T_{\text{ant}}=\frac{p_1}{p_4}\left(1-e^{-\tau (\upsilon )}\right).$$(B.7)

To get the excitation temperature, we mask all pixels having the opacity higher than 10 and lower than 0.2 to avoid degeneracies. The excitation temperature is derived assuming the Local Thermodynamical Equilibrium (LTE) case: $$T_{\text{ex}}=\frac{hv}{k}{\left\{\ln \left[1+\frac{hv}{k{T}_{\text{ant}}}\left(1-e^{-{\tau }_v}\right)\right]\right\}}^{-1}.$$(B.8)

We fit the CN emission in CUBE that uses an optimized version of the minimization method taken from the MINUIT system of CERN. To do so, we load the data into CUBE and run the commands \fit/minimize and add \hfs to specify we want to fit the hyperfine structure. Therefore, we also provide a file containing information about the components of the hyperfine multiplet, shown in Table 1. In the case of the hyperfine fitting, CUBE fits four parameters, p₁, p₂, p₃, and p₄, described in Eq. B.1, B.2, B.3, and B.4 respectively.

Before fitting, we define a two-dimensional mask based on the S/N of the CN integrated intensity. We select the threshold of 10 and fit two CN components in the area that goes inside the mask and one component outside the mask. By adding this step, we provide CUBE with additional information about the area where the CN emission is bright enough to observe double components in its spectrum. CUBE requires the initial values of the free parameters for each component we want to fit. We specify our initial guess for the centroid velocities of the components to be 10km/s and 7km/s. However, we do not specify initial conditions for the rest of the fitting parameters and let CUBE find the best possible values. We show the results of the CN emission fitting for both components in Fig. B.1.

Appendix B.2 Gaussian lines of the HCO⁺ emission

The HCO⁺ emission line can be described using the Gaussian function. Similarly, as in the case of the CN, by inspecting the HCO⁺ data cube, we notice the presence of two velocity components across the large portion of the map. In addition, we observe a region in which HCO⁺ spectra contain extended line wings, suggestive of the molecular outflow also seen in the CO emission, whose presence is also known from previous studies (for instance, see, Richer et al. 1992). Different velocity components we observe correspond to the primary and intermediate velocity layers seen in the CO emission and its isotopologues across NGC 2024 (Gaudel et al. 2023).

Before fitting the Gaussian line to the HCO⁺ spectra, similarly as in the case of CN, we want to locate regions where HCO⁺ the emission has more than one peak. Therefore, we decompose the HCO⁺ emission using SCOUSEPY. For the full description of the fitting procedures in SCOUSEPY, we refer to work by Henshaw et al. (2016, 2019). In the first step, SCOUSEPY divides the spectral cube into spectral averaging areas (SAAs). The spectrum of each SAA is the average spectrum of all pixels found within the SAA. Then, by assuming the shape of the spectral line, CUBE decomposes the spectrum of each SAA, identifies a number of velocity components, and suggests a model in a second step. In this step, the user can modify a model suggested by SCOUSEPY. Next, in the third step, SCOUSEPY fits emission in each pixel based on the fitting model of each SAA. The final step allows the user to check the fitting result within each pixel and, similarly to the second step, modify the model if needed.

We fit a Gaussian line profile to the HCO⁺ emission in SCOUSEPY. The output parameters describing a Gaussian function are each component’s peak temperature, centroid velocity, and FWHM. Additionally, SCOUSEPY computes the rms, S/N and the residuals. Based on these results we derive from SCOUSEPY, we define two two-dimensional masks. The first mask contains a region where SCOUSEPY identifies three velocity components. The second mask is a region within which SCOUSEPY finds two velocity components. We use these masks as the input to CUBE to mark regions where we observe three (first mask), two (second mask), and one velocity component (outside these two masks).

Similarly, as for the CN emission, we use the command \fit/minimize in CUBE and add \gaussian to specify the shape of the spectral line. In CUBE, the Gaussian line is described using three parameters: the area (in Kkm/s), the centroid velocity, and the FWHM. We also specify initial guesses for fitting parameters. In the last step, we select the brightest component to be the first, and then combine results of the second (fainter) component and the outflowing feature into one map. We show the results of this fitting procedure in Fig. B.2.

Appendix C Radiative transfer modeling of CN and HCO⁺ emission

Here, we show results from the radiative transfer modeling of CN and HCO⁺ excitation, presented in Sec. 4.3. For each region we study in this work (edges of the bubble, filament, and the overlap region), we show how our output parameters from RADEX (excitation temperature, opacity, and the peak temperature) vary as a function of input parameters (line width, column and volume density) in Fig. C.1, C.2, C.3 and C.4 for CN emission. In the case of HCO⁺ emission, we show these results in Fig. C.5, C.6, C.7 and C.8. We also include a χ² as a function of the input parameters at the bottom row in each figure.

Based on the results presented in Sec. 4.3, we computed model spectrum of CN and HCO⁺ for edges of the bubble, filament and the overlap area. We show the beam-averaged spectrum of CN and HCO⁺ of each region studied in this work and their model in Fig. C.9.

Fig. C.1 Parameter space used for the radiative transfer modeling of CN in the region at the western side of the bubble. The first row shows results for the excitation temperature. The color bar shows values of excitation temperature, Tex for the grid of volume densities (x-axis), line widths (y-axis) and column densities (each panel) for fixed electron fraction. The second row shows results for opacity, τ, and the third row shows the peak temperature, Tpeak. The last row shows the χ2 minimization of the modeled peak temperature and the observed value computed from the CN(1−0) spectra.

Fig. C.2 Same as in Fig. C.1, but for the eastern part of the bubble.

Fig. C.3 Same as in Fig. C.1, but for the filament.

Fig. C.4 Same as in Fig. C.1, but for the overlap region.

Fig. C.5 Same as in Fig. C.1, but for HCO+ (1−0).

Fig. C.6 Same as in Fig. C.5, but for the eastern side of the bubble.

Fig. C.7 Same as in Fig. C.5, but for the filament.

Fig. C.8 Same as in Fig. C.5, but for the overlap region.

Fig. C.9 Beam averaged spectrum of CN and HCO+ taken within four different regions shown in the left panel of Fig. 1 and in Fig. 2. Colored line represent the corresponding model spectrum infrerred from the non-LTE RADEX modeling for physical conditions presented in Sec. C and shown in Table 3. Dashed lines show the corresponding residuals.

Appendix D Measuring magnetic field strength in NGC 2024

In this part, we provide additional information on the sliding window technique in Sec. D.1, an alternative approach to measuring the angle dispersion using the histogram analysis in Sec. D.2, and results of magnetic field strength computed using different methods in Sec. D.3.

Appendix D.1 Sliding window

We investigate the change of the angle dispersion measured using the sliding window (Sec. 5.3) by varying its size. We used the size of the sliding window from 1 to 3.5 beam sizes, with a step of 0.5, which corresponds to sizes from 4 × 4 to 14 × 14 pixels, respectively. We show these results in Fig. D.1. The beam-sized sliding window is too small to remove any large-scale contribution of the magnetic field. We find that the rms of the angle changes by a factor of around 3 as we increase the size of the slidding window.

Appendix D.2 Histogram analysis

In this section, we describe an alternative approach to compute ${\widehat{\sigma }}_{\text{c}}(\phi )$. We consider all points within each circularly shaped region and create a histogram showing the distribution of magnetic field angles. Next, we compute the circular mean and fit the Gaussian function to the distribution of angles. We show the histograms for each beam-averaged area we used in this work, including the fitting parameters of each distribution in Fig. D.2. In most cases, the histogram corresponds to a single Gaussian. In Fig. D.3, we show comparison between ${\widehat{\sigma }}_{\text{c}}(\phi )$ from sliding window (y-axis) and histogram analysis (x-axis). The derived standard deviations from the sliding window (Sec. 5) are higher than those derived from the histogram analysis.

Fig. D.1 Angle rms as a function of the size of the sliding window.

Fig. D.2 Top left: Bubble on the west. Top right: Bubble on the eastern side of NGC 2024. Bottom left: Filament. Bottom right: Mixture of bubble and dusty filament.

Fig. D.3 Comparison of the angle dispersion using two different approaches, sliding window (y-axis, see Sec. 5) and histogram analysis (x-axis).

Fig. D.4 Magnetic field strength derived using different methods (Sec. 5.1, D.3) for computing the magnetic field strength and the angle dispersion (Sec. 5.3, D.2) as a function of the environment we analyze in this work.

Appendix D.3 B_POS using different approaches

We compare results on POS B-field strengths using several methods in Fig. D.4. As described in Sec. 5, results derived from the classical DCF are higher than those derived using ST, even when using the histogram analysis instead of the sliding window approach.

References

All Tables

All Figures

Fig. 1 Multiwavelength image of NGC 2024. Left panel: Color-composite image of the Flame nebula (NGC 2024) showing peak intensities of 12CO (1−0) (blue) emission and isotopologues, 13CO (1−0) (green) and C18O (1−0) (red) obtained by the IRAM 30-meter telescope (image credits: Pety et al. 2017). We overlaid the SOFIA HAWC+ field of view as a white dashed rectangle. Gray contours show a network of filaments presented in Orkisz et al. (2019). We label the regions we investigate in this work: Bubble East, Bubble West, Filament, and Bubble+Filament. Right panel: ESO-VISTA image (ESO/J. Emerson/VISTA, Cambridge Astronomical Survey Unit). We sketched the environments seen across NGC 2024: the H II region (white dashed circle, see Table 2 in Gaudel et al. 2023) and the filament (orange line, Orkisz et al. 2019). Light blue points are the positions of far-infrared sources in the background (Mezger et al. 1988, 1992). In addition, we labeled the edges of the expanding H II region and the direction of the stellar feedback driven by the radiation produced by recent star formation.

In the text

Fig. 2 SOFIA HAWC+ 154 μm (Band D, all panels) dust continuum maps at 13.6″ angular resolution. The maps are masked based on the measured S/N in measured Stokes intensity and polarization angle. Black lines in all panels show the orientation of the magnetic field for every fifth pixel. We show Band D dust polarization map across NGC 2024 (top panel). Bottom rows show the zoomed-in regions we analyse in this work: egde of the bubble on the western (middle left panel), eastern (middle right panel), in the filament (bottom left panel), and in the overlap region of the filament and bubble (bottom right panel). A circle in each of these zoomed-in panels shows the region we use to compute the dispersion in the mean angle of magnetic field, as well as number density and turbulent velocity dispersion.

In the text

Fig. 3 Moment maps of 12CN(1−0) and HCO+ (1−0) emission across NGC 2024. Top: maps of the CN emission: integrated intensity (the zeroth moment, left panel), the centroid velocity (the first moment, middle panel), and the FWHM (the second moment, right). Bottom: the same as in the top row, but for HCO+. We show the beam size for both molecular lines at the bottom left corner of each panel and the 0.5 pc scalebar at the bottom right corner. All pixels shown in these maps result from the masking technique presented in Einig et al. (2023). We refer the reader to the description of the production of these moment maps in Sect. 4.1.

In the text

Fig. 4 Direction of the magnetic field (left panel), the mean angle of the magnetic field (middle panel), and the rms of the magnetic field direction (right panel). Maps in the middle and right panels have been made using a 12 × 12 pixels sliding window with weight inversely proportional to the measured uncertainty variance. All maps have been regridded to match the grid size of the CN and HCO+ data. The gray contour indicates the mask we defined for dust polarization measurements, described in Sect. 2.1.

In the text

	Fig. 5 Ratio between measured line widths (right panels in Fig. 3) and the square root of the circular standard deviation of the magnetic field angle shown in the right panel of Fig. 4. The left panel shows the line width of CN divided by the square root of the rms of an angle, and the right panel shows the same as on the left panel, but here we use the HCO+ line width.
In the text

	Fig. A.1 Pixel-by-pixel comparison of the HAWC+ Band D and PACS 160 μm fluxes. The PACS data have been convolved to the spatial resolution of the HAWC+ image. Only high S/N pixels included in the polarization analysis are shown. As described in the text, the HAWC+ fluxes were scaled by a factor of 1.12.
In the text

	Fig. A.2 Comparison between Stokes parameters, I, Q, and U (each row) for two setups of data reduction (the one used in this work is shown in the left panel, and the other one in the middle panel; see Sec. 2.1). The right panels show the difference between the left and middle panels (computed for each Stokes parameter).
In the text

	Fig. A.3 Same as in Fig. A.2, but for the polarization angles (see Eq. 6 in Sec. 2.1).
In the text

	Fig. A.4 SOFIA HAWC+ 214 μm (Band E) dust continuum map at 18.2″ angular resolution corresponding to linear scales of ~0.037 pc. The map is masked the same as the Band D map (Fig. 2, see Sec. 2.1). Black lines in both panels show the orientation of the magnetic field for every fifth pixel.
In the text

	Fig. A.5 Difference in polarization angles measured from Band D and E data is shown in the left panel. The right panel shows the rms of the angle differences.
In the text

Fig. A.6 SOFIA HAWC+ dust continuum Band D and Band E measurements (left and right panels in the top row, respectively), but with overlayed positions of protostars and NIR polarization (dark purple circles and lines, Kandori et al. 2007). Dark blue contours represent the outflow observed in the HCO+ emission. Light blue points indicate the positions of FIR sources. Black dashed rectangles show the central area and the zoom-in panels on the bottom row. Here, we show magnetic field lines from the FIR dust polarization at 100 μm (Dotson et al. 2000) in dark magenta and from (sub)millimeter dust polarization at 850 μm (Matthews et al. 2002) in black.

In the text

	Fig. B.1 Fitting parameters for the hyperfine CN emission in NGC 2024 measured for both velocity components (top and bottom row each). The first parameter, p1 is shown in the left panel, followed by the centroid velocity, FWHM and the final parameter, the opacity (p4).
In the text

	Fig. B.2 Fitting parameters for the HCO+ emission assuming a Gaussian line profile. Each row shows parameters for each fitted component. The results in the first row represent the brightest HCO+ component. The second row shows the fainter component and the outflow observed close to the center of NGC 2024.
In the text

Fig. C.1 Parameter space used for the radiative transfer modeling of CN in the region at the western side of the bubble. The first row shows results for the excitation temperature. The color bar shows values of excitation temperature, Tex for the grid of volume densities (x-axis), line widths (y-axis) and column densities (each panel) for fixed electron fraction. The second row shows results for opacity, τ, and the third row shows the peak temperature, Tpeak. The last row shows the χ2 minimization of the modeled peak temperature and the observed value computed from the CN(1−0) spectra.

In the text

	Fig. C.2 Same as in Fig. C.1, but for the eastern part of the bubble.
In the text

	Fig. C.3 Same as in Fig. C.1, but for the filament.
In the text

	Fig. C.4 Same as in Fig. C.1, but for the overlap region.
In the text

	Fig. C.5 Same as in Fig. C.1, but for HCO+ (1−0).
In the text

	Fig. C.6 Same as in Fig. C.5, but for the eastern side of the bubble.
In the text

	Fig. C.7 Same as in Fig. C.5, but for the filament.
In the text

	Fig. C.8 Same as in Fig. C.5, but for the overlap region.
In the text

	Fig. C.9 Beam averaged spectrum of CN and HCO+ taken within four different regions shown in the left panel of Fig. 1 and in Fig. 2. Colored line represent the corresponding model spectrum infrerred from the non-LTE RADEX modeling for physical conditions presented in Sec. C and shown in Table 3. Dashed lines show the corresponding residuals.
In the text

	Fig. D.1 Angle rms as a function of the size of the sliding window.
In the text

	Fig. D.2 Top left: Bubble on the west. Top right: Bubble on the eastern side of NGC 2024. Bottom left: Filament. Bottom right: Mixture of bubble and dusty filament.
In the text

	Fig. D.3 Comparison of the angle dispersion using two different approaches, sliding window (y-axis, see Sec. 5) and histogram analysis (x-axis).
In the text

	Fig. D.4 Magnetic field strength derived using different methods (Sec. 5.1, D.3) for computing the magnetic field strength and the angle dispersion (Sec. 5.3, D.2) as a function of the environment we analyze in this work.
In the text