© The Authors 2026

1 Introduction

Magnetic fields are thought to influence the process of star formation by regulating the rate at which dense interstellar material collapses under its own gravity (see, e.g., Krumholz & Federrath 2019 and Pattle et al. 2022 for reviews). One method that has proven useful in inferring the morphology and strength of magnetic fields in star-forming regions is based on polarization measurements. As first pointed out by Hall (1949) and Hiltner (1949), the optical to near-infrared polarization of stellar radiation arises from magnetically aligned interstellar dust grains. The mechanisms responsible for the alignment continue to be an active topic of research today. Early explanations, such as alignment through paramagnetic dissipation (Davis & Greenstein 1951), have struggled to hold up against observational testing. Today, the most promising candidate for a theory of dust grain alignment is radiative torque (RAT) alignment (Dolginov & Mytrophanov 1976; Draine & Weingartner 1996; Lazarían & Hoang 2007; Andersson et al. 2015). Here, differences in the scattering efficiency of left- and right-circularly polarized light of an anisotropic radiation field cause dust grains with net helicity to spin up to suprathermal rates. If the material of the dust grain is paramagnetic, this rotation then leads to the magnetization of the grain through the Barnett effect (Barnett 1915) and consequently causes the grain to align preferentially with its short axis parallel to the magnetic field lines. At far-infrared to submillimeter wavelengths, this preferred orientation causes a net polarization of the thermally emitted light that is perpendicular to the plane-of-sky (POS) component of the magnetic field.

To accurately infer the strength and morphology of the underlying magnetic field, it is therefore necessary to first constrain the dust properties, as well as the efficiency of grain alignment across different environments.

In this work, we study the far-infrared polarization spectrum (Hildebrand et al. 1999) - namely, variations in the degree of linear polarization with wavelength. In particular, we study the molecular cloud core L43, as it contains a dense starless core in close proximity to a young stellar object (YSO), acting as an anisotropic radiation source. Previous work (e.g., Santos et al. 2019; Michail et al. 2021; Seifried et al. 2023; Tram et al. 2024) indicates that changes in the polarization spectrum are driven by temperature differences between dust species that align with the magnetic field (e.g., silicate grains) and those that are thought not to align (e.g., carbonaceous grains), as well as differences in the local temperatures along the line of sight (LOS). The latter effect is also termed the “heterogeneous cloud effect” (HCE; Michail et al. 2021). As local temperature differences are inherently coupled to the strength of the radiation field that heats the grains, these differences cause the efficiency of the RAT mechanism to vary as well. As polarization measurements are constrained to the POS component of the magnetic field, changes in the orientation of the field lines influence the total degree of polarization. Therefore, by measuring the polarization of thermal reemission radiation at multiple wavelengths, with each preferentially tracing grains at different temperatures and depths, variable orientations of the magnetic field along the LOS also influence the shape of the polarization spectrum. We have attempted to separate the effect of magnetic field orientation from HCE by comparing observations with corresponding simulations of a simple molecular cloud core, modeled with the 3D Monte Carlo radiative transfer code POLARIS (Reissl et al. 2016).

This paper is structured as follows: in Sect. 2, we give a brief description of the source, L43. Sect. 3 describes the observations, as well as the data reduction and preparation techniques used. In Sect. 4, we analyze the polarization data and explain how we derived temperature and column density maps, which we used to calculate magnetic field strengths. In Sect. 5, we explain how we employed the radiative transfer code POLARIS to construct a simple cloud core model for further analysis of the spatially resolved polarization spectrum. Finally, we give a brief discussion of our results in Sect. 6 and summarize our conclusions in Sect. 7.

2 Source description

L43 is a dark filamentary cloud with a Lynds opacity class of six (the highest opacity class in the Lynds catalog Lynds 1962). It contains a dense starless core to the east and a partially embedded Class I protostar (classified by Chen et al. 2009), commonly referred to as RNO 91 (“red nebulous object”; Cohen 1980), to the west. The cloud lies in Ophiuchus North (Oph N), the northern part of the Ophiuchus star-forming region (see, e.g., Hatchell et al. 2012; Alves et al. 2025 for an overview of the region). We assume that the distance to L43 is not significantly different from the rest of the cloud complex, for which de Geus et al. (1989) find d = 125 ± 25 pc.

The protostar RNO 91 is associated with a bipolar CO outflow (Bence et al. 1998; Lee et al. 2002; Lee & Ho 2005), seen in ¹²CO J = 1-0, J = 2-1 and J = 3-2 spectral line data. By comparing the ¹²CO data with I-band continuum data, Mathieu et al. (1988) pointed out that the southern lobe of the outflow has influenced the dust morphology of the cloud by clearing a U-shaped cavity around RNO 91. This influence was also noticed by Ward-Thompson et al. (2000), who studied the polarized flux of L43 with the SCUBA/SCUPOL polarimeter at the James Clerk Maxwell Telescope (JCMT) and noted that the outflow has likely influenced the morphology of the underlying magnetic field. More recent observations with its successor SCUBA-2/POL-2 by Karoly et al. (2023) support this as well. This possible interaction between the outflow and the magnetic field complicates the determination of magnetic field strengths through polarization measurements from thermal dust emission in this region. Nonetheless, Crutcher et al. (2004) apply the Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951; Chandrasekhar & Fermi 1953) to the densest region of L43, assuming this region is still unaffected by the outflow. They find POS field strengths of B_POS ≈ 160 μG for an area with an average volume density of n_H2 ≈ 3.8 × 10⁵cm⁻³. Karoly et al. (2023) divide L43 in three different subregions and find POS field strengths between B_POS ≈ (40-257) μG and volume densities of n_H2 ≈ (0.6-4.9) × 10⁵ cm⁻³.

3 Observations and data reduction

A summary of all the instruments used, along with their measuring wavelengths, pixel, and beam sizes, is given in Table 1.

3.1 SOFIA/HAWC+

SOFIA/HAWC+ (Harper et al. 2018) D-band observations (154 μm) of L43 were carried out on 11 June 2022 as part of observation cycle 9 (Proposal: 09_0159, Obs. ID: P_2022-06-11_HA_F886-09_0159_4-044), utilizing the on-the-fly mapping (OTFMAP) mode¹. The data were processed using the “HAWC DRP v3.0.0” pipeline² (see, e.g., Santos et al. 2019), resulting in “Level 4” (science quality) data. These include images of the total intensity (Stokes I), the linearly polarized intensities (Stokes Q and Stokes U), as well as the degree of polarization p′, the debiased degree of polarization p (see Sect. 3.5), the polarization angle θ, and their corresponding errors.

3.2 SCUBA-2/POL-2

L43 was observed at 450 μm and 850 μm simultaneously with the Submillimetre Common-User Bolometer Array 2 (SCUBA-2; Holland et al. 2013) and its polarimeter POL-2 on the James Clerk Maxwell Telescope (JCMT) as part of the B fields In STar Forming Regions Observations-3 (BISTRO-3) large survey program (Ward-Thompson et al. 2017, Project ID: M20AL018). Twenty-seven repeated observations with an integration time of 30 minutes each were carried out in the time frame between February 2020 and March 2021. One observation is flagged as questionable and omitted from the data reduction process for both wavelengths. Additionally, data from one observation at 450 μm are incomplete and were also omitted. For the data reduction, we used the Sub-Millimetre User Reduction Facility (SMURF³; Chapin et al. 2013), which is part of the Starlink software package (Currie et al. 2014). SMURF contains an iterative reduction routine for POL-2 data, called pol2map, which is used to clear the raw bolometer time streams from atmospheric and instrumental noise and create separate maps of Stokes I, Q, U, and their corresponding errors. We closely follow the reduction steps laid out by Karoly et al. (2023), resulting in one reduction with a pixel size of 4″ that is used in the calculation of column density and temperature maps (see Sect. 4.4) and a second reduction with a pixel size of 8″ that is used to create a polarization vector catalog, which itself is binned to 12″.

The reduced flux data from SCUBA-2, given in units of pW, are converted to Jy beam⁻¹ by applying a flux conversion factor (FCF), calculated from calibration sources with well-known fluxes. We apply FCFs of (472 ± 76) Jy beam⁻¹ pW⁻¹ and (495 ± 32) Jy beam⁻¹ pW⁻¹ multiplied by additional factors of 1.96 and 1.35 for 450 μm and 850 μm, respectively (Mairs et al. 2021). The latter factors are additional correction factors applied to cover flux losses induced by POL-2.

