11institutetext: LERMA, Observatoire de Paris, PSL Research University, CNRS,
Sorbonne Universités, 75014 Paris, France 22institutetext: Department of
Earth, Environment, and Physics, Worcester State University, Worcester, MA
01602, USA 33institutetext: Center for Astrophysics $|$ Harvard & Smithsonian,
60 Garden Street, Cambridge, MA 02138, USA 44institutetext: Jet Propulsion
Laboratory, California Institute of Technology, 4800 Oak Grove Drive,
Pasadena, CA 91109, USA 55institutetext: IRAM, 300 rue de la Piscine, 38406
Saint Martin d’Hères, France 66institutetext: Université de Toulon, Aix
Marseille Univ, CNRS, IM2NP, Toulon, France 77institutetext: IRAP, CNRS,
Université de Toulouse, 9 avenue du Colonel Roche, 31028 Toulouse Cedex 4,
France 88institutetext: Université Grenoble Alpes, CNRS, GIPSA-lab, Grenoble
INP, Grenoble, 38000, France 99institutetext: Instituto de Física Fundamental
(CSIC). Calle Serrano 121-123, 28006, Madrid, Spain 1010institutetext:
Laboratoire de Physique de l’ENS, ENS, Université PSL, CNRS, Sorbonne
Université, Université Paris Cité, Observatoire de Paris, F-75005 Paris,
France 1111institutetext: Cardiff Hub for Astrophysics Research $\&$
Technology, School of Physics $\&$ Astronomy, Cardiff University, Queens
Buildings, The parade, Cardiff CF24 3AA, UK 1212institutetext: Université
Paris-Saclay, CNRS, Institut d’Astrophysique Spatiale, 91405, Orsay, France

# The magnetic field in the Flame nebula

I. Bešlić The magnetic field in the Flame nebulaThe magnetic field in the
Flame nebula 0000-0003-0583-7363 S. Coudé The magnetic field in the Flame
nebulaThe magnetic field in the Flame nebulaThe magnetic field in the Flame
nebulaThe magnetic field in the Flame nebula 0000-0002-0859-0805 D. C. Lis The
magnetic field in the Flame nebulaThe magnetic field in the Flame nebula
0000-0002-0500-4700 M. Gerin The magnetic field in the Flame nebulaThe
magnetic field in the Flame nebula 0000-0002-2418-7952 P. F. Goldsmith The
magnetic field in the Flame nebulaThe magnetic field in the Flame nebula
0000-0002-6622-8396 J. Pety The magnetic field in the Flame nebulaThe magnetic
field in the Flame nebulaThe magnetic field in the Flame nebulaThe magnetic
field in the Flame nebula 0000-0003-3061-6546 A. Roueff The magnetic field in
the Flame nebulaThe magnetic field in the Flame nebula K. Demyk The magnetic
field in the Flame nebulaThe magnetic field in the Flame nebula
0000-0002-5019-8700 C. D. Dowell The magnetic field in the Flame nebulaThe
magnetic field in the Flame nebula L. Einig The magnetic field in the Flame
nebulaThe magnetic field in the Flame nebulaThe magnetic field in the Flame
nebulaThe magnetic field in the Flame nebula 0000-0003-4250-7638 J. R.
Goicoechea The magnetic field in the Flame nebulaThe magnetic field in the
Flame nebula 0000-0001-7046-4319 F. Levrier The magnetic field in the Flame
nebulaThe magnetic field in the Flame nebula 0000-0002-3065-9944 J. Orkisz The
magnetic field in the Flame nebulaThe magnetic field in the Flame nebula
0000-0002-3382-9208 N. Peretto The magnetic field in the Flame nebulaThe
magnetic field in the Flame nebula 0000-0002-6893-602X M. G. Santa-Maria The
magnetic field in the Flame nebulaThe magnetic field in the Flame nebula
0000-0002-3941-0360 N. Ysard The magnetic field in the Flame nebulaThe
magnetic field in the Flame nebulaThe magnetic field in the Flame nebulaThe
magnetic field in the Flame nebula 0000-0003-1037-4121 A. Zakardjian The
magnetic field in the Flame nebulaThe magnetic field in the Flame nebula
0000-0002-4240-6012

(Received 24 October 2023; accepted 20 January 2024)

###### Abstract

Context. Star formation drives the evolution of galaxies and the cycling of
matter between different phases of the interstellar medium and stars. The
support of interstellar clouds against gravitational collapse by magnetic
fields has been proposed as a possible explanation for the low observed star
formation efficiency in galaxies and the Milky Way. The Planck satellite
provided the first all-sky map of the magnetic field geometry in the diffuse
interstellar medium on angular scales of 5-15′. However, higher spatial
resolution observations are required to understand the transition from
diffuse, subcritical gas to dense, gravitationally unstable filaments.

Aims. NGC 2024, also known as the Flame Nebula, is located in the nearby Orion
B molecular cloud. It contains a young, expanding HII region and a dense
supercritical filament. This filament harbors embedded protostellar objects
and is likely not supported by the magnetic field against gravitational
collapse. Therefore, NGC 2024 provides an excellent opportunity to study the
role of magnetic fields in the formation, evolution, and collapse of dense
filaments, the dynamics of young HII regions, and the effects of mechanical
and radiative feedback from massive stars on the surrounding molecular gas.

Methods. We combine new 154 and 216 µm dust polarization measurements carried
out using the HAWC+ instrument aboard SOFIA with molecular line observations
of 12CN(1–0) and HCO+(1–0) from the IRAM 30-meter telescope to determine the
magnetic field geometry, and to estimate the plane of the sky magnetic field
strength across the NGC 2024 HII region and the surrounding molecular cloud.

Results. The HAWC+ observations show an ordered magnetic field geometry in NGC
2024 that follows the morphology of the expanding HII region and the direction
of the main dense filament. The derived plane of the sky magnetic field
strength is moderate, ranging from 30 to 80 µG. The strongest magnetic field
is found at the northern west edge of the HII region, characterized by the
highest gas densities and molecular line widths. In contrast, the weakest
field is found toward the main, dense filament in NGC 2024.

Conclusions. We find that the magnetic field has a non negligible influence on
the gas stability at the edges of the expanding HII shell (gas impacted by the
stellar feedback) and the filament (site of the current star formation).

###### Key Words.:

ISM: magnetic fields – ISM: molecules – Galaxies: star formation

##  1 Introduction

The question of what controls 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⊙/year (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
supernovae explosions (Dobbs et al., 2011; Walch & Naab, 2015), or due to
radiative and mechanical stellar feedback from high-mass stars 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 the efforts over the past 40 years 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 will have 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, e.g., Pattle et al., 2019, for a recent
review).

