# Appendix A: HAWC+ dust continuum and polarization

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

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

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

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

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

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

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

To compare dust polarization angles derived from Bands D and E observations, we first convolve Band D to match the lower spatial resolution of the Band E data. Next, we compute a difference between the two angles, $\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. A.5. The right panel in

<figure><img src="image_6.png" /></figure>

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

Fig. A.5 shows the corresponding rms of the map shown in the left panel.

## Appendix A.3: Comparison with the literature data

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

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

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

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

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

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

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

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

To compare dust polarization angles derived from Bands D and E observations, we first convolve Band D to match the lower spatial resolution of the Band E data. Next, we compute a difference between the two angles, $\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. A.5. The right panel in

<figure><img src="image_7.png" /></figure>

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

Fig. A.5 shows the corresponding rms of the map shown in the left panel.

## Appendix A.3: Comparison with the literature data

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

**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).
<figure><img src="image_8.png" /></figure>

**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).
<figure><img src="image_8.png" /><figcaption>**Fig. A.3.** Same as in Fig. A.2, but for the polarization angles (see Eq. 6 in Sec. 2.1).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. A.3.** Same as in Fig. A.2, but for the polarization angles (see Eq. 6 in Sec. 2.1).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. A.5.** Difference in polarization angles measured from Band D and E data is shown in the left panel. The right panel shows the rms of the angle differences.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. A.5.** Difference in polarization angles measured from Band D and E data is shown in the left panel. The right panel shows the rms of the angle differences.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

**Fig. A.6.** SOFIA HAWC+ dust continuum Band D and Band E measurements (left and right panels in the top row, respectively), but with overlayed positions of protostars and NIR polarization (dark purple circles and lines, Kandori et al. 2007). Dark blue contours represent the outflow observed in the HCO$^+$ emission. Light blue points indicate the positions of FIR sources. Black dashed rectangles show the central area and the zoom-in panels on the bottom row. Here, we show magnetic field lines from the FIR dust polarization at 100 µm (Dotson et al. 2000) in dark magenta and from (sub)millimeter dust polarization at 850 µm (Matthews et al. 2002) in black.
<figure><img src="image_12.png" /><figcaption>**Fig. A.6.** SOFIA HAWC+ dust continuum Band D and Band E measurements (left and right panels in the top row, respectively), but with overlayed positions of protostars and NIR polarization (dark purple circles and lines, Kandori et al. 2007). Dark blue contours represent the outflow observed in the HCO$^+$ emission. Light blue points indicate the positions of FIR sources. Black dashed rectangles show the central area and the zoom-in panels on the bottom row. Here, we show magnetic field lines from the FIR dust polarization at 100 µm (Dotson et al. 2000) in dark magenta and from (sub)millimeter dust polarization at 850 µm (Matthews et al. 2002) in black.</figcaption></figure>
## Appendix B: Spectral line fitting

<figure><img src="image_12.png" /><figcaption>**Fig. B.1.** Fitting parameters for the hyperfine CN emission in NGC 2024 measured for both velocity components (top and bottom row each). The first parameter, $p_1$ is shown in the left panel, followed by the centroid velocity, FWHM and the final parameter, the opacity ($p_4$).</figcaption></figure>

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

## Appendix B.1: Hyperfine structure of the CN line

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

Four parameters describe the hyperfine structure model: $p_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_{\text{ant}} \cdot \tau \qquad (B.1)$$

The second parameter is the centroid velocity of the main hyperfine component (the $F = 5/2 - 3/2$ transition):

$$p_2 = v_{0,F=5/2-3/2} \qquad (B.2)$$

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

$$p_3 = \text{FWHM}, \qquad (B.3)$$

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

$$p_4 = \tau. \qquad (B.4)$$

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

$$\tau_i(\nu) = \tau_i \cdot \exp \left[ -4 \ln 2 \left( \frac{\nu - \nu_{0,i}}{p_3} \right)^2 \right]. \qquad (B.5)$$

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

$$\tau = \sum_{i=1}^{N} \tau_i \qquad (B.6)$$

where $\tau_i$ is the opacity of the i-th component.
## Appendix B: Spectral line fitting

<figure><img src="image_13.png" /><figcaption>**Fig. B.1.** Fitting parameters for the hyperfine CN emission in NGC 2024 measured for both velocity components (top and bottom row each). The first parameter, $p_1$ is shown in the left panel, followed by the centroid velocity, FWHM and the final parameter, the opacity ($p_4$).</figcaption></figure>

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

## Appendix B.1: Hyperfine structure of the CN line

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

Four parameters describe the hyperfine structure model: $p_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_{\text{ant}} \cdot \tau \qquad (B.1)$$