The ¹²CO (3-2) molecular line emission is known to contaminate the 850 μm data obtained with SCUBA-2 (Drabek et al. 2012), leading to an overestimation of the continuum flux. To correct for this overestimation, we follow the instructions on the East Asian Observatory website⁴, by introducing ¹²CO (3-2) line emission data from the HARP/ACSIS instrument (Buckle et al. 2009) (Proposal ID: M07AU11) as a negative fakemap⁵ to the reduction with pol2map. The percentage flux originating from CO contamination amounts to 10%-15% to the west of RNO 91 and less than 10% to the east.

3.3 Herschel PACS and SPIRE

In order to calculate column density and temperature maps (see Sect. 4.4), we obtained data in the 160 μm band of the Photometer Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) and in the 250 μm, 350 μm and 500 μm band of the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) Instruments on board the Herschel Space Observatory (Pilbratt et al. 2010), which are publicly available at the Herschel Science Archive⁶. The Herschel data (Obs. ID: OT1_jhatchel_1) were acquired in the time frame between 17 February 2013 and 17 March 2013, utilizing the PACS/SPIRE Parallel Mode. Observations made in the PACS/SPIRE Parallel Mode are processed through the SPIRE Photometer Large Map pipeline and stored as a FITS file. In contrast to the reduction of the PACS/SPIRE data, the SCUBA-2 reduction pipeline needs to account for atmospheric noise. This leads to differences in the amount of recovered large-scale emission between the instruments. To account for this discrepancy, we follow the method described by Sadavoy et al. (2013), where the Herschel flux maps are inserted into the SCUBA-2 data reduction pipeline as a fakemap. The resulting flux maps are sensitive to the same emission as the SCUBA-2 data. We adopt flux calibration uncertainties of 20% for the PACS data and 10% for the SPIRE data (Chuss et al. 2019).

3.4 Nobeyama Radio Observatory/FOREST

The calculation of magnetic field strengths (see Sect. 4.5) requires information about the turbulent gas kinematics of L43. For this purpose, we made use of four N₂H⁺ (1-0) line emission observations taken with the superconductor-insulator-superconductor receiver FOREST (Minamidani et al. 2016) at the 45-m telescope operated by the Nobeyama Radio Observatory (NRO). The observations (Obs. ID’s: 20180504020849, 20180505015815, 20180506020433, 20190401033107) were taken between 4 May 2018 and 1 April 2019 and are publicly available at the Nobeyama-45m/ASTE Science Data Archive⁷. The data is stored in position-position-velocity cubes with a spectral resolution of 0.098 km s⁻¹. We combined the observations utilizing the reprojection package of Spectral-Cube⁸. The resulting spectra have an average rms noise of ~0.5 K, which allows us to only consider data exceeding an S/N-threshold of 3.

3.5 Data preparation

For an accurate comparison of the flux and polarization maps obtained with the different instruments, their pixel- and beam sizes are standardized. For the polarization and spectral line data, the command WCSALIGN⁹ from Starlink’s KAPPA package is applied to the SOFIA/HAWC+ and NRO data. This command takes the SCUBA-2/POL-2 observations as a reference, aligns the World Coordinate System (WCS) and adjusts the pixel size of the data to match the pixel size of the vector catalog (12″). The pixel sizes used for the column density and temperature fits (see Sect. 4.4) are standardized to 4″, by applying the adaptive anti-aliasing algorithm reproject_adaptive¹⁰ to the Herschel PACS and SPIRE data.

To standardize the beam sizes, we follow Chuss et al. (2019) and Michail et al. (2021), who smooth the images with a smaller beam size with a Gaussian kernel. The FWHM of these smoothing kernels are determined via $$\mathrm{F}\mathrm{W}\mathrm{H}{\mathrm{M}}_{\mathrm{s}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}}=\sqrt{\mathrm{F}\mathrm{W}\mathrm{H}{\mathrm{M}}_{\mathrm{r}\mathrm{e}\mathrm{f}}^2-\mathrm{F}\mathrm{W}\mathrm{H}{\mathrm{M}}_{\lambda }^2},$$(1)

where FWHM_ref describes the FWHM of the reference beam (the largest FWHM in the dataset) and FWHM_λ the original FWHM of the instrument beam at wavelength λ.

We apply this correction to the intensity data and the polarization data independently. The intensity data are smoothed to a resolution of 36.6″, corresponding to the Herschel/SPIRE beam size at 500 μm, while the polarization data are smoothed to a resolution of 14.6″, corresponding to the effective beam size of SCUBA-2/POL-2 at 850 μm (see Table 1).

Unlike the SOFIA/HAWC+ data, the degree of linear polarization p′ and angle θ are not precalculated for the SCUBA-2/POL-2 data. These values are determined from the Stokes I, Q and U data and calculated from the equations $$p^{\mathrm{\prime }}=\frac{1}{I}\sqrt{Q^2+U^2},\theta =\frac{12}{}{\tan }^{-1}\left(\frac{U}{Q}\right).$$(2)

Since the degree of polarization is a positive-definite quantity (see, e.g., Vaillancourt 2006), we attempt to correct for this bias by following Wardle & Kronberg (1974) in introducing the debiased polarization fraction p, calculated from the Stokes parameters I, Q and U and their respective uncertainties δQ and δU $$p=\sqrt{p^{\mathrm{\prime }}-\delta p^{\mathrm{\prime }}}=\frac{1}{I}\sqrt{Q^2+U^2-\frac{Q^2\delta Q^2+U^2\delta U^2}{Q^2+U^2}},$$(3)

where δp′ is the corresponding uncertainty of p′. Propagating the uncertainty of p yields $$\delta p=\sqrt{\frac{Q^2\delta Q^2+U^2\delta U^2}{I^2(Q^2+U^2)}+\frac{\delta I^2(Q^2+U^2)}{I^4}},$$(4)

where δI describes the uncertainty for the total intensity I. For the uncertainty on θ follows $$\delta \theta =\frac{12}{}\frac{\sqrt{Q^2\delta U^2+U^2\delta Q^2}}{Q^2+U^2}\times \frac{{180}^{\circ }}{\pi }.$$(5)

The following criteria are applied in the selection of the polarization vectors: $$\frac{I}{\delta I}>15,\frac{p}{\delta p}>2.$$(6)

Fig. 1 Dust continuum maps of the L43 region. Black lines indicate the polarization vectors. The vectors are rotated by 90° to indicate the magnetic field direction and are binned to the same 12″ grid. Black contour lines mark regions where the flux has doubled. White dashed ellipses show the two regions of interest that are referred to throughout this work. The effective beam size of each instrument is indicated with a red circle. Left: 154 μm continuum emission (SOFIA/HAWC+). Top right: 450 μm continuum emission (SCUBA-2/POL-2). Bottom right: 850 μm continuum emission (SCUBA-2/POL-2).

4 Observational results

4.1 Emission structure

Fig. 1 shows the dust continuum emission at 154 μm (left), 450 μm (top right), and 850 μm (bottom right), along with the rotated polarization vectors that fall above the S/N-cutoff defined in Eq. (6). The impact of the irradiation by the YSO RNO 91 is clearly seen at shorter wavelengths, whereas for the 850 μm observation, the flux decreases to a similar peak level as in L43E. These differences in the emission structure are likely to be attributed to variations in the local dust temperature. Note that the SOFIA/HAWC+ observation does not recover the extended emission seen in the northeastern part of the SCUBA-2/POL-2 observations in the L43E region.

4.2 Polarization properties

4.2.1 Magnetic field morphology

The magnetic field structure inferred from the polarization vectors in Fig. 1 shows a clear wavelength dependence. As Karoly et al. (2023) points out, the inferred magnetic field direction in the southern region of L43E coincides well with the cavity walls created by the molecular outflow associated with the YSO, suggesting that the magnetic field has been influenced by the outflow. This influence is also seen in the 450 μm emission, but is absent in the 154 μm emission. Fig. A.1 shows the distribution of the degree of polarization, as well as the distribution of the rotated polarization angles. The estimated outflow direction of 150° ± 10° (Karoly et al. 2023) coincides well with the circular mean of the 450 μm and 850 μm data in L43E, which are 151° and 153°, respectively. The stark difference in magnetic field morphology as well as emission structure between the SOFIA/HAWC+ data and the SCUBA-2/POL-2 data suggests that we might be observing different layers of the cloud, caused by differences in temperatures along the LOS. While the 154 μm emission traces the warmer, diffuse regions of L43, the 450 μm and 850 μm emission trace the inner, colder regions.

4.2.2 The polarization hole in L43

