The magnetic field in the Flame nebula

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 H II 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 H II 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 12 CN(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 H II 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 H II 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 H II 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 H II shell (gas impacted by the stellar feedback) and the filament (site of the current star formation).


Introduction
The question of what controls the star formation efficiency in molecular clouds has long been at the center of star formation research.Early studies (Zuckerman & Evans 1974) showed that if all the gas within dense interstellar clouds were to collapse freely under self-gravity, the star formation rate in the Milky Way would be two orders of magnitude higher than the observed rate of 2 M ⊙ /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 lowdensity parts of the cloud, which contain up to 90% of the cloud mass (Battisti & Heyer 2014), disperse over the scales of free fall time as a result of turbulence generated by galactic shear or by energy input from supernova explosions (Dobbs et al. 2011;Walch & Naab 2015), or due to radiative and mechanical stellar feedback from high-mass stars that formed early on during the evolution of the cloud (Murray 2011; Colín et al. 2013;Pabst et al. 2020;Chevance et al. 2022;Suin et al. 2022).
In addition to preventing the gas from collapsing, the turbulence and the magnetic field can also bolster the star formation processes.For example, the coupling between the magnetic field and the neutral gas can allow parts of low-density clouds to fragment and initiate star formation (Fiedler & Mouschovias 1993).In this scenario, while neutral gas is collapsing, a part of the magnetic flux is removed, which further contributes to the gravitational instability (Mouschovias & Ciolek 1999;Lazarian et al. 2012;Priestley & Whitworth 2022;Tritsis et al. 2023).
However, it has been difficult to test these star formation models observationally.Despite all of the effort made over the past 40 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 has energetic properties (at least in terms of the total kinetic energy) similar to those of an identical cloud supported against gravity by turbulence.In addition, measurements of magnetic fields are observationally extremely challenging, as they typically involve relatively weak signals, as is the case for Zeeman gas line splitting (Crutcher 2012) or (sub)millimeter/far-infrared (FIR) dust polarization measurements (see, Pattle et al. 2019, for a recent review).
Interstellar dust thermal emission is polarized (Hall & Mikesell 1949;Hiltner 1949) due to the presence of B fields and nonspherical dust grains in the interstellar medium (ISM).The explanation for this phenomenon was firstly proposed by Davis & Greenstein (1951) as a paramagnetic alignment with the magnetic field.Later on, this process was described by the radiative alignment torques (RAT) theory (Hoang & Lazarian 2014;Andersson et al. 2015, and references therein), in which the minor axis of dust grains is aligned parallel to the direction of the B field.Consequently, the (sub)millimeter dust continuum emission is polarized perpendicular to the direction of the component of the B field in the plane of the sky (POS).
Regardless of the observational challenges, technological advances in the past decade are now allowing for polarized dust emission to be measured over increasingly extended regions.The Planck satellite provided an all-sky map of the magnetic field geometry in the diffuse ISM, albeit at a low angular resolution of 5-15 ′ .These observations revealed that the galactic B field is intertwined with the filamentary structure of the ISM.In particular, the POS orientation of the B field and the filamentary structures is correlated with the observed column density, with the field largely parallel to diffuse structures (striations) while perpendicular to dense filaments (Planck Collaboration 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 (HII bubble shells, Chapman et al. 2011).
The Balloon-borne Large Aperture Submillimeter Telescope for Polarimetry (BLASTPol) has also observed a similar B-field morphology at a higher angular resolution (2.5 ′ , Fissel et al. 2019), which has been interpreted as evidence of the magnetic field dominating the energy balance in the diffuse gas.In contrast, gravity is dominant in the dense, star-forming regions.This interpretation has recently been reinforced by high angular resolution (14 ′′ ) ground-based observations using the POL-2 instru-ment on the James Clerk Maxwell Telescope (JCMT; Pattle et al. 2017;Liu et al. 2019;Wang et al. 2019), which have shown that field lines bend along dense star-forming filaments, possibly as a result of gravitationally driven flows (Goldsmith et al. 2008;Chapman et al. 2011;Hwang et al. 2023).
A transition from magnetically dominated to gravitationally dominated gas requires a redistribution of magnetic flux (Tritsis et al. 2022).In principle, such a transition occurs before the cores form (Ching et al. 2022), which implies that cores should be generally super-critical.Nevertheless, a change in the gravitational stability of the gas should be accompanied by a corresponding change in its kinematic properties, arguing for the necessity of combining high angular resolution dust polarization and velocity-resolved molecular emission imaging.Such studies have only recently started and cover relatively small isolated fields.A recent study (Tang et al. 2019) has indeed investigated the relation between dense gas velocity gradients (as imaged in N 2 H + ) and the magnetic field direction in massive star-forming infrared dark clouds and concluded that the two are strongly correlated as a result of gravity dragging matter toward the center of the massive ridge.Extending such studies to larger fields is of crucial importance and the High-resolution Airborne Wideband Camera-Plus (HAWC+, Harper et al. 2018) instrument on the Stratospheric Observatory for Infrared Astronomy (SOFIA) offered exceptional capabilities for mapping the magnetic field geometry, as demonstrated, for example, by the observations of the OMC-1 region at 53, 89, 154, and 214 µm (Chuss et al. 2019).The capabilities of HAWC+ for large-area magnetic field mapping have further improved with the commissioning of the on-the-fly map (OTFMAP) polarimetric mode for fast wide-field polarimetric imaging.
When estimating the strength of the POS B field, it is critical to know the level of turbulence in the gas because the small-scale turbulence causes local deviations from the mean direction of the magnetic field.The Davis-Chandrasekhar-Fermi (DCF) method (Davis & Greenstein 1951;Chandrasekhar & Fermi 1953) provides a way of calculating the POS B-field strength using information about gas density, turbulence, and changes in the direction of the magnetic field.This method assumes linear geometry and sub-Alfvénic (magnetically dominated) turbulence.Therefore, it is necessary to obtain information about the turbulence of the gas and dust polarization to use this method.Several variations of the DCF have been developed over the last few decades.The main modification lies in taking into account the number of turbulent cells present along the line of sight and those captured within the beam by including a correction factor Q, which takes a value between zero and one (see Eq. 16 in Ostriker et al. 2001).In another variation, the deviation in the mean direction of the B field is substituted by the ratio between the ordered and turbulent magnetic field component (Houde et al. 2009, and references therein).The DCF method assumes that imcompressible motions cause the dispersions of the observed polarization angles, which is not always applicable within the ISM.Skalidis & Tassis (2021) provided an alternative method for deriving the magnetic field strength from dust polarization measurementsthe Skalidis-Tassis (ST) method.In addition to Alfvénic modes, this study accounts for the presence of compressible motions in the gas without discarding the anisotropic nature of the turbulence, providing a physically motivated approach to estimating the magnetic field strength (Skalidis et al. 2021).For this work, we used the ST method to investigate highly structured regions in the Flame nebula in the Orion B complex, while also presenting the results from the modified DCF method (for example,

Expanding bubble
Expanding bubble   (Mezger et al. 1988(Mezger et al. , 1992)).In addition, we labeled the edges of the expanding HII region and the direction of the stellar feedback driven by the radiation produced by recent star formation.Ostriker et al. 2001;Crutcher 2004;Lyo et al. 2021) for comparison.
Although OMC-1 is the closest and best-studied high-mass star-forming region, it is not a typical cloud to study the role of the magnetic field in star formation.The reason is that OMC-1 is particularly dense (10 5 cm −3 ) and very active in terms of star formation, affected by strong shocks (including an explosive outflow resulting from the merger of three massive protostars, see, Bally 2008).These shocks exhibit exceptionally intense UV illumination, with enhancement factors of G 0 ∼10 4 , up to 10 5 relative to the standard interstellar radiation field (ISRF) (in the vicinity of the Trapezium cluster, Goicoechea et al. 2015Goicoechea et al. , 2019)).The presence of a strong radiation field implies that the thermal pressure is high, and the magnetic field plays a limited role in the photon-dominated region (PDR) gas dynamics.
Fortunately, the vast quantity of molecular data already obtained in Orion B makes this region a particularly well-suited location for studying the correlation between changes in the magnetic field geometry and associated gas velocity gradients.The Orion B giant molecular cloud (GMC) hosts NGC 2024, one of the closest high-mass star-forming regions (at a distance of d = 410 pc; Cao et al. 2023), and it is overall more representative of a standard GMC in our galaxy and normal galaxies, with the far-ultraviolet (FUV) radiation field G 0 in the range of ∼ 10 3 − 10 4 (Santa-Maria et al. 2023).
In a recent study, Orkisz et al. (2017) used 12 CO(1-0) and 13 CO(1-0) lines to characterize observationally the ratio of compressive versus solenoidal motions in the turbulent flow and to relate this to the star formation efficiency in various regions of Orion B. Orkisz et al. (2019) accurately analyzed the dynamics of the filamentary network using C 18 O(1-0).Most identified filaments in Orion B are low-density, thermally subcritical structures, not collapsing to form stars.Only about 1% of the Orion B cloud mass can be found in super-critical, star-forming filaments, consistent with the low overall star formation efficiency of the region (Orkisz et al. 2019).
NGC 2024 (Fig. 1) is located east of Alnitak (ζ Ori) in the Orion B complex (Meyer et al. 2008).This region contains a massive and young HII region (age 2 • 10 5 yr, Tremblin et al. 2014, and private communication), deeply embedded in dust, located in the foreground, and extending to the south (Barnes et al. 1989, see gray contrours showing the network of filaments from Orkisz et al. (2019) in the left panel and labels in the right panel in Fig. 1).The north-south filament in NGC 2024 is super-critical (Orkisz et al. 2019) and the site of ongoing star formation as recently confirmed in the southern part of this filament observed by NOEMA (Shimajiri et al. 2023).
The dust bar observed along the line of sight is visible in the ESO-Vista image (right panel in Fig. 1) as the dark extinction pattern across the image, leading to an apparent cooler dust temperature derived from Herschel observations (Lombardi et al. 2014).The central part of NGC 2024 contains warm dust and gas, heated by the HII region, as well as embedded protostellar objects (Mezger et al. 1988(Mezger et al. , 1992;;Lis et al. 1991) located in the background (light blue points in the right panel in Fig. 1).The young HII region is expanding, strongly impacting its parental cloud, and creating sharp ionization fronts toward the south, which makes it a good example demonstrating how such sys-Article number, page 3 of 36 tems can efficiently exert stellar feedback.The edges of the bubble are seen toward the west and east of the center of NGC 2024 (labeled as Eastern and Western Loop in Fig. 1).

Dust polarization measurements using SOFIA HAWC+
Our work employs the dust polarization measurements acquired using HAWC+ on SOFIA in September 2021 for program 09_0015 (PI: D. Lis).Specifically, we observed NGC 2024 at 154 µm (Band D) and 214 µm (Band E) with HAWC+ in polarization mode on flights F779 (8 September 2021), F780 (9 September 2021), F782 (11 September 2021), F783 (14 September 2021), and F784 (15 September 2021).The map obtained at each wavelength comprises ten 20 ′ × 7 ′ strips in a weave pattern: five vertical and five horizontal.These strips were designed using the On-The-Fly scan mode of the HAWC+ camera, which covers the requested area using Lissajous patterns on the sky (Harper et al. 2018).Each strip was repeated at least twice.
The level 0 data for each flight were downloaded from the NASA/IPAC Infrared Science Archive (IRSA) and reduced in July 2023 using the SOFIA Data Reduction software (SOFIA Redux version 1.3.0;Clarke & Vliet 2023).The HAWC+ scan mode reduction package from SOFIA Redux was initially built in Python from the Java-based CRUSH data reduction software (Kovács 2008).All available files at a given wavelength (104 each for Band D and Band E) were loaded in SOFIA Redux to be reduced using the default parameters of the software except for two options at the Compute Scan Map step of the pipeline.Specifically, the fixjumps option was set to True to filter out flux jumps in individual detectors during observations, and the rounds option was set to 40 to improve the recovery of astronomical large scale flux from the background subtraction.Appendix A.1 shows the comparison of the Stokes I values from the resulting data to archival Herschel fluxes and the improvement relative to the Level 4 data available on IRSA at the time of writing.
The pixel scale of the HAWC+ polarization maps and the effective beam size is 3.4 ′′ and 14.0 ′′ , respectively, in Band D, and 4.6 ′′ and 18.7 ′′ in Band E. The final data products for both Band D and E contain the I, Q and U Stokes parameters, the polarization fraction P and angle θ, 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 is calculated using the following equation: The uncertainty δI ′ p on the polarized intensity I ′ p is then where δQ and δU are the uncertainties on Stokes Q and U, respectively.
The polarized intensity I ′ p has to be corrected for the bias created by the quadratic addition of the noise in the Stokes Q and U maps (Wardle & Kronberg 1974;Naghizadeh-Khouei & Clarke 1993).This de-biased polarized intensity I p (although this is not the only way to do so, Montier et al. 2015), used in this work, is calculated using with uncertainty δI p = δI ′ p .The polarization fraction P is then obtained from with the uncertainty δP given by The polarization angle θ p is defined using the Stokes parameters Q and U as with the uncertainty δθ p given by Since the thermal emission from interstellar dust grains is preferentially polarized perpendicular to the plane of the sky magnetic field lines (Hoang & Lazarian 2014;Andersson et al. 2015, and references therein), the direction of the magnetic field in the plane of the sky can be obtained by adding π/2 to Eq. 6.
Prior to the data analysis, we filter out low signal-to-noise ratio (S/N) data points.After masking the data, we keep pixels that show an S/N ≥ 50 in total intensity, ≤ 30 o in polarization angle uncertainty, and ≤ 30% in polarization fraction.Our HAWC+ observations filter out a non-negligible fraction of low-level extended emission.We estimate the amount of missing flux in the SOFIA observations by comparing our data with PACS measurements in Appendix A.1.The data set used in this study was reduced using a larger number (40) of iterations, or rounds, than the default (15) used for the original Level 4 data available on IRSA.A larger number of iterations during data reduction typically improves the recovery of diffuse large-scale emission for scan mode maps.We compare our reduced maps of each Stokes parameter and the polarization angle with those derived from the standard pipeline setup and show them in App.A.1 in Figs.A.2 and A.3.We find an overall good agreement between these two data sets.

IRAM 30-meter observations of Orion B
In our work, we make use of information of the 12 CN(1-0) and HCO + (1-0) emission from the ongoing IRAM-30m ORION-B Large Program (PIs: M. Gerin & J. Pety, see the left panel in Fig. 1 taken from Pety et al. 2017).ORION B images a 5 square degree field (∼20 pc across) in the Orion B molecular cloud, at an angular resolution of 26 ′′ (10 4 AU, 0.05 pc) in at least 30 molecular lines in the (71-79) and (84-116) GHz range with a spectral resolution ∼0.6 km s −1 .These include common molecular tracers such as 12 CO(1-0), HCO + (1-0), HCN(1-0), and CS(2-1), as well as their optically thin isotopologues, which have narrow line widths and are most sensitive to kinematic variations.The resulting wide-field hyper-spectral data cube is genuinely unique in terms of its massive information content (∼820,000 Article number, page 4 of 36 Bešlić: The magnetic field in the Flame nebula Table 1.Components of the hyperfine structure of the 12 CN(1-0) transition, their relative offsets, and intensities relative to the sum of all components (Milam et al. 2009) pixels, ∼240,000 spectral channels per pixel), enabling an unprecedented characterization of the physical structure, chemistry, and dynamics of a GMC, and their connection to its star formation activity.
The J = 1 − 0 lines were observed in 2013-2020 in the context of the ORION-B Large Program, using a combination of the EMIR receivers and FTS spectrometers.The data reduction is described in Pety et al. (2017).It uses the standard methods provided by the GILDAS1 /CLASS software.The data cubes were reprojected on the same astrometric grid as the SOFIA data.We use the data cubes at their original spatial resolution of 30 ′′ and 23 ′′ for the HCO + and 12 CN lines, respectively.The velocity spacing is 0.5 km/s.The data are calibrated on the main beam temperature scale.The achieved noise levels are 0.26 and 0.58 K (see Tab. 2 in Gratier et al. 2017), respectively, over the studied field of view.
We use spectroscopic data of cyanide ( 12 CN) and formyl cation (HCO + ) (J = 1 − 0 transition line) to trace the UV-irradiated gas (Bron et al. 2018).CN is present in UVilluminated 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 the opacity (Milam et al. 2009) and excitation temperature of the (1-0) transition to be derived, enabling its excitation to be constrained accurately and for limits on the gas density to be obtained.
HCO + , similarly to CN, traces photon-dominated regions (Young Owl et al. 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 (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 C 18 O molecule, which gets photodissociated.Therefore, using CN and HCO + emission in our work, together with the FUV heated dust traced by SOFIA HAWC+ observations, we can directly assess the gas impacted by the radiative stellar feedback in NGC 2024.

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 ∼ 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 circularshaped, 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 The dust continuum emission is strongest toward the center of NGC 2024.A significant emission is observed along the filament that is in the front of the ionization region caused by the young massive star IRS 2b (the central part of our map, Bik et al. 2003), and along the ionization fronts on the southeastern and southwestern part of NGC 2024 (labeled as Bubble West and Bubble East, and shown in the middle row in Fig. 2).
In addition, we compare our magnetic field directions to the direction inferred from the near-infrared (NIR) polarization of young stars (Kandori et al. 2007).The NIR polarization directly traces the direction of the magnetic field (that is, there is no π/2 difference as for the FIR dust polarization -Sec.2.1).We show this comparison in the top panels of  The area labeled with dark blue contours in the top panels in Fig. A.6 is particularly interesting because we observe an outflowing feature in the HCO + emission at velocities of ∼ 14 km s −1 .The outflow originates from the FIR 5 source located in the dense molecular cloud behind the expanding HII region (Richer et al. 1992;Greaves et al. 2001;Choi et al. 2015).Dust continuum emission at 154 µm possibly traces the outflow in this region (see also the bottom middle panel in Fig. 3), as we note that magnetic field lines follow the direction of the outflow.We have not observed this behavior at 214 µm.
The magnetic flux is frozen within the molecular gas, and as a consequence, the B field will trace its morphology.Therefore, the local environmental conditions that shape the distribution of molecular gas will also impact the magnetic field.On the one hand, the magnetic field is highly ordered in some regions of NGC 2024; for instance, B field follows the dusty filament to the south of NGC 2024 (see, for instance, the bottom left panel of Fig. 2), except for the very southern part of the filament (close to the bottom left corner of the top panel in Fig. 2), where the magnetic field is perpendicular to the filament.At the northern part of NGC 2024, the magnetic field direction varies from following the filament to being perpendicular.
We observe the nearly horizontal direction of the B-field at the edges of the expanding HII region (middle panels in Fig. 2).Gas affected by the stellar feedback (for example, 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 (for example, Tahani et al. 2023).On the other hand, the magnetic field appears chaotic in the central and densest parts of NGC 2024.This area shows great complexity, as it results from the mixture of the filament located in the front, the HII region and the dense molecular gas in the background (see Fig. 8 in Matthews et al. 2002, or Fig. 5 in Roshi et al. 2014).One possible explanation for this observed morphology of the magnetic field is that the magnetic field changes its orientation relative to the filament depending on the column density contrast between the filament and the background emission, from being parallel to perpendicular to the molecular gas distribution (∆N H 2 > 10 20 cm −2 , for example, 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. 2020Doi et al. , 2021)), which can also explain the magnetic field direction at the northern part of NGC 2024 and the very southern part of the filament.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.

Characterizing the turbulence and density structures in NGC 2024
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: where σ obs,mol is the measured FWHM of a molecular line, and σ TH,mol is the line width of the thermal component.We calculate the thermal broadening using the following equation: where k is Boltzman's constant, T k is the kinetic temperature obtained using the CO(1-0) measurements (Orkisz et al. 2017), and m mol is the mass of a molecule.A final correction that we apply in our analysis is for opacity broadening of molecular lines, using the following equation: where the factor β is a function of the optical depth (τ 0 ) at the line center, defined as (Phillips et al. 1979;Hacar et al. 2016;Orkisz et al. 2017): Article number, page 8 of 36 Bešlić: The magnetic field in the Flame nebula The equations presented above are crucial at subparsec scales since these effects significantly contribute to the measured line width at these scales.For example, the optical depth effect could broaden the line by a few tens of percent, which is the case of CN and HCO + .This is discussed in the following section.
The FWHM of a molecular line, corrected for the contribution mentioned above of the thermal and opacity broadening, can be derived in multiple ways.For example, the FWHM can be derived directly from measuring the second order moment map (right panel in Fig. 3, 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.

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 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 kin ) derived from the 12 CO(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)•10 14 cm −2 for CN and (1-4)•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 H 2 is comprised of para-and ortho-H 2 , assuming the temperature dependence of the ortho-to-para ratio: 9 • e −170/T kin (Mandy & Martin 1993).Here, we assume a fixed electron fraction (electron-to-H 2 density ratios), f e = 1.4 • 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 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 H 2 ratio and electron fraction, we run a grid in three independent variables: line width, column density, and H 2 volume density.
We consider all relevant collision partners for CN and HCO + in our input files.For the CN emission, we take into account its hyperfine structure (Müller et al. 2005), and include collisions with ortho-and para-H 2 (Kalugina et al. 2012), as well as electrons (Harrison et al. 2013;Santa-Maria et al. 2023).The input file with collision partners comes from the combination of two data files from the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005) and Excitation of Molecules and Atoms for Astrophysics (EMAA; EMAA 2021) to include all relevant collisional partners and take into account the hyperfine structure of CN.In the case of HCO + , we consider collisions with ortho-, para-H 2 (Denis-Alpizar et al. 2020) and electrons (Fuente et al. 2008) as well.
The output parameters from the RADEX modeling are the excitation temperature, T ex and opacity, τ.We additionally compute the peak temperature of the spectrum generated from each model using SpectralRadex.In the case of CN, we calculate excitation temperature and the opacity of each component of the multiplet, which is scaled using information about relative intensities (Tab.1).
We compute model CN(1-0) and HCO + (1-0) spectra assuming the optical depth to be a Gaussian function of frequency, τ(ν): The best-fit model is found based on the modeled spectrum for each combination of (n H 2 , N mol , ∆υ).We perform a χ 2 minimization for the modeled and observed peak temperature, T 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. C.9 in App. C.

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: where ρ is the gas volume density, σ NT is the nonthermal velocity dispersion (σ NT = ∆υ/ √ 8 ln 2), where ∆υ is the measured FWHM of the line, corrected for the line broadening and opacity effects -Eq.8 and 10 -Sec.4.2, and σ c (φ) 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 (Lyo et al. 2021), Eq. 13 can be expressed as: We have included in the above a factor of 0.5 for overestimation of the magnetic field strength due to line of sight integration effects (see, for instance, Ostriker et al. 2001).
Article number, page 9 of 36 A&A proofs: manuscript no.NGC2024_final Table 2. Derived physical conditions of molecular gas traced by CN and HCO + emission for the 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 modeling using RADEX: volume number density, and measured turbulent line widths.

Skalidis-Tassis method
The DCF method is widely used in the literature to derive the strength of the magnetic field (see, Pattle et al. 2019).However, it assumes that isotropic turbulence and Alfvénic waves in an incompressible medium cause the observed dispersion in polarization angle.Other mechanisms such as magneto-hydrodynamic waves (Heyvaerts & Priest 1983) where ρ, σ NT , and σ c (φ) are the same as in Eq. 13.
Similarly as in Eq. 14, after substituting all constants, the above equation becomes: 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.

Sliding window
We present the dispersion of the mean direction of the magnetic field in Fig. 4. We produced this map by using a "sliding window" to remove the impact of gradients of the large-scale magnetic field ("unsharp-masking", Dharmawansa et al. 2009;Pattle et al. 2017).The "sliding window" technique is based on computing the mean and standard deviation of the magnetic field orientation within the window that is three times bigger than the beam and has 12 × 12 pixels.The size of the sliding window is made to ensure that we remove a large-scale magnetic field contribution without losing information on the small-scale perturbation of the magnetic field (Pattle et al. 2017).As shown in Mardia & Jupp (1999), a suitable way to estimate the mean and standard deviation from the 12×12 = 144 directions φ(l) of each window is to compute: Article number, page 10 of 36 Bešlić: The magnetic field in the Flame nebula where L is the number of pixels within the sliding window.Defining a and b such that z = a + i • b, one computes the circular mean using and the circular standard deviation using Moreover, to consider the uncertainty associated with the different directions φ(l), we use a weighted mean to compute z with weight inversely proportional to the variance.The error of the circular standard deviation is computed as follows: We show maps of the magnetic field angle, the mean angle computed using Eq.18 and 12 × 12 window, and the circular standard deviation (Eq.19) in Fig. 4. As seen in the left panel of Fig. 2, the magnetic field direction differs along the edges of the expanding shell and the filament.A similar behavior we observe in the mean angle of the POS magnetic field is shown in the middle panel of Fig. 4. The standard deviation of the angle that is caused by the small-scale turbulence also shows different values within the borders of the HII region and the filament.In particular, the measured σ c is higher along the filament (more than 10 degrees), whereas it is a few degrees at the edges of the bubble.This result indicates that the magnetic field is possibly stronger at the edges of a bubble than in the filament.This result could also be a consequence of a complex geometry of NGC 2024: near the edge of a bubble, we observe a coherent limb-brightened structure, whereas, near the center of NGC 2024, we find a superposition of several components.
The ratio between the line width and the angle dispersion in Eq. 14 is shown in Fig. 5.The two panels show the measured ratios between the line width of the CN line (left) and HCO + emission (right) to the dispersion of the angle of the magnetic field.We used the line width information derived from the second moment map shown in the right panels in Fig. 3.This approach gives a robust estimate of the range of the magnetic field strength assuming a specific number density of the gas because the second moment map provides information on the line width of all emission along the line of sight without considering possible impacts (see 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.

Magnetic field strength in NGC 2024
We calculate the strength of magnetic field across the beamaveraged 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: where ∆n, ∆(∆υ), and ∆( σ c (φ)) are measured uncertainties of n, ∆υ and σ c (φ) (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 1 4 in front of the last term in the equation.
We measure the uncertainties of n, σ NT from the χ 2 minimization (see the bottom panel in Fig. C.1-C.8).We select the range in parameter space within which the χ 2 reaches 0.5 to derive corresponding n, and σ NT .Thus, we found ∆n to be 10 3 cm −3 and ∆σ 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. D.4 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 POS strength is in the range of ∼ 20 − 90 µG.The magnetic field strength inferred from CN is overall smaller than that derived using HCO + , mainly due to HCO + showing slightly larger line widths than the CN in some regions, which directly impacts the strength of the magnetic field.Due to ambipolar diffusion effects (Zweibel 2004;Tritsis et al. 2023), it is expected that the emission lines of neutral molecular species (CN in our case) are wider than those of the ion molecular species (such as HCO + ) (Li & Houde 2008;Yin et al. 2021).However, in our case, we observe that HCO + is broader than the CN spectral lines, even after corrections for the opacity broadening.This result has two implications.Firstly, our observations did not spectrally resolve HCO + emission, which potentially results in larger line widths caused by the blending of velocity components.Secondly, it is possible that our measurements do not probe ambipolar diffusion scale (∼ 10 −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 toward the eastern side of the bubble than on the western edge.We found lower line widths of CN and HCO + at the western side of the bubble.The weakest magnetic field is derived toward the filament, although we note the B POS derived from HCO + in the western side of the bubble and filament to be comparable given the uncertainties.

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 Article number, page 11 of 36 Table 3. Measured angle dispersion from HAWC+ data and results on the 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 brand 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 modeled magnetic field in NGC 2024 with two components.The first component is related to the dense and dust obscured part, where magnetic field lines follow locations of the FIR sources (right panel in Fig. 1).The second component of magnetic field is linked to gas affected by the stellar feedback.In our case, we observe regions where these two features are dominating: Band D dust polarization traces the filament in the central and southern part of NGC 2024, whereas it is impacted by the stellar feedback at eastern and the western parts at ionization front (Eastern and Western Loop in the right panel in Fig. 1).
We overlay dust polarization vectors at 100 µm from Hildebrand et al. (1995) and at 850 µm from Matthews et al. (2002) in the bottom row of Fig. A.6 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.
Article number, page 12 of 36 Bešlić: The magnetic field in the Flame nebula The magnetic field lines are mainly parallel to the filament in NGC 2024, which is in the agreement with results presented in other studies (such as, Planck Collaboration 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).

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.

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 ∼ 20−40 µG on the west side, and ∼ 82 µG on the east side.The good agreement between the CN and HCO + B POS strengths in these regions, particularly at the eastern side, could indicate that these two molecular lines are tracing similar gas.Nevertheless, the factor of almost 2 difference in the B POS measured in the west and the east could be due to different gas densities at these sides of the bubble.As the eastern side of the bubble is denser than the western side (see, Meyer et al. 2008), we expect that the gas here is less impacted by the incoming radiation.In addition, the opacity of the HCO + emission is notably higher (τ > 4) on the western side than on the eastern side (τ = 2).
As previously presented in 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 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 ∼ 750 µG and ∼ 250 µG respectively (Brogan & Troland 2001), and ∼ (1000−1700) µG using dust polarization data (Hoang et al. 2022).The magnetic field strength of the PDR in the Horsehead nebula (SMM1, Hwang et al. 2023) is a few tens of µG, and comparable to our results in NGC 2024.That study used C 18 O to derive the line width and dust column density and effective radius to estimate the gas volume density.Although the edges of the bubble in our work and the Horsehead nebula have comparable densities (a few × 10 3 cm −3 ), measured line widths are different.It is worth pointing out that C 18 O and CN and HCO + are tracing different gas (Philipp et al. 2006;Bron et al. 2018).C 18 O traces more compact structures than the CN and HCO + , and it gets destructed by the stellar feedback, contrary to CN and HCO + , whose emission becomes enhanced in these regions (Bron et al. 2018).In addition, Hwang et al. (2023) used a modified DCF method that measures the ratio of the ordered and turbulent component of the B-field (Hildebrand et al. 2009) to derive the magnetic field strength.
Other PDR regions impacted by super stellar clusters have stronger magnetic field to those reported in our work.For example, Pattle et al. (2018) reported strong magnetic fields of a few hundred µG in the Pillars of M 16.This work used the DCF method to derive the magnetic field strength.Gas density of the Pillars are somewhat higher than those in our work 5 • 10 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 POS .We show how B POS varies with the method we selected in Fig. D.4.

Filament
The filament going across NGC 2024 is super-critical (Orkisz et al. 2019), which means it is gravitationally unstable and a potential site for star formation.We measure the lowest values of the POS magnetic field (∼ 30 − 50 µG) in this region, suggesting that the magnetic field cannot dominate gas dynamics.We return to this point in 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 C 18 O emission to measure the turbulence and number density and different method to derive the B POS .Similarly, Ching et al. (2018) investigated magnetic field toward the DRL1 filament and found B POS = 600 µG.

Magnetic field support in NGC 2024
Next, we investigate the role of the magnetic field at the edges of the bubble, in the filament, and the overlap region (bubble and Article number, page 13 of 36 Table 4.Estimated mass-to-flux ratio µ Φ , Alfvénic Mach number (M) A and plasma-beta β using B 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.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 M Φ and compare it to the critical mass-to-flux ratio.The critical mass-to-flux ratio M Φ crit depends on the assumed geometry.For instance, for a uniform disk, M Nakano & Nakamura 1978;Joos et al. 2012;Hanawa et al. 2019).In addition, for a spherical geometry, the critical mass-toflux ratio will be expressed as & Spitzer 1976).The µ Φ parameter is the ratio between estimated mass-to-flux and M Φ crit and can be computed as: 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 M Φ crit and the above equation: Therefore, the µ Φ will vary by a factor of 0.67 that comes from the assumed geometry.It is necessary to point out that there is not a critical mass-to-flux ratio for lateral contraction of a filament when it threaded by magnetic field parallel to its axis of symmetry (Mouschovias & Morton 1992).In our work, we use Eq.22 and assume the cylindrical geometry, but we note that this factor can vary depending on the assumed geometry of a system.Moreover, in case of the filament, the interpretation and physical meaning of mass-to-flux ratio in not straightforward, and it should be threated with causion.
Next, we compute the Alfvénic Mach number, M A : where σ υ,NT is the non-thermal velocity dispersion reported in Tab. 2. The √ 3 factor comes from the assumption of isotropic turbulence (Crutcher 2004;Stewart & Federrath 2022).υ A is the Alfvén velocity defined as where B is the total magnetic field strength in units of G, and ρ 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, β B , that gives information about the ratio of thermal and magnetic pressure: 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 µ Φ < 1 in each region.This result could imply the presence of magnetically supported gas (Pattle et al. 2019).In general, it is worth noting that we see a difference of one and two orders of magnitude in the µ Φ between the edges of the bubble and the filament and the overlap region.
We find M A greater than 1 in all regions, which suggests that the magnetic field does not govern the gas motions.Moreover, we measure a plasma-beta parameter lower than 1 in almost all environments, which agrees with the lower values of µ Φ , particularly at the edges of the bubble.A plasma beta derived from CN measurements of magnetic field at the western side of the bubble is slightly higher than one.Low values of plasma-beta suggest that the magnetic energy dominates the thermal motions of the gas.
The results reported in 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 the gravitational stability of the gas in these regions requires systematic analysis of the magnetic field.This also includes both POS and LOS components in these regions and constraints on the geometry of NGC 2024, which is beyond the scope of this study.
The properties of the gas we study in our work vary with the environment and this gas shows different characteristics.Our results thus highlight the non-negligible role of the magnetic field Article number, page 14 of 36 Bešlić: The magnetic field in the Flame nebula 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 µ Φ > 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.M 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.

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 combined these measurements with molecular data observations from the ORION-B Large Program on the IRAM 30-meter telescope.We summarize our findings obtained using these observations as follows: 1. Our results focus on a subset of environments found in NGC 2024, particularly on the shell near the edge of the 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. 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. We use HAWC+ Band D at 154 µm to characterize the geometry of the magnetic field across NGC 2024.We find that the structure of the magnetic field is ordered and follows the morphology of the expanding HII region and the direction of the filament. 4. Using the CN(1-0) and HCO + (1-0) molecular emission obtained using the IRAM 30-meter telescope, we characterize physical conditions (turbulence and gas density) in four specific regions in NGC 2024: the edges of the expanding 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. We derived the POS magnetic field strength using the ST (and classical DCF) method, with values ranging from 20 to 90 µG (50 -240 µG) in the environments mentioned above.
The strongest magnetic field is found in the region composed of the dusty filament and the edge of the expanding 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. Magnetic field strengths derived using CN(1-0) and HCO + (1-0) are comparable within the uncertainties at the edges of the bubble, especially in the eastern side of the bubble.The magnetic field strength calculated from HCO + (1-0) is generally higher than that inferred from CN(1-0).We note that B POS measured from CN(1-0) and HCO + (1-0) is larger on the eastern side of the bubble than in the western edge.This observed difference may result from changes in the magnetic field direction indicated in previous studies.In addition, the western side of the bubble is more impacted by stellar radiation due to its lower density.7. By analyzing the mass-to-flux ratio, Alfvénic Mach number, and plasma-beta parameter, we find that the edges of the bubble and the filament show different properties.Gas impacted by the stellar feedback and traced by the CN(1-0) and HCO + (1-0) emission seems better supported against gravitational collapse by the magnetic field than the gas in the filament, which represents a location where star formation can take place.However, we find the Alfvénic Mach number higher than one in all regions, which suggests that the magnetic field does not control the gas motions.Our results should be treated as upper limits given the use of the POS component of the B field.In addition, the estimated mass-toflux ratio has a large uncertainty due to the specific geometry of NGC 2024.We also note that these regions of NGC 2024 could be in the transition phase regarding gravitational stability.
Our research emphasizes the importance of utilizing dust polarization measurements to characterize the structure of the magnetic field and the need to combine these measurements with molecular data to accurately infer the magnetic field's strength.Our results demonstrate that the magnetic field plays a critical role and that its contribution and other factors cannot be ignored.We also demonstrate the significant impact that different methods used to derive the POS magnetic field can have on our results and how it affects the interpretation of the gravitational stability of the gas based on the corresponding parameters.Moreover, by removing the contribution of large-scale gradients of the magnetic field direction using techniques such as the sliding window, we can consider only the local magnetic field.In addition, we highlight the necessity of conducting a careful analysis of the Article number, page 15 of 36 line radiative transfer and using appropriate rate coefficients for inelastic collisions with H 2 and electrons, as well as line opacity corrections.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: The second parameter is the centroid velocity of the main hyperfine component (the F = 5/2 − 3/2 transition): The next parameter is the full width at half maximum (FWHM) of the hyperfine components: and the last parameter is the opacity of all components of the studied multiplet: (B.4) We assume that the optical depth of each component can be described using a Gaussian function of velocity: The opacity of the multiplet is then calculated as a sum of opacities of all hyperfine components: where τ i is the opacity of the i-th component.
Finally, the antenna temperature, T ant is thus derived as the ratio between the first and the fourth fitting parameters, p 1 , and p 4 : υ) .(B.7) To get the excitation temperature, we mask all pixels having the opacity higher than 10 and lower than 0.2 to avoid degeneracies.The excitation temperature is derived assuming the Local Thermodynamical Equilibrium (LTE) case: We fit the CN emission in CUBE that uses an optimized version of the minimization method taken from the MINUIT system of CERN.To do so, we load the data into CUBE and run the commands \fit/minimize and add \hfs to specify we want to fit the hyperfine structure.Therefore, we also provide a file containing information about the components of the hyperfine multiplet, shown in 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.B.1, B.2, B.3, and B.4 respectively.
Before fitting, we define a two-dimensional mask based on the S/N of the CN integrated intensity.We select the threshold of 10 and fit two CN components in the area that goes inside the mask and one component outside the mask.By adding this step, we provide CUBE with additional information about the area where the CN emission is bright enough to observe double components in its spectrum.CUBE requires the initial values of the free parameters for each component we want to fit.We specify our initial guess for the centroid velocities of the components to be 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. 1 .
Fig. 1.Multiwavelength image of NGC 2024.(Left panel) Color-composite image of the Flame nebula (NGC 2024) showing peak intensities of 12 CO (1 − 0) (blue) emission and isotopologues, 13 CO (1 − 0) (green) and C 18 O (1 − 0) (red) obtained by the IRAM 30-meter telescope (image credits: Pety et al. 2017).We overlaid the SOFIA HAWC+ field of view as a white dashed rectangle.Gray contours show a network of filaments presented in Orkisz et al. (2019).We label the regions we investigate in this work: Bubble East, Bubble West, Filament, and Bubble+Filament.(Right panel) ESO-VISTA image (ESO/J.Emerson/VISTA, Cambridge Astronomical Survey Unit).We sketched the environments seen across NGC 2024: the 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(Mezger et al. , 1992)).In addition, we labeled the edges of the expanding HII region and the direction of the stellar feedback driven by the radiation produced by recent star formation.
Fig. A.4 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. A.5 in App.A.2.While the rms is only a few degrees (∼ 2 deg) across a large part of NGC 2024, we observe higher values (> 10 deg) toward the center of NGC 2024.Nevertheless, we consider only Band D data in this work because of an overall good correspondence (measured as a low rms) between dust polarization angles measured from Band D and E across the regions we further analyze in this work.
Fig. A.6  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 (seeFig.A.6 in App.A.2).In addition, we compare the magnetic field direction with the 100 µm(Dotson et al. 2000) and submillimeter (850 µm) polarization measurements(Matthews et al.

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

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

Fig. 4 .
Fig. 4. Direction of the magnetic field (left panel), the mean angle of the magnetic field (middle panel), and the rms of the magnetic field direction (right panel).Maps in the middle and right panels have been made using a 12 × 12 pixels sliding window with weight inversely proportional to the measured uncertainty variance.All maps have been regridded to match the grid size of the CN and HCO + data.The gray contour indicates the mask we defined for dust polarization measurements, described in Sec.2.1.
and entropy modes (Lithwick & Goldreich 2001) also cause fluctuations in polarization angle.Therefore, to derive the magnetic field strength in NGC 2024, it is important to acknowledge the contribution of non-Alfvénic motions and the compressible nature of the ISM.In this work, we derive the B POS using the prescription presented in Skalidis & Tassis (2021); Skalidis et al. (2021):

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

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

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

Fig. C. 1 .
Fig. C.1.Parameter space used for the radiative transfer modeling of CN in the region at the western side of the bubble.The first row shows results for the excitation temperature.The color bar shows values of excitation temperature, T 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, τ, and the third row shows the peak temperature, T peak .The last row shows the χ 2 minimization of the modeled peak temperature and the observed value computed from the CN(1-0) spectra.

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

SOFIA Field of View SOFIA Field of View Alnitak Alnitak Stellar feedback Stellar feedback Dust filament Dust filament ESO VISTA ESO VISTA
.