The second parameter is the centroid velocity of the main hyperfine component (the $F = 5/2 - 3/2$ transition):

$$p_2 = v_{0,F=5/2-3/2} \qquad (B.2)$$

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

$$p_3 = \text{FWHM}, \qquad (B.3)$$

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

$$p_4 = \tau. \qquad (B.4)$$

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

$$\tau_i(\nu) = \tau_i \cdot \exp\left[-4 \ln 2 \left(\frac{\nu - \nu_{0,i}}{p_3}\right)^2\right]. \qquad (B.5)$$

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

$$\tau = \sum_{i=1}^{N} \tau_i \qquad (B.6)$$

where $\tau_i$ is the opacity of the i-th component.
<figure><img src="image_13.png" /><figcaption>Fig. B.2. Fitting parameters for the HCO+ emission assuming a Gaussian line profile. Each row shows parameters for each fitted component. The results in the first row represent the brightest HCO+ component. The second row shows the fainter component and the outflow observed close to the center of NGC 2024.</figcaption></figure>

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$:

$$T_{ant} = \frac{p_1}{p_4} (1 - e^{-\tau(v)}). \qquad (B.7)$$

To get the excitation temperature, we mask all pixels having the opacity higher than 10 and lower than 0.2 to avoid degeneracies. The excitation temperature is derived assuming the Local Thermodynamical Equilibrium (LTE) case:

$$T_{ex} = \frac{h\nu}{k} \left\{ \ln \left[ 1 + \frac{h\nu}{kT_{ant}} (1 - e^{-\tau_{v}}) \right] \right\}^{-1}. \qquad (B.8)$$

We fit the CN emission in CUBE that uses an optimized version of the minimization method taken from the MINUIT system of CERN. To do so, we load the data into CUBE and run the commands `\fit/minimize` and add `\hfs` to specify we want to fit

the hyperfine structure. Therefore, we also provide a file containing information about the components of the hyperfine multiplet, shown in Table 1. In the case of the hyperfine fitting, CUBE fits four parameters, $p_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. B.1.
<figure><img src="image_14.png" /><figcaption>Fig. B.2. Fitting parameters for the HCO+ emission assuming a Gaussian line profile. Each row shows parameters for each fitted component. The results in the first row represent the brightest HCO+ component. The second row shows the fainter component and the outflow observed close to the center of NGC 2024.</figcaption></figure>

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$:

$$T_{ant} = \frac{p_1}{p_4} (1 - e^{-\tau(v)}). \qquad (B.7)$$

To get the excitation temperature, we mask all pixels having the opacity higher than 10 and lower than 0.2 to avoid degeneracies. The excitation temperature is derived assuming the Local Thermodynamical Equilibrium (LTE) case:

$$T_{ex} = \frac{h\nu}{k} \left\{ \ln \left[ 1 + \frac{h\nu}{kT_{ant}} (1 - e^{-\tau_{v}}) \right] \right\}^{-1}. \qquad (B.8)$$

We fit the CN emission in CUBE that uses an optimized version of the minimization method taken from the MINUIT system of CERN. To do so, we load the data into CUBE and run the commands `\fit/minimize` and add `\hfs` to specify we want to fit

the hyperfine structure. Therefore, we also provide a file containing information about the components of the hyperfine multiplet, shown in Table 1. In the case of the hyperfine fitting, CUBE fits four parameters, $p_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. B.1.
## Appendix B.2: Gaussian lines of the HCO$^+$ emission

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

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

We fit a Gaussian line profile to the HCO$^+$ emission in SCOUSEPY. The output parameters describing a Gaussian function are each component's peak temperature, centroid velocity, and FWHM. Additionally, SCOUSEPY computes the rms, S/N and the residuals. Based on these results we derive from SCOUSEPY, we define two two-dimensional masks. The first mask contains

a region where SCOUSEPY identifies three velocity components. The second mask is a region within which SCOUSEPY finds two velocity components. We use these masks as the input to CUBE to mark regions where we observe three (first mask), two (second mask), and one velocity component (outside these two masks).