Previous analyses of near-infrared polarization properties in molecular clouds have shown that there exists a correlation between polarization efficiency, denoted as the degree of polarization due to extinction, per optical depth p_ext,λ/τ_λ, with visual extinction A_V (Whittet et al. 2008). This correlation, often referred to as the term “polarization hole”, is usually parameterized in the form of a power law $$p_{\mathrm{e}\mathrm{x}\mathrm{t},\lambda }/{\tau }_{\lambda }=a_1{A}_V^{-a_2}.$$(7)

Here, a₂ ≈ 0 indicates that alignment is retained throughout the denser regions of a cloud, while a₂ ≈ 1 is interpreted as polarized emission originating only from the outermost layer of dust (Andersson et al. 2015). However, additional effects, such as a decrease in polarization due to changes in the direction of the magnetic field lines on scales smaller than the beam size or the non-Gaussian nature of the noise on p (Vaillancourt 2006; Pattle et al. 2019) potentially lead to an overestimation of a₂.

In the context of far-infrared and submillimeter polarimetry, it is common to use the total intensity I as a proxy for the visual extinction A_V ¹¹. Together with the relation p_ext,λ = − τ_λp_emitλ¹² (Andersson et al. 2015), the dependency that is fitted to the polarization data becomes $$p_{\mathrm{e}\mathrm{m}\mathrm{i}\mathrm{t},\lambda }=a_1{\left(\frac{I}{I_{\mathrm{m}\mathrm{a}\mathrm{x}}}\right)}^{-a_2},$$(8)

where I_max denotes the intensity maximum in the considered region.

Fig. 2 shows the fit function from Eq. (8) applied to the data obtained at all three wavelengths in our regions of interest. Here, the fitting is done via the Markov chain Monte Carlo method. Every dataset is well fit¹³, indicated by values of ${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2<1$, with the only exception being the SOFIA/HAWC+ data in the RNO 91 region, where we find ${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=1.79$. The exponent a₂ shows little variation between regions and wavelengths, with values ranging from 0.81-1.03 in the L43E region and 0.82-1.04 in the RNO 91 region. These findings indicate that the polarization efficiency in the dense regions of L43 might be heavily attenuated. To assess whether the positive bias in p is still affecting our results, we repeat the fits with a stricter S/N-cutoff of p/δp > 3 and find no significant deviation from the previous values of a₂. One notable result is that we still see a strong depolarization occurring in the RNO 91 region, even though the presence of an embedded radiation source is expected to aid the alignment. To check whether this result still holds true in the direct vicinity of RNO 91, we repeat the fit in a circle with a 20″ radius around the intensity peak of the RNO 91 region with the SOFIA/HAWC+ data and find $a_2={0.65}_{-0.06}^{0.06}$. Due to the relatively low number of data points (14) in this region, the fit quality is slightly worse (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=2.53$).

Fig. 2 Degree of polarization as a function of the normalized intensity in the RNO 91 region (top panel) and the L43E region (bottom panel). The black lines indicate the best fit of Eq. (8) to each dataset. In the L43E region, we find $a_2={1.03}_{-0.10}^{0.10}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=0.63$) at 154 μm, $a_2={0.81}_{-0.07}^{0.08}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=0.73$) at 450 μm and$a_2={0.97}_{-0.05}^{0.05}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=0.99$) at 850 μm. In the RNO 91 region, we find $a_2={0.82}_{-0.02}^{0.02}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=1.79$) at 154 μm, $a_2={1.04}_{-0.10}^{0.10}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=0.73$) at 450 μm and $a_2={1.02}_{-0.09}^{0.09}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=0.69$) at 850 μm.

4.3 The far-infrared polarization spectrum of L43

Given a singular dust species (e.g., silicate or carbonaceous grains), we do not expect variations in the degree of polarization at far-infrared wavelengths for any common dust species in the interstellar medium (ISM; Hildebrand et al. 1999). However, many previous observations have shown that such variations exist (e.g., Vaillancourt & Matthews 2012; Zeng et al. 2013; Santos et al. 2019; Michail et al. 2021; Cox et al. 2025). To construct the far-infrared polarization spectrum, we restrict our analysis to those data points, where measurements at all three wavelengths are available, that fulfill the criteria defined in Eq. (6). Previous studies have further constrained the selection of data points by limiting the analysis to points where the difference in polarization angles between wavelengths is below a certain threshold (Vaillancourt 2002). This is done in an attempt to minimize the effects that a changing magnetic field along the LOS has on the polarization spectrum. Owing to the strong variability of polarization angles between wavelengths found in L43, we cannot apply this additional constraint and have to assume that the polarization spectrum is shaped in part by different inclinations of the magnetic field lines along the LOS.

Fig. 3 shows the global polarization spectrum for both considered subregions (as defined in Fig. 1), as well as for the entire region. Data points indicate the median degree of polarization and error bars the median absolute deviation (MAD, see also Table 2). Adhering to common practice, we normalize the spectra to the wavelength closest to 350 μm, which in our case is 450 μm. We find a decreasing spectrum in both subregions, with the only increase observed between 154 μm and 450 μm in the RNO 91 region. As Seifried et al. (2023) mention in their semi-analytical modeling for a homogeneous cloud consisting of silicate and carbonaceous grains, we expect a monotonically increasing spectrum caused by higher average temperatures of carbonaceous grains (compared to silicate grains) that are believed not to align with the external magnetic field. When a given LOS contains two or more dust phases with different temperatures and alignment properties, the resulting spectrum is expected to become V-shaped, with a local minimum at around 350 μm. The fact that our observed spectral slope is still negative between 450 μm and 850 μm implies that either the alignment of the dust in the cold phase is substantially impaired, or the inclination of the magnetic field is much higher in the cold phase (see Appendix B for details). Given the results obtained in Sect. 4.2.2, we can see evidence of a reduced grain alignment efficiency in both subregions. However, a magnetic field influenced by a potentially heavily inclined outflow, for which Lee & Ho (2005) find an inclination angle i ≈ (60 ± 15)° with respect to the POS, is likely contributing to the observed shape of the spectrum. Separating the contributions of these effects from observations alone is difficult. In Sect. 5, we investigate the effect of changing dust phases along the LOS.

Fig. 3 Global polarization spectra for the entire region (left panel), the L43E region (center panel), and the RNO 91 region (right panel). Points indicate the median, and error bars the MAD at a given wavelength. The spectra are normalized at λ0 = 450 μm. A black dashed line indicates unity.

4.4 Column density and temperature maps

Inferring the column density, LOS temperature, and optical depth of a molecular cloud requires a set of continuum measurements that sample the spectral energy distribution (SED) of the cloud. To do this, we use the measurements from the PACS and SPIRE instruments on board Herschel at 160 μm, 250 μm, 350 μm, and 500 μm as described in Sect. 3.3 in combination with the 850 μm measurement from SCUBA-2/POL-2 (Sect. 3.2). The beam and pixel sizes of all instruments were standardized according to Sect. 3.5 and the modified blackbody function (Hildebrand 1983) $$I_{\nu }=(1-\exp \left({\tau }_{\nu }\right))B_{\nu }(T_{\mathrm{d}})$$(9)

is fit to every pixel utilizing the Levenberg-Marquardt method. Here, B_v(T_d) refers to the Planck function at dust temperature T_d and τ_ν to the optical depth, described by the expression $${\tau }_{\nu }={\mu }_{{\mathrm{H}}_2}{m}_{\mathrm{H}}{N}_{{\mathrm{H}}_2}{\kappa }_{{\nu }_0}{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta },$$(10)

where μ_H2 is the mean molecular weight of molecular hydrogen, m_H the atomic mass of hydrogen, N_H2 the column density, κ_v0 the dust opacity per unit mass for the reference frequency v₀ and β the dust emissivity index. Keeping in line with other works (e.g., Sadavoy et al. 2013; Karoly et al. 2023), we adopt μ_H2 = 2.8, m_H = 1.00784u, κ_v0 = 0.1 cm²g⁻¹, v₀ = 10¹²Hz and β = 1.8. We simultaneously fit column density and temperature, as shown in Fig. 4. Pixels for which the error in T_d exceeds 5 K, or the error in the column density is larger than the column density itself, are excluded.