Figure 1: (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 IRAM 30-meter
telescope (image credits: Pety et al., 2017). We overlay 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
sketch the environments seen across NGC 2024: the HII region (white dashed
circle, see Tab. 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 label
the edges of the expanding HII region and the direction of the stellar
feedback driven by the radiation produced by recent star formation.

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 is
firstly proposed by Davis & Greenstein (1951) as a paramagnetic alignment with
the magnetic field. Later on, this process is described by the Radiative
Alignment Torques theory (RAT; 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 will be 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 measuring polarized dust emission 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 plane-of-the-sky
orientation of the $B$-field and the filamentary structures is correlated with
the observed column density, with the field largely parallel to diffuse
structures (i.e., striations) while perpendicular to dense filaments (Planck
Collaboration et al., 2016a, b). 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 (e.g.,
HII bubble shells, Chapman et al., 2011).

BLASTPol has also observed similar $B$-field morphology at 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^{\prime\prime}$) ground-based observations using the POL-2 instrument on
the James Clerk Maxwell Telescope (JCMT; e.g., 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 gravitationally-dominated gas
requires a redistribution of magnetic flux (e.g., 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. Indeed, a recent study
(Tang et al., 2019) has investigated the relation between dense gas velocity
gradients (as imaged in N2H+) 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 towards the center of the
massive ridge. Extending such studies to larger fields is of crucial
importance and the High-resolution Airborne Wide-band 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, e.g., 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.

For 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 method (DCF; 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 0 and 1
(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
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 magnetic field strength from dust polarization
measurements (Skalidis-Tassis method; hereafter ST). 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). In this work, we use 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 (e.g., 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 (105 cm-3) and very
active in terms of star formation, affected by strong shocks (including an
explosive outflow resulting from the merger of 3 massive protostars, see,
e.g., Bally, 2008). These shocks exhibit exceptionally intense UV
illumination, with enhancement factors of $G_{0}$ $\sim$104, up to 105
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 is overall more representative of a
standard GMC in our galaxy and normal galaxies, with FUV radiation field
$G_{0}$ in the range of $\sim 10^{3}-10^{4}$ (Santa-Maria et al., 2023).

In a recent study, Orkisz et al. (2017) used 12CO(1–0) and 13CO(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 C18O(1–0). Most
identified filaments in Orion B are low-density, thermally sub-critical
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 ($\zeta$ Ori) in the Orion B
complex (e.g., Meyer et al., 2008). This region contains a massive and young
HII region (age $2\cdot 10^{5}$ yr, Tremblin et al., 2014, and private
communication), deeply embedded in dust, located in the foreground, and
extending to the south (e.g., 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 HII region, as well as embedded protostellar
objects (Mezger et al., 1988, 1992; Lis et al., 1991) located in the
background (light blue points on the right panel in Fig. 1). The young HII
region is expanding, strongly impacting its parental cloud, and creating sharp
ionization fronts towards the south, which makes it a good example
demonstrating how such systems can efficiently exert stellar feedback. The
edges of the bubble are seen towards the west and east of the center of NGC
2024 (labeled as Eastern and Western Loop in Fig. 1).

##  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_0154 (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 (2021/09/08), F780 (2021/09/09), F782
(2021/09/11), F783 (2021/09/14), and F784 (2021/09/15). The map obtained at
each wavelength comprises ten $20^{\prime}\times 7^{\prime}$ 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 & 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^{\prime\prime}$ , respectively, in Band D, and
$4.6^{\prime\prime}$ and $18.7^{\prime\prime}$ in Band E. The final data
products for both Band D and E contain the $I$, $Q$ and $U$ Stokes parameters,
the polarization fraction $P$ and angle $\theta$, and their uncertainties
(Gordon et al., 2018). Clarke & 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}^{\prime}$ is calculated using the following equation:

$I_{p}^{\prime}=\sqrt{Q^{2}+U^{2}}.$ (1)

The uncertainty $\delta I_{p}^{\prime}$ on the polarized intensity
$I_{p}^{\prime}$ is then

$\delta I_{p}^{\prime}=\dfrac{\sqrt{(Q\,\delta Q)^{2}+(U\,\delta
U)^{2}}}{I_{p}^{\prime}},$ (2)

where $\delta Q$ and $\delta U$ are the uncertainties on Stokes $Q$ and $U$,
respectively.

The polarized intensity $I_{p}^{\prime}$ 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, e.g.,
Montier et al., 2015), used in this work, is calculated using

$I_{p}=\sqrt{I_{p}^{\prime 2}-\delta I_{p}^{\prime 2}},$ (3)

with uncertainty $\delta I_{p}=\delta I_{p}^{\prime}$.

The polarization fraction $P$ is then obtained from

$P=100\,\frac{I_{p}}{I},$ (4)

with the uncertainty $\delta P$ given by

$\delta P=P\,\sqrt{\left(\frac{\delta
I_{p}^{\prime}}{I_{p}^{\prime}}\right)^{2}+\left(\frac{\delta
I}{I}\right)^{2}}.$ (5)

The polarization angle $\theta_{p}$ is defined using the Stokes parameters $Q$
and $U$ as

$\theta_{p}=\dfrac{1}{2}\,\arctan{\dfrac{U}{Q}},$ (6)

with the uncertainty $\delta\theta_{p}$ given by

$\delta\theta_{p}=\dfrac{1}{2}\,\dfrac{\sqrt{(Q\,\delta U)^{2}+(U\,\delta
Q)^{2}}}{I_{p}^{\prime 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
$\pi/2$ to Eq. 6.

Prior to the data analysis, we filter out low signal-to-noise (S/N) data
points. After masking the data, we keep pixels that show S/N ratio $\geq 50$
in total intensity, $\leq 30^{o}$ in polarization angle uncertainty, and $\leq
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 App. A.1 in Figs. 7 and 8. We find an
overall good agreement between these two data sets.

###  2.2 IRAM 30-meter observations of Orion B

Table 1: Components of the hyperfine structure of the 12CN(1–0) transition, their relative offsets and intensities relative to the sum of all components (Milam et al., 2009). Transition F | Rest frequency [GHz] | Relative offset [km/s] | Relative intensity
---|---|---|---
3/2 $-$ 1/2 | 113.488142 | 7.51 | 0.1235
5/2 $-$ 3/2 | 113.490985 | 0 | 0.3333
1/2 $-$ 1/2 | 113.499643 | -22.88 | 0.0988
3/2 $-$ 3/2 | 113.508934 | -47.45 | 0.0988
1/2 $-$ 3/2 | 113.520414 | -77.8 | 0.0123

In our work, we make use of information of the 12CN(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 ($\sim$20 pc across) in the Orion B molecular
cloud, at an angular resolution of 26′′ ($10^{4}$ AU, $0.05$ pc) in at least
30 molecular lines in the (71–79) and (84–116) GHz range with a spectral
resolution $\sim$0.6 km s-1. These include common molecular tracers such as
12CO(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 ($\sim$820,000
pixels, $\sim$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 GILDAS111See
http://www.iram.fr/IRAMFR/GILDAS for more information about the GILDAS
software (Pety, 2005)./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^{\prime\prime}$ and $23^{\prime\prime}$ for
the HCO+ and 12CN 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 Tab. 2 in Gratier et al., 2017), respectively,
over the studied field of view.

We use spectroscopic data of cyanide (12CN) 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 Tab. 1, which allow derivation
of the opacity (Milam et al., 2009) and excitation temperature of the (1–0)
transition, enabling accurately constraining its excitation and obtaining
limits on the gas density.

HCO+ , similarly to CN, traces photon-dominated regions (Young Owl et al.,
2000a; 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., 2000b; van der Werf et al., 1996)
or CH+ molecule (e.g., Goicoechea et al., 2016). Recent studies detected HCO+
emission the edge of the PDR in the Orion bar (Goicoechea et al., 2016), and
the Horsehead 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 C18O
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.

Figure 2: SOFIA HAWC+ $154\,$µm (Band D, all panels) dust continuum maps at
$13.6^{\prime\prime}$ angular resolution. The maps are masked based on the
measured S/N ratio 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.

##  3 Dust polarization

In this section, we discuss our SOFIA HAWC+ dust polarization measurements
introduced in Sec. 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 $\sim 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. 9 in App. A.2. 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. 10 in App. A.2. While the rms is only a few
degrees ($\sim 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 towards 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 (i.e., there is no $\pi/2$ difference as for the FIR dust polarization -
Sec. 2.1). We show this comparison in the top panels of Fig. 11 in App. A.2.
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. 11 in App. A.2). In addition,
we compare the magnetic field direction with the 100 µm (Dotson et al., 2000)
and sub-millimeter (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. 11) and show zoomed-in panels in the bottom row of Fig. 11.

The area labeled with dark blue contours in the top panels in Fig. 11 is
particularly interesting because we observe an outflowing feature in the HCO+
emission at velocities of $\sim 14$ km s-1. The outflow originates from the
FIR 5 source located in the dense molecular cloud behind the expanding HII
region (e.g., 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, e.g., 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 HII region (middle panels in Fig. 2). Gas affected by the
stellar feedback (e.g., UV radiation, stellar winds) is pushed outward from
the HII region. Consequently, magnetic field lines become parallel to the
edges of such an expanding shell (e.g., 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 HII region and the dense molecular gas
in the background (see, e.g., 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 ($\Delta N_{\mathrm{H_{2}}}>10^{20}$ cm-2, e.g.,
Planck Collaboration et al., 2016c; 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.

##  4 Characterizing the turbulence and density structures in NGC 2024

The DCF method (Sec. 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
12CN(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).

Figure 3: Moment maps of 12CN(1–0) and HCO+(1–0) emission across NGC 2024.
Top: Maps of the CN emission: integrated intensity (i.e., 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 Sec. 4.1.

###  4.1 Moment maps

Fig. 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 ratio; we have selected all voxels with S/N ratio 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 Sec. 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. 12 and 13 in App. 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 towards 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 13CO(1–0) and C18O(1–0) emission
(Gaudel et al., 2023). We discuss velocity components in App. B.1 and B.2.

The second moment map (the observed line width ($\sigma_{\mathrm{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,
e.g., Barnes & Crutcher, 1990), which additionally broadens the line. Third,
it is possible that HCO+, similarly to 13CO and C18O 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. 12). 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. 13) 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 12CN and HCO+
line emission profiles. For instance, the 12CN(1–0) has a hyperfine structure
(see Fig. 22 in App. B.1). 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 12CN(1–0) emission in App. B.1. The profile of the
HCO+ can be described by a Gaussian function (see Fig. 22 in App. B.2)
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_{\mathrm{NT,mol}}=\sqrt{\sigma_{\mathrm{obs,mol}}^{2}-\sigma_{\mathrm{TH,mol}}^{2}},$
(8)

where $\sigma_{\mathrm{obs,mol}}$ is the measured FWHM of a molecular line,
and $\sigma_{\mathrm{TH,mol}}$ is the line width of the thermal component. We
calculate the thermal broadening using the following equation:

$\sigma_{\mathrm{TH,mol}}=\sqrt{\dfrac{kT_{\mathrm{k}}}{m_{\mathrm{mol}}}},$
(9)

where $k$ is Boltzman’s constant, $T_{\mathrm{k}}$ is the kinetic temperature
obtained using the CO(1–0) measurements (Orkisz et al., 2017), and
$m_{\mathrm{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:

$\Delta\upsilon=\dfrac{\sigma_{\mathrm{NT,mol}}}{\beta},$ (10)

where the factor $\beta$ is a function of the optical depth ($\tau_{0}$) at
the line center, defined as (Phillips et al., 1979; Hacar et al., 2016; Orkisz
et al., 2017):

$\beta=\dfrac{1}{\sqrt{\ln{2}}}\cdot\left[\ln{\left\\{\dfrac{\tau_{0}}{\ln{\left(\dfrac{2}{\exp{(-\tau_{0})}+1}\right)}}\right\\}}\right]^{\dfrac{1}{2}}.$
(11)

The equations presented above are crucial at sub-pc 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, Sec. 4.1), or from the line fitting by
using, for example, the spectral decomposition (see Apps. 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, 2021) 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_{\mathrm{H_{2}}}$ and line widths across four
regions in NGC 2024: two edges of the expanding HII 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_{\mathrm{kin}}$) derived from the 12CO(1–0)
peak intensity (Orkisz et al., 2017), and shown in Tab. 2. We create a grid of
column densities in the range of (1-4)$\cdot 10^{14}$ cm-2 for CN and
(1-4)$\cdot 10^{13}$ cm-2 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^{2}$ to
$10^{5}$ cm-3. The density of H2 is comprised of para- and ortho-H2, assuming
the temperature dependence of the ortho-to-para ratio: $9\cdot
e^{-170/T_{\mathrm{kin}}}$ (Mandy & Martin, 1993). Here, we assume a fixed
electron fraction (electron-to-H2 density ratios), $f_{\mathrm{e}}=1.4\cdot
10^{-4}$, 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 $f_{\mathrm{e}}$ 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 H2 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-H2
(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-H2 (Denis-Alpizar et al., 2020) and
electrons (Fuente et al., 2008) as well.

The output parameters from the RADEX modelling are the excitation temperature,
$T_{\mathrm{ex}}$ and opacity, $\tau$. 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 (Tab. 1).

We compute model CN(1–0) and HCO+(1–0) spectra assuming the optical depth to
be a Gaussian function of frequency, $\tau(\nu)$:

$T_{mb}(\nu)=(T_{\mathrm{ex}}-T_{\mathrm{bg}})\cdot(1-e^{-\tau(\nu)}).$ (12)

The best-fit model is found based on the modeled spectrum for each combination
of $(n_{\mathrm{H_{2}}},N_{\mathrm{mol}},\Delta\upsilon)$. We perform a
$\chi^{2}$ minimization for the modeled and observed peak temperature,
$T_{\mathrm{peak}}$. Then, we calculate the turbulent line width following the
prescription in Sec. 4.2. We report our turbulent line width and gas number
density results in Tab. 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.
22 in App. C.

Table 2: Derived physical conditions of molecular gas traced by CN and HCO+
emission for four regions we investigate in NGC 2024, which are shown in Fig.
2. We tabulate position of each region, and kinetic temperature inferred from
the CO data. Next, we tabulate results from the non-LTE radiative transfer
modelling using RADEX: volume number density, and measured turbulent line
widths.

Region | RA [deg] | Dec [deg] |  $T_{\mathrm{kin}}$ [K] |  $\Delta\upsilon_{\mathrm{CN}}$ [km/s] |  $n_{\mathrm{H_{2}}}\mathrm{(CN)}$ [$10^{3}$cm-3] |  $\Delta\upsilon_{\mathrm{HCO^{+}}}$ [km/s] |  $n_{\mathrm{H_{2}}}\mathrm{(HCO^{+})}$ [$10^{3}$cm-3]
---|---|---|---|---|---|---|---
Bubble West | 85.4018 | -1.9642 | 53.04 | 0.56 | 2.27 | 1.10 | 2.93
Bubble East | 85.4768 | -1.9449 | 32.63 | 1.24 | 4.87 | 1.24 | 4.93
Filament | 85.4658 | -1.9864 | 27.18 | 0.62 | 4.72 | 0.96 | 2.31
Filament+Bubble | 85.3490 | -1.7587 | 30.80 | 0.94 | 13.37 | 1.07 | 5.39

Figure 4: The direction of 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\times 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 grey contour indicates the mask we defined
for dust polarization measurements, described in Sec. 2.1.

##  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_{\mathrm{pos}}=\sqrt{4\pi\rho}\dfrac{\sigma_{\mathrm{NT}}}{\widehat{\sigma}_{\mathrm{c}}(\varphi)},$
(13)

where $\rho$ is the gas volume density, $\sigma_{\mathrm{NT}}$ is the non-
thermal velocity dispersion
($\sigma_{\mathrm{NT}}=\Delta\upsilon/\sqrt{8\ln{2}}$), where $\Delta\upsilon$
is the measured FWHM of the line, corrected for the line broadening and
opacity effects - Eq. 8 and 10 \- Sec. 4.2, and
$\widehat{\sigma}_{\mathrm{c}}(\varphi)$ 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-3, km
s-1, and deg, respectively (e.g., Lyo et al., 2021), Eq. 13 can be expressed
as:

$B_{\mathrm{pos}}\approx
9.3\cdot\sqrt{n_{\mathrm{H_{2}}}}\cdot\dfrac{\Delta\upsilon}{\widehat{\sigma}_{\mathrm{c}}(\varphi)}\,[\rm\upmu
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, e.g.,
Ostriker et al., 2001).

###  5.2 Skalidis-Tassis method

The DCF method is widely used in the literature to derive the strength of the
magnetic field (see, e.g., 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_{\mathrm{POS}}$ using the prescription
presented in Skalidis & Tassis (2021); Skalidis et al. (2021):

$B_{\mathrm{pos}}=\sqrt{2\pi\rho}\dfrac{\sigma_{\mathrm{NT}}}{\sqrt{\widehat{\sigma}_{\mathrm{c}}(\varphi)}},$
(15)

where $\rho$, $\sigma_{\mathrm{NT}}$, and
$\widehat{\sigma}_{\mathrm{c}}(\varphi)$ are the same as in Eq. 13.

Similarly as in Eq. 14, after substituting all constants, the above equation
becomes:

$B_{\mathrm{pos}}\approx
1.8\cdot\sqrt{n_{\mathrm{H_{2}}}}\cdot\dfrac{\Delta\upsilon}{\sqrt{\widehat{\sigma}_{\mathrm{c}}(\varphi)}}\,[\rm\upmu
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 (e.g., ”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\times 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\times 12=144$ directions $\varphi(l)$ of each
window is to compute:

$z=\dfrac{1}{L}\sum_{l=1}^{L}e^{2i\varphi(l)},$ (17)

where $L$ is the number of pixels within the sliding window. Defining $a$ and
$b$ such that $z=a+i\cdot b$, one computes the circular mean using

$\widehat{m}_{\mathrm{c}}(\varphi)=\frac{1}{2}\arctan\dfrac{b}{a},$ (18)

and the circular standard deviation using

$\widehat{\sigma}_{\mathrm{c}}(\varphi)=\sqrt{\dfrac{1}{2}\left(1-\left|z\right|\right)}.$
(19)

Moreover, to consider the uncertainty associated with the different directions
$\varphi(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:

$\Delta(\widehat{\sigma}_{\mathrm{c}}(\varphi))=\dfrac{\widehat{\sigma}_{\mathrm{c}}(\varphi)}{\sqrt{2\cdot
L-2}}.$ (20)

We show maps of the magnetic field angle, the mean angle computed using Eq. 18
and $12\times 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
HII region and the filament. In particular, the measured
$\widehat{\sigma}_{\mathrm{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.

Figure 5: The ratio between measured line widths (right panels in Fig. 3) and
the square root of 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 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.

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 Secs. 4.1 and 4.2).

By assuming a gas number density of, for instance, $10^{4}$ cm-3, using Eq.
14, the ratio of 0.5 km s-1 deg-1 in Fig. 5 corresponds to $90$ µG, whereas a
ratio of 2 km s-1 deg-1 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.

###  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 HII region (Filament+Bubble) (Fig. 2, also see
Sec. 4). We measure the uncertainty of magnetic field strength derived using
the DCF method (Sec. 5.1) as follows:

$\Delta
B_{\mathrm{POS}}=B_{\mathrm{POS}}\cdot\sqrt{\dfrac{1}{4}\cdot\dfrac{(\Delta
n)^{2}}{n^{2}}+\dfrac{(\Delta(\Delta\upsilon))^{2}}{(\Delta\upsilon)^{2}}+\dfrac{(\Delta(\widehat{\sigma}_{\mathrm{c}}(\varphi)))^{2}}{\widehat{\sigma}_{\mathrm{c}}(\varphi)^{2}}},$
(21)

where $\Delta n$, $\Delta(\Delta\upsilon)$, and
$\Delta(\widehat{\sigma}_{\mathrm{c}}(\varphi))$ are measured uncertainties of
$n$, $\Delta\upsilon$ and $\widehat{\sigma}_{\mathrm{c}}(\varphi)$ (Sec. 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$, $\sigma_{\mathrm{NT}}$ from the
$\chi^{2}$ minimization (see the bottom panel in Fig. 14-21). We select the
range in parameter space within which the $\chi^{2}$ reaches 0.5 to derive
corresponding $n$, and $\sigma_{\mathrm{NT}}$. Thus, we found $\Delta n$ to be
$10^{3}$ cm-3 and $\Delta\sigma_{\mathrm{NT}}=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 Tab. 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. 26 in App. D.3. 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 Sec. 5.2. The $B_{\mathrm{POS}}$ strength is in the range of $\sim 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 (e.g. 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 ($\sim 10^{-3}$ 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
HII 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 Sec. C) found in this work. The magnetic field strength measured
using both CN and HCO+ emission is stronger towards 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 towards
the filament, although we note the $B_{\mathrm{POS}}$ derived from HCO+ in the
western side of the bubble and filament to be comparable given the
uncertainties.

Table 3: Measured angle dispersion from HAWC+ data and results on measured
magnetic field strength and its uncertainty, using results from modeling the
CN and HCO+ excitation (App. C, Tab. 2). The POS magnetic field strength is
derived using the classical DCF method (third and fourth column) and the ST
method presented in Skalidis & Tassis (2021); Skalidis et al. (2021) (fifth
and sixth column).

Region |  $\sigma_{\theta}$ [deg] |  $B_{\mathrm{DCF,CN}}$ [µG] |  $B_{\mathrm{DCF,HCO^{+}}}$ [µG] |  $B_{\mathrm{ST,CN}}$ [µG] |  $B_{\mathrm{ST,HCO^{+}}}$ [µG]
---|---|---|---|---|---
Bubble West | 5.5 $\pm$ 0.3 | 46 $\pm$ 23 | 90 $\pm$ 26 | 20 $\pm$ 10 | 40 $\pm$ 11
Bubble East | 3.5 $\pm$ 0.2 | 235 $\pm$ 55 | 235 $\pm$ 55 | 82 $\pm$ 19 | 82 $\pm$ 19
Filament | 6.1 $\pm$ 0.4 | 66 $\pm$ 28 | 102 $\pm$ 35 | 30 $\pm$ 13 | 47 $\pm$ 16
Filament+Bubble | 6.3 $\pm$ 0.4 | 162 $\pm$ 45 | 185 $\pm$ 48 | 76 $\pm$ 21 | 87 $\pm$ 22

##  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 Sec. 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 HI 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 modelled 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. 11 in App.
A.3. 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 HII 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 (e.g.,
Planck Collaboration et al., 2016a; 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 HII region to be $\sim 20-40$ µG on the west side, and $\sim 82$ µG
on the east side. The good agreement between the CN and HCO+
$B_{\mathrm{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_{\mathrm{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 e.g. 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 ($\tau>4$) on the western side than on the eastern
side ($\tau=2$).

As previously presented in Sec. 5, the eastern edge of HII 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 (Tab. 2),
stellar feedback did not sweep up gas as far as on the other side (Barnes et
al., 1989).

The values of $B_{\mathrm{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 HI and OH in M 17 is $\sim
750\,$µG and $\sim 250\,$µG respectively (Brogan & Troland, 2001), and
$\sim(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 C18O 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
$\times$ $10^{3}$ cm-3), measured line widths are different. It is worth
pointing out that C18O and CN and HCO+ are tracing different gas (Philipp et
al., 2006; Bron et al., 2018). C18O 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\cdot
10^{4}$ cm-3 (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
Tab. 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_{\mathrm{POS}}$. We show how $B_{\mathrm{POS}}$ varies with the method we
selected in Fig. 26.

####  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 ($\sim
30-50$ µG) in this region, suggesting that the magnetic field can not dominate
gas dynamics. We return to this point in Sec. 6.3. These low values in the
magnetic field strength come from the largest dispersion angles and narrow
spectral lines computed for this region (Tab. 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 C18O emission to measure the turbulence and
number density and different method to derive the $B_{\mathrm{POS}}$.
Similarly, Ching et al. (2018) investigated magnetic field towards the DRL1
filament and found $B_{\mathrm{POS}}=600\,$µG.

###  6.3 Is the magnetic field sufficient to prevent star formation in NGC
2024?

Table 4: Estimated mass-to-flux ratio $\mu_{\Phi}$, Alfvénic Mach number
$\mathcal{(}M)_{\mathrm{A}}$ and plasma-beta $\beta$ using $B_{\mathrm{POS}}$
reported in Tab. 3. These parameters indicate role of radiation, gravitational
and magnetic field and their impact on the gas in four regions in NGC 2024.

Region | $\mu_{\Phi}(\mathrm{CN})$ | $\mu_{\Phi}(\mathrm{HCO^{+}})$ | $\mathcal{M}_{\mathrm{A}}(\mathrm{CN})$ | $\mathcal{M}_{\mathrm{A}}(\mathrm{HCO^{+}})$ | $\beta_{B}(\mathrm{CN})$ | $\beta_{B}(\mathrm{HCO^{+}})$
---|---|---|---|---|---|---
Bubble West | 0.38 | 0.19 | 1.78 | 2.03 | 1.051 | 0.345
Bubble East | 0.09 | 0.09 | 1.42 | 1.43 | 0.083 | 0.084
Filament | 2.50 | 1.61 | 1.88 | 1.31 | 0.483 | 0.098
Filament+Bubble | 1.00 | 0.88 | 1.91 | 1.21 | 0.247 | 0.077

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^{-n}$.

We make an estimate on the mass-to-flux ratio of gas in NGC 2024
$\left(\frac{M}{\Phi}\right)$ and compare it to the critical mass-to-flux
ratio. The critical mass-to-flux ratio
$\left(\frac{M}{\Phi}\right)_{\mathrm{crit}}$ depends on the assumed geometry.
For instance, for a uniform disk,
$\left(\frac{M}{\Phi}\right)_{\mathrm{crit,cyl}}=\dfrac{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}{\Phi}\right)_{\mathrm{crit,sph}}=\dfrac{1}{3\pi\sqrt{\frac{G}{5}}}$
(Mouschovias & Spitzer, 1976). The $\mu_{\Phi}$ parameter is the ratio between
estimated mass-to-flux and $\left(\frac{M}{\Phi}\right)_{\mathrm{crit}}$ and
can be computed as:

$\left(\mu_{\Phi}\right)_{\mathrm{cyl}}=\dfrac{\left(M/\Phi\right)}{\left(M/\Phi\right)_{\mathrm{crit}}}=7.6\cdot
10^{-21}\dfrac{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-2, 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}{\Phi}\right)_{\mathrm{crit}}$ and the
above equation:

$\left(\mu_{\Phi}\right)_{\mathrm{sph}}=0.67\cdot\left(\mu_{\Phi}\right)_{\mathrm{cyl}}.$
(23)

Therefore, the $\mu_{\Phi}$ 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
threated with causion.

Next, we compute the Alfvénic Mach number, $\mathcal{M}_{\mathrm{A}}$:

$\mathcal{M}_{\mathrm{A}}=\dfrac{\sqrt{3}\cdot\sigma_{\upsilon,\mathrm{NT}}}{\upsilon_{\mathrm{A}}}.$
(24)

where $\sigma_{\upsilon,NT}$ is the non-thermal velocity dispersion reported
in Tab. 2. The $\sqrt{3}$ factor comes from the assumption of isotropic
turbulence (Crutcher, 2004; Stewart & Federrath, 2022).
$\upsilon_{\mathrm{A}}$ is the Alfvén velocity defined as

$\upsilon_{\mathrm{A}}=\dfrac{B}{\sqrt{4\pi\rho}},$ (25)

where $B$ is the total magnetic field strength in units of G, and $\rho$ is
the gas mass volume density in units of g cm-3. Eq. 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, $\beta_{B}$, that gives
information about the ratio of thermal and magnetic pressure:

$\beta_{B}=\dfrac{nk_{\mathrm{B}}T}{B^{2}/(8\pi)},$ (26)

where $n$ is the number density, $T$ is the gas temperature (Tab. 2), and $B$
magnetic field strength (Tab. 3).

We report our results in Tab. 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 $\mu_{\Phi}<1$ in each region.
This result could imply the presence of magnetically supported gas (see e.g.
Pattle et al., 2019). In general, it is worth noting that we see a difference
of one and two orders of magnitude in the $\mu_{\Phi}$ between the edges of
the bubble and the filament and the overlap region.

We find $M_{\mathrm{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 $\mu_{\Phi}$, 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 Tab. 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 Tab. 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 gravitational
stability of the gas in these regions require 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 HII region is not always the case.
For example, the expansion of the HII region can also trigger star formation,
as seen in Galactic HII regions RCW 82 (Pomarès et al., 2009) and RCW 120
(Figueira et al., 2017). Since the HII 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 HII region.

The role of the magnetic field is somewhat different in the filament. High
$\mu_{\Phi}>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.
$\mathcal{M}_{\mathrm{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 Tab. 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
HII region and associated molecular cloud in the Orion B molecular cloud. We
combine 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. 1.

Our results focus on a subset of environments found in NGC 2024, particularly
on the shell near the edge of the HII 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. 2.

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

  3. 3.

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

  4. 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 HII
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. 5.

We derive 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 HII 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. 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. Magnetic field strength calculated from HCO+(1–0) is
generally higher than that inferred from CN(1–0). $B_{\mathrm{POS}}$ measured
from CN(1–0) and HCO+(1–0) is larger on the eastern side of the bubble than on
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. 7.

By analyzing 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 1 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 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 (emaa.osug.fr and EMAA reference). 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.

## References

  * Alina et al. (2019) Alina, D., Ristorcelli, I., Montier, L., et al. 2019, MNRAS, 485, 2825
  * Andersson et al. (2015) Andersson, B. G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501
  * Ballesteros-Paredes et al. (2011) Ballesteros-Paredes, J., Hartmann, L. W., Vázquez-Semadeni, E., Heitsch, F., & Zamora-Avilés, M. A. 2011, MNRAS, 411, 65
  * Bally (2008) Bally, J. 2008, in Handbook of Star Forming Regions, Volume I, ed. B. Reipurth, Vol. 4, 459
  * Barnes & Crutcher (1990) Barnes, P. J. & Crutcher, R. M. 1990, ApJ, 351, 176
  * Barnes et al. (1989) Barnes, P. J., Crutcher, R. M., Bieging, J. H., Storey, J. W. V., & Willner, S. P. 1989, ApJ, 342, 883
  * Basu (2000) Basu, S. 2000, ApJ, 540, L103
  * Battisti & Heyer (2014) Battisti, A. J. & Heyer, M. H. 2014, ApJ, 780, 173
  * Bik et al. (2003) Bik, A., Lenorzer, A., Kaper, L., et al. 2003, A&A, 404, 249
  * Bonnell et al. (2011) Bonnell, I. A., Smith, R. J., Clark, P. C., & Bate, M. R. 2011, MNRAS, 410, 2339
  * Brogan & Troland (2001) Brogan, C. L. & Troland, T. H. 2001, ApJ, 560, 821
  * Bron et al. (2018) Bron, E., Daudon, C., Pety, J., et al. 2018, A&A, 610, A12
  * Bron et al. (2021) Bron, E., Roueff, E., Gerin, M., et al. 2021, A&A, 645, A28
  * Cao et al. (2023) Cao, Z., Jiang, B., Zhao, H., & Sun, M. 2023, ApJ, 945, 132
  * Chandrasekhar & Fermi (1953) Chandrasekhar, S. & Fermi, E. 1953, ApJ, 118, 113
  * Chapman et al. (2011) Chapman, N. L., Goldsmith, P. F., Pineda, J. L., et al. 2011, ApJ, 741, 21
  * Chevance et al. (2022) Chevance, M., Kruijssen, J. M. D., Krumholz, M. R., et al. 2022, MNRAS, 509, 272
  * Ching et al. (2018) Ching, T.-C., Lai, S.-P., Zhang, Q., et al. 2018, ApJ, 865, 110
  * Ching et al. (2022) Ching, T. C., Li, D., Heiles, C., et al. 2022, Nature, 601, 49
  * Choi et al. (2015) Choi, M., Kang, M., & Lee, J.-E. 2015, AJ, 150, 29
  * Chuss et al. (2019) Chuss, D. T., Andersson, B. G., Bally, J., et al. 2019, ApJ, 872, 187
  * Clark & Hensley (2019) Clark, S. E. & Hensley, B. S. 2019, ApJ, 887, 136
  * Clarke & Vliet (2023) Clarke, M. & Vliet, R. V. 2023, SOFIA-USRA/sofia_redux: v1.3.2
  * Colín et al. (2013) Colín, P., Vázquez-Semadeni, E., & Gómez, G. C. 2013, MNRAS, 435, 1701
  * Crutcher et al. (2009) Crutcher, R., Roberts, D., Troland, T., & Goss, a. 2009, The Astrophysical Journal, 515, 275
  * Crutcher (2004) Crutcher, R. M. 2004, in The Magnetized Interstellar Medium, ed. B. Uyaniker, W. Reich, & R. Wielebinski, 123–132
  * Crutcher (2012) Crutcher, R. M. 2012, ARA&A, 50, 29
  * Crutcher et al. (1999) Crutcher, R. M., Roberts, D. A., Troland, T. H., & Goss, W. M. 1999, ApJ, 515, 275
  * Davis & Greenstein (1951) Davis, Leverett, J. & Greenstein, J. L. 1951, ApJ, 114, 206
  * Denis-Alpizar et al. (2020) Denis-Alpizar, O., Stoecklin, T., Dutrey, A., & Guilloteau, S. 2020, MNRAS, 497, 4276
  * Dharmawansa et al. (2009) Dharmawansa, P., Rajatheva, N., & Tellambura, C. 2009, IEEE Transactions on Communications, 57, 915
  * Dobbs et al. (2011) Dobbs, C. L., Burkert, A., & Pringle, J. E. 2011, MNRAS, 413, 2935
  * Doi et al. (2020) Doi, Y., Hasegawa, T., Furuya, R. S., et al. 2020, ApJ, 899, 28
  * Doi et al. (2021) Doi, Y., Tomisaka, K., Hasegawa, T., et al. 2021, ApJ, 923, L9
  * Dotson et al. (2000) Dotson, J. L., Davidson, J., Dowell, C. D., Schleuning, D. A., & Hildebrand, R. H. 2000, ApJS, 128, 335
  * Einig et al. (2023) Einig, L., Pety, J., Roueff, A., et al. 2023, A&A, 677, A158
  * EMAA (2021) EMAA. 2021
  * Fiedler & Mouschovias (1993) Fiedler, R. A. & Mouschovias, T. C. 1993, ApJ, 415, 680
  * Figueira et al. (2017) Figueira, M., Zavagno, A., Deharveng, L., et al. 2017, A&A, 600, A93
  * Fissel et al. (2019) Fissel, L. M., Ade, P. A. R., Angilè, F. E., et al. 2019, ApJ, 878, 110
  * Fuente et al. (2008) Fuente, A., García-Burillo, S., Usero, A., et al. 2008, A&A, 492, 675
  * Gaudel et al. (2023) Gaudel, M., Orkisz, J. H., Gerin, M., et al. 2023, A&A, 670, A59
  * Goicoechea et al. (2016) Goicoechea, J. R., Pety, J., Cuadrado, S., et al. 2016, Nature, 537, 207
  * Goicoechea et al. (2019) Goicoechea, J. R., Santa-Maria, M. G., Bron, E., et al. 2019, A&A, 622, A91
  * Goicoechea et al. (2015) Goicoechea, J. R., Teyssier, D., Etxaluze, M., et al. 2015, ApJ, 812, 75
  * Goldsmith et al. (2008) Goldsmith, P. F., Heyer, M., Narayanan, G., et al. 2008, ApJ, 680, 428
  * Gordon et al. (2018) Gordon, M. S., Lopez-Rodriguez, E., Andersson, B. G., et al. 2018, arXiv e-prints, arXiv:1811.03100
  * Graf et al. (2012) Graf, U. U., Simon, R., Stutzki, J., et al. 2012, A&A, 542, L16
  * Gratier et al. (2017) Gratier, P., Bron, E., Gerin, M., et al. 2017, A&A, 599, A100
  * Greaves et al. (2001) Greaves, J. S., Holland, W. S., & Ward-Thompson, D. 2001, ApJ, 546, L53
  * Guerra et al. (2021) Guerra, J. A., Chuss, D. T., Dowell, C. D., et al. 2021, ApJ, 908, 98
  * Hacar et al. (2016) Hacar, A., Alves, J., Burkert, A., & Goldsmith, P. 2016, A&A, 591, A104
  * Hall & Mikesell (1949) Hall, J. S. & Mikesell, A. H. 1949, AJ, 54, 187
  * Hanawa et al. (2019) Hanawa, T., Kudoh, T., & Tomisaka, K. 2019, ApJ, 881, 97
  * Harper et al. (2018) Harper, D. A., Runyan, M. C., Dowell, C. D., et al. 2018, Journal of Astronomical Instrumentation, 7, 1840008
  * Harrison et al. (2013) Harrison, S., Faure, A., & Tennyson, J. 2013, Monthly Notices of the Royal Astronomical Society, 435, 3541
  * Hennebelle & Falgarone (2012) Hennebelle, P. & Falgarone, E. 2012, A&A Rev., 20, 55
  * Henshaw et al. (2019) Henshaw, J. D., Ginsburg, A., Haworth, T. J., et al. 2019, MNRAS, 485, 2457
  * Henshaw et al. (2016) Henshaw, J. D., Longmore, S. N., Kruijssen, J. M. D., et al. 2016, SCOUSE: Semi-automated multi-COmponent Universal Spectral-line fitting Engine, Astrophysics Source Code Library, record ascl:1601.003
  * Hernández-Vera et al. (2023) Hernández-Vera, C., Guzmán, V. V., Goicoechea, J. R., et al. 2023, A&A, 677, A152
  * Heyvaerts & Priest (1983) Heyvaerts, J. & Priest, E. R. 1983, A&A, 117, 220
  * Hildebrand et al. (1995) Hildebrand, R. H., Dotson, J. L., Dowell, C. D., et al. 1995, in Astronomical Society of the Pacific Conference Series, Vol. 73, From Gas to Stars to Dust, ed. M. R. Haas, J. A. Davidson, & E. F. Erickson, 97–104
  * Hildebrand et al. (2009) Hildebrand, R. H., Kirby, L., Dotson, J. L., Houde, M., & Vaillancourt, J. E. 2009, ApJ, 696, 567
  * Hiltner (1949) Hiltner, W. A. 1949, Nature, 163, 283
  * Hoang & Lazarian (2014) Hoang, T. & Lazarian, A. 2014, MNRAS, 438, 680
  * Hoang et al. (2022) Hoang, T. D., Ngoc, N. B., Diep, P. N., et al. 2022, ApJ, 929, 27
  * Hogerheijde et al. (1995) Hogerheijde, M. R., de Geus, E. J., Spaans, M., van Langevelde, H. J., & van Dishoeck, E. F. 1995, ApJ, 441, L93
  * Holdship (2021) Holdship, J. 2021
  * Houde et al. (2009) Houde, M., Vaillancourt, J. E., Hildebrand, R. H., Chitsazzadeh, S., & Kirby, L. 2009, ApJ, 706, 1504
  * Hwang et al. (2023) Hwang, J., Pattle, K., Parsons, H., Go, M., & Kim, J. 2023, AJ, 165, 198
  * Joos et al. (2012) Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128
  * Kalugina et al. (2012) Kalugina, Y., Lique, F., & Kłos, J. 2012, Monthly Notices of the Royal Astronomical Society, 422, 812
  * Kandori et al. (2007) Kandori, R., Tamura, M., Kusakabe, N., et al. 2007, PASJ, 59, 487
  * Kovács (2008) Kovács, A. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7020, Millimeter and Submillimeter Detectors and Instrumentation for Astronomy IV, ed. W. D. Duncan, W. S. Holland, S. Withington, & J. Zmuidzinas, 70201S
  * Kovács (2008) Kovács, A. 2008, in Millimeter and Submillimeter Detectors and Instrumentation for Astronomy IV, ed. W. D. Duncan, W. S. Holland, S. Withington, & J. Zmuidzinas, Vol. 7020, International Society for Optics and Photonics (SPIE), 70201S
  * Lai et al. (2002) Lai, S.-P., Crutcher, R. M., Girart, J. M., & Rao, R. 2002, ApJ, 566, 925
  * Larson (1981) Larson, R. B. 1981, MNRAS, 194, 809
  * Lazarian et al. (2012) Lazarian, A., Esquivel, A., & Crutcher, R. 2012, ApJ, 757, 154
  * Li & Houde (2008) Li, H.-b. & Houde, M. 2008, ApJ, 677, 1151
  * Lis et al. (1991) Lis, D. C., Carlstrom, J. E., & Phillips, T. G. 1991, ApJ, 370, 583
  * Lis & Schilke (2003) Lis, D. C. & Schilke, P. 2003, ApJ, 597, L145
  * Lithwick & Goldreich (2001) Lithwick, Y. & Goldreich, P. 2001, ApJ, 562, 279
  * Liu et al. (2019) Liu, J., Qiu, K., Berry, D., et al. 2019, ApJ, 877, 43
  * Lombardi et al. (2014) Lombardi, M., Bouy, H., Alves, J., & Lada, C. J. 2014, A&A, 566, A45
  * Lyo et al. (2021) Lyo, A. R., Kim, J., Sadavoy, S., et al. 2021, ApJ, 918, 85
  * Mandy & Martin (1993) Mandy, M. E. & Martin, P. G. 1993, ApJS, 86, 199
  * Mardia & Jupp (1999) Mardia, K. & Jupp, P. E. 1999, Directional Statistics, 2nd edn. (New York: John Wiley and Sons)
  * Matthews et al. (2002) Matthews, B. C., Fiege, J. D., & Moriarty-Schieven, G. 2002, ApJ, 569, 304
  * McClure-Griffiths et al. (2006) McClure-Griffiths, N. M., Dickey, J. M., Gaensler, B. M., Green, A. J., & Haverkorn, M. 2006, ApJ, 652, 1339
  * Meyer et al. (2008) Meyer, M. R., Flaherty, K., Levine, J. L., et al. 2008, in Handbook of Star Forming Regions, Volume I, ed. B. Reipurth, Vol. 4, 662
  * Mezger et al. (1988) Mezger, P. G., Chini, R., Kreysa, E., Wink, J. E., & Salter, C. J. 1988, A&A, 191, 44
  * Mezger et al. (1992) Mezger, P. G., Sievers, A., Zylka, R., et al. 1992, A&A, 265, 743
  * Milam et al. (2009) Milam, S. N., Woolf, N. J., & Ziurys, L. M. 2009, ApJ, 690, 837
  * Montier et al. (2015) Montier, L., Plaszczynski, S., Levrier, F., et al. 2015, A&A, 574, A136
  * Mouschovias & Ciolek (1999) Mouschovias, T. C. & Ciolek, G. E. 1999, in NATO Advanced Study Institute (ASI) Series C, Vol. 540, The Origin of Stars and Planetary Systems, ed. C. J. Lada & N. D. Kylafis, 305
  * Mouschovias & Morton (1992) Mouschovias, T. C. & Morton, S. A. 1992, ApJ, 390, 166
  * Mouschovias & Spitzer (1976) Mouschovias, T. C. & Spitzer, L., J. 1976, ApJ, 210, 326
  * Müller et al. (2005) Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, Journal of Molecular Structure, 742, 215
  * Murray (2011) Murray, N. 2011, ApJ, 729, 133
  * Naghizadeh-Khouei & Clarke (1993) Naghizadeh-Khouei, J. & Clarke, D. 1993, A&A, 274, 968
  * Nakano & Nakamura (1978) Nakano, T. & Nakamura, T. 1978, PASJ, 30, 671
  * Orkisz et al. (2019) Orkisz, J. H., Peretto, N., Pety, J., et al. 2019, A&A, 624, A113
  * Orkisz et al. (2017) Orkisz, J. H., Pety, J., Gerin, M., et al. 2017, A&A, 599, A99
  * Ostriker et al. (2001) Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980
  * Pabst et al. (2020) Pabst, C. H. M., Goicoechea, J. R., Teyssier, D., et al. 2020, A&A, 639, A2
  * Panopoulou et al. (2016) Panopoulou, G. V., Psaradaki, I., & Tassis, K. 2016, MNRAS, 462, 1517
  * Pattle et al. (2019) Pattle, K., Lai, S.-P., Hasegawa, T., et al. 2019, ApJ, 880, 27
  * Pattle et al. (2017) Pattle, K., Ward-Thompson, D., Berry, D., et al. 2017, ApJ, 846, 122
  * Pattle et al. (2018) Pattle, K., Ward-Thompson, D., Hasegawa, T., et al. 2018, ApJ, 860, L6
  * Penzias et al. (1974) Penzias, A. A., Wilson, R. W., & Jefferts, K. B. 1974, Phys. Rev. Lett., 32, 701
  * Peretto et al. (2012) Peretto, N., André, P., Könyves, V., et al. 2012, A&A, 541, A63
  * Pety (2005) Pety, J. 2005, in SF2A-2005: Semaine de l’Astrophysique Francaise, ed. F. Casoli, T. Contini, J. M. Hameury, & L. Pagani, 721
  * Pety et al. (2017) Pety, J., Guzmán, V. V., Orkisz, J. H., et al. 2017, A&A, 599, A98
  * Philipp et al. (2006) Philipp, S. D., Lis, D. C., Güsten, R., et al. 2006, A&A, 454, 213
  * Phillips et al. (1979) Phillips, T. G., Huggins, P. J., Wannier, P. G., & Scoville, N. Z. 1979, ApJ, 231, 720
  * Pillai (2009) Pillai, T. 2009, in American Astronomical Society Meeting Abstracts, Vol. 213, American Astronomical Society Meeting Abstracts #213, 324.05
  * Pillai et al. (2015) Pillai, T., Kauffmann, J., Tan, J. C., et al. 2015, ApJ, 799, 74
  * Planck Collaboration et al. (2016a) Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016a, A&A, 586, A135
  * Planck Collaboration et al. (2016b) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016b, A&A, 586, A138
  * Planck Collaboration et al. (2016c) Planck Collaboration, Aghanim, N., Alves, M. I. R., et al. 2016c, A&A, 586, A137
  * Pomarès et al. (2009) Pomarès, M., Zavagno, A., Deharveng, L., et al. 2009, A&A, 494, 987
  * Priestley & Whitworth (2022) Priestley, F. D. & Whitworth, A. P. 2022, MNRAS, 512, 1407
  * Richer et al. (1992) Richer, J. S., Hills, R. E., & Padman, R. 1992, MNRAS, 254, 525
  * Robitaille & Whitney (2010) Robitaille, T. P. & Whitney, B. A. 2010, ApJ, 710, L11
  * Roshi et al. (2014) Roshi, D. A., Goss, W. M., & Jeyakumar, S. 2014, ApJ, 793, 83
  * Ryutov et al. (2005) Ryutov, D. D., Kane, J. O., Mizuta, A., Pound, M. W., & Remington, B. A. 2005, Ap&SS, 298, 183
  * Santa-Maria et al. (2023) Santa-Maria, M. G., Goicoechea, J. R., Pety, J., et al. 2023, arXiv e-prints, arXiv:2309.03186
  * Santos et al. (2016) Santos, F. P., Ade, P., Angilè, F. E., et al. 2016, in American Astronomical Society Meeting Abstracts, Vol. 228, American Astronomical Society Meeting Abstracts #228, 405.03
  * Schöier et al. (2005) Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369
  * Shimajiri et al. (2023) Shimajiri, Y., André, P., Peretto, N., et al. 2023, A&A, 672, A133
  * Skalidis et al. (2021) Skalidis, R., Sternberg, J., Beattie, J. R., Pavlidou, V., & Tassis, K. 2021, A&A, 656, A118
  * Skalidis & Tassis (2021) Skalidis, R. & Tassis, K. 2021, A&A, 647, A186
  * Sofia et al. (2004) Sofia, U. J., Lauroesch, J. T., Meyer, D. M., & Cartledge, S. I. B. 2004, ApJ, 605, 272
  * Stewart & Federrath (2022) Stewart, M. & Federrath, C. 2022, MNRAS, 509, 5237
  * Suin et al. (2022) Suin, P., Shore, S. N., & Pavlík, V. 2022, A&A, 667, A69
  * Tahani et al. (2023) Tahani, M., Bastien, P., Furuya, R. S., et al. 2023, ApJ, 944, 139
  * Tan et al. (2006) Tan, J. C., Krumholz, M. R., & McKee, C. F. 2006, ApJ, 641, L121
  * Tang et al. (2019) Tang, Y.-W., Koch, P. M., Peretto, N., et al. 2019, ApJ, 878, 10
  * Tremblin et al. (2014) Tremblin, P., Schneider, N., Minier, V., et al. 2014, A&A, 564, A106
  * Tritsis et al. (2023) Tritsis, A., Basu, S., & Federrath, C. 2023, MNRAS, 521, 5087
  * Tritsis et al. (2022) Tritsis, A., Federrath, C., Willacy, K., & Tassis, K. 2022, MNRAS, 510, 4420
  * van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627
  * van der Werf et al. (1996) van der Werf, P. P., Stutzki, J., Sternberg, A., & Krabbe, A. 1996, A&A, 313, 633
  * Walch & Naab (2015) Walch, S. & Naab, T. 2015, MNRAS, 451, 2757
  * Wang et al. (2019) Wang, J.-W., Lai, S.-P., Eswaraiah, C., et al. 2019, ApJ, 876, 42
  * Ward-Thompson et al. (2017) Ward-Thompson, D., Pattle, K., Bastien, P., et al. 2017, ApJ, 842, 66
  * Wardle & Kronberg (1974) Wardle, J. F. C. & Kronberg, P. P. 1974, ApJ, 194, 249
  * White et al. (1999) White, G. J., Nelson, R. P., Holland, W. S., et al. 1999, A&A, 342, 233
  * Yin et al. (2021) Yin, C., Priestley, F. D., & Wurster, J. 2021, MNRAS, 504, 2381
  * Young Owl et al. (2000a) Young Owl, R. C., Meixner, M. M., Wolfire, M., Tielens, A. G. G. M., & Tauber, J. 2000a, ApJ, 540, 886
  * Young Owl et al. (2000b) Young Owl, R. C., Meixner, M. M., Wolfire, M., Tielens, A. G. G. M., & Tauber, J. 2000b, ApJ, 540, 886
  * Zuckerman & Evans (1974) Zuckerman, B. & Evans, N. J., I. 1974, ApJ, 192, L149
  * Zweibel (2004) Zweibel, E. G. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 323, Star Formation in the Interstellar Medium: In Honor of David Hollenbach, ed. D. Johnstone, F. C. Adams, D. N. C. Lin, D. A. Neufeeld, & E. C. Ostriker, 97

##  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 Sec. A.3.

###  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 ratio
pixels are included in the dust polarization analysis (see Sec. 2.1). Fig. 6
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-1 (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-1. 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 ($\pm
5\%$) intervals centered at PACS fluxes of 30,000 and 10,000 MJy sr-1, marked
with vertical dotted lines in Fig. 6. 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.

Figure 6: A 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 ratio pixels included in the polarization analysis
are shown. As described in the text, the HAWC+ fluxes were scaled by a factor
of 1.12.

Figure 7: Comparison between Stokes parameters, $I$, $Q$, and $U$ (each row)
for two setups of data reduction (the one used in this work is shown on 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). Figure 8: Same as in Fig. 7, but for the polarization
angles (see Eq. 6 in Sec. 2.1).

###  A.2 SOFIA HAWC+ Band E dust continuum

We show SOFIA HAWC+ Band E dust continuum map in Fig. 9 at its native
resolution of $18.7^{\prime\prime}$. 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,
$\theta_{D}-\theta_{E}$. Finally, we measure the rms of this difference, as
described in Sec. 5.3, Eq. 19. We show a map of $\theta_{D}-\theta_{E}$ in the
left panel of Fig. 10. The right panel in Fig. 10 shows the corresponding rms
of the map shown in the left panel.

Figure 9: SOFIA HAWC+ $214\,$µm (Band E) dust continuum map at
$18.2^{\prime\prime}$ angular resolution corresponding to linear scales of
$\sim 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. Figure 10: Difference in polarization angles
measured from Band D and E data is shown on the left panel. The right panel
shows the rms of the angle differences.

###  A.3 Comparison with the literature data

In Fig. 11, 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. 11. 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).

Figure 11: 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

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

###  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 Tab. 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_{1},p_{2},p_{3}$
and $p_{4}$. The first parameter is the antenna temperature multiplied by the
optical depth of the CN emission:

$p_{1}=T_{\mathrm{ant}}\cdot\tau$ (27)

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}.$ (28)

The next parameter is the full width at half maximum (FWHM) of the hyperfine
components:

$p_{3}=FWHM,$ (29)

and the last parameter is the opacity of all components of the studied
multiplet:

$p_{4}=\tau.$ (30)

We assume that the optical depth of each component can be described using a
Gaussian function of velocity:

$\tau_{\mathrm{i}}(\upsilon)=\tau_{\mathrm{i}}\cdot\exp{\left[-4\ln{2}\left(\dfrac{\upsilon-\upsilon_{\mathrm{0,i}}}{p_{3}}\right)^{2}\right]}.$
(31)

The opacity of the multiplet is then calculated as a sum of opacities of all
hyperfine components:

$\tau=\sum_{\mathrm{i}=1}^{N}\tau_{i}$ (32)

where $\tau_{\mathrm{i}}$ is the opacity of the $i$-th component.

Finally, the antenna temperature, $T_{\mathrm{ant}}$ is thus derived as the
ratio between the first and the fourth fitting parameters, $p_{1}$, and
$p_{4}$:

$T_{\mathrm{ant}}=\dfrac{p_{1}}{p_{4}}\left(1-e^{-\tau(\upsilon)}\right).$
(33)

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_{\mathrm{ex}}=\dfrac{h\nu}{k}\left\\{\ln{\left[1+\dfrac{h\nu}{kT_{\mathrm{ant}}}\left(1-e^{-\tau_{\nu}}\right)\right]}\right\\}^{-1}.$
(34)

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 $\backslash$fit/minimize and add
$\backslash$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 Tab. 1. In the case of the hyperfine fitting,
CUBE fits four parameters, $p_{1},p_{2},p_{3},$ and $p_{4}$, described in Eq.
27, 28, 29, and 30 respectively.

Before fitting, we define a two-dimensional mask based on the S/N ratio 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 10 km/s and 7 km/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. 12.

Figure 12: Fitting parameters for the hyperfine CN emission in NGC 2024
measured for both velocity components (top and bottom row each). The first
parameter, $p_{1}$ is shown in the left panel, followed by the centroid
velocity, FWHM and the final parameter, the opacity ($p_{4}$).  Figure 13:
Fitting parameters for the HCO+ emission assuming 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.

###  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 (e.g., 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 $\backslash$fit/minimize
in CUBE and add $\backslash$gaussian to specify the shape of the spectral
line. In CUBE, the Gaussian line is described using three parameters: the area
(in K km/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. 13.

##  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. 14, 15, 16 and 17 for CN emission. In the case of HCO+
emission, we show these results in Fig. 18, 19, 20and 21. We also include a
$\chi^{2}$ 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. 22.

Figure 14: 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 colour bar shows values of excitation
temperature, $T_{\mathrm{ex}}$ 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, $\tau$, and the third row shows the
peak temperature, $T_{\mathrm{peak}}$. The last row shows the $\chi^{2}$
minimization of the modelled peak temperature and the observed value computed
from the CN(1–0) spectra.

Figure 15: The same as in Fig. 14, but for the eastern part of the bubble.

Figure 16: The same as in Fig. 14, but for the filament.

Figure 17: The same as in Fig. 14, but for the overlap region.

Figure 18: The same as in Fig. 14, but for HCO+(1–0).

Figure 19: The same as in Fig. 18, but for the eastern side of the bubble.

Figure 20: The same as in Fig. 18, but for the filament.

Figure 21: The same as in Fig. 18, but for the overlap region.

Figure 22: A beam averaged spectrum of CN and HCO+ taken within four different
region shown on the left panel of Fig. 1 and in Fig. 2. Coloured line
represent the corresponding model spectrum infrerred from the non-LTE RADEX
modelling for physical conditions presented in Sec. C and shown in Tab. 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.

###  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\times 4$ to $14\times 14$ pixels, respectively. We show these results in
Fig. 23. 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.

Figure 23: The angle rms as a function of the size of the sliding window.

###  D.2 Histogram analysis

In this section, we describe an alternative approach to compute
$\widehat{\sigma}_{\mathrm{c}}(\varphi)$. 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. 24. In most cases, the histogram corresponds to a single
Gaussian. In Fig. 25, we show comparison between
$\widehat{\sigma}_{\mathrm{c}}(\varphi)$ 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.

Table 5: Results on magnetic field strength from histogram analysis.

Region |  $\widehat{\sigma}_{\mathrm{c}}(\varphi)$ [deg] |  $B_{\mathrm{DCF,CN}}$ [µG] |  $B_{\mathrm{DCF,HCO^{+}}}$ [µG] |  $B_{\mathrm{ST,CN}}$ [µG] |  $B_{\mathrm{ST,HCO^{+}}}$ [µG]
---|---|---|---|---|---
Bubble West | 4.7 $\pm$ 0.2 | 53 $\pm$ 27 | 106 $\pm$ 30 | 22 $\pm$ 11 | 43 $\pm$ 12
Bubble East | 2.56 $\pm$ 0.04 | 317 $\pm$ 72 | 318 $\pm$ 72 | 95 $\pm$ 22 | 95 $\pm$ 22
Filament | 4.4 $\pm$ 0.2 | 91 $\pm$ 38 | 141 $\pm$ 48 | 36 $\pm$ 15 | 55 $\pm$ 19
Filament+Bubble | 6.5 $\pm$ 0.5 | 158 $\pm$ 44 | 180 $\pm$ 47 | 75 $\pm$ 21 | 86 $\pm$ 22

Figure 24: 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. Figure 25: Comparison of angle dispersion using two different
approaches, sliding window (y-axis, see Sec. 5) and histogram analysis
(x-axis).

###  D.3 $B_{\mathrm{POS}}$ using different approaches

We compare results on POS $B$-field strengths using several methods in Fig.
26. 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.

Figure 26: The 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.