Similarly, as for the CN emission, we use the command \fit/minimize in CUBE and add \gaussian to specify the shape of the spectral line. In CUBE, the Gaussian line is described using three parameters: the area (in 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. B.2.

## Appendix C: Radiative transfer modeling of CN and HCO$^+$ emission

Here, we show results from the radiative transfer modeling of CN and HCO$^+$ excitation, presented in Sec. 4.3. For each region we study in this work (edges of the bubble, filament, and the overlap region), we show how our output parameters from RADEX (excitation temperature, opacity, and the peak temperature) vary as a function of input parameters (line width, column and volume density) in Fig. C.1, C.2, C.3 and C.4 for CN emission. In the case of HCO$^+$ emission, we show these results in Fig. C.5, C.6, C.7 and C.8. We also include a $\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. C.9.
## Appendix B.2: Gaussian lines of the HCO$^+$ emission

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

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

We fit a Gaussian line profile to the HCO$^+$ emission in SCOUSEPY. The output parameters describing a Gaussian function are each component's peak temperature, centroid velocity, and FWHM. Additionally, SCOUSEPY computes the rms, S/N and the residuals. Based on these results we derive from SCOUSEPY, we define two two-dimensional masks. The first mask contains

a region where SCOUSEPY identifies three velocity components. The second mask is a region within which SCOUSEPY finds two velocity components. We use these masks as the input to CUBE to mark regions where we observe three (first mask), two (second mask), and one velocity component (outside these two masks).

Similarly, as for the CN emission, we use the command \fit/minimize in CUBE and add \gaussian to specify the shape of the spectral line. In CUBE, the Gaussian line is described using three parameters: the area (in 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. B.2.

## Appendix C: Radiative transfer modeling of CN and HCO$^+$ emission

Here, we show results from the radiative transfer modeling of CN and HCO$^+$ excitation, presented in Sec. 4.3. For each region we study in this work (edges of the bubble, filament, and the overlap region), we show how our output parameters from RADEX (excitation temperature, opacity, and the peak temperature) vary as a function of input parameters (line width, column and volume density) in Fig. C.1, C.2, C.3 and C.4 for CN emission. In the case of HCO$^+$ emission, we show these results in Fig. C.5, C.6, C.7 and C.8. We also include a $\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. C.9.
<figure><img src="image_14.png" /><figcaption>**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*<sub>ex</sub> 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*<sub>peak</sub>. The last row shows the χ<sup>2</sup> minimization of the modeled peak temperature and the observed value computed from the CN(1–0) spectra.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**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*<sub>ex</sub> 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*<sub>peak</sub>. The last row shows the χ<sup>2</sup> minimization of the modeled peak temperature and the observed value computed from the CN(1–0) spectra.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. C.2.** Same as in Fig. C.1, but for the eastern part of the bubble.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. C.2. Same as in Fig. C.1, but for the eastern part of the bubble.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. C.3. Same as in Fig. C.1, but for the filament.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. C.3. Same as in Fig. C.1, but for the filament.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. C.4.** Same as in Fig. C.1, but for the overlap region.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. C.4.** Same as in Fig. C.1, but for the overlap region.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. C.5. Same as in Fig. C.1, but for HCO⁺(1–0).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. C.5. Same as in Fig. C.1, but for HCO⁺(1–0).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. C.6.** Same as in Fig. C.5, but for the eastern side of the bubble.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. C.6.** Same as in Fig. C.5, but for the eastern side of the bubble.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. C.7. Same as in Fig. C.5, but for the filament.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. C.7. Same as in Fig. C.5, but for the filament.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. C.8.** Same as in Fig. C.5, but for the overlap region.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. C.8.** Same as in Fig. C.5, but for the overlap region.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**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 inferred from the non-LTE RADEX modeling for physical conditions presented in Sec. C and shown in Table 3. Dashed lines show the corresponding residuals.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**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 inferred from the non-LTE RADEX modeling for physical conditions presented in Sec. C and shown in Table 3. Dashed lines show the corresponding residuals.</figcaption></figure>
# Appendix D: Measuring magnetic field strength in NGC 2024

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

## Appendix D.1: Sliding window

We investigate the change of the angle dispersion measured using the sliding window (Sec. 5.3) by varying its size. We used the size of the sliding window from 1 to 3.5 beam sizes, with

a step of 0.5, which corresponds to sizes from 4 × 4 to 14 × 14 pixels, respectively. We show these results in Fig. D.1. The beam-sized sliding window is too small to remove any large-scale contribution of the magnetic field. We find that the rms of the angle changes by a factor of around 3 as we increase the size of the sliding window.

## Appendix D.2: Histogram analysis

In this section, we describe an alternative approach to compute $\hat{\sigma}_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

<figure><img src="image_23.png" /><figcaption>**Fig. D.1.** Angle rms as a function of the size of the sliding window.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. D.2.** Top left: Bubble on the west. Top right: Bubble on the eastern side of NGC 2024. Bottom left: Filament. Bottom right: Mixture of bubble and dusty filament.</figcaption></figure>
# Appendix D: Measuring magnetic field strength in NGC 2024

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

## Appendix D.1: Sliding window

We investigate the change of the angle dispersion measured using the sliding window (Sec. 5.3) by varying its size. We used the size of the sliding window from 1 to 3.5 beam sizes, with

a step of 0.5, which corresponds to sizes from 4 × 4 to 14 × 14 pixels, respectively. We show these results in Fig. D.1. The beam-sized sliding window is too small to remove any large-scale contribution of the magnetic field. We find that the rms of the angle changes by a factor of around 3 as we increase the size of the sliding window.

## Appendix D.2: Histogram analysis

In this section, we describe an alternative approach to compute $\hat{\sigma}_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

<figure><img src="image_24.png" /><figcaption>**Fig. D.1.** Angle rms as a function of the size of the sliding window.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. D.2.** Top left: Bubble on the west. Top right: Bubble on the eastern side of NGC 2024. Bottom left: Filament. Bottom right: Mixture of bubble and dusty filament.</figcaption></figure>
<table><caption>**Table D.1.** Results on the magnetic field strength from histogram analysis.</caption><thead><tr><th>Region</th><th>σ̂<sub>c</sub>(φ) [deg]</th><th>B<sub>DCF,CN</sub> [μG]</th><th>B<sub>DCF,HCO+</sub> [μG]</th><th>B<sub>ST,CN</sub> [μG]</th><th>B<sub>ST,HCO+</sub> [μG]</th></tr></thead><tbody><tr><td>Bubble West</td><td>4.7 ± 0.2</td><td>53 ± 27</td><td>106 ± 30</td><td>22 ± 11</td><td>43 ± 12</td></tr><tr><td>Bubble East</td><td>2.56 ± 0.04</td><td>317 ± 72</td><td>318 ± 72</td><td>95 ± 22</td><td>95 ± 22</td></tr><tr><td>Filament</td><td>4.4 ± 0.2</td><td>91 ± 38</td><td>141 ± 48</td><td>36 ± 15</td><td>55 ± 19</td></tr><tr><td>Filament+Bubble</td><td>6.5 ± 0.5</td><td>158 ± 44</td><td>180 ± 47</td><td>75 ± 21</td><td>86 ± 22</td></tr></tbody></table>

<figure><img src="image_25.png" /><figcaption>**Fig. D.3.** Comparison of the angle dispersion using two different approaches, sliding window (y-axis, see Sec. 5) and histogram analysis (x-axis).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. D.4.** Magnetic field strength derived using different methods (Sec. 5.1, D.3) for computing the magnetic field strength and the angle dispersion (Sec. 5.3, D.2) as a function of the environment we analyze in this work.</figcaption></figure>

histograms for each beam-averaged area we used in this work, including the fitting parameters of each distribution in Fig. D.2. In most cases, the histogram corresponds to a single Gaussian. In Fig. D.3, we show comparison between $\hat{\sigma}_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.

## Appendix D.3: $B_{POS}$ using different approaches

We compare results on POS B-field strengths using several methods in Fig. D.4. As described in Sec. 5, results derived from the classical DCF are higher than those derived using ST, even when using the histogram analysis instead of the sliding window approach.
<table><caption>**Table D.1.** Results on the magnetic field strength from histogram analysis.</caption><thead><tr><th>Region</th><th>σ̂<sub>c</sub>(φ) [deg]</th><th>B<sub>DCF,CN</sub> [μG]</th><th>B<sub>DCF,HCO+</sub> [μG]</th><th>B<sub>ST,CN</sub> [μG]</th><th>B<sub>ST,HCO+</sub> [μG]</th></tr></thead><tbody><tr><td>Bubble West</td><td>4.7 ± 0.2</td><td>53 ± 27</td><td>106 ± 30</td><td>22 ± 11</td><td>43 ± 12</td></tr><tr><td>Bubble East</td><td>2.56 ± 0.04</td><td>317 ± 72</td><td>318 ± 72</td><td>95 ± 22</td><td>95 ± 22</td></tr><tr><td>Filament</td><td>4.4 ± 0.2</td><td>91 ± 38</td><td>141 ± 48</td><td>36 ± 15</td><td>55 ± 19</td></tr><tr><td>Filament+Bubble</td><td>6.5 ± 0.5</td><td>158 ± 44</td><td>180 ± 47</td><td>75 ± 21</td><td>86 ± 22</td></tr></tbody></table>

<figure><img src="image_26.png" /><figcaption>**Fig. D.3.** Comparison of the angle dispersion using two different approaches, sliding window (y-axis, see Sec. 5) and histogram analysis (x-axis).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. D.4.** Magnetic field strength derived using different methods (Sec. 5.1, D.3) for computing the magnetic field strength and the angle dispersion (Sec. 5.3, D.2) as a function of the environment we analyze in this work.</figcaption></figure>

histograms for each beam-averaged area we used in this work, including the fitting parameters of each distribution in Fig. D.2. In most cases, the histogram corresponds to a single Gaussian. In Fig. D.3, we show comparison between $\hat{\sigma}_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.

## Appendix D.3: $B_{POS}$ using different approaches

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