For the L43E region, we find an average column density of $\overline{N_{{\mathrm{H}}_2}}$ = (1.38 ± 0.34) 10²² cm⁻² with a peak value of N_H2,max = (4.72 ± 1.25) 10²² cm⁻² and a mean Temperature of T_d = (12.5 ± 0.4) K with a peak value of T_d,max = (16.6 ± 0.4) K. In the RNO 91 region, we find $\overline{N_{{\mathrm{H}}_2}}$ = (0.86 ± 0.27) 10²² cm⁻² and N_H2,max = (3.83 ± 1.55) 10²² cm⁻² with temperatures T_d = (15.2 ± 0.9) K and T_d,max = (23.5 ± 3.9) K. The optical depth can now be calculated via Eq. (10). For our shortest wavelength considered in this work, at 154 μm, we find τ_154μm,max ≈ 0.03.

Fig. 4 Spatially resolved column density (top panel) and LOS temperature (bottom panel) distributions, acquired by fitting Eq. (9) to Herschel PACS/SPIRE and JCMT/SCUBA-2 data at observing wavelengths of 160 μm, 250 μm, 350 μm, 500 μm, and 850 μm. Points where ∆Td > 5 K and NH2/ΔNH2 < 1 are excluded, where ∆Td and ∆NH2 refer to the errors on the respective fit parameter. Overlaid are contour lines of the 850 μm continuum emission from Fig. 1.

4.5 Magnetic field strength

4.5.1 Differential measure analysis

We apply the recently developed differential measure analysis (DMA; Lazarian et al. 2022) technique to our polarization and spectral line data to infer magnetic field strengths for the two considered subregions. Following their approach, the POS magnetic field strength in the most general form is given by $$B_{\mathrm{P}\mathrm{O}\mathrm{S}}\approx f\sqrt{4\pi ⟨\rho ⟩}\frac{D_{2D}^{1/2}\{C\}(\ell )}{D_{2D}^{1/2}\{\theta \}(\ell )},$$(11)

where f is a correction factor and 〈ρ〉 is the mean density of the gas. The expressions D_2D{C}(ℓ) and D_2D{θ}(ℓ) describe the second-order structure functions of the mean velocity and the polarization angle, respectively, measured at distance ℓ. For some quantity Q(X) at POS coordinate X, the second-order structure function D_2D{Q}(ℓ) is given by the expression $$D_{2\mathrm{D}}\{Q\}(\ell )=⟨[Q(\mathbf{X}+\ell )-Q(\mathbf{X}){]}^2{⟩}_{\mathbf{X}}.$$(12)

Similar to the regularly employed DCF method (Davis 1951; Chandrasekhar & Fermi 1953) and its more recent modifications (e.g., Hildebrand et al. 2009; Houde et al. 2009; Skalidis & Tassis 2021), with DMA, one compares the dispersion of velocities and magnetic field directions. This relation stems from the Alfvénic coupling between perturbations in turbulent velocity δv_turb and magnetic field δB, given by δv_turb = δB/√4πρ. However, Lazarian et al. (2022) argue that the use of global measures, as is the case with DCF, is only applicable when the sampling size for the cloud is greater than the scale L_inj, at which turbulent motions are injected into the system¹⁴. If this is not the case, DCF will overestimate the magnetic field strength.

DMA also allows one to account for the anisotropy of MHD turbulence, by encoding the dependence of B_POS on the Alfvénic Mach number M_A, the angle between the mean magnetic field and the POS γ, and the ratio between gas pressure and magnetic pressure β_gas/mag in the correction factor f = f (M_A,γ,β_gas/mag). They derive the correction factor for the low β_gas/mag-regime through numerical simulations and find $$f(\gamma )\approx \{\begin{array}{l}0.5{\mathcal{M}}_{\mathrm{A}}^{1/2},\mathrm{i}\mathrm{f}\gamma =\pi /2\\ 0.83{\mathcal{M}}_{\mathrm{A}}^{1/2},\mathrm{i}\mathrm{f}\gamma =\pi /4\end{array}$$(13)

(see their Figs. 15 and 20). Inferring the angle γ is not straightforward. Therefore, we will use these values as our upper and lower bounds for the estimation of B_POS.

4.5.2 Velocity structure

DMA makes use of velocity centroids, available from spectroscopic observations of Doppler-shifted molecular line transitions, as a proxy for the spatial distribution of the mean velocity. Given the Intensity I(X, v_i) of the molecular line in the i-th velocity channel v_i and at POS coordinate X, one can find the velocity centroid C(X) by computing $$C(\mathbf{X})=\sum _i{v}_{i}I(\mathbf{X},v_i)/\sum _{i}I(\mathbf{X},v_i)$$(14)

over a chosen molecular line. One caveat of this approach is that the line profile is not solely influenced by velocity fluctuations along the LOS, but also by fluctuations in density. We aim to mitigate this issue by decomposing the channel intensity into a contribution that is correlated with density (which we disregard for our analysis) and one that is not (for details, we refer to Yuen et al. 2021).

We compute density-corrected velocity centroids of the N₂H⁺(101-012) hyperfine transition from the data described in Sect. 3.4 over a range of 3 km s⁻¹. N₂H⁺ has proven to be a useful tracer of gas kinematics inside dense cloud cores, as it is less affected by condensation on dust surfaces compared to CO and remains optically thin even at high volume densities (Tafalla et al. 2002; Lippok et al. 2013). The resulting centroid map is then regridded to a resolution of 12″, according to the procedure described in Sect. 3.5. Given the larger beam size of the NRO data (see Table 1), for this section, we reapply Eq. (1) to the polarization data, as we need to make sure the resolution of all datasets remains comparable¹⁵.

The study of velocity centroids also allows us to find an independent measure for the Alfvénic Mach number M_A through the application of the gradient technique (GT, González-Casanova & Lazarian 2017; Lazarian & Yuen 2018; Lazarian et al. 2024), where gradients of velocity centroids are used to infer properties of the underlying magnetic field. Lazarian et al. (2024) derive a connection between the gradient covariance tensor along the directions i and j, i.e., ${\sigma }_{ij}\equiv ⟨{\mathrm{\nabla }}_{X^i}C(\mathbf{X}){\mathrm{\nabla }}_{X^j}C(\mathbf{X})⟩$ and M_A, given by the expression $$\frac{({\sigma }_{xx}-{\sigma }_{yy}{)}^2+4{\sigma }_{xy}}{({\sigma }_{xx}+{\sigma }_{yy}{)}^2}\approx \exp (-2({\mathcal{M}}_{\mathrm{A}}/\sin \gamma {)}^2).$$(15)

4.5.3 Application of DMA

We estimate the mean density 〈ρ〉 in the two subregions assuming that all dust and gas along the LOS is located inside an ellipsoid with semimajor and semiminor POS half-axes a_HA and b_HA and the depth of the ellipsoid being taken as the geometric mean $c_{\mathrm{H}\mathrm{A}}=\sqrt{a_{\mathrm{H}\mathrm{A}}{b}_{\mathrm{H}\mathrm{A}}}$ of these half-axes (Karoly et al. 2023). The mean volume density for each subregion is then given by $$⟨\rho ⟩=\mu m_{\mathrm{H}}{n}_{{\mathrm{H}}_2}=\frac{3\mu m_{\mathrm{H}}}{4}\frac{N_{{\mathrm{H}}_2}}{c_{\mathrm{H}\mathrm{A}}}[\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}].$$(16)

By combining Eqs. (13) and (15), one obtains the value of f for the two considered inclinations of the magnetic field. To obtain an estimate of B_POS, we calculate the structure functions $D_{2D}^{1/2}\{C\}(\ell )$ and $D_{2D}^{1/2}\{\theta \}(\ell )$ at a separation distance ℓ = 24″. This distance ensures the independence of the data pairs, as it is larger than the largest beam size, and likely also fulfills the second criterion for the application of DMA, which is that ℓ is required to be smaller than the LOS extent of the cloud (Lazarian et al. 2022). Fig. 5 shows the values of $B_{\mathrm{P}\mathrm{O}\mathrm{S},\gamma =\pi /4}$, calculated from Eq. (11), for each subregion and wavelength as a function of the distance ℓ between the data pairs.

Given the POS magnetic field strength and the column density (see Sect. 4.4), one can use the mass-to-flux ratio λ_M/Φ (Nakano & Nakamura 1978; Crutcher et al. 2004) $${\lambda }_{M/\mathrm{\Phi }}=\frac{(M/\mathrm{\Phi }{)}_{\mathrm{o}\mathrm{b}\mathrm{s}}}{(M/\mathrm{\Phi }{)}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}=7.6\times {10}^{-21}\frac{N_{{\mathrm{H}}_2}}{B_{\mathrm{P}\mathrm{O}\mathrm{S}}}\left[\frac{\mathrm{c}{\mathrm{m}}^{-2}}{\mu \mathrm{G}}\right]$$(17)

to obtain a measure for the required magnetic flux Φ that would prevent the gravitational collapse of a structure with mass M. If λ_M/Φ < 1, the magnetic field is strong enough to prevent the structure from collapsing (magnetically subcritical regime), while λ_M/Φ > 1 indicates that the magnetic field is unable to prevent the collapse by itself (magnetically super-critical regime). The determination of λ_M/Φ through this method suffers from geometrical biases, depending on the LOS and the orientation of the cloud, as explained in Crutcher et al. (2004). A statistical correction for this bias yields λ_M/Φ,cor = λ_M/Φ/3. A summary of all results obtained in this section is given in Table 3. A discussion of the results can be found in Sect. 6.2.

Fig. 5 POS magnetic field strength $B_{\mathrm{P}\mathrm{O}\mathrm{S},\gamma =\pi /4}$ as a function of the separation distance ℓ for the L43E (solid lines) and RNO 91 (dashed lines) regions. Polarization measurements are taken at 154 μm (blue), 450 μm (green) and 850 μm (red).

4.5.4 Comparison with DCF

The DCF technique has been the primary method for inferring magnetic field strengths in the ISM for many years. To test the prediction made by DMA, namely that DCF will overestimate the magnetic field strength for ℓ ≪ L_inj, we also applied the modified DCF approach of Hildebrand et al. (2009) to our data. We find that for both subregions, the DCF estimates are higher by a factor of ~ 1.5-2.5, in agreement with the prediction of Lazarian et al. (2022).

5 Modeling the polarization spectrum of a molecular cloud core with POLARIS

5.1 The model

In this Section, we employ the 3D Monte Carlo radiative transfer code POLARIS (Reissl et al. 2016) to analyze how variations in dust alignment properties and temperatures along the LOS affect the slope of the polarization spectrum under the radiative torque paradigm. For this purpose, we construct a simple model of a molecular cloud core that is irradiated by the interstellar radiation field (ISRF), as well as a nearby star. We assume a Bonnor-Ebert density distribution of the form $$\rho (r)=\{\begin{array}{ll}{\rho }_0,\mathrm{i}\mathrm{f}& r\le r_0,\\ {\rho }_0{r}^{-2}/r_0,\mathrm{i}\mathrm{f}& r_0<r\le r_{\mathrm{o}\mathrm{u}\mathrm{t}},\\ 0,\mathrm{i}\mathrm{f}& r_{\mathrm{o}\mathrm{u}\mathrm{t}}<r,\end{array}$$(18)

which has already been used in other studies to model molecular cloud cores (e.g., Brauer et al. 2016; Zielinski et al. 2021; Giang et al. 2023). Here, r₀ is the truncation radius, at which the density distribution begins to decrease with r⁻², until the outer radius r_out is reached. The model parameters (see Table 4) are chosen to resemble the environment of L43. Based on the emission profiles in Fig. 1 and the distance estimate of de Geus et al. (1989), we assume a truncation radius of 3000 au, an outer radius of 3 × 10⁴ au, and a distance between the YSO and the density peak of L43E of 10⁴ au. The total mass is taken to be 3 M_⊙, matching the previously derived column density in L43 (see Fig. 4). We consider two different setups. In Setup one, we only consider the ISRF, where photon packages start propagating from the outside of the model space. The total energy density of the ISRF is taken to be u_ISRF = 8.64 × 10⁻¹³ erg cm⁻³, which is equal to the total energy density of the solar neighborhood (Mathis et al. 1983). Setup two considers the ISRF, as well as a stellar source, located inside the model space at a displacement radius of 10⁴ au. We apply a constant, homogeneous magnetic field B = 100 μG, aligned perpendicular to the LOS and rotated by 150° with respect to the axis defined by the center of the model space and the location of the star in Setup two, to roughly match the large-scale field of L43 seen by Planck (Ade et al. 2015). The energy of the photon packages propagating from the star is given by its SED, which is modeled to resemble RNO 91. We assume a surface temperature of 4000 K, motivated by the average effective temperatures of young stars in the region (Doppmann et al. 2005; McClure et al. 2010) with a bolometric luminosity of L_* = 2.5 L_⊙ (Chen et al. 2009). We consider oblate dust grains with effective radii a ranging from 5 nm to 500 nm and an axis ratio of 1:2, distributed in the form of a power law (n(a) ∝ a^α), with α = −3.5. The dust consists of a mixture of 62.5% silicate and 37.5% carbonaceous grains, whose optical properties are calculated with DDSCAT (Draine & Flatau 1994) using the optical properties given in Draine (2003). As already pointed out by Seifried et al. (2023), differences in temperatures between aligned and unaligned dust species are a key contributor to the observed shape of the polarization spectrum. For a given density distribution consisting of a composition of dust species, POLARIS calculates a singular dust temperature for each cell by default. Thus, we consider two distinct but spatially overlapping dust density distributions for each dust species. Guidance on the implementation of this feature is given in Appendix A of Seifried et al. (2023). Size-dependent dust temperatures are not considered. Instead, we estimate their impact in Appendix B. We assume a constant gas-to-dust mass ratio of 100 : 1.

5.2 Alignment mechanisms

In the following, we assume that the alignment of the grains along their minor axis and the direction of the external magnetic field is given by RAT alignment. In POLARIS, the torque Γ_rad exerted on a grain is calculated using the expression (Hoang & Lazarian 2014; Reissl et al. 2020) $${\mathrm{\Gamma }}_{\mathrm{r}\mathrm{a}\mathrm{d}}=\pi a^2\int \left(\frac{\lambda }{2\pi }\right){\gamma }_{\lambda }\cos (\psi )Q_{\mathrm{\Gamma }}(a,\lambda )u_{\lambda }\mathrm{d}\lambda .$$(19)

Here, γ_λ is a factor quantifying the degree of anisotropy of the radiation field. A value of γ_λ = 1 indicates a unidirectional radiation field, while γ_λ = 0 occurs if the radiation field is isotropic. The quantity ψ refers to the angle between the direction of the incident radiation and that of the magnetic field, Q_Γ(a, λ) is a factor determining the RAT efficiency, and u_λ is the spectral energy density of the radiation field. Dust grains are assumed to be aligned if the condition $$\frac{J_{\mathrm{r}\mathrm{a}\mathrm{d}}}{J_{\mathrm{t}\mathrm{h}}}=\frac{{\tau }_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}{\mathrm{\Gamma }}_{\mathrm{r}\mathrm{a}\mathrm{d}}}{J_{\mathrm{t}\mathrm{h}}}\ge 3$$(20)

is met, where J_rad and J_th refer to the angular momentum caused by radiative torques and thermal excitations, respectively. The quantity τ_drag refers to the damping time, after which the grain’s direction of rotation is randomized through gas collisions.

We also consider imperfect internal alignment. Here, internal alignment refers to the ability of the grain to align its angular momentum parallel to the short axis of the grain. This condition depends on the angular momentum a grain receives under the influence of RATs. Lazarian & Hoang (2007) found that the angular momentum of a grain settles at either a high-J or low-J attractor point. Grains at high-J are less prone to processes that randomize their orientation and therefore contribute more to the polarized emission than grains sitting at a low-J attractor point. The ratio f_high-J of grains that settle at the high-J attractor point depends on the size and magnetic properties of the grain (Hoang & Lazarian 2016; Giang et al. 2023), but is still poorly constrained. Therefore, we assume a conservative estimate of f_high-J = 0.25.

5.3 Fitting the spatially resolved polarization spectrum

Polarization maps with a pixel resolution of 255 × 255, corresponding to 3.14″px⁻¹, are calculated for the observing wavelengths of 154 μm, 450 μm and 850 μm. The maps are convolved using a Gaussian kernel with a FWHM of 14.6″, corresponding to the effective beam size of the set of polarization observations used in Sect. 4, after Eq. (1) is applied.

Keeping in line with previous works (e.g., Gandilo et al. 2016; Michail et al. 2021; Lee et al. 2024; Cox et al. 2025), we investigate changes in the slope of the polarization spectrum across different LOSs by fitting a function of the form $$p(\lambda )/p_{{\lambda }_0}=a_l(b_l[\lambda -{\lambda }_0]+1),$$(21)

to each pixel of the model space. Here, similar to Sect. 4.3, the reference wavelength λ₀ is taken to be 450 μm. Given that a_l ≈ 1¹⁶, we investigate correlations between the effective slope b_l and cloud properties, specifically the column density and LOS temperature.

Fig. 6 2D histograms of the effective slope bl of the normalized polarization spectrum, as a function of LOS temperature (left panels) and column density (right panels), respectively, colored by the count of data points in each bin. Top panels show results obtained with the ISRF as the only considered radiation source (Setup one), whereas bottom panels show results where the ISRF as well as an offset stellar source are considered (Setup two, see Sect. 5.1 for details). Pearson r, as well as two-tailed p values (ptt), are calculated from 50 segmented median values in each distribution.

5.4 Results

In Fig. 6, the value of the effective slope b_l is shown as a function of LOS temperature and column density. The distribution is segmented into 50 chunks of equal size. The median and MAD of each chunk are displayed as a point and corresponding error bar, respectively. From these medians, we infer correlations by computing the Pearson r correlation coefficient along with a twotailed p value, denoted as p_tt ¹⁷, indicating the probability that the null hypothesis of no correlation is incorrect.

For Setup one, where the ISRF acts as the only radiation source (top panels of Fig. 6), we find LOS temperatures ranging from 13 K in the outer regions of the model space to 11 K in the center regions. This decline is expected, since the radiation responsible for heating the grains is increasingly attenuated as it traverses through the model space. Therefore, due to the symmetry of this setup, column density and LOS temperature are strongly anticorrelated. With Setup one, we find spectra with positive slopes throughout the model space that decline steeply for column densities ≳ 10²² cm⁻² or LOS temperatures of ≲ 12 K. Similar trends have been observed in observational studies (e.g., Michail et al. 2021; Cox et al. 2025). However, more diffuse LOSs in the outer regions of the model space exhibit the opposite of this trend, where b_l first increases with increasing column density (or decreasing LOS temperature), but this effect is not of any physical nature. Rather, it is the result of fitting a linear function to a normalized, nonlinear distribution. In Appendix C, we further investigate this effect by applying the fits to the data before normalization.

In Setup two (bottom rows in Fig. 6), the influence of the stellar radiation results in higher LOS temperatures. In contrast to Setup one, after an initial increase, bl decreases with increasing LOS temperatures. Data points on the far side of the stellar source are less affected by the stellar source and thus show a similar dependency to those in Setup one, leading to the initial positive correlation. The data points closest to the stellar source exhibit spectra with the shallowest slopes. To explain these differences in correlation between external and embedded radiation sources, we employ another measure, that is, the standard deviation of the LOS temperature profile σ_T (Lee et al. 2024). This measure aims to quantify the degree of heterogeneity along a given LOS, which is thought to be the cause of the observed negative or V-shaped spectra (Santos et al. 2019; Michail et al. 2021; Seifried et al. 2023). Fig. 7 shows the correlation between b_l and σ_T. In contrast to LOS temperature and column density, b_l decreases with σ_T independently of the radiation source considered, indicating that heterogeneity along the LOS governs the slope of the polarization spectrum.

Fig. 7 2D histograms of the effective slope bl of the normalized polarization spectrum, as a function of the standard deviation of the LOS temperature profile σT, colored by the count of data points in each bin. Top panels show results obtained with the ISRF as the only considered radiation source (Setup one), whereas bottom panels show results where the ISRF as well as an offset stellar source are considered (Setup two, see Sect. 5.1 for details). Pearson r, as well as two-tailed p values (ptt), are calculated from 50 segmented median values in each distribution.

6 Discussion

6.1 Polarization efficiency

The RAT mechanism provides a natural explanation for the observed trend between the degree of polarization and the normalized intensity (Sect. 4.2.2), in which dust that is shielded from radiation exhibits a lower polarization efficiency compared to dust in a radiation-rich environment, such as the outer layers of molecular clouds. The fact that the polarization efficiency shows a similar decrease in the RNO 91 region, which contains an embedded radiation source, challenges this interpretation. However, when considering only the innermost region around the YSO, we observe a decrease in the value of the fit parameter a₂, indicative of a slight but significant enhancement of the polarization efficiency in this region compared to the rest of the cloud.

Although some alternative explanations for the occurrence of depolarization toward density peaks, such as the influence of dichroic extinction (Brauer et al. 2016) can be ruled out due to the low optical depths at all observing wavelengths (see Sect. 4.4), other mechanisms, such as beam-averaging effects due to changes in the magnetic field smaller than the beam size or changes in grain composition cannot.

6.2 Gravitational stability

In Sects. 4.4 and 4.5, we have determined column densities and magnetic field strengths to constrain the ability of the magnetic field to inhibit the gravitational collapse in the starless core L43E. We find peak column densities of ~5 × 10²² cm⁻², which are approximately six times lower than what Karoly et al. (2023) find, who used a similar dataset and method, but agrees well with the findings of Chen et al. (2009). This discrepancy is likely caused by the fact that Karoly et al. (2023) fit temperature and column density in two steps, whereas in our method, we fit them simultaneously. Furthermore, the poorly constrained dust properties used in the fit cause additional uncertainties of up to 50% (Roy et al. 2014; Zielinski & Wolf 2024).

The differences in polarization angles observed at various positions, as well as their variations between observing wavelengths, indicate a complex magnetic field structure in the region. The alignment of the inferred magnetic field lines with the cavity walls of the CO outflow seen at 450 μm and 850 μm indicates that the magnetic field in the colder and denser regions of L43E has been influenced by the outflow, further complicating the determination of magnetic field strengths. The absence of this influence in the 154 μm data suggests that we are measuring different parts of the cloud along the LOS, caused by variations in local temperatures. The mass-to-flux ratio λ_M/Φ,cor (see Eq. (17)), corrected for geometrical biases (Crutcher et al. 2004), is greater than unity in L43E (see Table 3), suggesting that the cloud is likely gravitationally unstable. L43E is slightly sub-Alfvénic, indicating that the dynamics are dominated to a greater extent by the magnetic field, rather than by turbulence. For the RNO 91 region, we do not provide mass-to-flux ratios because the region already contains a YSO. The Alfvénic Mach number in this region slightly exceeds unity when assuming that the mean magnetic field direction lies perpendicular to the LOS. The low dispersion seen especially in the 154 μm data suggests that the dynamics are dominated by the magnetic field, though this might not hold true for all regions within the cloud, as suggested by the high dispersion measured at 850 μm.

6.3 Influences on the shape of the polarization spectrum

The simple model we employed to study the polarization spectrum of a molecular cloud core under a homogeneous magnetic field (see Sect. 5) shows that an increase in heterogeneity along a given LOS significantly decreases the observed slope of the polarization spectrum - a phenomenon known as the heterogeneous cloud effect (HCE, Michail et al. 2021).

In the diffuse regions of the model, the spectrum steadily increases, resulting from differences in temperatures between silicate and carbonaceous grains, as predicted by Seifried et al. (2023). Under the assumption that the dust properties, radiation field strength, and column densities that we used in our model are indeed similar to those found in L43, additional influences, such as differences in the inclination of the magnetic field lines along the LOS, are needed to explain the observed negative slopes in Fig. 3. The existence of such differences in inclination is plausible, considering that the inclination of the CO outflow is thought to be ~60° with respect to the POS (Lee & Ho 2005). Given an accurate characterization of the dust properties in a molecular cloud, one can infer the underlying 3D magnetic field structure through multiwavelength polarimetry (e.g., Tram et al. 2024).

7 Conclusions

In this paper, we have analyzed the polarization signal in the molecular cloud core L43, arising from magnetically aligned dust grains in the 154-850 μm wavelength range and compared our findings with radiative transfer simulations of a simple cloud core model. Our findings are summarized as follows:

1. SOFIA/HAWC+ measurements at 154 μm, as well as SCUBA-2/POL-2 measurements at 450 μm and 850 μm reveal wavelength-dependent changes in emission structure, polarization angles and degrees in both considered subregions of the molecular cloud core L43, indicating a complex magnetic field morphology;

2. The CO outflow associated with the YSO RNO 91 seems to have influenced the polarization signal at 450 μm and 850 μm, indicating the alignment of the magnetic field with the outflow direction. This influence is not seen in the 154 μm data;

3. The relationship between the degree of polarization and the normalized intensity suggests that the polarization efficiency is heavily attenuated in the inner parts of both subregions. However, analysis constrained to the 20″ region around the YSO suggests that alignment is significantly more retained here, in line with predictions that RATs cause more efficient alignment in the presence of an anisotropic radiation source;

4. We find peak column densities of N_H2,max = (4.72 ± 1.25) 10²² cm⁻² in the main starless core L43e, which is considerably lower than what Karoly et al. (2023) find using a similar method, but agrees well with previous estimates of Chen et al. (2009). Temperatures in L43E range from 1015 K. The peak temperature is found in the southern part of the RNO 91 region, giving T_d,max = (23.5 ± 3.9) K;

5. We apply the differential measure analysis method (Lazarian et al. 2022) to infer magnetic field strengths in both subregions at all three wavelengths. In the L43E region, we find POS field strengths ranging from 13-24 μG, and in the RNO 91 region, we find field strengths ranging from 1960 μG. These estimates are subject to several uncertainties and should be interpreted with care. We find mass-to-flux ratios λ_M/Φ,cor exceeding unity in L43E, suggesting a gravitationally unstable core;

6. We find polarization spectra with negative slopes throughout the cloud. A comparison with radiative transfer simulations suggests that HCE does have a significant influence on the spectral slope, but additional influences, such as varying inclinations of the magnetic field lines along the LOS, are likely needed to fully explain observations.

Acknowledgements

The authors thank the anonymous referee for helpful comments and suggestions. We also thank Dr. Janik Karoly for helpful insights on the data reduction with SMURF and Dr. Moritz Lietzow-Sinjen, as well as Dr. Stefan Reissl, for their helpful comments on working with POLARIS. M.S. and S.W. thank the DFG for financial support under grant WO 857/251. This research is based (in part) on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 2002 to the University of Stuttgart. These observations were obtained (in part) by the James Clerk Maxwell Telescope, operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; the National Astronomical Research Institute of Thailand; Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities and organizations in the United Kingdom and Canada. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation. The Nobeyama 45-m radio telescope is operated by Nobeyama Radio Observatory, a branch of the National Astronomical Observatory of Japan. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. Facilities: SOFIA (HAWC+), JCMT (SCUBA-2, POL-2), Herschel (PACS, SPIRE) Software: numpy (van der Walt et al. 2011), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), astropy (Astropy Collaboration 2013, 2018), emcee (Foreman-Mackey et al. 2013), Starlink (Currie et al. 2014).

References

Appendix A Histograms of polarization degree and angle

Fig. A.1 Histograms of the degree of polarization (top panels) and the by 90° rotated polarization angles (bottom panels) for the three observing wavelengths at 154 μm (blue), 450 μm (green) and 850 μm (red). The dashed black line indicates the mean, and the gray shaded area the standard deviation of the distribution. Left panels: L43E region. Right panels: RNO 91 region.

Appendix B The impact of grain-size-dependent dust temperatures on the polarization spectrum

In this Section, we estimate the impact of grain-size-dependent dust temperatures on the shape of the polarization spectrum. In contrast to Sect. 5, we do not consider the full radiative transfer of individual photon packages. Instead, we calculate the temperature T_d(a) of a dust grain when subjected to a radiation field with a known spectral energy density u_λ, assuming local thermal equilibrium. Equating the heating and cooling rates of a grain yields (Draine & Lee 1984) $$\frac{4\pi }{c}{\int }_0^{\mathrm{\infty }}{Q}_{\mathrm{a}\mathrm{b}\mathrm{s}}(a,\lambda )B_{\lambda }(T_{\mathrm{d}}(a))\mathrm{d}\lambda ={\int }_0^{\mathrm{\infty }}{Q}_{\mathrm{a}\mathrm{b}\mathrm{s}}(a,\lambda )u_{\lambda }\mathrm{d}\lambda ,$$(B.1)

where Q_abs(a, λ) describes the absorption efficiency for a given grain size a and wavelength λ. The left side of this equation characterizes the cooling of the grain by the thermal radiation it emits, whereas the right side describes the heating through the absorption of light from the incident radiation field. Given the precalculated absorption efficiencies used in the modeling in Sect. 5.1, combined with the spectral energy density of the ISRF in the solar neighborhood (Mathis et al. 1983), we can calculate the grain-size-dependent temperatures for both of the dust species we consider in Sect. 5.

As the radiation field interacts with dust in a molecular cloud, back-scattering and absorption events gradually attenuate the radiation field, causing the spectral energy density to redden. The grain-size-dependent temperature distribution of an ensemble of dust grains exposed to a reddened radiation field, therefore, changes its shape as well. In Fig. B.1, we use Eq. B.1 to calculate the grain-size-dependent temperatures for silicate grains (left) and carbonaceous grains (right), as well as their mass-weighted average temperatures, for varying optical depths at a reference wavelength of 1 μm, given the optical properties of the dust used in Sect. 5. It should be noted that this approach, unlike the modeling with POLARIS in Sect. 5, neglects the influence of heating through thermal reemission of dust grains, as well as changes in the shape of the local spectral energy density due to scattering. However, the dust temperatures shown in Fig. B.1 correspond to grain-size-dependent temperatures calculated with POLARIS within 0.1 K.

Fig. B.1 Grain-size-dependent temperatures for differing levels of attenuation (solid lines) along with their mass-weighed averages (dashed lines), calculated with Eq. B.1. Left: Silicate grains. Right: Carbonaceous grains.

To investigate the impact of grain-size-dependent dust temperatures on the polarization spectrum, we employ the semi-analytical model from Seifried et al. (2023). In this model, changes in the local dust temperature along a given LOS are approximated by two independent phases, each phase contributing to the polarized-, as well as the total intensity at a given wavelength. The total intensity of the cold phase is given by $$I_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}=\sum _i{f}_i\times {\int }_{a_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{a_{\mathrm{m}\mathrm{a}\mathrm{x}}}{B}_{\lambda }(T_{\mathrm{d},i}(a))\times \pi a^2\times Q_{\mathrm{a}\mathrm{b}\mathrm{s},i}(a,\lambda )\times a^{\alpha }\mathrm{d}a,$$(B.2)

where i denotes the dust species with mass fraction f_i, and πa² × Q_abs(a, λ) the absorption cross-section. We assume that only silicate grains of sizes within a_{alig,cold,min} and a_{alig,cold,max} contribute to the polarized emission in the cold phase. Therefore, the polarized intensity of this phase is given by $$I_{\mathrm{p}\mathrm{o}\mathrm{l},\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}=f_{\mathrm{s}\mathrm{i}\mathrm{l}}\times {\int }_{a_{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g},\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d},\mathrm{m}\mathrm{i}\mathrm{n}}}^{a_{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g},\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d},\mathrm{m}\mathrm{a}\mathrm{x}}}{B}_{\lambda }(T_{\mathrm{d},i}(a))\times \pi a^2\times Q_{\mathrm{a}\mathrm{b}\mathrm{s},i}(a,\lambda )\times a^{\alpha }\mathrm{d}a.$$(B.3)

Similarly, we model the total and polarized intensities from the hot phase as $$I_{\mathrm{h}\mathrm{o}\mathrm{t}}=f_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}}\times \sum _i{f}_i\times {\int }_{a_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{a_{\mathrm{m}\mathrm{a}\mathrm{x}}}{B}_{\lambda }(T_{\mathrm{d},i}^{\mathrm{\prime }}(a))\times \pi a^2\times Q_{\mathrm{a}\mathrm{b}\mathrm{s},i}(a,\lambda )\times a^{\alpha }\mathrm{d}a$$(B.4)

and $$I_{\mathrm{p}\mathrm{o}\mathrm{l},\mathrm{h}\mathrm{o}\mathrm{t}}=f_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}}\times f_{\mathrm{s}\mathrm{i}\mathrm{l}}\times {\int }_{a_{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g},\mathrm{h}\mathrm{o}\mathrm{t},\mathrm{m}\mathrm{i}\mathrm{n}}}^{a_{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g},\mathrm{h}\mathrm{o}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}}{B}_{\lambda }(T_{\mathrm{d},i}^{\mathrm{\prime }}(a))\times \pi a^2\times Q_{\mathrm{a}\mathrm{b}\mathrm{s},i}(a,\lambda )\times a^{\alpha }\mathrm{d}a,$$(B.5)

where we include a factor f_scale that determines the relative abundance of the material, i.e., gas and dust, in the hot phase in relation to that in the cold phase (Seifried et al. 2023). Here, T′_di(a) refers to either the un-attenuated temperature distribution T_{τ(1 μm)=0.0}(a) or its mass-weighted average T_{const,τ<1 μm)=0.0} (see Fig. B.1).

The wavelength-dependent polarization fraction then becomes $$p(\lambda )=\frac{I_{\mathrm{p}\mathrm{o}\mathrm{l},\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}+I_{\mathrm{p}\mathrm{o}\mathrm{l},\mathrm{h}\mathrm{o}\mathrm{t}}}{I_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}+I_{\mathrm{h}\mathrm{o}\mathrm{t}}}.$$(B.6)

For all cases considered, we assume that the upper size limit for alignment is higher than the maximum grain size, that is, a_alig,hot,max = a_{alig,cold,max} = a_max and that the grains in the hot phase have a lower size limit for alignment of a_alig,hot,min = 50 nm, which is similar to values found for dust in the diffuse ISM (Andersson et al. 2015). For the cold phase, we consider a case where alignment is moderately well retained, setting a_{alig,cold,min} = 250 nm and another case where alignment is lost, resulting in a_{alig,cold,min} = a_max = 500 nm. In Fig. B.2 we compare the polarization spectrum obtained with the grain-size-dependent temperatures from Fig. B.1 with the spectra obtained through their mass-weighted average temperatures at moderate levels of attenuation (τ_cold(1 μm) = 5), as well as high levels of attenuation (τ_cold(1 μm) = 50). Furthermore, we examine how the variation of f_scale influences this comparison, by considering f_scale = 1 and f_scale = 10⁻². In the case where alignment is still retained (Fig. B.2, left), we generally see increasing or V-shaped spectra. Compared to the case where alignment in the colder regions is lost (Fig. B.2, right), we observe that the spectra that were previously V-shaped are now continuously decreasing, even at wavelengths λ > 1 mm. This shows that at longer wavelengths, the shape of the spectrum is defined by the contribution of the aligned population in the cold phase. In general, the difference between spectra where we consider grain-size-dependent temperatures and those where constant temperatures are applied is negligible compared to other uncertainties, such as differences in dust properties, temperatures, or the direction of the magnetic field lines with respect to the observer. The highest deviations are observed for the cases where f_scale = 10⁻² and τ_cold(1 μm) = 5, regardless of whether the alignment is retained or not.

Appendix C The effect of normalization on the polarization spectrum slope

To investigate whether the initial increase in the slope of the polarization spectrum b_l, seen in Figs. 6 and 7, can be attributed solely to the use of normalization, we fit Eq. 21 to the same data we previously analyzed in Sect. 5.

Fig. B.2 Semi-analytical polarization spectra for varying levels of relative abundance and attenuation between the hot and the cold phases. Exemplary polarization spectra, where the cold phase retains moderate levels of alignment (left panel), and where the alignment in the cold phase is lost (right panel).

Fig. C.1 Same as Fig. 6, but with the true slope albl and without the use of normalization.

Here, we cannot apply the assumption of a_l ≈ 1 anymore and refer to the product a_lb_l as the slope instead. In Fig. C.1, the Pearson r correlation coefficient between the slope of the spectrum and the LOS temperature (left), as well as the column density (right), is calculated for the two setups, similar to Sect. 5.4. The initial increase of the spectrum slope for more heterogeneous regions is no longer observed, resulting in a strong correlation with cloud properties. Fig. C.2 shows the level of correlation between the spectral slope and the standard deviation of the LOS temperature distribution, similar to Fig. 7. For Setup one, where the ISRF acts as the only radiation source (top panel), we again see a strong correlation between these quantities. In the case where we consider a stellar source as well, we find a weaker correlation. This behavior arises even though the presence of either radiation source causes the spectrum slope to decrease. When calculating the correlation coefficient through the medians of chunked data points, differences in the sensitivity to σ_T can lead to the median value increasing with higher levels of heterogeneity.

We note that an increase in the spectral slope with column density can alternatively arise through effects related to the optical depth. As Lee et al. (2024) mention, if the optical depth is high, polarized emission from within the cloud is partially absorbed before exiting, resulting in a reduction of the degree of polarization at shorter wavelengths and, in turn, a more positive slope of the spectrum. Since in our case τ ≪ 1 holds for all wavelengths considered, this effect is negligible.

Fig. C.2 Same as Fig. 7, but with the true slope albl and without the use of normalization.

All Tables

All Figures

Fig. 1 Dust continuum maps of the L43 region. Black lines indicate the polarization vectors. The vectors are rotated by 90° to indicate the magnetic field direction and are binned to the same 12″ grid. Black contour lines mark regions where the flux has doubled. White dashed ellipses show the two regions of interest that are referred to throughout this work. The effective beam size of each instrument is indicated with a red circle. Left: 154 μm continuum emission (SOFIA/HAWC+). Top right: 450 μm continuum emission (SCUBA-2/POL-2). Bottom right: 850 μm continuum emission (SCUBA-2/POL-2).

In the text

Fig. 2 Degree of polarization as a function of the normalized intensity in the RNO 91 region (top panel) and the L43E region (bottom panel). The black lines indicate the best fit of Eq. (8) to each dataset. In the L43E region, we find $a_2={1.03}_{-0.10}^{0.10}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=0.63$) at 154 μm, $a_2={0.81}_{-0.07}^{0.08}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=0.73$) at 450 μm and$a_2={0.97}_{-0.05}^{0.05}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=0.99$) at 850 μm. In the RNO 91 region, we find $a_2={0.82}_{-0.02}^{0.02}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=1.79$) at 154 μm, $a_2={1.04}_{-0.10}^{0.10}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=0.73$) at 450 μm and $a_2={1.02}_{-0.09}^{0.09}$ (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2=0.69$) at 850 μm.

In the text

	Fig. 3 Global polarization spectra for the entire region (left panel), the L43E region (center panel), and the RNO 91 region (right panel). Points indicate the median, and error bars the MAD at a given wavelength. The spectra are normalized at λ0 = 450 μm. A black dashed line indicates unity.
In the text

Fig. 4 Spatially resolved column density (top panel) and LOS temperature (bottom panel) distributions, acquired by fitting Eq. (9) to Herschel PACS/SPIRE and JCMT/SCUBA-2 data at observing wavelengths of 160 μm, 250 μm, 350 μm, 500 μm, and 850 μm. Points where ∆Td > 5 K and NH2/ΔNH2 < 1 are excluded, where ∆Td and ∆NH2 refer to the errors on the respective fit parameter. Overlaid are contour lines of the 850 μm continuum emission from Fig. 1.

In the text

	Fig. 5 POS magnetic field strength $B_{\mathrm{P}\mathrm{O}\mathrm{S},\gamma =\pi /4}$ as a function of the separation distance ℓ for the L43E (solid lines) and RNO 91 (dashed lines) regions. Polarization measurements are taken at 154 μm (blue), 450 μm (green) and 850 μm (red).
In the text

Fig. 6 2D histograms of the effective slope bl of the normalized polarization spectrum, as a function of LOS temperature (left panels) and column density (right panels), respectively, colored by the count of data points in each bin. Top panels show results obtained with the ISRF as the only considered radiation source (Setup one), whereas bottom panels show results where the ISRF as well as an offset stellar source are considered (Setup two, see Sect. 5.1 for details). Pearson r, as well as two-tailed p values (ptt), are calculated from 50 segmented median values in each distribution.

In the text

Fig. 7 2D histograms of the effective slope bl of the normalized polarization spectrum, as a function of the standard deviation of the LOS temperature profile σT, colored by the count of data points in each bin. Top panels show results obtained with the ISRF as the only considered radiation source (Setup one), whereas bottom panels show results where the ISRF as well as an offset stellar source are considered (Setup two, see Sect. 5.1 for details). Pearson r, as well as two-tailed p values (ptt), are calculated from 50 segmented median values in each distribution.

In the text

	Fig. A.1 Histograms of the degree of polarization (top panels) and the by 90° rotated polarization angles (bottom panels) for the three observing wavelengths at 154 μm (blue), 450 μm (green) and 850 μm (red). The dashed black line indicates the mean, and the gray shaded area the standard deviation of the distribution. Left panels: L43E region. Right panels: RNO 91 region.
In the text

	Fig. B.1 Grain-size-dependent temperatures for differing levels of attenuation (solid lines) along with their mass-weighed averages (dashed lines), calculated with Eq. B.1. Left: Silicate grains. Right: Carbonaceous grains.
In the text

	Fig. B.2 Semi-analytical polarization spectra for varying levels of relative abundance and attenuation between the hot and the cold phases. Exemplary polarization spectra, where the cold phase retains moderate levels of alignment (left panel), and where the alignment in the cold phase is lost (right panel).
In the text

	Fig. C.1 Same as Fig. 6, but with the true slope albl and without the use of normalization.
In the text

	Fig. C.2 Same as Fig. 7, but with the true slope albl and without the use of normalization.
In the text