# A three-dimensional reconstruction of the interstellar magnetic field toward a star-forming region

K. Ferrière, L. Montier, J.-S. Carrière

► To cite this version:

K. Ferrière, L. Montier, J.-S. Carrière. A three-dimensional reconstruction of the interstellar magnetic field toward a star-forming region. Astronomy & Astrophysics - A&A, 2026, 705, pp.A26. 10.1051/0004-6361/202451034. [hal-05521168](https://hal.science/hal-05521168v1)

HAL Id: hal-05521168

https://hal.science/hal-05521168v1

Submitted on 20 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A three-dimensional reconstruction of the interstellar magnetic field toward a star-forming region

K. Ferrière*, L. Montier, and J.-S. Carrière

IRAP, Université de Toulouse, CNRS, 9 avenue du Colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France

Received 7 June 2024 / Accepted 29 September 2025

## ABSTRACT

**Context.** The polarized thermal emission from interstellar dust offers a valuable tool for probing both the dust and the magnetic field in the interstellar medium (ISM). However, existing observations only yield the total amount of dust emission along the line of sight (LoS), with no information on its LoS distribution.

**Aims.** We present a new method designed to give access to the LoS distribution of the dust emission, both in terms of intensity and polarization.

**Methods.** We relied on three kinematic gas tracers (H<sub>I</sub>, <sup>12</sup>CO, and <sup>13</sup>CO emission lines) to identify the different clouds present along the LoS. We decomposed the measured intensity of the dust emission, I<sub>d</sub>, into separate contributions from these clouds. We performed a similar decomposition of the measured Stokes parameters for linear polarization, Q<sub>d</sub> and U<sub>d</sub>, to derive the polarization parameters of the different clouds, and from this we inferred the clouds' magnetic field orientations.

**Results.** We applied our method to a 3 deg<sup>2</sup> region of the sky, centered on (l, b) = (139°30', -3°16') and exhibiting signs of star formation activity. We found this region to be dominated by an extended and bright cloud with nearly horizontal magnetic field, as expected from the nearly vertical polarization angles measured by Planck. More importantly, we detected the presence of two smaller, depolarizing molecular clouds with very different magnetic field orientations in the plane of the sky (≈65° and ≈45° from the horizontal). This is a novel and viable result, which cannot be directly read off the Planck polarization maps.

**Conclusions.** The application of our method to the G139 region convincingly demonstrates the need to complement 2D polarization maps with 3D kinematic information when looking for reliable estimates of magnetic field orientations.

**Key words.** ISM: clouds – dust, extinction – ISM: kinematics and dynamics – ISM: magnetic fields

## 1. Introduction

The sky offers a 2D view of the cosmos and one of the main challenges of observational astrophysics is to gain access to the line-of-sight (LoS) dimension. The latter is crucial to uncover the true physical properties and understand the exact inner workings of our cosmic environment. For instance, in the case of our Galaxy, access to the LoS dimension makes it possible to retrieve the 3D distribution of interstellar matter (Lallement et al. 2019, 2022; Green et al. 2019; Chen et al. 2019; Leike & Enßlin 2019; Leike et al. 2020; Hottier et al. 2020; Zucker et al. 2021; O'Neill et al. 2024); in turn, this paves the way to studying the dynamics of the interstellar medium (ISM), the formation and evolution of interstellar structures, and, ultimately, the whole cycle of matter between stars and the ISM. It is evident that 2D images alone are insufficient to reach that goal, as they afford only a partial, and possibly misleading, perspective on the observed medium. In general, integration along the LoS results in a loss of information on the spatial variations of physical properties such as temperature, volume density, emissivity, and so on. The loss of information is more severe in the case of vectorial quantities, such as polarization vectors, which can add up either constructively or destructively. As a result, LoS integration can cause depolarization of an intrinsically polarized signal; for instance, when the magnetic field orientation varies along the LoS (Planck Collaboration 2015b, 2017; Clark 2018; Pelgrims et al. 2021) or in the presence of Faraday rotation (Sokoloff et al. 1998;

Beck 2001). Another issue is that projecting onto the plane of the sky (PoS) can distort our perception of certain classes of objects (cores, clumps, filaments, etc.) as well as our estimation of their geometric characteristics (size, aspect ratio, relative orientation angle, etc.) and, accordingly, bias the related statistical analyses (e.g., Planck Collaboration 2016a; Padoan et al. 2023).

There are a range of tools available to help probe the LoS dimension. For instance, turning to spectroscopy, we can infer, from the profile of a spectral line, the brightness distribution as a function of radial velocity (RV) and then we can use the RV as a proxy for LoS coordinate. The direct output is a 3D map of the line brightness as a function of position in the sky and RV, known as a spectral cube. To convert from RV to LoS coordinate (and, thus, obtain a physical cube), we can rely on objects whose distances can be measured directly (e.g., via the parallax) or indirectly (e.g., via bracketing or standard candles).

Observations of the Galactic magnetic field are no exception, as they are also plagued by the lack of information along the LoS. The classical observational methods based on Faraday rotation, synchrotron emission, and dust polarization provide only LoS-integrated quantities, with no information on how the integrant (Faraday rotation rate, synchrotron emissivity, or dust emissivity) varies along the LoS. Several avenues have been proposed to (partially) overcome this problem and, thus, gain (partial) access to the 3D structure of the Galactic magnetic field. Here, we note three possible approaches, based on polarization of starlight, Faraday rotation of pulsar signals, and Faraday tomography, respectively.

<footnote>* Corresponding author: katia.ferriere@irap.omp.eu</footnote>
The polarization of starlight presumably results from its interaction with interstellar dust grains that are aligned by the interstellar magnetic field. The measured polarization orientation directly gives the orientation of $B_{\perp}^1$, while the measured polarization fraction can, under certain conditions, provide an estimate of the inclination of $B$ to the PoS. If $B$ varies along the LoS, the inferred orientation and inclination angles are dust-weighted averages between the considered star and the observer. By considering a large number of stars distributed in space and combining their measured polarization parameters with their measured distances (e.g., from Gaia), we can reconstruct the orientation of $B_{\perp}$ (or possibly $B$) in 3D. This type of tomographic mapping of $B_{\perp}$ was performed toward the Perseus molecular cloud (Doi et al. 2021), toward the Southern Coalsack dark nebula (Versteeg et al. 2024) and toward a portion of the Sagittarius spiral arm (Doi et al. 2024). For the larger datasets from the optical polarization survey PASIPHAE, an automated, Bayesian LoS-inversion method was developed by Pelgrims et al. (2023) and later applied by Pelgrims et al. (2024) to a 4 deg² region of the sky, centered on $(l, b) = (103.3^{\circ}, 22.3^{\circ})$.

Faraday rotation of the linearly polarized radio waves emitted by Galactic pulsars occurs in ionized regions of the ISM, i.e., mainly in the warm ionized medium (WIM), where radio waves interact with the thermal electrons of the medium. The angle by which the polarization orientation rotates is equal to the observing wavelength squared, $\lambda^2$, times the rotation measure, $RM \propto \int_0^L n_e B_\parallel ds$, where $n_e$ is the thermal-electron density and $L$ the path length from the pulsar to the observer. Aside from the RM of a pulsar, we can also measure its dispersion measure, $DM = \int_0^L n_e ds$. The ratio of RM to DM directly yields the $n_e$-weighted average value of $B_\parallel$ between the pulsar and the observer. By combining the average values of $B_\parallel$ toward all pulsars with measured RM and DM (roughly 1500 at the present time) with their measured distances, we can map out the large-scale 3D distribution of $B_\parallel$ (e.g., Han et al. 2018; Sobey et al. 2019).

Faraday tomography also relies on Faraday rotation, but instead of considering the Faraday rotation of the linearly polarized radiation from a background radio source, we can exploit the $\lambda^2$-dependent Faraday rotation of the synchrotron radiation from the Galaxy itself (Burn 1966; Brentjens & de Bruyn 2005). More precisely, we can measure the Galactic polarized intensity at many different radio wavelengths and then convert its variation with $\lambda^2$ into a variation with LoS coordinate. The standard output of Faraday tomography is a so-called Faraday cube, namely, a 3D map of the synchrotron polarized emission as a function of position in the sky and Faraday depth, which is the equivalent of physical depth measured in terms of Faraday rotation. In practice, Faraday tomography can be used to separate synchrotron-emitting regions located at different Faraday depths and to estimate their respective synchrotron polarized intensities, which, in turn, can lead to the strength and the orientation of their $B_\perp$. Faraday tomography can also be used to uncover intervening Faraday screens and to estimate their Faraday thicknesses, which, in turn, can lead to their $B_\parallel$. The method is particularly interesting when the uncovered Faraday screens can be identified with known gaseous structures because it then offers a new way of probing their magnetic fields. Faraday tomography was successfully applied to several regions of the sky, including small fields centered on the nearby galaxy IC 342 (Van Eck et al. 2017)

and the extra-galactic point source 3C 196 (Turić et al. 2021), as well as a much larger area toward the high-latitude outer Galaxy (Erceg et al. 2022, 2024).

Aside from these three classical methods, Hu & Lazarian (2023) proposed a more indirect approach to map out the 3D distribution of the orientation and strength of $B_\perp$, based on the application of the velocity gradient and two Mach number techniques to HI spectroscopic observations.

In this paper, we propose a new 3D polarimetric approach, which combines two of the well-proven observational tools described above, namely, 3D spectral cubes and dust polarization. Compared to the starlight polarimetric approach, which utilizes the polarization of stars with measured distances, we rely on the polarization of the thermal emission from interstellar dust, which we connect to kinematic gas tracers with measured spectral cubes. Thus, our method combines two complementary datasets: spectral cubes of atomic (HI) and molecular ($^{12}$CO and $^{13}$CO) gas tracers and polarization maps of the dust emission. The former contain the kinematic information needed to locate the dust-emitting structures along the LoS, albeit in terms of RV rather than physical distance, and the latter provide the polarization information needed to reconstruct the magnetic field orientation.

In Sect. 2, we lay out the general method. In Sect. 3, we apply our method to a $104' \times 104'$ area of the sky containing a star-forming region. In Sect. 4, we discuss our results and present our conclusions.

## 2. General method

In this section, we introduce the general equations needed in our study and we explain how we can identify the different dust-emitting clouds along the LoS and estimate their magnetic field orientations. In Sect. 2.1, we present the basic equations describing the polarized dust emission. In Sect. 2.2, we introduce our atomic (HI) and molecular ($^{12}$CO and $^{13}$CO) kinematic gas tracers, combine the latter into a single molecular (CO) tracer, and connect the dust emission to the HI and CO tracers via conversion factors. In Sect. 2.3, we explain how the measured intensity of the dust emission can be decomposed along the LoS into the contributions from different clouds identified with the help of the kinematic gas tracers. In Sect. 2.4, we show how the polarization parameters of the different clouds can be fitted to the observed polarization maps of the dust emission and used to infer the orientations of their internal magnetic fields, which we assume to be uniform. Along the way, we make it clear that our LoS decomposition has degeneracies, with implications for the validity of our derived magnetic field orientations. We also propose a convergence test to identify the magnetic field orientations that are truly reliable.

### 2.1. Polarized dust emission

The polarized thermal emission from interstellar dust, integrated along the LoS, can be described by three quantities: the intensity,

$$I_d = \int_0^\infty \varepsilon_d(r) dr, \qquad (1)$$

and the two Stokes parameters for linear polarization,

$$Q_d = \int_0^\infty p_d^{\text{loc}}(r) \varepsilon_d(r) \cos(2\psi_d^{\text{loc}}(r)) dr \qquad (2)$$

<footnote n="1">¹ In this paper, the magnetic field vector is denoted by $B$, its component in the plane of the sky (PoS) by $B_\perp$, and its component along the line of sight (LoS) by $B_\parallel$.</footnote>
and

$$U_d = \int_0^\infty p_d^{\text{loc}}(r) \mathcal{E}_d(r) \sin(2\psi_d^{\text{loc}}(r)) dr, \quad (3)$$

where r denotes the LoS distance from the observer, subscript d stands for dust, $\mathcal{E}_d(r)$ is the local emissivity of the dust thermal emission, $p_d^{\text{loc}}(r)$ the local polarization fraction, and $\psi_d^{\text{loc}}(r)$ the local polarization angle (increasing counterclockwise from Galactic north). Equations (2) and (3) can be combined into a single equation for the complex polarized intensity,

$$\mathcal{P}_d \equiv Q_d + i U_d = \int_0^\infty p_d^{\text{loc}}(r) \mathcal{E}_d(r) e^{2i\psi_d^{\text{loc}}(r)} dr. \quad (4)$$

The norm of the complex polarized intensity is generally referred to as the (real) polarized intensity,

$$P_d = |\mathcal{P}_d| = \sqrt{Q_d^2 + U_d^2}. \quad (5)$$

The local polarization fraction can be written as

$$p_d^{\text{loc}} = (p_d^{\text{loc}})_{\text{max}} \cos^2 \gamma_B, \quad (6)$$

where $\gamma_B$ is the inclination angle of the local magnetic field to the PoS, $(p_d^{\text{loc}})_{\text{max}} = (p_{\text{intr}} R)$ is the theoretical maximum polarization fraction that can be achieved locally, $p_{\text{intr}}$ is the intrinsic polarization fraction, which depends on the shape, elongation, and material of dust grains, and $R$ is the Rayleigh reduction factor, which accounts for the imperfect alignment of dust grains<fnref n="2" /> (see, e.g., Planck Collaboration 2015a, and references therein). The local polarization angle is related to the orientation angle of the local magnetic field in the PoS, $\psi_B$, via

$$\psi_d^{\text{loc}} = \psi_B \pm 90^{\circ}. \quad (7)$$

The ±90° term in Eq. (7) arises because both orientation angles are defined within a 180° range, for instance, in the range [-90°, +90°].

The LoS-averaged polarization fraction is given by

$$p_d^{\text{los}} = \frac{P_d}{I_d} = \frac{\sqrt{Q_d^2 + U_d^2}}{I_d}, \quad (8)$$

and the LoS-averaged polarization angle by

$$\psi_d^{\text{los}} = \frac{1}{2} \arctan \left( \frac{U_d}{Q_d} \right), \quad (9)$$

with arctan the two-argument arctangent function defined from -180° to +180°. Equations (8) and (9) are equivalent to the pair of equations

$$Q_d = p_d^{\text{los}} I_d \cos(2\psi_d^{\text{los}}) \quad (10)$$

and

$$U_d = p_d^{\text{los}} I_d \sin(2\psi_d^{\text{los}}), \quad (11)$$

which can be rewritten in complex form as

$$\mathcal{P}_d \equiv Q_d + i U_d = p_d^{\text{los}} I_d e^{2i\psi_d^{\text{los}}}. \quad (12)$$

The link between the LoS-averaged polarization fraction and angle and their local counterparts is easily obtained by equating Eq. (12) to Eq. (4), while keeping Eq. (1) in mind,

$$p_d^{\text{los}} = \frac{\left| \int_0^\infty p_d^{\text{loc}}(r) \mathcal{E}_d(r) e^{2i\psi_d^{\text{loc}}(r)} dr \right|}{\int_0^\infty \mathcal{E}_d(r) dr} \quad (13)$$

and

$$e^{2i\psi_d^{\text{los}}} = \frac{\int_0^\infty p_d^{\text{loc}}(r) \mathcal{E}_d(r) e^{2i\psi_d^{\text{loc}}(r)} dr}{\left| \int_0^\infty p_d^{\text{loc}}(r) \mathcal{E}_d(r) e^{2i\psi_d^{\text{loc}}(r)} dr \right|} \quad (14)$$

The physical meaning of Eqs. (13) and (14) is pretty straightforward: the LoS-averaged polarization fraction, $p_d^{\text{los}}$, is a dust-weighted LoS average of the local polarization fraction, $p_d^{\text{loc}}(r)$, reduced by a LoS depolarization factor, $e^{2i\psi_d^{\text{loc}}(r)}$, due to fluctuations in the local polarization angle; the LoS-averaged polarization orientation, defined by the LoS-averaged polarization angle, $\psi_d^{\text{los}}$, is a dust-weighted LoS average of the local polarization orientation, defined by the local polarization angle, $\psi_d^{\text{loc}}(r)$. In reality, observations do not strictly capture a single LoS, but rather a whole telescope beam. Therefore, in practice, the LoS integrals in Eqs. (13) and (14) are actually integrals over a telescope beam, and the associated averaging and depolarization actually occur both along the LoS and across the telescope beam. Observational values of $p_d^{\text{los}}$ are discussed in Appendix B.

## 2.2. Kinematic gas tracers and conversion factors

### 2.2.1. Kinematic gas tracers

It is often implicitly assumed that the local polarization fraction and angle are uniform along the LoS through (most of) the dust-emitting region. This assumption makes it possible to infer their values from the observed LoS-averaged polarization fraction and angle, respectively, $p_d^{\text{loc}}(r) = p_d^{\text{los}}$ and $\psi_d^{\text{loc}}(r) = \psi_d^{\text{los}}$. In reality, however, the LoS is likely to intersect several clouds with different polarization properties. In that case, the polarization fraction and angle of the different clouds cannot be immediately inferred from the observed dust emission. Other tracers, such as kinematic gas tracers, are needed to estimate them separately.

Here, we rely on spectral cubes, namely, 3D maps in $(l, b, v)$ space<fnref n="3" />, of the brightness temperature, $T_b$, of the HI 21 cm, $^{12}$CO (1–0) 2.6 mm, and $^{13}$CO (1–0) 2.7 mm emission lines, based on the notion that the HI line traces atomic gas, while the $^{12}$CO and $^{13}$CO lines together trace molecular gas. In principle, we should also include a spectral cube of the ionized gas, but we assume that the contribution from ionized gas is negligible in the region of interest.

The brightness temperature, $T_b$, is affected by optical depth effects, undergoing gradual saturation with increasing opacity. Correcting for opacity saturation is a difficult and subtle

<footnote n="2">² There is no consensus in the community on the exact definition of the intrinsic polarization fraction, $p_{\text{intr}}$. For some authors, $p_{\text{intr}}$ is the local polarization fraction achieved when **B** lies in the PoS (the quantity that we denote $(p_d^{\text{loc}})_{\text{max}}$). For others, $p_{\text{intr}}$ is the local polarization fraction that would be achieved if, in addition to **B** lying in the PoS, dust grains were perfectly aligned (the definition adopted here).</footnote>

<footnote n="3">³ Throughout the paper, $l$ = Galactic longitude; $b$ = Galactic latitude; $v$ = radial velocity (RV); $d$ = distance.</footnote>
task, which we address in Appendix A. There, we define an opacity-corrected brightness temperature, $T_{b*}$, and we derive the equation relating $T_{b*}$ to $T_b$ (Eq. (A.5)). Since this equation involves the excitation temperature, $T$, we discuss the choice of the value of $T$. In brief, we argue in favor of choosing $T = 80$ K for HI (see paragraph preceding Eq. (A.13)) and $T = (T_b^{12CO})_{max}$ for CO (see paragraph preceding Eq. (A.20)). Once the values of $T$ are set, we can convert the spectral cubes of $T_b^{HI}$, $T_b^{12CO}$, and $T_b^{13CO}$ into spectral cubes of $T_{b*}^{HI}$, $(T_{b*}^{12CO})_1$, and $(T_{b*}^{12CO})_2$, respectively, where $(T_{b*}^{12CO})_1$ and $(T_{b*}^{12CO})_2$ are two complementary estimates of $T_{b*}^{12CO}$. We can then combine the $(T_{b*}^{12CO})_1$ and $(T_{b*}^{12CO})_2$ spectral cubes into the spectral cube of a best-estimate $T_{b*}^{12CO}$ (Eq. (A.29)). From now on, superscript 12 in $^{12}$CO is dropped for notational simplicity.

The spectral cubes of $T_{b*}^{HI}$ and $T_{b*}^{CO}$ generally have different angular and spectral resolutions, so they first need to be brought to a common angular resolution, $\delta\theta$, a common spectral resolution, $\delta v$, and a common $(l, b, v)$ grid with, say, $n_1 \times n_2$ pixels. The maps of $I_d$, $Q_d$, and $U_d$ (see Sect. 2.1) also need to be brought to the angular resolution $\delta\theta$ and the $n_1 \times n_2$ pixel grid.

We assume that our $T_{b*}^{HI}$ and $T_{b*}^{CO}$ cubes are complementary, in the sense that together they account for all the gas along the LoS, with no omission and no overlap. It then follows that the hydrogen column density, $N_H$, can be decomposed into the contributions from the gas components probed by HI and CO,

$$N_H = N_H^{HI} + N_H^{CO}. \qquad (15)$$

The intensity of the dust emission can be decomposed in a similar way,

$$I_d = I_d^{HI} + I_d^{CO}. \qquad (16)$$

We make the standard assumption that for each gas tracer A (A = HI or CO), the intensity of the dust emission associated with that tracer, $I_d^A$, is simply proportional to the hydrogen column density probed by that tracer, $N_H^A$ (e.g., Hildebrand 1983). If we denote the conversion factor from hydrogen column density to dust intensity by $X_{I_d/N_H}^A$, we can then write

$$I_d^A = X_{I_d/N_H}^A N_H^A. \qquad (17)$$

We note that HI and CO are expected to have slightly different conversion factors, mostly because dust grains have different properties (composition, size, emissivity) in the atomic and molecular media. We further assume that for each tracer A, the hydrogen column density, $N_H^A$, is proportional to the opacity-corrected brightness temperature, $T_{b*}^A(v)$, integrated over the RV (see Eq. (A.7) in Appendix A),

$$N_H^A = X_{N_H/T_b}^A \int T_{b*}^A(v) dv, \qquad (18)$$

where $X_{N_H/T_b}^A$ is the conversion factor from velocity-integrated opacity-corrected brightness temperature to hydrogen column density. Combining Eqs. (17) and (18) leads to

$$I_d^A = X_{I_d/T_b}^A \int T_{b*}^A(v) dv, \qquad (19)$$

with $X_{I_d/T_b}^A = X_{I_d/N_H}^A X_{N_H/T_b}^A$. If we now insert Eq. (19) into Eq. (16), we finally obtain for the dust intensity

$$I_d = X_{I_d/T_b}^{HI} \int T_{b*}^{HI} dv + X_{I_d/T_b}^{CO} \int T_{b*}^{CO} dv. \qquad (20)$$

### 2.2.2. Conversion factors

The observables in Eq. (20) are the dust intensity, $I_d$, and the opacity-corrected brightness temperatures, $T_{b*}^{HI}$ and $T_{b*}^{CO}$. The conversion factors, $X_{I_d/T_b}^{HI}$ and $X_{I_d/T_b}^{CO}$, are not directly observable, and we did not find any values for them in the literature. We did find separate estimates for some of the intermediate conversion factors, $X_{I_d/N_H}$ and $X_{N_H/T_b}$ (see Sect. 3.2), but they were obtained for restricted regions of the sky and they tend to have wide scatter, so they cannot be directly applied to our region of interest.

Here, we treat the conversion factors, $X_{I_d/T_b}^{HI}$ and $X_{I_d/T_b}^{CO}$, as free parameters, and we derive their best-fit values in the considered region by minimizing the reduced $\chi^2$, defined by

$$\chi_r^2 = \frac{1}{(n_{\text{dat}} - n_{\text{par}})} \sum_{n_{\text{pix}}} \frac{[I_d^{\text{obs}} - (I_d^{\text{HI}} + I_d^{\text{CO}})]^2}{\sigma_I^2}, \qquad (21)$$

where $I_d^{\text{obs}}$ is the observed dust intensity, $I_d^{HI}$ and $I_d^{CO}$ are the dust intensities associated with HI and CO, respectively (Eq. (19) with A = HI and CO), $\sigma_I$ is the total observational uncertainty, $n_{\text{dat}} = n_{\text{pix}} = n_1 \times n_2$ is the total number of data points, $n_{\text{par}} = 2$ is the number of free parameters, and the sum runs over the $n_{\text{pix}}$ pixels. The total observational uncertainty is the quadratic sum of the measurement errors in $I_d^{\text{obs}}$, $I_d^{HI}$, and $I_d^{CO}$,

$$\sigma_I^2 = \sigma^2(I_d^{\text{obs}}) + \sigma^2(I_d^{\text{HI}}) + \sigma^2(I_d^{\text{CO}}), \qquad (22)$$

where

$$\sigma^2(I_d^A) = (X_{I_d/T_b}^A)^2 \left( \sum_v \sigma^2(T_{b*}^A) \right) \delta v^2 \qquad (23)$$

for A = HI and CO, $\sigma^2(T_{b*}^{HI})$ and $\sigma^2(T_{b*}^{CO})$ are given by Eqs. (A.13) and (A.30), respectively, $\delta v$ is the spectral resolution, and the sum runs over all the velocity bins.

Minimization of $\chi_r^2$ is performed through Markov chain Monte Carlo (MCMC) simulations. Since MCMC simulations alone consider only measurement errors, which in the case at hand are potentially dominated by modeling errors (namely, errors in Eq. (20) with constant conversion factors), they are likely to underestimate the uncertainties in the best-fit parameters. To obtain meaningful uncertainties, $\sigma(X_{I_d/T_b}^{HI})$ and $\sigma(X_{I_d/T_b}^{CO})$, we resort to parametric bootstrap sampling, with each bootstrap sample involving MCMC simulations.

### 2.3. LoS decomposition into clouds

Consider a given LoS and assume that the dust emission measured on that LoS arises from $n_{\text{cl}}$ distinct clouds. Here, the term “cloud” is to be understood in a broad sense, which may include intercloud regions. The intensity of the dust emission can then be modeled as the sum of the contributions from the $n_{\text{cl}}$ clouds,

$$I_d^{\text{mod}} = \sum_{i=1}^{n_{\text{cl}}} I_{d,i}, \qquad (24)$$

and similarly for the two Stokes parameters for linear polarization,

$$Q_d^{\text{mod}} = \sum_{i=1}^{n_{\text{cl}}} Q_{d,i}, \qquad (25)$$
and

$$U_d^{\text{mod}} = \sum_{i=1}^{n_{\text{cl}}} U_{d,i}, \qquad (26)$$

where the subscript i refers to cloud i. The contributions $I_{d,i}$, $Q_{d,i}$, and $U_{d,i}$ from cloud i have the same expressions as the total $I_d$ (Eq. (1)), $Q_d$ (Eq. (2)), and $U_d$ (Eq. (3)), respectively, with the LoS integral over an infinite path length replaced by a LoS integral over the path length through cloud i, $L_i$.

Next, we use the spectral cubes of the HI and CO opacity-corrected brightness temperatures, $T_{b\star}^{\text{HI}}$ and $T_{b\star}^{\text{CO}}$, introduced in Sect. 2.2 to identify the different clouds along the LoS.

### 2.3.1. ROHSA decomposition

In a first step, we apply the algorithm ROHSA (Regularized Optimization for Hyper-Spectral Analysis) developed by Marchal et al. (2019) to our $T_{b\star}^{\text{HI}}$ and $T_{b\star}^{\text{CO}}$ spectral cubes separately. This algorithm is designed to decompose a spectral data cube, say, a cube of $T_{b\star}(l, b, v)$, into several spatially coherent Gaussian kinematic components. The Gaussian decomposition is optimized over the entire data cube at once, with the requirement that the solution must be spatially smooth.

In mathematical terms, the observed $T_{b\star}(l, b, v)$ is approximated by a modeled $\tilde{T}_{b\star}(l, b, v)$, equal to the sum of $n_G$ Gaussian components,

$$\tilde{T}_{b\star}(l, b, v) = \sum_{j=1}^{n_G} \tilde{T}_{b\star,j}(l, b, v), \qquad (27)$$

where each component j is described by an amplitude $A_j(l, b)$, a mean RV, $\bar{v}_j(l, b)$, and a standard deviation, $\sigma_j(l, b)$,

$$\tilde{T}_{b\star,j}(l, b, v) = A_j(l, b) \exp \left[ - \frac{(v - \bar{v}_j(l, b))^2}{2 \sigma_j^2(l, b)} \right]. \qquad (28)$$

The best-fit values of the $3n_G$ Gaussian parameters, $A_j$, $\bar{v}_j$, and $\sigma_j$, are derived through minimization of a cost function that includes the standard $\chi^2$ term plus a regularization term meant to ensure a spatially smooth solution.

ROHSA has several free parameters (six in the initial version described in Marchal et al. (2019), and more in the present online version). We keep the default values of these parameters, except for the three hyper-parameters entering the regularization term, $\lambda_A$, $\lambda_{\bar{v}}$, and $\lambda_{\sigma}$. Since our ultimate purpose is to construct smooth, coherent clouds, we need to impose strong constraints of spatial coherence on $A_j$, $\bar{v}_j$, and $\sigma_j$, which is done by choosing large values for $\lambda_A$, $\lambda_{\bar{v}}$, and $\lambda_{\sigma}$. After verifying that the exact values are not critical, we adopt $\lambda_A = \lambda_{\sigma} = 1000$ and $\lambda_{\bar{v}} = 100$ for our application in Sect. 3. Furthermore, since we want to retain only truly physical components, we discard the extracted components whose velocity-integrated $\tilde{T}_{b\star,j}$ fall below twice the noise level at every pixel.

### 2.3.2. Cloud reconstruction

The second step is to reconstruct the different clouds along the LoS. To that end, we collect all the Gaussian components from both HI and CO, and we group together the components that have similar velocity profiles. Since the very definition and the exact outline of an interstellar cloud are moot points, the criterion we choose to rely on is necessarily subjective.

To start with, there is a total of $n_G^{\text{tot}} = n_G^{\text{HI}} + n_G^{\text{CO}}$ Gaussian components, which we consider two by two. For every pair of components j and j' ($j \neq j'$), we compute a velocity correlation coefficient at each pixel, $(l, b)$,

$$C_{jj'}(l, b) = \frac{\sum_{v} \tilde{T}_{b\star,j} \tilde{T}_{b\star,j'}}{\sqrt{\sum_{v} \tilde{T}_{b\star,j}^2} \sqrt{\sum_{v} \tilde{T}_{b\star,j'}^2}}, \qquad (29)$$

where the sums run over all the velocity bins. Clearly, $C_{jj'}$ depends only on the velocity profiles (mean velocities, $\bar{v}_j$ and $\bar{v}_{j'}$, and standard deviations, $\sigma_j$ and $\sigma_{j'}$) of both components, not on their amplitudes ($A_j$ and $A_{j'}$). We then compute the average value of $C_{jj'}(l, b)$ over all the pixels, weighted by the sum of the dust intensities of components j and j'.

$$\langle C_{jj'} \rangle = \frac{\sum_{n_{\text{pix}}} (\tilde{I}_{d,j} + \tilde{I}_{d,j'}) C_{jj'}}{\sum_{n_{\text{pix}}} (\tilde{I}_{d,j} + \tilde{I}_{d,j'})}, \qquad (30)$$

where $\tilde{I}_{d,j}$ is the dust intensity of component j (defined by Eq. (33) below) and the sum runs over the $n_{\text{pix}}$ pixels.

We consider that components j and j' are correlated, and thus parts of a same cloud, if $\langle C_{jj'} \rangle$ lies above a certain threshold, $C_{\text{th}}$:

$$\langle C_{jj'} \rangle \geq C_{\text{th}}. \qquad (31)$$

It is important to realize that this is a purely kinematic criterion, independent of any possible spatial correlation. This choice is motivated by the fact that while two components of a given cloud are expected to have similar velocity profiles, they do not need to be correlated in space; for instance, they could very well be co-spatial (e.g., a region containing both HI and CO) or adjacent (e.g., an HI envelope around a CO core).

A critical issue concerns the choice of the value of $C_{\text{th}}$. For our application in Sect. 3, we tested all the integer values of $C_{\text{th}}$ between 0 and 100%. Based on the results (summarized in Table 3), we decided to present a detailed analysis in the case $C_{\text{th}} = 50\%$ (Sects. 3.3–3.4), considering that this value strikes a good balance between preserving physical clouds and separating distinct clouds along the LoS. The results obtained for other values of $C_{\text{th}}$ as well as their sensitivity to the exact value of $C_{\text{th}}$ are discussed in Sect. 3.5.

Pairs satisfying Eq. (31) are themselves grouped together into one multicomponent cloud if they possess a component in common. For instance, if pairs jk and kl satisfy $\langle C_{jk} \rangle \geq C_{\text{th}}$ and $\langle C_{kl} \rangle \geq C_{\text{th}}$, they are grouped into one cloud enclosing components j, k, and l, even if $\langle C_{jl} \rangle < C_{\text{th}}$.

All the remaining Gaussian components, i.e., the components that are not correlated with any other component (in the sense of satisfying Eq. (31)), are considered to each form a separate cloud.

Altogether, we end up with $n_{\text{cl}}$ distinct clouds made up of one, two, or more Gaussian components. The dust intensity of cloud i, $I_{d,i}$, can be written as a sum over the $n_{G,i}$ Gaussian components of cloud i:

$$I_{d,i} = \sum_{j} \tilde{I}_{d,j}, \qquad (32)$$
where the dust intensity of component j, $\tilde{I}_{d,j}$, is related to its opacity-corrected brightness temperature, $\tilde{T}_{b*,j}$ (Eq. (28)), through an equation similar to Eq. (19):

$$ \tilde{I}_{d,j} = X_{I_d/T_b}^A \int \tilde{T}_{b*,j}(v) dv. \quad (33) $$

Superscript A in the conversion factor $X_{I_d/T_b}^A$ refers to the gas tracer (HI or CO) associated with component j. As a reminder, the best-fit value of $X_{I_d/T_b}^A$ and the attendant uncertainty, $\sigma(X_{I_d/T_b}^A)$, are determined as explained in Sect. 2.2.

## 2.4. Derivation of the magnetic field orientation in each cloud

At this point, we have identified $n_{cl}$ clouds along the LoS, and we have derived their contributions $I_{d,i}$ (Eqs. (32) and (33)) to the dust intensity, $I_d^{\text{mod}}$ (Eq. (24)). We now turn to their contributions $Q_{d,i}$ and $U_{d,i}$ to the two Stokes parameters for linear polarization, $Q_d^{\text{mod}}$ and $U_d^{\text{mod}}$ (Eqs. (25) and (26)). As mentioned below Eq. (26), $I_{d,i}$, $Q_{d,i}$, and $U_{d,i}$ have the same expressions as $I_d$, $Q_d$, and $U_d$ (Eqs. (1), (2), and (3)), respectively, with the LoS integral reduced to the path length through cloud i, $L_i$. By analogy with Eqs. (10)–(14), we can then directly write

$$ Q_{d,i} = p_{d,i} I_{d,i} \cos(2\psi_{d,i}) \quad (34) $$

and

$$ U_{d,i} = p_{d,i} I_{d,i} \sin(2\psi_{d,i}), \quad (35) $$

with

$$ p_{d,i} = \frac{\left| \int_{L_i} p_d^{\text{loc}}(r) \mathcal{E}_d(r) e^{2i\psi_d^{\text{loc}}(r)} dr \right|}{\int_{L_i} \mathcal{E}_d(r) dr} \quad (36) $$

and

$$ e^{2i\psi_{d,i}} = \frac{\int_{L_i} p_d^{\text{loc}}(r) \mathcal{E}_d(r) e^{2i\psi_d^{\text{loc}}(r)} dr}{\left| \int_{L_i} p_d^{\text{loc}}(r) \mathcal{E}_d(r) e^{2i\psi_d^{\text{loc}}(r)} dr \right|} \quad (37) $$

The physical meaning of $p_{d,i}$ and $\psi_{d,i}$ for cloud i is similar to that of $p_d^{\text{los}}$ and $\psi_d^{\text{los}}$ given below Eq. (14) for the entire LoS.

To proceed, we consider that $p_{d,i}$ and $\psi_{d,i}$ are uniform across the PoS surface of cloud i, in contrast to $I_{d,i}$, which generally varies. Strictly speaking, this is unlikely to be correct, as magnetic field lines have probably been distorted by internal motions. However, we are a priori entitled to take this approach if we are only interested in dust-weighted average values of $p_{d,i}$ and $\psi_{d,i}$ at the cloud scale, which we denote with an overscore$^4$. Thus, the pair of equations we work here with is obtained by inserting Eqs. (34) and (35), with $p_{d,i}$ and $\psi_{d,i}$ overscored, into Eqs. (25) and (26):

$$ Q_d^{\text{mod}} = \sum_{i=1}^{n_{\text{cl}}} \bar{p}_{d,i} I_{d,i} \cos(2\bar{\psi}_{d,i}) \quad (38) $$

and

$$ U_d^{\text{mod}} = \sum_{i=1}^{n_{\text{cl}}} \bar{p}_{d,i} I_{d,i} \sin(2\bar{\psi}_{d,i}). \quad (39) $$

In the above equations, the dust intensity of cloud i, $I_{d,i}$, was derived in Sect. 2.3 (Eqs. (32) and (33)), while its average polarization fraction and angle, $\bar{p}_{d,i}$ and $\bar{\psi}_{d,i}$, are treated as free parameters. This leaves us with two free parameters per cloud and hence a total of $2n_{cl}$ free parameters.

The best-fit values of $\bar{p}_{d,i}$ and $\bar{\psi}_{d,i}$ are the values that minimize

$$ \chi_r^2 = \frac{1}{(n_{\text{dat}} - n_{\text{par}})} \sum_{n_{\text{pix}}} \left[ \frac{(Q_d^{\text{obs}} - Q_d^{\text{mod}})^2}{\sigma_Q^2} + \frac{(U_d^{\text{obs}} - U_d^{\text{mod}})^2}{\sigma_U^2} \right], \quad (40) $$

where $Q_d^{\text{obs}}$ and $U_d^{\text{obs}}$ are the observed Stokes parameters, $Q_d^{\text{mod}}$ and $U_d^{\text{mod}}$ are the modeled Stokes parameters given by Eqs. (38) and (39), respectively, $\sigma_Q$ and $\sigma_U$ are the total “observational uncertainties” associated with $Q_d$ and $U_d$, respectively, $n_{\text{dat}} = 2n_{\text{pix}}$ is the total number of data points, $n_{\text{par}} = 2n_{\text{cl}}$ is the number of free parameters, and the sum runs over the $n_{\text{pix}}$ pixels. The total “observational uncertainties” are the quadratic sums of the measurement errors in the observed Stokes parameters, $\sigma(Q_d^{\text{obs}})$ and $\sigma(U_d^{\text{obs}})$, and the “decomposition errors” in the modeled Stokes parameters arising from decomposition errors in the modeled cloud intensities, $I_{d,i}$. The contributions from the different clouds cannot be derived separately, but we may reasonably consider that the “decomposition errors” in $Q_d^{\text{mod}}$ and $U_d^{\text{mod}}$ are both equal to the error in the modeled intensity, $I_d^{\text{mod}}$, times the observed LoS-averaged polarization fraction, $p_d^{\text{los,obs}}$. We may further approximate the error in $I_d^{\text{mod}}$ by the residual $\Delta(I_d) = I_d^{\text{obs}} - I_d^{\text{mod}}$. Altogether, we have

$$ \sigma_Q^2 = \sigma^2(Q_d^{\text{obs}}) + (p_d^{\text{los,obs}})^2 \Delta^2(I_d) \quad (41) $$

and

$$ \sigma_U^2 = \sigma^2(U_d^{\text{obs}}) + (p_d^{\text{los,obs}})^2 \Delta^2(I_d). \quad (42) $$

Minimization of $\chi_r^2$ is performed through MCMC simulations, leading to the best-fit values of $\bar{p}_{d,i}$ and $\bar{\psi}_{d,i}$, together with the attendant uncertainties, $\sigma(\bar{p}_{d,i})$ and $\sigma(\bar{\psi}_{d,i})$.

Once the polarization fraction, $\bar{p}_{d,i}$, and the polarization angle, $\bar{\psi}_{d,i}$, of cloud i have been determined, it is possible to estimate the orientation of its internal magnetic field, $\mathbf{B}_i$. Here, too, we are referring to dust-weighted average values, denoted with an overscore. The orientation angle of $\mathbf{B}_i$ in the PoS, $\bar{\psi}_{B,i}$, can be directly inferred from the polarization angle, $\bar{\psi}_{d,i}$, with the help of Eq. (7) applied to cloud i:

$$ \bar{\psi}_{B,i} = \bar{\psi}_{d,i} \pm 90^{\circ}. \quad (43) $$

The inclination angle of $\mathbf{B}_i$ to the PoS, $\bar{\gamma}_{B,i}$, can in principle be inferred from the polarization fraction, $\bar{p}_{d,i}$, using Eq. (6) applied to cloud i,

$$ \cos^2 \bar{\gamma}_{B,i} = \frac{\bar{p}_{d,i}}{(\bar{p}_{d,i})_{\text{max}}}, \quad (44) $$

<footnote>$^4$ The exact expressions of the dust-weighted average values of $p_{d,i}$ and $\psi_{d,i}$ obtained through minimization of $\chi_r^2$ (Eq. (40)) are derived in Appendix C.</footnote>
<table><caption>**Table 1.** Characteristics of the input maps of the G139 region displayed in the top rows of Figs. 1 and 2.</caption><thead><tr><th>Input map</th><th>δθ<sup>a</sup></th><th>δv (km s<sup>-1</sup>)<sup>b</sup></th><th>Δv (km s<sup>-1</sup>)<sup>c</sup></th><th>Survey/instrument</th><th>Reference</th></tr></thead><tbody><tr><td>Dust 353 GHz</td><td>7'</td><td></td><td></td><td>Planck</td><td>1</td></tr><tr><td>HI 21 cm</td><td>10.8'</td><td>1.44</td><td>[-600, +600]</td><td>EBHIS</td><td>2</td></tr><tr><td><sup>12</sup>CO 2.6 mm</td><td>50"</td><td>0.16</td><td>[-95, +25]</td><td>MWISP</td><td>3</td></tr><tr><td><sup>13</sup>CO 2.7 mm</td><td>50"</td><td>0.17</td><td>[-95, +25]</td><td>MWISP</td><td>4</td></tr></tbody></table>

**Notes.** (a) Angular resolution. (b) Spectral resolution. (c) Spectral extent. **References.** (1) Planck Collaboration (2020b); (2) Winkel et al. (2016); (3) Su et al. (2019); Yuan et al. (2021); (4) Yuan et al. (2022).

together with an adopted value of the maximum polarization fraction averaged at the cloud scale, $(\bar{p}_{d,i})_{\text{max}}$. Following Planck Collaboration (2015a), we could, for instance, take $(\bar{p}_{d,i})_{\text{max}} = 23\%$ (see Appendix B). Because of the uncertainty in $(\bar{p}_{d,i})_{\text{max}}$, the derived value of $\bar{\gamma}_{B,i}$ is much less reliable than the derived value of $\bar{\psi}_{B,i}$. Moreover, the existence of two opposite-signed solutions for $\bar{\gamma}_{B,i}$ implies that the orientation of $\bar{B}_i$ can only be determined with a mirror ambiguity with respect to the PoS. Zee-man observations, which are sensitive to $B_\parallel$, would be needed to set the sign of $\bar{\gamma}_{B,i}$ for the dominant cloud.

### 3. Application to the G139 region

To illustrate the method presented in Sect. 2, we now apply it to a small region of the sky encompassing the Herschel G139 field, which is one of the 116 Galactic fields observed with Herschel as part of the Herschel Galactic cold core (GCC) key-program (Juvela et al. 2010, 2012). This small region, which we refer to as the G139 region, is a 104' × 104' square centered on $(l, b) = (139^\circ 30', -3^\circ 16')$. The Herschel maps reveal a long filamentary structure with signs of star formation activity, including a number of embedded cores and young stellar object candidates (Montillaud et al. 2015). This filamentary structure is surrounded by a more diffuse and extended emission.

In Sect. 3.1, we present the relevant available data. In Sect. 3.2, we derive the conversion factors from gas tracers to dust emission. In Sect. 3.3, we decompose the measured intensity of the dust emission into the contributions from seven separate clouds along the LoS. In Sect. 3.4, we derive the magnetic field orientation in each cloud. In Sect. 3.5, we examine alternative configurations of clouds, involving an increasing number of clouds. The input maps of the G139 region used in our study are listed in Table 1, along with their angular resolution, $\delta\theta$, and (when relevant) their spectral resolution, $\delta v$, and their spectral extent, $\Delta v$.

#### 3.1. Available data for G139

The polarization information needed to infer the magnetic field orientations toward G139 can be extracted from the all-sky maps of the polarized dust emission measured by Planck at 353 GHz (Planck Collaboration 2015b, 2020b). The three panels in the top row of Fig. 1 show the maps of the intensity, $I_d$ (Eq. (1)), and of the two Stokes parameters for linear polarization, $Q_d$ (Eq. (2)) and $U_d$ (Eq. (3)), of the 353 GHz dust emission toward G139.

The kinematic information needed to separate the different clouds along the LoS is provided by the spectral cubes of the brightness temperature, $T_b$, of the three gas tracers introduced in Sect. 2.2, namely, the HI 21 cm, $^{12}$CO 2.6 mm, and $^{13}$CO 2.7 mm emission lines. Here, we resort to the following cubes: the HI

cube from the Effelsberg-Bonn HI Survey (EBHIS) of the whole northern sky, with $\delta \ge -5^\circ$ (Winkel et al. 2016), and the $^{12}$CO and $^{13}$CO cubes from the Milky Way Imaging Scroll Painting (MWISP) survey of the Galactic plane, with $-10^\circ \le l \le +250^\circ$ and $|b| \le 5.2^\circ$ (Su et al. 2019; Yuan et al. 2021, 2022). The corresponding maps of the velocity-integrated $T_b$ toward G139 are displayed in the three panels of the top row of Fig. 2. We see that HI (left panel) is quite uniformly distributed, with a weak northward gradient, while $^{12}$CO (middle panel) and $^{13}$CO (right panel) have more structured distributions, with a prominent peak at the position of the bright core in the dust intensity map (top-left panel of Fig. 1). Interestingly, the $^{12}$CO peak does not stand out as prominently as the $^{13}$CO peak, which is most likely because $^{12}$CO emission from the underlying emitting region is optically very thick.

Following the procedure explained in detail in Appendix A, we correct $T_b$ for opacity saturation, thereby obtaining an opacity-corrected brightness temperature, $T_{b*}$. For HI, we derive $T_{b*}^{\text{HI}}$ using Eq. (A.5) with $T = 80$ K. For CO, we derive two complementary estimates of $T_{b*}^{\text{CO}}$ assuming $T = (T_b^{\text{CO}})_{\text{max}}: (T_b^{\text{CO}})_1$ based on the $^{12}$CO cube (Eq. (A.18)) and $(T_b^{\text{CO}})_2$ based on the rescaled $^{13}$CO cube (Eq. (A.19)); we then combine $(T_b^{\text{CO}})_1$ and $(T_b^{\text{CO}})_2$ to obtain a best-estimate $T_{b*}^{\text{CO}}$ (Eq. (A.29)), from now on referred to as $T_{b*}^{\text{CO}}$. We find that $T_{b*}^{\text{CO}}$ is nearly equal to the estimate with the smaller uncertainty, which is $(T_b^{\text{CO}})_1$ away from the CO peak and $(T_b^{\text{CO}})_2$ toward the CO peak. This dichotomy results from the huge jump in the opacity correction to $T_b^{\text{CO}}$ toward the CO peak, which renders $(T_b^{\text{CO}})_1$ extremely uncertain.

The maps of the velocity-integrated $T_{b*}^{\text{HI}}$ and $T_{b*}^{\text{CO}}$ are displayed in the second row of Fig. 2. The opacity-corrected HI map (left panel) is similar to, but more contrasted than, its observed counterpart (left panel in the first row). It also differs by the emergence of a slight over-intensity along the lower-left boundary, which probably indicates the presence of a small, HI-bright cloud. The opacity-corrected $^{12}$CO and $^{13}$CO maps (not shown) are quite similar, except toward the CO peak, where only $T_{b*}^{\text{CO}}$ is reliable. Altogether, the combined opacity-corrected CO map (right panel in the second row) is very close to the opacity-corrected $^{12}$CO map, with only the CO peak taken from the rescaled opacity-corrected $^{13}$CO map.

Before jointly exploiting the $T_{b*}^{\text{HI}}$ and $T_{b*}^{\text{CO}}$ spectral cubes, we bring them to common angular and spectral resolutions and to a common $(l, b, v)$ grid. We also bring the Planck maps of $I_d$, $Q_d$, and $U_d$ to the common angular resolution and common $(l, b)$ grid. Since the HI cube has the lowest angular resolution (see Table 1), we adopt its angular resolution, $\delta\theta = 10.8'$, together with a pixel size at the (rounded) Nyquist limit, $\delta\theta_{\text{pix}} = 4'$. Accordingly, the $(l, b)$ grid needed to cover our $104' \times 104'$ G139
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Maps of the intensity, *I*<sub>d</sub> (left), and of the two Stokes parameters for linear polarization, *Q*<sub>d</sub> (middle) and *U*<sub>d</sub> (right), of the polarized dust emission at 353 GHz toward the G139 region defined in the first paragraph of Sect. 3. Top row: observational maps from Planck (Planck Collaboration 2020b). Middle row: same maps resampled to the common 26 × 26 pixel grid of the gas tracers (pixel size = 4′). Bottom row: statistical uncertainties in the Planck maps resampled to the 26 × 26 pixel grid. The total uncertainties are equal to the quadratic sums of the statistical uncertainties and the photometric calibration uncertainties, which in turn are given by 0.0078 *I*<sub>d</sub> for *I*<sub>d</sub>(Planck Collaboration 2016b, 2020b) and 0.015 *P*<sub>d</sub> for *Q*<sub>d</sub> and *U*<sub>d</sub> (Planck Collaboration 2020b,c).</figcaption></figure>

region possesses 26 × 26 pixels. Clearly, this grid undersamples
the CO and dust maps. Along the v-axis, we retain a total velocity
range [−120, +20] km s⁻¹, which is broad enough to encompass
all the true emission from HI and CO, and we adopt a spec-
tral resolution δv = 0.3 km s⁻¹, which is much finer than the
spectral resolution of the HI cube and roughly twice the spec-
tral resolution of the CO cube. This spectral oversampling of
T<sub>b★</sub><sup>HI</sup> makes it easier to extract meaningful Gaussian kinematic
components with ROHSA. The re-gridded maps of I<sub>d</sub>, Q<sub>d</sub>, and
U<sub>d</sub> and those of the velocity-integrated T<sub>b★</sub><sup>HI</sup> and T<sub>b★</sub><sup>CO</sup> are dis-
played in the second row of Fig. 1 and the third row of Fig. 2,
respectively.

**3.2. Derivation of the conversion factors**

For each gas tracer A (A = HI or CO), we need to determine
the conversion factor from velocity-integrated opacity-corrected
brightness temperature to dust intensity at 353 GHz, X<sup>A</sup><sub>I<sub>d</sub>/T<sub>b</sub></sub> (see
Eq. (19)). As explained in Sect. 2.2, the best-fit values of X<sup>HI</sup><sub>I<sub>d</sub>/T<sub>b</sub></sub>
and X<sup>CO</sup><sub>I<sub>d</sub>/T<sub>b</sub></sub>, together with their uncertainties, are obtained by min-
imizing χ<sub>r</sub><sup>2</sup> in Eq. (21) through bootstrap sampling + MCMC
simulations.

We find that the best fit has $\chi_r = 1.63$. This small value of $\chi_r$ indicates that our model is globally satisfactory, in the sense that the errors in the reconstruction of the dust intensity with Eq. (20) are only slightly larger than the total observational uncertainty, $\sigma_I$ (Eq. (22)). The latter is dominated by $\sigma(I_d^{HI})$, which generally exceeds $\sigma(I_d^{obs})$ and $\sigma(I_d^{CO})$ by a factor ~2–10, except toward the CO peak, where $\sigma(I_d^{CO}) \sim 2 \sigma(I_d^{HI}) \sim 3.5 \sigma(I_d^{obs})$. The dominant contribution to $\sigma(I_d^{HI})$, in turn, comes by far from the uncertainty in the spin temperature of HI (see second term in the r.h.s. of Eq. (A.13)), while the dominant contribution to $\sigma(I_d^{CO})$ comes from measurement errors (see first terms in the r.h.s. of Eqs. (A.20) and (A.21)).

Displayed in Fig. 3 is a corner plot of the marginal and
joint probability density functions of $X_{I_d/T_b}^{HI}$ and $X_{I_d/T_b}^{CO}$. Their
best-fit values and standard deviations, written above the cor-
responding 1D histograms (upper-left and lower-right panels),
are $X_{I_d/T_b}^{HI} = (0.00081 \pm 0.000006) (\text{MJy sr}^{-1}) (\text{K km s}^{-1})^{-1}$ and
$X_{I_d/T_b}^{CO} = (0.0685 \pm 0.0019) (\text{MJy sr}^{-1}) (\text{K km s}^{-1})^{-1}$, respectively.
The shape of the joint distribution (lower-left panel) indicates
that $X_{I_d/T_b}^{HI}$ and $X_{I_d/T_b}^{CO}$ are strongly anticorrelated. This result can
be explained by the fact that HI and CO both have widespread
spatial distributions, with significant overlap (see third row of
Fig. 2), such that an increase in the dust emission associated
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Top row: observational maps of the velocity-integrated brightness temperatures, $T_b$, of the HI 21 cm (**left**), $^{12}$CO 2.6 mm (**middle**), and $^{13}$CO 2.7 mm (**right**) emission lines toward the G139 region; these maps are from Winkel et al. (2016), Yuan et al. (2021), and Yuan et al. (2022), respectively. Second row: maps of the velocity-integrated opacity-corrected brightness temperatures, $T_{b*}$, of the HI (**left**) and $^{12}$CO (**right**) lines, where $T_{b*}^{12}$CO is derived based on the combined $^{12}$CO and $^{13}$CO data. Third row: same maps resampled to the common 26 × 26 pixel grid (pixel size = 4 arcsec). Fourth row: reconstructed maps after application of the Gaussian decomposition algorithm ROHSA. Bottom row: maps of the residuals obtained by subtracting the reconstructed ROHSA maps from the resampled input maps. The color bars for HI and CO (superscript 12 dropped) can be rescaled to dust intensity at 353 GHz using the best-fit values of the conversion factors derived in Sect. 3.2 (see Fig. 3). The resulting ranges in the third and fourth rows are ≈[0, 3.5 MJy sr⁻¹] for HI and ≈[0, 2.5 MJy sr⁻¹] for CO.</figcaption></figure>

with one tracer must be accompanied by a decrease in the dust emission associated with the other tracer.

The best-fit values of the conversion factors enable us to rescale the HI and CO maps in the third row of Fig. 2 to dust intensity at 353 GHz and thus obtain the best-fit maps of $I_d^{HI}$ and $I_d^{CO}$ as well as the best-fit map of the reconstructed dust intensity, $I_d^{HI} + I_d^{CO}$. The latter is displayed in the top-middle panel of Fig. 4, where it can be compared to the observational *Planck* map, $I_d^{obs}$, in the left panel. Both maps look similar, and the bright regions in the *Planck* map are well reproduced. The HI medium provides the general dust-emission background, including its northward

gradient, while the CO medium is responsible for the bright core to the right and for the weaker dust enhancement that runs obliquely north of the southeast-northwest diagonal. Globally, the HI and CO media account for 84.3% and 15.7%, respectively, of the dust emission from the G139 region.

The map of the residuals, $[I_d^{obs} - (I_d^{HI} + I_d^{CO})]$, is shown in the top-right panel of Fig. 4. The largest residuals are positive and arise on the northwest and southeast sides of the bright core (which itself is almost residual-free). These positive residuals result from the bright core being intrinsically more extended in the observational *Planck* map than in the opacity-corrected
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Corner plot of the marginal (1D) and joint (2D) probability density functions of the conversion factors from velocity-integrated opacity-corrected brightness temperature to dust intensity at 353 GHz, $X_{I_d/T_b}^{HI}$ and $X_{I_d/T_b}^{CO}$ [in (MJy sr$^{-1}$) (K km s$^{-1}$)$^{-1}$], for the two gas tracers, HI and CO, in the G139 region.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Maps of the intensity, $I_d$, of the dust emission at 353 GHz toward the G139 region. Left: observational map from Planck. Middle: best-fit maps reconstructed with the two gas tracers, HI and CO, before (top) and after (bottom) application of the Gaussian decomposition algorithm ROHSA to each gas tracer. Right: maps of the residuals obtained by subtracting the respective reconstructed maps from the observational Planck map. The color code refers to the absolute residuals, whereas the contour lines follow the relative residuals.</figcaption></figure>
CO map. They possibly reveal the presence of so-called dark gas, namely, gas that is undetected in HI and CO (e.g., Grenier et al. 2005; Planck Collaboration 2011b). Large positive and negative residuals also arise in the upper-left and lower-right corners, respectively, i.e., in two regions with little CO emission and where the northward gradient of the opacity-corrected HI emission remains too shallow to reproduce the observed dust-emission gradient. The negative residual along the lower-left boundary coincides with the slight over-intensity appearing in the HI map after opacity correction (see Sect. 3.1), which suggests that the HI emission could have been over-corrected. Other residuals could also result from imperfect opacity corrections, following a poor estimation of the excitation temperature. Alternatively, residuals could indicate that dust emission does not exactly follow the gas distribution – in other words, that the conversion factors are not perfectly uniform.

It would be interesting to compare our best-fit values of the conversion factors $X_{I_d/T_b}^{HI}$ and $X_{I_d/T_b}^{CO}$ to previous estimates, but we did not find such estimates in the literature. However, we found a number of estimates for the intermediate conversion factors $X_{I_d/N_H}^{HI}$, $X_{N_H/T_b}^{HI}$, $X_{I_d/N_H}^{12CO}$, and $X_{N_H/T_b}^{12CO}$ defined in Sect. 2.2.

Regarding HI, the value of the conversion factor from velocity-integrated brightness temperature to hydrogen column density is generally taken as $X_{N_H/T_b}^{HI} = 0.0182 (10^{20} \text{ cm}^{-2}) (\text{K km s}^{-1})^{-1}$ (from Wilson et al. 2013), which is strictly valid in the optically thin case. This value was used by HI4PI Collaboration (2016) and by Planck Collaboration (2011b,a), with, in the latter study, an opacity correction equivalent to that in our Eq. (A.5) with $T = 80$ K. The two Planck papers also discussed the value of the conversion factor from hydrogen column density to dust intensity at 353 GHz. Planck Collaboration (2011b) obtained $X_{I_d/N_H}^{HI} = (0.0484 \pm 0.0010) (\text{MJy sr}^{-1}) (10^{20} \text{ cm}^{-2})^{-1}$ in their reference region defined by $|b| > 20^\circ$ and $N_{HI} < 1.2 \times 10^{21} \text{ cm}^{-2}$, while Planck Collaboration (2011a) derived values in the range $[0.021, 0.071] (\text{MJy sr}^{-1}) (10^{20} \text{ cm}^{-2})^{-1}$ for local (low-velocity) clouds toward 14 high-latitude fields covering $\approx 825 \text{ deg}^2$ on the sky. Combining the above values of $X_{N_H/T_b}^{HI}$ and $X_{I_d/N_H}^{HI}$ gives for the conversion factor from velocity-integrated (opacity-corrected in the second case) brightness temperature to dust intensity at 353 GHz $X_{I_d/T_b}^{HI} = X_{I_d/N_H}^{HI} X_{N_H/T_b}^{HI} = (0.00088 \pm 0.000018) (\text{MJy sr}^{-1}) (\text{K km s}^{-1})^{-1}$ and $[0.00038, 0.00129] (\text{MJy sr}^{-1}) (\text{K km s}^{-1})^{-1}$, respectively, where the uncertainties include only uncertainties in $X_{I_d/N_H}^{HI}$, not uncertainties in $X_{N_H/T_b}^{HI}$ associated with, for instance, opacity saturation. Both estimates are consistent. Our best-fit value, $X_{I_d/T_b}^{HI} = (0.00081 \pm 0.000006) (\text{MJy sr}^{-1}) (\text{K km s}^{-1})^{-1}$, is slightly smaller than the value from Planck Collaboration (2011b), and it falls right within the range from Planck Collaboration (2011a).

Regarding CO, Remy et al. (2017) derived values of $X_{CO}$ and $(\tau_{353}/N_H)^{CO}$ for six nearby molecular clouds between longitude $l = 139^\circ$ and $191^\circ$ and between latitude $b = -3^\circ$ and $-56^\circ$. Relying on 353 GHz dust emission data and on [0.4, 100] GeV $\gamma$-ray emission data, they obtained values of $X_{CO}$ in the range $[0.86, 1.73] (10^{20} \text{ cm}^{-2}) (\text{K km s}^{-1})^{-1}$ and $[0.36, 1.14] (10^{20} \text{ cm}^{-2}) (\text{K km s}^{-1})^{-1}$, respectively. They further obtained values of $(\tau_{353}/N_H)^{CO}$ in the range $[12.6, 44] (10^{27} \text{ cm}^{-2})^{-1}$. Noting that $X_{N_H/T_b}^{CO} = 2 X_{CO}$ and $I_d = (4.76 \times 10^4 \text{ MJy sr}^{-1}) \tau_{353}$ for a dust temperature $T_d = 19.7$ K (Planck Collaboration 2014), we find that the corresponding ranges of $X_{I_d/T_b}^{CO} =$

$X_{I_d/N_H}^{CO}$ and $X_{N_H/T_b}^{CO}$ are $[0.10, 0.73] (\text{MJy sr}^{-1}) (\text{K km s}^{-1})^{-1}$ and $[0.043, 0.48] (\text{MJy sr}^{-1}) (\text{K km s}^{-1})^{-1}$, respectively. These ranges include no opacity correction for the CO line, which might explain why they are so broad. Our best-fit value, $X_{I_d/T_b}^{CO} = (0.0685 \pm 0.0019) (\text{MJy sr}^{-1}) (\text{K km s}^{-1})^{-1}$ falls somewhat below the former range and within, though close to the lower end, of the latter range. Finding a lower value of $X_{I_d/T_b}^{CO}$ here is not surprising given that our opacity correction increases the brightness temperatures, and this increase must be offset by a decrease in the conversion factor.

### 3.3. LoS decomposition into clouds

#### 3.3.1. ROHSA decomposition

Next, we apply the Gaussian decomposition algorithm ROHSA (presented in Sect. 2.3) to the observed spectral cubes of the opacity-corrected brightness temperatures, $T_{b*}$, of our two gas tracers, HI and CO, separately. The procedure yields 19 and 7 Gaussian kinematic components, respectively. The spectra of these components, averaged over the $26 \times 26$ pixels of the common $(l, b)$ grid, are plotted in the top and bottom panels, respectively, of Fig. 5. For simplicity, the components of each tracer are ordered by increasing mean velocity. Also plotted in Fig. 5 are the total reconstructed $T_{b*}^{HI}$ and $T_{b*}^{CO}$ spectra, obtained by summing the averaged spectra of their respective kinematic components (black solid lines), as well as the corresponding observed spectra (black dashed lines).

Comparing the reconstructed and observed spectra indicates that ROHSA does on average a very good job at reconstructing the observed spectra. Since we discarded the extracted components falling everywhere below twice the noise level, the reconstructed spectra are automatically free of the measurement noise apparent in the observational $^{12}$CO and $^{13}$CO cubes.

The HI Gaussian decomposition (top panel of Fig. 5) results in 19 components, which peak at various velocities between $\approx -90$ and $+2$ km s$^{-1}$ and which, together, cover almost the entire velocity range $[-120, +20]$ km s$^{-1}$. The HI spectrum is dominated by a cluster of 15 components (HI [5] - HI [19]) peaking between $\approx -15$ and $+2$ km s$^{-1}$. Also prominent in the HI spectrum is a strong and broad component (HI [2]) centered at $v \approx -40$ km s$^{-1}$.

The CO decomposition (bottom panel of Fig. 5) leads to 7 components peaking between $\approx -32$ and $0$ km s$^{-1}$. The four dominant components (CO [3] - CO [6]) are clustered around $v \sim -10$ km s$^{-1}$, with peak velocities between $\approx -16$ and $0$ km s$^{-1}$; each of these components could easily be related to one of the 15 clustered HI components. The CO spectrum also contains a very weak component (CO [7]) centered at $v \approx 0$ km s$^{-1}$ as well as two very close components (CO [1] and CO [2]) centered at $v \approx -32$ km s$^{-1}$. The latter probably form a single physical entity, which was artificially split by the sharp transition between our two estimates of $T_{b*}^{12CO}$; they will naturally be recombined at the cloud-reconstruction step. All the CO components are narrower than the HI components.

The velocity-integrated maps of all the kinematic components of our two gas tracers are plotted in the top and bottom parts, respectively, of Fig. 6, along with the total reconstructed maps of both tracers (leftmost column), obtained by superposing their respective kinematic components. The brightest HI component is clearly HI [2], which, we noted earlier, is also prominent in the average HI spectrum (top panel of Fig. 5). The globally brightest CO component is CO [3], which leads to the
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Spectra of the opacity-corrected brightness temperatures, $T_{b*}^{\text{HI}}$ (top) and $T_{b*}^{\text{CO}}$ (bottom), averaged over the 26 × 26 pixels of the common $(l, b)$ grid, toward the G139 region. The observed spectra are plotted in black dashed line, the reconstructed spectra obtained after application of the Gaussian decomposition algorithm ROHSA are plotted in black solid line, and the spectra of the individual Gaussian kinematic components, which are ordered (for each tracer) by increasing mean velocity, are plotted in color. The $T_{b*}^{\text{HI}}$ and $T_{b*}^{\text{CO}}$ spectra can be rescaled to dust intensity per unit velocity at 353 GHz using the best-fit values of the conversion factors derived in Sect. 3.2 (see Fig. 3).</figcaption></figure>

highest peak in the average CO spectrum (bottom panel of Fig. 5). The locally brightest CO feature is the bright core to the right, which, we now see, can be attributed to the CO [1] - CO [2] pair at $v \approx -32 \text{ km s}^{-1}$.

The total reconstructed maps of both gas tracers are also plotted in the fourth row of Fig. 2, where they can be compared to the corresponding re-gridded observational maps in the third row. The comparison confirms the generally very good quality of the ROHSA reconstruction, which we already noted when discussing the average spectra in Fig. 5. Furthermore, in the same way as the re-gridded observational maps of both gas tracers were rescaled (see last two sentences in caption of Fig. 2) and combined into a map of the 353 GHz dust intensity (top-middle panel of Fig. 4), their reconstructed ROHSA maps can be rescaled and combined into a reconstructed ROHSA map of the 353 GHz dust intensity (bottom-middle panel of Fig. 4). Here, too, the agreement with the pre-ROHSA dust intensity map is very good.

It is interesting to compare our Gaussian decomposition of the CO cube to previous decompositions. The G139 region is part of the $^{12}$CO cube constructed by Dame et al. (2001) from a composite CO survey of the entire Galactic plane with the CfA & Cerro Tololo 1.2 m telescopes. Straižys & Laugalys (2007)

divided the G139 region of this cube into three layers with different RVs (as defined earlier by Digel et al. 1996): the Gould Belt layer, with $v \approx [-5, +10] \text{ km s}^{-1}$, the Camelopardalis (Cam) OB1 association layer, with $v \approx [-20, -5] \text{ km s}^{-1}$, and the Perseus arm, with $v \approx [-60, -30] \text{ km s}^{-1}$. Using the Galactic rotation curve to convert RVs to distances, they estimated the distances to the clouds of the Gould Belt layer, the Cam OB1 layer, and the Perseus arm at $d \approx [150, 300] \text{ pc}$, $[800, 900] \text{ pc}$, and $[2, 3] \text{ kpc}$, respectively. Later, Montillaud et al. (2015) identified three kinematic components, with RVs $v \approx -34 \text{ km s}^{-1}$, $-16 \text{ km s}^{-1}$, and $-9 \text{ km s}^{-1}$, in the same $^{12}$CO data from Dame et al. (2001) as those studied by Straižys & Laugalys (2007). The $v \approx -34 \text{ km s}^{-1}$ component clearly peaks at the position of the bright core in the dust intensity map (top-left panel of Fig. 1), the $v \approx -16 \text{ km s}^{-1}$ component has a more extended and more diffuse emission, and the $v \approx -9 \text{ km s}^{-1}$ component is fainter. Referring to the work of Straižys & Laugalys (2007), they concluded that the $v \approx -34 \text{ km s}^{-1}$ component must be part of the Perseus arm, while the $v \approx -16 \text{ km s}^{-1}$ and $v \approx -9 \text{ km s}^{-1}$ components must belong to the Cam OB1 layer. In addition, their identification of the $v \approx -34 \text{ km s}^{-1}$ component with the bright core led them to place the latter at a distance $d = 2.5 \pm 0.5 \text{ kpc}$. To make
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Reconstructed maps of the opacity-corrected brightness temperatures, $T_{b*}^{HI}$ (top) and $T_{b*}^{CO}$ (bottom), integrated over velocity, toward the G139 region, after application of the Gaussian decomposition algorithm ROHSA (leftmost column), together with the contributions from the individual Gaussian kinematic components (columns 2–6). The maps can be rescaled to dust intensity at 353 GHz using the best-fit values of the conversion factors derived in Sect. 3.2 (see Fig. 3).</figcaption></figure>

the link with our own study, the $v \approx -34 \text{ km s}^{-1}$, $-16 \text{ km s}^{-1}$, and $-9 \text{ km s}^{-1}$ components of Montillaud et al. (2015) presumably correspond, in our decomposition, to the CO [1] - CO [2] pair (centered at $v \approx -32 \text{ km s}^{-1}$) and to the two peaks formed by the clustered CO [3] - CO [6] components (at $v \approx -16$ and $-7 \text{ km s}^{-1}$).

### 3.3.2. Cloud reconstruction

Moving on to the second step of the procedure described in Sect. 2.3, we examine the $n_G^{\text{tot}} = 19 + 7 = 26$ Gaussian kinematic components identified with ROHSA and seek to group together, i.e., assign to a same cloud, those with similar velocity profiles. To that end, we consider every possible pair of components *j* and *j'* (*j* ≠ *j'*), compute the correlation coefficient $C_{jj'}$ (Eq. (29)) at each of the 26 × 26 pixels of our grid, retain the weighted average value of $C_{jj'}$ over all the pixels, $\langle C_{jj'} \rangle$ (Eq. (30)), and compare $\langle C_{jj'} \rangle$ to the velocity-coherence threshold, $C_{\text{th}}$, entering Eq. (31). Here, for the purpose of illustration, we adopt $C_{\text{th}} = 50\%$. The results obtained for this value of $C_{\text{th}}$ are presented in Sects. 3.3–3.4, while an overview of the results obtained for all the integer values of $C_{\text{th}}$ between 0 and 100% is provided in Sect. 3.5. When $\langle C_{jj'} \rangle \ge C_{\text{th}}$, we combine components *j* and *j'* and assign them to a same cloud. The results of our correlation analysis for the 325 possible pairs of components are reported in Fig. 7, with mini-maps of $C_{jj'}$ displayed in the upper-right half and the derived values of $\langle C_{jj'} \rangle$ indicated in the lower-left half.

It emerges from Fig. 7 that 83 pairs satisfy the condition $\langle C_{jj'} \rangle \ge C_{\text{th}}$; for better visibility, the corresponding small squares in the lower-left half of the figure are shaded in red (with increasing level of red as $\langle C_{jj'} \rangle$ increases), as opposed to light blue for the other pairs. Amongst the pairs with $\langle C_{jj'} \rangle \ge C_{\text{th}}$, those having a component in common are further grouped together into a same multicomponent cloud. The end result is a set of seven clouds, which we name $C_1, C_2, C_3, C_4, C_5, C_6$, and $C_7$, and which enclose 1, 2, 1, 4, 15, 2, and 1 components, respectively. The different components of each cloud can be retrieved from the labels with the cloud's name along the diagonal.

Using the best-fit values of the conversion factors derived in Sect. 3.2 (see Fig. 3), we can now rescale the $T_{b*}$ spectral cubes of all the kinematic components to dust emission at 353 GHz. This common dust scale enables us to combine the different components of each cloud and thus obtain its dust emission spectral cube. In the left part of Fig. 8, we plot the dust emission spectra of the seven clouds, averaged over the 26 × 26 pixels (black solid lines), together with the spectra of their individual components (color lines). The dust intensity maps of the seven clouds are displayed in the right column.

$C_1, C_2$, and $C_3$ are three purely atomic clouds, containing 1, 2, and 1 HI components, respectively. $C_1$ is faint, with an enhancement in the southwest corner; $C_2$ is much brighter, especially in the northwest part; $C_3$ is faint, knotty, and mostly confined to the northeast corner. $C_1$ and $C_2$ both cover a rather broad velocity range, which could perhaps indicate that they are quite extended along the LoS.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Correlation coefficient, $C_{jj'}$ (Eq. (29)), expressed in percent, for the 325 possible pairs of Gaussian kinematic components *j* and *j'* ($j \neq j'$) extracted with R0HSA in the HI and CO spectral cubes of the G139 region. The names of the 26 components, together with the mini-maps of their dust intensities, are shown along the top row and left column. Mini-maps of $C_{jj'}$ are displayed in the upper-right half of the figure. The weighted average values of $C_{jj'}$ over the 26 × 26 pixels, $\langle C_{jj'} \rangle$ (Eq. (30)), are indicated in the lower-left half, with red or light-blue shading according to whether $\langle C_{jj'} \rangle \ge 50\%$ or $\langle C_{jj'} \rangle < 50\%$ and with a level of shading increasing with increasing $\langle C_{jj'} \rangle$. Pairs for which $\langle C_{jj'} \rangle \ge 50\%$ have their two components *j* and *j'* combined and assigned to a same cloud. The 26 components are thus grouped into seven different clouds, $C_1, ..., C_7$, as indicated along the diagonal.</figcaption></figure>

$C_4$ and $C_5$ are two mixed atomic-molecular clouds. $C_4$ encloses an oblique, elongated CO structure (2 CO components), partially surrounded by an extended HI envelope (2 HI components). $C_5$ is the richest cloud, with 13 HI and 2 CO components, which together fill the entire region. In both clouds, HI appears to surround CO in velocity.

$C_6$ and $C_7$ are two purely molecular clouds, with two very close CO components for $C_6$ and a single CO component for $C_7$. $C_6$ is bright and localized both in the sky and in velocity; it clearly corresponds to the bright core in the observational dust intensity map (top-left panel of Fig. 1). $C_7$ is much fainter, confined to the northeast corner, and also localized in velocity.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Reconstructed spectra (**left**) and maps (**right**) of the 353 GHz dust intensities, $I_{d,i}$, of the seven clouds, $C_i$ ($i$ = 1, ..., 7), identified in the G139 region. The spectra are averaged over the 26 × 26 pixels of the (l, b) grid. The total spectra are plotted in black solid line, while the spectra of the individual kinematic components are plotted in color, using the same colors as in Fig. 5 and two different linestyles to distinguish the two gas tracers (solid for HI and dashed for CO).</figcaption></figure>

Most of the dust emission from the G139 region arises in $C_5$ (≈42%), $C_2$ (≈30%), and $C_4$ (≈18%), with a small, but locally high contribution from $C_6$ (≈1.4%). $C_5$ and $C_2$ account for most of the HI emission (third panel in the left column of Fig. 2), while $C_4$ and $C_6$ account for most of the CO emission (right column): $C_4$ produces the oblique band north of the southeast-northwest diagonal with the two bright spots at its north end, and $C_6$ is responsible for the high CO peak to the right.

### 3.4. Derivation of the magnetic field orientation in each cloud

Our LoS decomposition of the 353 GHz Planck dust intensity map of the G139 region in Sect. 3.3 led to the identification of seven clouds, $C_1, ..., C_7$. For each cloud $C_i$, we obtained a map of the dust intensity, $I_{d,i}$ (right column of Fig. 8). It now remains to determine the magnetic field orientation in each cloud, or, equivalently, its polarization fraction, $\bar{p}_{d,i}$, and polarization angle, $\bar{\psi}_{d,i}$. Following the procedure described in Sect. 2.4, we derive the best-fit values of $\bar{p}_{d,i}$ and $\bar{\psi}_{d,i}$, together with their uncertainties, $\sigma(\bar{p}_{d,i})$ and $\sigma(\bar{\psi}_{d,i})$, by minimizing $\chi_r^2$ in Eq. (40) through MCMC simulations. As a prior, we require $\bar{p}_{d,i} \le 23\%$, as suggested by Planck Collaboration (2015a) (see Appendix B).

The best fit has $\chi_r = 3.17$. This fairly large value indicates that our model composed of seven clouds with uniform polarization parameters ($\bar{p}_{d,i}$ and $\bar{\psi}_{d,i}$) does not reproduce the Planck polarization maps within the total “observational uncertainties”, $\sigma_Q$ and $\sigma_U$, defined below Eq. (40). These uncertainties are dominated by “decomposition errors” in $Q_d^{\text{mod}}$ and $U_d^{\text{mod}}$ arising

from decomposition errors in the cloud intensities (second term in the r.h.s of Eqs. (41)–(42)), which generally exceed measurement errors in $Q_d^{\text{obs}}$ and $U_d^{\text{obs}}$ (first term) by a factor ~2–3. The main reason for the fairly large value of $\chi_r$ is that the polarization parameters of our seven clouds are actually not uniform. Allowing for large-scale variations in these parameters would almost certainly improve the situation. It could also be that some of our clouds do not form magnetically coherent structures, in the sense that they are actually composed of two or more magnetically distinct regions, with different polarization parameters.

The best-fit values of the free parameters, $\bar{p}_{d,i}$ and $\bar{\psi}_{d,i}$, together with their standard deviations, are listed in Table 2. It appears that the polarization fractions vary widely, from $\bar{p}_{d,4} \approx 2\%$ to $\bar{p}_{d,1}$, $\bar{p}_{d,7}$ hitting our imposed upper limit of 23%. Clearly, $C_7$ is trying to make up for the missing matter in the northeast corner (see positive residual in the bottom-right panel of Fig. 4), while $C_1$ is trying to contribute to the enhanced LoS-averaged polarization fraction in the south (see left panel in the fourth row of Fig. 9 below). The polarization angles also vary widely, from $|\bar{\psi}_{d,1}|, |\bar{\psi}_{d,5}| \lesssim 5^\circ$ to $|\bar{\psi}_{d,3}| \approx 90^\circ$. Yet most of the matter is contained within clouds with small polarization angles, consistent with the small values of the LoS-averaged polarization angle measured by Planck (left panel in the bottom row of Fig. 9).

We can now use the best-fit maps of $I_{d,i}$ (right column of Fig. 8) in conjunction with the best-fit values of $\bar{p}_{d,i}$ and $\bar{\psi}_{d,i}$ (third and fourth columns of Table 2) to compute the Stokes parameters $Q_{d,i}$ and $U_{d,i}$ of each cloud $C_i$ (Eqs. (34) and (35))
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Maps of the two Stokes parameters for linear polarization, $Q_d$ and $U_d$, the polarized intensity, $P_d$, the LoS-averaged polarization fraction, $p_d^{\text{los}}$, and the LoS-averaged polarization angle, $\psi_d^{\text{los}}$, of the dust emission at 353 GHz toward the G139 region. Left: observational maps from Planck. Middle: best-fit maps reconstructed with our seven clouds, $C_1, ..., C_7$. Right: maps of the residuals obtained by subtracting the respective reconstructed maps from the observational Planck maps.</figcaption></figure>

<table><caption>**Table 2.** List of the seven clouds, $C_i$ ($i = 1, ..., 7$), identified toward the G139 region.</caption><thead><tr><th>Cloud</th><th>Kinematic components</th><th>p&#773;<sub>d,i</sub> [%]<sup>a</sup></th><th>&psi;&#773;<sub>d,i</sub> [deg]<sup>b</sup></th></tr></thead><tbody><tr><td>C<sub>1</sub></td><td>HI [1]</td><td>22.9 &plusmn; 0.1</td><td>5.2 &plusmn; 1.3</td></tr><tr><td>C<sub>2</sub></td><td>HI [2-3]</td><td>6.4 &plusmn; 0.2</td><td>16.8 &plusmn; 0.8</td></tr><tr><td>C<sub>3</sub></td><td>HI [4]</td><td>14.2 &plusmn; 0.4</td><td>89.4 &plusmn; 0.9</td></tr><tr><td>C<sub>4</sub></td><td>HI [5-6], CO [3-4]</td><td>2.2 &plusmn; 0.2</td><td>46.3 &plusmn; 1.9</td></tr><tr><td>C<sub>5</sub></td><td>HI [7-19], CO [5-6]</td><td>17.9 &plusmn; 0.1</td><td>-3.4 &plusmn; 0.2</td></tr><tr><td>C<sub>6</sub></td><td>CO [1-2]</td><td><b>5.8 &plusmn; 0.3</b></td><td><b>66.1 &plusmn; 1.7</b></td></tr><tr><td>C<sub>7</sub></td><td>CO [7]</td><td>22.5 &plusmn; 0.5</td><td>-21.2 &plusmn; 3.4</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Best-fit value and standard deviation of the polarization fraction of cloud $C_i$. <sup>(b)</sup>Best-fit value and standard deviation of the polarization angle of cloud $C_i$. The boxed values pertain to the two most interesting molecular clouds: the bright CO core, $C_6$ (bold style), and the oblique, elongated cloud, $C_4$ (normal style).

as well as the resulting Stokes parameters, $Q_d^{\text{mod}}$ and $U_d^{\text{mod}}$ (Eqs. (25) and (26)), and the associated polarized intensity, $P_d^{\text{mod}}$ (Eq. (5)), LoS-averaged polarization fraction, $p_d^{\text{los,mod}}$ (Eq. (8)), and LoS-averaged polarization angle, $\psi_d^{\text{los,mod}}$ (Eq. (9)). The best-fit maps of $Q_d^{\text{mod}}$, $U_d^{\text{mod}}$, $P_d^{\text{mod}}$, $p_d^{\text{los,mod}}$, and $\psi_d^{\text{los,mod}}$, are displayed in the middle column of Fig. 9, where they can be compared to the observational Planck maps in the left column. The observational LoS-averaged polarization fraction, $p_d^{\text{los,obs}}$, was debiased using the modified asymptotic (MAS) estimator (Planck Collaboration 2020a). The maps of the residuals, i.e., the differences between the observational and reconstructed maps, are shown in the right column.

The map of $P_d^{\text{obs}}$ (left panel in the third row of Fig. 9) is globally more uniform than the map of $I_d^{\text{obs}}$ (left panel of Fig. 4). This is because brighter [fainter] regions are generally less [more] polarized, as can be seen by comparing the maps of $I_d^{\text{obs}}$ (left
panel of Fig. 4) and $p_d^{\text{los,obs}}$ (left panel in the fourth row of Fig. 9). The map of $P_d^{\text{obs}}$ actually exhibits two faint spots along a central horizontal band. The east faint spot and the westernmost part of the west faint spot correspond to regions that are both fainter and less polarized, whereas the easternmost part of the west faint spot is bright and weakly polarized. Our model does a good job at reproducing the map of $P_d^{\text{obs}}$, including the two faint spots (see middle panel in the third row of Fig. 9). The similarity between the maps of $P_d^{\text{mod}}$ and $I_{d,5}$ (fifth panel in the right column of Fig. 8) suggest that the bright and strongly polarized ($\bar{\psi}_{d,5} \approx 18\%$) cloud $\mathbb{C}_5$ provides a dominant contribution to $P_d^{\text{mod}}$.

The map of $p_d^{\text{los,obs}}$ (left panel in the fourth row of Fig. 9) is roughly divided into a strongly polarized region south and a more weakly polarized region north, with a distinct low-polarization spot (blue region). All these features are nicely reproduced in our model (middle panel), which also provides a natural explanation for the low-polarization spot. This spot arises precisely at the location of $\mathbb{C}_6$, which is very bright and has a polarization angle very different from those of the other bright clouds, $\mathbb{C}_5$ and $\mathbb{C}_2$, detected along its LoS (see fourth column of Table 2). As a result, $\mathbb{C}_6$ contributes constructively to $I_d^{\text{mod}}$, leading to a bright spot in the map of $I_d^{\text{mod}}$ (Fig. 4), and contributes destructively to $P_d^{\text{mod}}$, leading to a dip in the map of $P_d^{\text{mod}}$ and a low-polarization spot in the map of $p_d^{\text{los,mod}}$ (Fig. 9).

The map of $\psi_d^{\text{los,obs}}$ (left panel in the fifth row of Fig. 9) is fairly homogeneous, with $|\psi_d^{\text{los,obs}}| \lesssim 10^\circ$ everywhere. The four large-$\psi_d^{\text{los,obs}}$ spots (orange regions) near the middle are associated with low $p_d^{\text{los,obs}}$, which most likely indicates the presence of depolarizing clouds along the LoS. The reason why the low-polarization spot in the map of $p_d^{\text{los,obs}}$ is not associated with a particularly large $\psi_d^{\text{los,obs}}$ (it actually sits between two large-$\psi_d^{\text{los,obs}}$ spots) is probably because the putative depolarizing cloud has a polarization angle relatively close to $90^\circ$. The model map of $\psi_d^{\text{los,mod}}$ (middle panel) is even more homogeneous than its observational counterpart, except toward $\mathbb{C}_6$, where $\psi_d^{\text{los,mod}}$ reaches $\approx 16^\circ$. The homogeneous, low-$|\psi_d^{\text{los,mod}}|$ background can be mostly attributed to the widespread, bright, and strongly polarized cloud $\mathbb{C}_5$, which dominates the polarized emission and has a small polarization angle ($\bar{\psi}_{d,5} \approx -3^\circ$). The jump in $\psi_d^{\text{los,mod}}$ toward $\mathbb{C}_6$ can be explained by this very bright

and moderately polarized cloud having a much larger polarization angle ($\bar{\psi}_{d,6} \approx 66^\circ$). The discrepancies between $\psi_d^{\text{los,obs}}$ and $\psi_d^{\text{los,mod}}$ are mostly due to our assumption that each cloud has a uniform polarization angle. In particular, the behavior of $\psi_d^{\text{los,obs}}$ in the vicinity of $\mathbb{C}_6$ suggests that the true polarization angle of $\mathbb{C}_6$ has values closer to $90^\circ$ in the bright center of the cloud (hence a small impact on $\psi_d^{\text{los,obs}}$) and values closer to $45^\circ$ on either side (hence a more significant increase of $\psi_d^{\text{los,obs}}$).

Once we have derived the best-fit values of the polarization fraction, $\bar{p}_{d,i}$, and polarization angle, $\bar{\psi}_{d,i}$, of every cloud $\mathbb{C}_i$, we can, in principle, obtain the best-fit orientation of its internal magnetic field, $\bar{B}_i$, defined by the orientation angle of $\bar{B}_i$ in the PoS, $\bar{\psi}_{B,i}$ (Eq. (43)), and the inclination angle of $\bar{B}_i$ to the PoS, $\bar{\gamma}_{B,i}$ (Eq. (44)).

The small values of $|\bar{\psi}_{d,1}|$ and $|\bar{\psi}_{d,5}|$ indicate that the faint $\mathbb{C}_1$ and the dominant $\mathbb{C}_5$ both have nearly horizontal $\bar{B}_\perp$. At the other extreme, $|\bar{\psi}_{d,3}| \approx 90^\circ$ corresponds to a vertical $\bar{B}_\perp$; however, the faint $\mathbb{C}_3$, which appears to be mostly confined to the northeast corner, is probably part of a larger cloud that extends beyond the G139 region, so its derived polarization parameters might not be representative of this larger cloud. In between, the bright $\mathbb{C}_2$ and the faint $\mathbb{C}_7$ have slightly tilted $\bar{B}_\perp$, whereas the two most interesting molecular clouds, $\mathbb{C}_6$ and $\mathbb{C}_4$, have significantly tilted $\bar{B}_\perp$. Our previous argument regarding $\psi_d^{\text{los,obs}}$ suggests that $\bar{B}_\perp,6$ is more strongly tilted in the bright center of $\mathbb{C}_6$ than around it. Also noteworthy is that $\bar{B}_\perp,4$ appears to be nearly aligned with the elongated main segment of $\mathbb{C}_4$.

The difficulty with $\bar{\gamma}_{B,i}$ is that the maximum polarization fraction, $(\bar{p}_{d,i})_{\text{max}}$, entering Eq. (44) is unknown. Below Eq. (44), we suggested adopting $(\bar{p}_{d,i})_{\text{max}} = 23\%$. Evidently, this is not a reasonable choice for $\mathbb{C}_1$ and $\mathbb{C}_7$, whose polarization fractions would have exceeded 23% had they not been constrained by our imposed upper limit of 23% (see third column of Table 2). All we can reasonably conclude regarding the magnetic field inclinations to the PoS is that $\bar{B}_1$ and $\bar{B}_7$ are probably close to the PoS; $\bar{B}_2, \bar{B}_3, \bar{B}_5$, and $\bar{B}_6$ are probably moderately inclined to the PoS; and $\bar{B}_4$ is probably close to the LoS. An additional conclusion can potentially be drawn regarding $\mathbb{C}_4$, for which we just noted that $\bar{B}_\perp,4$ is nearly aligned with the elongated main segment of the cloud. If this alignment also exists in 3D, the main segment of $\mathbb{C}_4$ must also be close to the LoS, with the implication that it is actually much more elongated in 3D than in the PoS.

$$p_d^{\text{los}} = \frac{\sqrt{P_{d,0}^2 + P_{d,1}^2 + 2 P_{d,0} P_{d,1} \cos(2\psi_{d,1})}}{I_{d,0} + I_{d,1}} \quad (45)$$

and

$$\tan(2\psi_d^{\text{los}}) = \frac{P_{d,1} \sin(2\psi_{d,1})}{P_{d,0} + P_{d,1} \cos(2\psi_{d,1})}, \quad (46)$$

with $P_{d,0} = (I_{d,0} p_{d,0})$ and $P_{d,1} = (I_{d,1} p_{d,1})$. The maximum rotation of the polarization orientation occurs for

$$\cos(2\psi_{d,1}) = \frac{P_{d,1}}{P_{d,0}}, \quad (47)$$

which, in the case $P_{d,1} \ll P_{d,0}$, reduces to $|\psi_{d,1}| \approx 45^\circ$.

### 3.5. Results obtained for other values of the velocity-coherence threshold, $C_{th}$

The 26 Gaussian kinematic components extracted with ROHSA each have a well-determined dust intensity, $\tilde{I}_{d,j}$, and the total dust intensity, $I_d^{\text{mod}}$ (Eqs. (24) and (32)), is just the sum of the 26 $\tilde{I}_{d,j}$, independent of how the 26 components are subsequently grouped into clouds. In contrast, the polarization fraction and angle of a given component are not determined by ROHSA, but are those of the cloud that this component is assigned to. As a result, the Stokes parameters, $Q_d^{\text{mod}}$ and $U_d^{\text{mod}}$ (Eqs. (38) and (39)), as well as the associated polarized intensity, $P_d^{\text{mod}}$ (Eq. (5)), LoS-averaged polarization fraction, $p_d^{\text{los,mod}}$ (Eq. (8)), and LoS-averaged polarization angle, $\psi_d^{\text{los,mod}}$ (Eq. (9)), are sensitive to the exact distribution of the 26 components between different clouds. This distribution, in turn, is governed by the

<footnote n="5">5 A cloud with parameters $I_{d,1}$, $p_{d,1}$, $\psi_{d,1}$ against a brighter emission background with parameters $I_{d,0}$, $p_{d,0}$, $\psi_{d,0} = 0^\circ$ causes depolarization and rotation of the polarization orientation. The resulting LoS-averaged polarization fraction and angle are given by</footnote>
<table><caption>**Table 3.** Results obtained for the first eight cloud configurations identified in the G139 region.</caption><thead><tr><th>Configuration</th><th>Range of C<sub>th</sub> [%]<sup>a</sup></th><th>Clouds</th><th>Kinematic components</th><th>p&#773;<sub>d,i</sub> (%)<sup>b</sup></th><th>&psi;&#773;<sub>d,i</sub> (deg)<sup>c</sup></th><th>&chi;<sub>r</sub><sup>d</sup></th></tr></thead><tbody><tr><td>1 cloud</td><td>[0, 1]</td><td>C<sub>1</sub></td><td>H1 [1–19], CO [1–7]</td><td>8.7 &pm; 0.1</td><td>2.7 &pm; 0.1</td><td>4.66</td></tr><tr><td rowspan="2">2 clouds</td><td rowspan="2">[2, 14]</td><td>C<sub>1</sub></td><td>H1 [1]</td><td>22.9 &pm; 0.1</td><td>-1.5 &pm; 0.8</td><td rowspan="2">4.53</td></tr><tr><td>C<sub>2</sub></td><td>H1 [2–19], CO [1–7]</td><td>8.2 &pm; 0.1</td><td>3.2 &pm; 0.1</td></tr><tr><td rowspan="3">3 clouds</td><td rowspan="3">[15, 27]</td><td>C<sub>1</sub></td><td>H1 [1]</td><td>22.9 &pm; 0.1</td><td>-1.5 &pm; 0.8</td><td rowspan="3">4.30</td></tr><tr><td>C<sub>2</sub></td><td>H1 [2–4], CO [1–2]</td><td>4.1 &pm; 0.1</td><td>25.2 &pm; 0.8</td></tr><tr><td>C<sub>3</sub></td><td>H1 [5–19], CO [3–7]</td><td>11.5 &pm; 0.1</td><td>-1.1 &pm; 0.2</td></tr><tr><td rowspan="4">4 clouds</td><td rowspan="4">[28, 31]</td><td>C<sub>1</sub></td><td>H1 [1]</td><td>22.9 &pm; 0.1</td><td>-1.6 &pm; 0.8</td><td rowspan="4">4.29</td></tr><tr><td>C<sub>2</sub></td><td>H1 [2–4], CO [1–2]</td><td>4.1 &pm; 0.1</td><td>25.2 &pm; 0.8</td></tr><tr><td>C<sub>3</sub></td><td>H1 [5–19], CO [3–6]</td><td>11.6 &pm; 0.1</td><td>-1.0 &pm; 0.2</td></tr><tr><td>C<sub>4</sub></td><td>CO [7]</td><td>22.5 &pm; 0.5</td><td>-63.7 &pm; 3.6</td></tr><tr><td rowspan="5">5 clouds</td><td rowspan="5">[32, 41]</td><td>C<sub>1</sub></td><td>H1 [1]</td><td>22.9 &pm; 0.1</td><td>-1.6 &pm; 0.8</td><td rowspan="5">4.29</td></tr><tr><td>C<sub>2</sub></td><td>H1 [2–4]</td><td>4.5 &pm; 0.2</td><td>21.5 &pm; 1.0</td></tr><tr><td>C<sub>3</sub></td><td>H1 [5–19], CO [3–6]</td><td>11.2 &pm; 0.1</td><td>-0.9 &pm; 0.2</td></tr><tr><td>C<sub>4</sub></td><td>CO [1–2]</td><td>3.5 &pm; 0.3</td><td>40.8 &pm; 2.8</td></tr><tr><td>C<sub>5</sub></td><td>CO [7]</td><td>22.5 &pm; 0.5</td><td>-64.4 &pm; 3.6</td></tr><tr><td rowspan="6">6 clouds</td><td rowspan="6">[42, 43]</td><td>C<sub>1</sub></td><td>H1 [1]</td><td>22.9 &pm; 0.1</td><td>-1.8 &pm; 0.8</td><td rowspan="6">4.14</td></tr><tr><td>C<sub>2</sub></td><td>H1 [2–3]</td><td>7.2 &pm; 0.2</td><td>14.9 &pm; 0.7</td></tr><tr><td>C<sub>3</sub></td><td>H1 [4]</td><td>13.0 &pm; 0.4</td><td>89.3 &pm; 1.0</td></tr><tr><td>C<sub>4</sub></td><td>H1 [5–19], CO [3–6]</td><td>11.2 &pm; 0.1</td><td>-0.9 &pm; 0.2</td></tr><tr><td>C<sub>5</sub></td><td>CO [1–2]</td><td>3.4 &pm; 0.3</td><td>43.5 &pm; 2.9</td></tr><tr><td>C<sub>6</sub></td><td>CO [7]</td><td>22.2 &pm; 0.8</td><td>-26.4 &pm; 4.2</td></tr><tr><td rowspan="7">7 clouds</td><td rowspan="7">[44, 51]</td><td>C<sub>1</sub></td><td>H1 [1]</td><td>22.9 &pm; 0.1</td><td>5.2 &pm; 1.3</td><td rowspan="7">3.17</td></tr><tr><td>C<sub>2</sub></td><td>H1 [2–3]</td><td>6.4 &pm; 0.2</td><td>16.8 &pm; 0.8</td></tr><tr><td>C<sub>3</sub></td><td>H1 [4]</td><td>14.2 &pm; 0.4</td><td>89.4 &pm; 0.9</td></tr><tr><td>C<sub>4</sub></td><td>H1 [5–6], CO [3–4]</td><td>2.2 &pm; 0.2</td><td>46.3 &pm; 1.9</td></tr><tr><td>C<sub>5</sub></td><td>H1 [7–19], CO [5–6]</td><td>17.9 &pm; 0.1</td><td>-3.4 &pm; 0.2</td></tr><tr><td>C<sub>6</sub></td><td>CO [1–2]</td><td>5.8 &pm; 0.3</td><td>66.1 &pm; 1.7</td></tr><tr><td>C<sub>7</sub></td><td>CO [7]</td><td>22.5 &pm; 0.5</td><td>-21.2 &pm; 3.4</td></tr><tr><td rowspan="8">8 clouds</td><td rowspan="8">[52, 58]</td><td>C<sub>1</sub></td><td>H1 [1]</td><td>22.9 &pm; 0.1</td><td>4.8 &pm; 1.3</td><td rowspan="8">3.14</td></tr><tr><td>C<sub>2</sub></td><td>H1 [2–3]</td><td>6.4 &pm; 0.2</td><td>17.2 &pm; 0.8</td></tr><tr><td>C<sub>3</sub></td><td>H1 [4]</td><td>13.8 &pm; 0.4</td><td>89.5 &pm; 1.1</td></tr><tr><td>C<sub>4</sub></td><td>H1 [5–6], CO [3–4]</td><td>2.1 &pm; 0.2</td><td>42.5 &pm; 2.0</td></tr><tr><td>C<sub>5</sub></td><td>H1 [7–18], CO [5–6]</td><td>17.2 &pm; 0.1</td><td>-3.1 &pm; 0.2</td></tr><tr><td>C<sub>6</sub></td><td>H1 [19]</td><td>22.9 &pm; 0.1</td><td>-5.4 &pm; 1.0</td></tr><tr><td>C<sub>7</sub></td><td>CO [1–2]</td><td>6.0 &pm; 0.3</td><td>66.2 &pm; 1.7</td></tr><tr><td>C<sub>8</sub></td><td>CO [7]</td><td>22.5 &pm; 0.5</td><td>-21.1 &pm; 3.8</td></tr></tbody></table>

Notes. <sup>(a)</sup> $C_{th}$ is the threshold imposed on the weighted-averaged correlation coefficient, $\langle C_{ij'} \rangle$ (Eq. (30)), as a condition for components j and j' to be assigned to a same cloud (Eq. (31)). <sup>(b)</sup> Best-fit values and standard deviations of the polarization fractions of the different clouds $C_i$ found in every configuration. <sup>(c)</sup> Best-fit values and standard deviations of the polarization angles of the different clouds $C_i$ found in every configuration. <sup>(d)</sup> Best-fit reduced $\chi$. The boxed values pertain to the two most interesting molecular clouds: the bright CO core (bold style) and the oblique, elongated cloud (normal style).

value adopted for the velocity-coherence threshold, $C_{th}$, entering Eq. (31).

The purpose of this section is to examine the impact of $C_{th}$ on the polarization results. By considering all the integer values of $C_{th}$ in the range [0, 100%], we are led to identify different cloud configurations with increasing numbers of clouds, from $n_{cl} = 1$ to 26. In Table 3 and Fig. 10, we present the results obtained for the first eight configurations, with $n_{cl} = 1$ to 8, corresponding to $C_{th} = 0$ to 58%. For every configuration, Table 3 lists the relevant range of $C_{th}$, the ROHSA kinematic components of all the clouds, the best-fit values and standard deviations of their

polarization parameters, and the best-fit $\chi_r$, while Fig. 10 displays the spatially-averaged spectra of all the clouds, their dust intensity maps together with their polarization half-vectors, and the reconstructed maps of $P_d^{\text{mod}}$, $P_d^{\text{los,mod}}$, and $\psi_d^{\text{los,mod}}$.

The single cloud formed at small values of $C_{th}$ loses its weak, high-velocity H1 [1] component as soon as $C_{th}$ reaches 2%. It then breaks up into two mixed atomic-molecular clouds when $C_{th}$ reaches 15%: the localized cloud $C_2$, which includes the bright CO core, and the pervasive cloud $C_3$, which contains the rest of the CO gas. $C_3$ loses its very weak, low-velocity CO [7] component at $C_{th} = 28\%$, while $C_2$ successively loses the bright
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Results obtained for the eight cloud configurations described in Table 3. **Left:** spatially-averaged spectra of the 353 GHz dust intensity of the $n_{\text{el}}$ clouds, $C_i$ ($i = 1, ..., n_{\text{el}}$), all in the same (arbitrary) units. **Center:** maps of the 353 GHz dust intensity (in MJ sr$^{-1}$) of the $n_{\text{el}}$ clouds. The polarization fractions, $\bar{P}_{d,i}$, and polarization angles, $\bar{\psi}_{d,i}$, of the clouds are indicated by the length and orientation of the overplotted white headless vectors. To help make the connection between spectra and maps, each map is topped with a thick band having the same color as the associated spectrum. **Right:** maps of the resulting polarized intensity, $P_d^{\text{mod}}$, LoS-averaged polarization fraction, $P_d^{\text{los,mod}}$, and LoS-averaged polarization angle, $\psi_d^{\text{los,mod}}$, all in the same units as in Fig. 9. For reference, the corresponding observational Planck maps are shown in the top row.</figcaption></figure>
CO core (CO [1–2]) at $C_{\text{th}} = 32\%$ and the knotty HI [4] compo-
nent at $C_{\text{th}} = 42\%$. At $C_{\text{th}} = 44\%$, an oblique, elongated chunk
(HI [5–6], CO [3–4]) splits off from the dominant cloud, leav-
ing behind a mostly atomic cloud (HI [7–19], CO [5–6]). Further
fissions successively occur above $C_{\text{th}} = 51\%$ until every cloud
reduces to a single kinematic component.

The general quality of the fit, measured by the value of $\chi_r$, gradually improves as $C_{\text{th}}$ increases from 0 to 42% – and, accordingly, $n_{\text{cl}}$ increases from 1 to 6. This is because at each division the two new clouds are allowed to take on more representative polarization parameters, presumably closer to reality. When $C_{\text{th}}$ reaches 44%, $\chi_r$ drops by 23% and the fit suddenly looks much better. What happens is that the oblique, elongated cloud (HI [5–6], CO [3–4]) breaks away with a polarization angle very different from that of its parent cloud. As a result, it acts as a depolarizing cloud, which manages to reproduce a large portion of the region of lower polarization in the map of $p_d^{\text{los,obs}}$ (second-to-last map in the top row of Fig. 10) – in the same way as the bright CO core was argued in Sect. 3.4 to explain the localized low-polarization spot in the map of $p_d^{\text{los,obs}}$. In the rest of this subsection, these two depolarizing clouds are referred to as $C_4^{[7]}$ and $C_6^{[7]}$, respectively, based on their names in the seven-cloud configuration.

As $C_{\text{th}}$ keeps increasing above 51%, the emergence of new
clouds only leads to marginal improvements in the reconstructed
polarization maps, with no significant decrease of $\chi_r$. Most
importantly, the polarization parameters of $C_4^{[7]}$ and $C_6^{[7]}$ remain
stable until these clouds themselves break apart, which we take
as evidence that they have nearly reached their true values. Phys-
ically, when $C_4^{[7]}$ becomes a separate cloud, it directly takes on
the polarization parameters that make it possible to reproduce
not only the observed Stokes parameters in its own direction, but
also the general polarized background against which $C_6^{[7]}$ acts
as a small depolarizing cloud. $C_6^{[7]}$ then automatically adjusts its
polarization parameters to reproduce the low-polarization spot
observed at its location.

To sum up, the status of $C_4^{[7]}$ and $C_6^{[7]}$ as depolarizing clouds
appears to be robust, and the best-fit values of their polarization
parameters can be considered to be trustworthy.

**4. Discussion and conclusions**

In this paper, we present a new method designed (1) to identify
along the LoS the different clouds that contribute to the observed
dust emission and (2) to estimate the orientations of their inter-
nal magnetic fields. The cloud identification is performed with
the help of three kinematic gas tracers: the HI 21 cm, ¹²CO
2.6 mm, and ¹³CO 2.7 mm emission lines. The 3D spectral cubes
of these three tracers are corrected for opacity saturation, the
corrected ¹²CO and ¹³CO cubes are combined into a single CO
cube, and the HI and combined CO cubes are each decomposed
with the algorithm ROHSA into several spatially coherent Gaus-
sian kinematic components. All the kinematic components from
both tracers are rescaled to dust intensity, and those (from either
tracer) with similar velocity profiles are grouped together into
clouds. The result is a set of $n_{\text{cl}}$ clouds, $C_i$ ($i = 1, \dots, n_{\text{cl}}$), with
given dust intensities, $I_{d,i}(l, b)$. The estimation of their magnetic
field orientations rests on the linear polarization of the observed
dust emission, described by the two Stokes parameters, $Q_d$ and
$U_d$. By decomposing the latter in the basis formed by the 2D dust
intensities of the $n_{\text{cl}}$ clouds (Eqs. (38) and (39)), we can obtain

the polarization parameters (polarization fraction, $\bar{p}_{d,i}$, and polar-
ization angle, $\bar{\psi}_{d,i}$) of each cloud, $C_i$, which, in turn, lead to its
magnetic field orientation.

As an illustration of our method, we proposed a first applica-
tion to the G139 region, for which we had access to 353 GHz
Planck maps of the polarized dust emission as well as spec-
tral cubes of the HI 21 cm, ¹²CO 2.6 mm, and ¹³CO 2.7 mm
emission lines. The exact number of clouds present along the
LoS is poorly constrained, such that too much reality should
not be ascribed to the exact values of their individual polariza-
tion parameters (listed in Table 3). These parameters show great
variability: the polarization fractions span a broad range from
≈2% up to our imposed upper limit of 23% and the polariza-
tion angles span most of the range [−90°, +90°]. The dominant
cloud ($C_5$ in the seven-cloud configuration) consistently has a
very small polarization angle, which accounts for most of the
observed low-$\psi_d^{\text{los,obs}}$ background. The two molecular clouds that
stand out against this background in the seven-cloud configura-
tion ($C_6$ and $C_4$) have much larger polarization angles (≈65° and
≈45°, respectively), such that they act as depolarizing clouds and
produce the localized low-polarization spot and a large portion
of the more extended low-polarization region, respectively, in the
map of $p_d^{\text{los,obs}}$.

The broad ranges obtained for the clouds’ polarization
parameters imply that their magnetic fields are variously inclined
to the PoS, from nearly perpendicular to nearly parallel, and vari-
ously tilted to the Galactic plane, from nearly horizontal to nearly
vertical. The magnetic field of the dominant cloud is consistently
found to be nearly horizontal, in accordance with the general ori-
entation of the interstellar magnetic field in the Galactic disk.
The magnetic fields of the two depolarizing molecular clouds
are significantly tilted to the Galactic plane, which again is not
surprising given that the formation of molecular clouds is often
accompanied by rotation and shearing of magnetic field lines.
$C_6$ shows evidence that its bright center has rotated more than
its periphery, while $C_4$ appears to have its magnetic field nearly
aligned with its elongated main segment (at least in the PoS).

Our study highlights the need to incorporate the LoS dimen-
sion and to understand how different clouds along the LoS
contribute to the observed polarized dust emission. Had we only
relied on 2D sky maps of the polarized dust emission and tried to
read polarization angles directly off the map of $\psi_d^{\text{los,obs}}$, we would
have mistakenly concluded that the bright CO core causes no
more than a small distortion in a nearly horizontal magnetic field,
and we would have completely missed the contribution from the
oblique, elongated molecular cloud near the center of the G139
region. Despite the degeneracies in the LoS decomposition, the
values derived here for the polarization parameters of these two
clouds appear to be robust, as they remain stable from the seven-
cloud configuration where the oblique, elongated cloud becomes
a separate entity to the higher-$n_{\text{cl}}$ configuration where each cloud
breaks apart. These values also admit a plausible, simple physi-
cal interpretation. In principle, our method can be applied to any
region of the sky for which polarization and kinematic data are
available.

Below, we discuss the main assumptions underlying our
method and comment on its limitations.

First, our study relies on two kinds of tracers: polarized
thermal emission from dust and spectral lines of HI and CO
(including ¹²CO and ¹³CO). Our first important assumption
is that these two kinds of tracers are linearly related through
Eq. (20). This, in turn, assumes that (1) although dust accounts
for only ~1% in mass of the interstellar matter, it is well mixed in
with the gas, and its column density is proportional to that of the gas; (2) the emissivity of dust grains is uniform across each of the atomic and molecular media; (3) HI and CO are complementary tracers of the gas, which together detect all the gas along the LoS, with no omission and no overlap; and (4) the brightness temperatures of the considered spectral lines are properly corrected for opacity saturation. The validity of this multiple assumption can be assessed with the help of Fig. 4, which compares the dust intensity map reconstructed with the two gas tracers (top-middle panel) to the observational Planck map (left panel). The map of the residuals (top-right panel) indicates that the reconstruction is generally good to within a few percent, although three regions show residuals of up to ~20%. Residuals around the bright CO core suggest that our HI and CO tracers combined miss a small fraction of the gas along the LoS (the so-called dark gas). Other residuals are probably linked to our opacity corrections, which were shown to be quite sensitive to the excitation temperatures and their uncertainties. The adopted values of $T^{\text{HI}}$ and $\sigma(T^{\text{HI}})$ are particularly critical, affecting most of the G139 region. The values of $T^{\text{CO}}$ and $\sigma(T^{\text{CO}})$ could potentially have an even greater impact in saturated $^{12}\text{CO}$ regions, but this impact is considerably reduced by our using a combination of $^{12}\text{CO}$ and $^{13}\text{CO}$ data. Another, more minor source of residuals could be that dust emission does not perfectly trace the gas distribution (nonuniform conversion factors).

Second, the matter associated with the observed dust emission is decomposed, on purely kinematic grounds into several clouds along the LoS. This LoS decomposition is by no means unique: both the Gaussian kinematic decomposition of each tracer with ROHSA and the grouping of all kinematic components into separate clouds involve free parameters whose values are chosen a bit arbitrarily. This is unavoidable, as the very definition of an interstellar cloud, the criteria used to locate its physical boundaries, and the extent to which adjacent structures can be grouped together into a single cloud are all somewhat subjective. The ROHSA decomposition of each tracer turns out to be very good, as the spatially-averaged spectrum and the velocity-integrated map reconstructed with all the kinematic components of each tracer are very close to their observational counterparts. The component grouping into clouds has no impact on the reconstructed dust intensity, but it does affect the reconstructed Stokes parameters for linear polarization. This is why we looked into different cloud configurations having from $n_{\text{cl}} = 1$ to 26 clouds. We found that the reconstruction of the Stokes parameters improves as $n_{\text{cl}}$ increases from 1 to 7, then remains quite stable as $n_{\text{cl}}$ increases even further.

Our method has a wide range of applications. One of them is to derive the relative orientation angles between filamentary structures and the ambient magnetic field in star-forming regions, with the aim of gaining insight into the role played by magnetic fields in the processes of filament formation and star formation. With that perspective in mind, we intend to apply our recently developed method for Filament Detection and Reconstruction at Multiple Scales (FilDReaMS; Carrière et al. 2022b,a) to a large sample of star-forming regions. This kind of application would benefit from higher angular-resolution data for both the gas distribution and dust polarization.

Our proposed new method requires a lot more testing and validation than we have provided in the present paper. Systematic tests should be conducted to examine the degeneracies in the LoS decomposition and their impact on the derived cloud parameters, $I_{d,i}$, $\bar{p}_{d,i}$, and $\bar{\psi}_{d,i}$. More generally, it would be very useful to identify the conditions under which the method produces reliable results and to assess the degree of accuracy of the results.

We could, for instance, start with simple toy models, involving an increasing number of partly overlapping homogeneous clouds. The method could then be applied to mock data and to simulations, with an increasing number of more realistic clouds. In each case, it would be possible to quantify the accuracy with which the cloud parameters are recovered and to determine how the validity of the reconstruction depends on various factors, such as the relative (total and polarized) intensities of the clouds, their spatial distribution, and the presence of depolarizing clouds.

Along the same lines, it would be very interesting to compare the results of our LoS decomposition to the results of tomographic decomposition with stellar polarization data (see, e.g., Pelgrims et al. 2024). Both methods have their intrinsic limitations, and comparing them would make it possible not only to improve each of them, but also to connect the RVs of the dust-emitting clouds to their distances.

Our results could also be confronted with complementary magnetic field observations, such as Zeeman measurements. This could help to lift (at least part of) the degeneracies and improve the robustness of our results. Moreover, this would lead to a more complete view of the magnetic fields of the dominant clouds.

**Acknowledgements.** We wish to extend our deepest thanks to Isabelle Ristorcelli, with whom we had many interesting and lively discussions related to our study. We also thank Mika Juvela, Dana Alina, Julien Montillaud, and Tie Liu, as well as the anonymous referee. This work was supported by the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. We made use of the data from the Milky Way Imaging Scroll Painting (MWISP) project, which is a multiline survey in $^{12}\text{CO}/^{13}\text{CO}/\text{C}^{18}\text{O}$ along the northern galactic plane with the PMO-13.7m telescope. We are grateful to all the members of the MWISP working group, particularly the staff members at the PMO-13.7m telescope, for their long-term support. MWISP was sponsored by National Key R&D Program of China with grants 2023YFA1608000 & 2017YFA0402701 and by CAS Key Research Program of Frontier Sciences with grant QYZDJ-SSW-SLH047.

## References

Alves, F. O., Frau, P., Girart, J. M., et al. 2014, A&A, 569, L1

Beck, R. 2001, Space Sci. Rev., 99, 243

Brentjens, M. A., & de Bruyn, A. G. 2005, A&A, 441, 1217

Burn, B. J. 1966, MNRAS, 133, 67

Carrière, J. S., Ferrière, K., Ristorcelli, I., & Montier, L. 2022a, A&A, 668, A42

Carrière, J. S., Montier, L., Ferrière, K., & Ristorcelli, I. 2022b, A&A, 668, A41

Chen, B. Q., Huang, Y., Yuan, H. B., et al. 2019, MNRAS, 483, 4277

Clark, S. E. 2018, ApJ, 857, L10

Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792

Digel, S. W., Lyder, D. A., Philbrick, A. J., Puche, D., & Thaddeus, P. 1996, ApJ, 458, 561

Doi, Y., Hasegawa, T., Bastien, P., et al. 2021, ApJ, 914, 122

Doi, Y., Nakamura, K., Kawabata, K. S., et al. 2024, ApJ, 961, 13

Draine, B. T., & Weingartner, J. C. 1997, ApJ, 480, 633

Erceg, A., Jelić, V., Haverkorn, M., et al. 2022, A&A, 663, A7

Erceg, A., Jelić, V., Haverkorn, M., et al. 2024, A&A, 687, A23

Ferrière, K. 2020, Plasma Phys. Controlled Fusion, 62, 014014

Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93

Grenier, I. A., Casandjian, J.-M., & Terrier, R. 2005, Science, 307, 1292

Han, J. L., Manchester, R. N., van Straten, W., & Demorest, P. 2018, ApJS, 234, 11

HI4PI Collaboration (Ben Bekhti, N., et al.) 2016, A&A, 594, A116

Hildebrand, R. H. 1983, QJRAS, 24, 267

Ho, P. T. P., Moran, J. M., & Rodriguez, L. F. 1982, ApJ, 262, 619

Hoang, T., & Lazarian, A. 2008, MNRAS, 388, 117

Hottier, C., Babusiaux, C., & Arenou, F. 2020, A&A, 641, A79

Hu, Y., & Lazarian, A. 2023, MNRAS, 524, 2379

Juvela, M., Ristorcelli, I., Montier, L. A., et al. 2010, A&A, 518, L93

Juvela, M., Ristorcelli, I., Pagani, L., et al. 2012, A&A, 541, A12

Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135
Lallement, R., Vergely, J. L., Babusiaux, C., & Cox, N. L. J. 2022, A&A, 661, A147

Leike, R. H., & Enßlin, T. A. 2019, A&A, 631, A32

Leike, R. H., Glatzle, M., & Enßlin, T. A. 2020, A&A, 639, A138

Marchal, A., Miville-Deschênes, M.-A., Orieux, F., et al. 2019, A&A, 626, A101

Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wyckoff, S. 2005, ApJ, 634, 1126

Montillaud, J., Juvela, M., Rivera-Ingraham, A., et al. 2015, A&A, 584, A92

O'Neill, T. J., Zucker, C., Goodman, A. A., & Edenhofer, G. 2024, ApJ, 973, 136

Padoan, P., Pelkonen, V. M., Juvela, M., Haugbølle, T., & Nordlund, Å. 2023, MNRAS, 522, 3548

Pattle, K., & Fissel, L. 2019, Front. Astron. Space Sci., 6, 15

Pelgrims, V., Clark, S. E., Hensley, B. S., et al. 2021, A&A, 647, A16

Pelgrims, V., Panopoulou, G. V., Tassis, K., et al. 2023, A&A, 670, A164

Pelgrims, V., Mandarakas, N., Skalidis, R., et al. 2024, A&A, 684, A162

Planck Collaboration XXIV. 2011a, A&A, 536, A24

Planck Collaboration XIX., Ade, P. A. R., Aghanim, N., et al. 2011b A&A, 536, A19

Planck Collaboration XI. et al. 2014, A&A, 571, A11

Planck Collaboration CV. 2015a, A&A, 576, A105

Planck Collaboration CIV. et al. 2015b, A&A, 576, A104

Planck Collaboration CXXXV. 2016a, A&A, 586, A135

Planck Collaboration VIII. 2016b, A&A, 594, A8

Planck Collaboration XXVIII., Ade, P. A. R., Aghanim, N., et al. 2016c A&A, 594, A28

Planck Collaboration LI. 2017, A&A, 599, A51

Planck Collaboration XII., Aghanim, N., Akrami, Y., et al. 2020a A&A, 641, A12

Planck Collaboration III., Aghanim, N., Akrami, Y., et al. 2020b, A&A, 641, A3

Planck Collaboration XI. 2020c, A&A, 641, A11

Remy, Q., Grenier, I. A., Marshall, D. J., & Casandjian, J. M. 2017, A&A, 601, A78

Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (Mörlenbach: John Wiley & Sons)

Sobey, C., Bilous, A. V., Grießmeier, J. M., et al. 2019, MNRAS, 484, 3646

Sokoloff, D. D., Bykov, A. A., Shukurov, A., et al. 1998, MNRAS, 299, 189

Straižys, V., & Laugalys, V. 2007, Baltic Astron., 16, 167

Su, Y., Yang, J., Zhang, S., et al. 2019, ApJS, 240, 9

Turić, L., Jelić, V., Jaspers, R., et al. 2021, A&A, 654, A5

Van Eck, C. L., Haverkorn, M., Alves, M. I. R., et al. 2017, A&A, 597, A98

Versteeg, M. J. F., Angarita, Y., Magalhães, A. M., et al. 2024, AJ, 167, 177

Ward-Thompson, D., Kirk, J. M., Crutcher, R. M., et al. 2000, ApJ, 537, L135

Wilson, T. L., Rohlfs, K., & Hüttemeister, S. 2013, Tools of Radio Astronomy

Winkel, B., Kerp, J., Flöer, L., et al. 2016, A&A, 585, A41

Yuan, L., Yang, J., Du, F., et al. 2021, ApJS, 257, 51

Yuan, L., Yang, J., Du, F., et al. 2022, ApJS, 261, 37

Zucker, C., Goodman, A., Alves, J., et al. 2021, ApJ, 919, 35
# Appendix A: Opacity corrections

## Appendix A.1: Opacity corrections for H1

For notational simplicity, superscript HI is dropped throughout this subsection.

In a homogeneous medium in thermal equilibrium, the brightness temperature of the H1 21 cm line at LoS velocity $v$, $T_b(v)$, can be written in terms of the kinetic temperature, $T$, and the optical depth, $\tau_v$, as

$$T_b(v) = T (1 - e^{-\tau_v}). \qquad (A.1)$$

In the optically thin regime ($\tau_v \ll 1$), $T_b(v) = T \tau_v$. In the general case, we can correct $T_b(v)$ for opacity saturation and define an opacity-corrected or "de-saturated" brightness temperature,

$$T_{b*}(v) = T \tau_v. \qquad (A.2)$$

If the medium is out of thermal equilibrium, $T$ should be understood as an excitation temperature (the spin temperature in this case). By combining Eqs. (A.1) and (A.2), we can easily relate $T_{b*}$ to $T_b$, either in terms of $\tau_v$,

$$T_{b*} = T_b \frac{\tau_v}{1 - e^{-\tau_v}}, \qquad (A.3)$$

or in terms of $T$,

$$T_b = T \left( 1 - \exp\left(-\frac{T_{b*}}{T}\right) \right), \qquad (A.4)$$

which is equivalent to

$$T_{b*} = -T \ln\left(1 - \frac{T_b}{T}\right). \qquad (A.5)$$

The optical depth, $\tau_v$, is directly related to the hydrogen column density, $N_H$, and the temperature, $T$,

$$\tau_v = C \frac{N_H}{T} \phi(v), \qquad (A.6)$$

where $C$ is a constant for the considered spectral line and $\phi(v)$ is the so-called line profile function (e.g., Rybicki & Lightman 1979). Substituting Eq. (A.6) into Eq. (A.2) and integrating over $v$ leads to a relation between the velocity-integrated opacity-corrected brightness temperature and the hydrogen column density,

$$\int T_{b*}(v) dv = C N_H, \qquad (A.7)$$

when use is made of the normalization condition $\int \phi(v) dv = 1$.

In an inhomogeneous medium, Eq. (A.7) remains valid, but Eq. (A.5) strictly makes sense only if the spin temperature is uniform. If this is not the case, it might still be possible to use Eq. (A.5), with $T$ replaced by an appropriate effective temperature, $T_{eff}$. Remember that HI gas in the ISM can basically be found in two different phases: a cold phase with $T_c \approx 80$ K and a warm phase with $T_w \approx 8000$ K (e.g., Ferrière 2020), where subscripts c and w denote the cold and warm phases, respectively. The total optical depth can then be written as

$$\tau_v = \tau_{v,c} + \tau_{v,w}, \qquad (A.8)$$

with, in view of Eq. (A.6),

$$\tau_{v,i} = C \frac{N_{H,i}}{T_i} \phi_i(v) = C \frac{f_i(v)}{T_i} N_H \phi(v), \qquad (A.9)$$

where $f_i(v)$ is the mass fraction of phase i at velocity $v$, for $i = c, w$. Similarly, Eq. (A.2) becomes

$$T_{b*}(v) = T_c \tau_{v,c} + T_w \tau_{v,w}, \qquad (A.10)$$

which can also be written as

$$T_{b*}(v) = T_{eff}(v) \tau_v, \qquad (A.11)$$

where the effective temperature is the weighted harmonic mean of $T_c$ and $T_w$:

$$T_{eff}(v) = \left( \frac{f_c(v)}{T_c} + \frac{f_w(v)}{T_w} \right)^{-1}. \qquad (A.12)$$

Clearly, unless the mass fraction of the cold phase at velocity $v$ is much smaller than that of the warm phase, $T_{eff} \approx T_c/f_c$.

The similarity between Eq. (A.11) and Eq. (A.2) suggests that the most reasonable approach is to adopt $T = T_{eff}$ in Eq. (A.5). Incidentally, this approach would be exact if the cold and warm phases were well mixed along the LoS, with a uniform mass ratio, because then $T_{eff}$ would also be uniform, and we would formally be back to the case of a homogeneous medium. The problem is that $T_{eff}$ varies across the H1 spectrum. However, aside from the requirement $T_{eff} \ge T_c$, $T_{eff}$ does not need to be determined with accuracy across the entire spectrum. The only spectral regions where the exact value of $T_{eff}$ is critical are regions where $T_b$ approaches (or exceeds) $T_c$. But these regions are almost always dominated by cold gas, so they can be described by $T_{eff} = T_c$. Spectral regions not dominated by cold gas have $T_b$ well below $T_c$, and hence well below $T_{eff}$, so their $T_{b*}$ inferred from Eq. (A.5) with $T = T_{eff}$ is little sensitive to the adopted value of $T_{eff}$, which can again be set to $T_c$. Altogether, $T_{b*}$ can be approximated by Eq. (A.5) with $T = T_{eff} = T_c = 80$ K throughout the spectrum. This choice for $T$ results in a small overestimation of $N_H$.

The uncertainty in $T_{b*}$, $\sigma(T_{b*})$, is the quadratic sum of the uncertainties arising from $T_b$ and from $T$:

$$\sigma^2(T_{b*}) = [\sigma^2(T_{b*})]_{T_b} + [\sigma^2(T_{b*})]_{T}, \qquad (A.13)$$

The uncertainty arising from $T_b$, $[\sigma(T_{b*})]_{T_b}$, directly follows from the measurement error in $T_b$, $\sigma(T_b)$:

$$[\sigma(T_{b*})]_{T_b} = \frac{\partial T_{b*}}{\partial T_b} \sigma(T_b) = \frac{T}{T - T_b} \sigma(T_b), \qquad (A.14)$$

where the second identity is obtained with the help of Eq. (A.5). Similarly, the uncertainty arising from $T$, $[\sigma(T_{b*})]_{T}$, reads

$$[\sigma(T_{b*})]_{T} = \frac{\partial T_{b*}}{\partial T} \sigma(T) = - \left[ \ln\left(1 - \frac{T_b}{T}\right) + \frac{T_b}{T - T_b} \right] \sigma(T). \qquad (A.15)$$

Here, we adopt $\sigma(T) = 20$ K (e.g., Ferrière 2020).
## Appendix A.2: Opacity corrections for ¹²CO

In this subsection, superscript ¹²CO is dropped unless there is a possible ambiguity with ¹³CO, in which case superscripts ¹²CO and ¹³CO are reduced to ¹² and ¹³, respectively.

In principle, we can derive the opacity-corrected brightness temperature of the ¹²CO 2.6 mm line, $T_{b\star}^{12}$, in the same way as explained in Sect. A.1 for the HI 21 cm line. There are, however, a few important differences.

First, Eq. (A.1) is strictly valid only in the Rayleigh-Jeans regime ($h\nu \ll kT$) and in the absence of a significant background radiation, both of which are satisfied in the case of HI, but not in the case of ¹²CO. Indeed, the higher frequency of the ¹²CO 2.6 mm line and the lower temperature of molecular gas cause small deviations from the Rayleigh-Jeans law and require taking the 2.73 K background radiation into account. Moreover, the observed ¹²CO signal may be subject to beam dilution. Under these conditions, the correct expression of the brightness temperature reads (e.g., Ho et al. 1982),

$$T_b(v) = \eta_b (J(T_{ex}) - J(T_{bg})) (1 - e^{-\tau_v}), \quad (A.16)$$

where $\eta_b$ is the beam dilution factor,

$$J(T) = \frac{T_0}{\exp\left(\frac{T_0}{T}\right) - 1} \quad (A.17)$$

is the intensity in temperature units, $T_0 = h\nu/k$ the intrinsic temperature, $T_{ex}$ the excitation temperature, and $T_{bg} = 2.73$ K the blackbody background temperature. In the case of the ¹²CO 2.6 mm line, $T_0 = 5.53$ K, $J(T_{ex}=20$ K) = 17.4 K, and $J(T_{bg}) = 0.84$ K. Although the correct expression is now Eq. (A.16), we are still allowed to use Eq. (A.1) provided $T$ refers to the corrected excitation temperature, $\eta_b[J(T_{ex}) - J(T_{bg})]$. Eq. (A.5) then provides a first estimate of $T_{b\star}^{12}$,

$$(T_{b\star}^{12})_1 = -T \ln\left(1 - \frac{T_b^{12}}{T}\right). \quad (A.18)$$

Second, in contrast to the HI 21 cm line, the ¹²CO 2.6 mm line is often strongly saturated ($T_b^{12} \to T$), with the implication that $(T_{b\star}^{12})_1$ is highly uncertain near the spectral peaks. This leads us to consider a second estimate of $T_{b\star}^{12}$, based on the optically much thinner ¹³CO 2.7 mm line,

$$(T_{b\star}^{12})_2 = r T_{b\star}^{13} = r \left[ -T \ln \left( 1 - \frac{T_b^{13}}{T} \right) \right], \quad (A.19)$$

where $r$ is the ¹²CO/¹³CO abundance ratio. In writing Eq. (A.19), we implicitly assume that ¹³CO has the same excitation temperature as ¹²CO, and we neglect the slight differences in $J(T_{ex})$ and $J(T_{bg})$ (Eq. (A.16)) and in $C$ (Eq. (A.7)) arising from the slight frequency difference between the ¹²CO and ¹³CO lines. We further assume that the ¹²CO/¹³CO abundance ratio is equal to the ¹²C/¹³C isotopic ratio, for which we adopt the value $r = 68 \pm 15$ (Milam et al. 2005). The advantage of considering both estimates of $T_{b\star}^{12}$ together is that they nicely complement each other: as the CO column density increases, $T_b^{12}$ approaches $T$, and thus $(T_{b\star}^{12})_1$ becomes more uncertain, but at the same time $T_b^{13}$ increasingly rises above the noise level, and thus $(T_{b\star}^{12})_2$ becomes more reliable.

The temperature, $T$, is a priori unknown, while its exact value is critical for the derivation of $(T_{b\star}^{12})_1$ (Eq. (A.18)) near the ¹²CO spectral peaks. Here, we consider that the ¹²CO line is fully saturated $(T_b^{12} \approx T)$ at the 3D positions of the ¹³CO spectral peaks exceeding a certain threshold, $(T_b^{13})_{thr}$,<fnref n="6" /> and we retain the small subset of fully saturated $T_b^{12}$. We might be tempted to assign their average value to $T$, but this would automatically lead to a number of unphysical voxels with $T_b^{12} > T$. Instead, we assign the largest fully saturated $T_b^{12}$, $(T_b^{12})_{max}$, to $T$ and the rms deviation of the fully saturated $T_b^{12}$ from $(T_b^{12})_{max}$ to the uncertainty in $T$, $\sigma(T)$. Our estimates for $T$ and $\sigma(T)$ are obviously biased toward the brightest ¹³CO-emitting regions, which are also the regions where the ¹²CO line is the most saturated, and hence where the opacity correction is the most critical. We note, however, that our particular choice $T = (T_b^{12})_{max}$ leads to a systematic underestimation of the opacity correction in these regions. We also note that our subset of fully saturated $T_b^{12}$ might miss unresolved bright sources, for which the measured $T_b^{13}$ is artificially reduced by beam dilution (factor $\eta_b$ in Eq. (A.16)).

The uncertainties in $(T_{b\star}^{12})_1$ and $(T_{b\star}^{12})_2$ are given by equations similar to Eq. (A.13):

$$\sigma^2((T_{b\star}^{12})_1) = [\sigma^2((T_{b\star}^{12})_1)]_{T_b^{12}} + [\sigma^2((T_{b\star}^{12})_1)]_T \quad (A.20)$$

and

$$\sigma^2((T_{b\star}^{12})_2) = [\sigma^2((T_{b\star}^{12})_2)]_{T_b^{12}} + [\sigma^2((T_{b\star}^{12})_2)]_T + [\sigma^2((T_{b\star}^{12})_2)]_r, \quad (A.21)$$

where (see Eqs. (A.14) and A.15))

$$[\sigma((T_{b\star}^{12})_1)]_{T_b^{12}} = \frac{T}{T - T_b^{12}} \sigma(T_b^{12}), \quad (A.22)$$

$$[\sigma((T_{b\star}^{12})_1)]_T = -\left[\ln\left(1 - \frac{T_b^{12}}{T}\right) + \frac{T_b^{12}}{T - T_b^{12}}\right] \sigma(T), \quad (A.23)$$

$$[\sigma((T_{b\star}^{12})_2)]_{T_b^{12}} = r \frac{T}{T - T_b^{12}} \sigma(T_b^{12}), \quad (A.24)$$

$$[\sigma((T_{b\star}^{12})_2)]_T = -r \left[ \ln\left(1 - \frac{T_b^{12}}{T}\right) + \frac{T_b^{12}}{T - T_b^{12}} \right] \sigma(T), \quad (A.25)$$

and

$$[\sigma((T_{b\star}^{12})_2)]_r = -T \ln\left(1 - \frac{T_b^{12}}{T}\right) \sigma(r). \quad (A.26)$$

In the above equations, $\sigma(T_b^{12})$ and $\sigma(T_b^{13})$ are the measurement errors in $T_b^{12}$ and $T_b^{13}$, respectively, $\sigma(T)$ is the uncertainty in $T$

<footnote n="6">⁶ For the application presented in Sect. 3, we adopt $(T_b^{13})_{thr} = 4.25$ K, to include the $n$ highest peaks that are clearly fully saturated in ¹²CO, as indicated by their $T_b^{12}/T_b^{13}$ ratio being very close to that of the strongest peak. The corresponding uncertainty in $T$, rounded to the nearest whole number, is $\sigma(T) = 5$ K.</footnote>
(discussed in the paragraph preceding Eq. (A.20)), and $\sigma(\mathbf{r})$ is the uncertainty in $\mathbf{r}$, for which we adopt $\sigma(\mathbf{r}) = 15$ (Milam et al. 2005, see paragraph following Eq. (A.19)).

It is important to realize that Eqs. (A.22) and (A.23) are strictly valid only in the limit of small errors. When $T_b^{12}$ approaches $T$, the small-error approximation fails, and Eqs. (A.22) and (A.23) should be replaced by

$$[\sigma((T_{b\star}^{12})_1)]_{T_b^{12}} = ((T_{b\star}^{12})_1)(T_b^{12} + \sigma(T_b^{12})) - ((T_{b\star}^{12})_1)(T_b^{12}) \quad (A.27)$$

and

$$[\sigma((T_{b\star}^{12})_1)]_T = ((T_{b\star}^{12})_1)(T - \sigma(T)) - ((T_{b\star}^{12})_1)(T), \quad (A.28)$$

respectively. In case $T_b^{12} + \sigma(T_b^{12}) \ge T$, the r.h.s. of Eq. (A.27) is undetermined (see Eq. (A.18)), and we simply set $[\sigma((T_{b\star}^{12})_1)]_{T_b^{12}}$ to an arbitrarily large value. Similarly for $[\sigma((T_{b\star}^{12})_1)]_T$ in case $T - \sigma(T) \le T_b^{12}$.

The best combined estimate of the opacity-corrected brightness temperature of the $^{12}$CO line, $T_{b\star}^{12}$, is the value that minimizes the sum of the weighted (by the inverse-variance) deviations squared from our two estimates, $(T_{b\star}^{12})_1$ and $(T_{b\star}^{12})_2$ – or, equivalently, the weighted (by the inverse-variance) mean of $(T_{b\star}^{12})_1$ and $(T_{b\star}^{12})_2$,

$$(T_{b\star}^{12CO})_{\text{best}} = \frac{\frac{(T_{b\star}^{12})_1}{\sigma^2((T_{b\star}^{12})_1)} + \frac{(T_{b\star}^{12})_2}{\sigma^2((T_{b\star}^{12})_2)}}{\frac{1}{\sigma^2((T_{b\star}^{12})_1)} + \frac{1}{\sigma^2((T_{b\star}^{12})_2)}}. \quad (A.29)$$

As expected, $(T_{b\star}^{12CO})_{\text{best}}$ is closer to the better-constrained of our two estimates, $(T_{b\star}^{12})_1$ and $(T_{b\star}^{12})_2$. The uncertainty in $(T_{b\star}^{12CO})_{\text{best}}$, $\sigma((T_{b\star}^{12CO})_{\text{best}})$, is given by

$$\sigma^2((T_{b\star}^{12CO})_{\text{best}}) = \left( \frac{1}{\sigma^2((T_{b\star}^{12})_1)} + \frac{1}{\sigma^2((T_{b\star}^{12})_2)} \right)^{-1}, \quad (A.30)$$

as is easily verified in the case of Gaussian distributions. Thus, the uncertainty in $(T_{b\star}^{12CO})_{\text{best}}$ depends only on the uncertainties in $(T_{b\star}^{12})_1$ and $(T_{b\star}^{12})_2$, not on the difference between $(T_{b\star}^{12})_1$ and $(T_{b\star}^{12})_2$. Moreover, $\sigma((T_{b\star}^{12CO})_{\text{best}})$ is smaller than each of $\sigma((T_{b\star}^{12})_1)$ and $\sigma((T_{b\star}^{12})_2)$, which is in line with our expectation that combining two different estimates of $T_{b\star}^{12}$ reduces the final uncertainty.

## Appendix B: Observational values of the LoS-averaged polarization fraction

The all-sky map of the 353 GHz polarized emission from Galactic dust presented by Planck Collaboration (2015b) revealed a large scatter in $p_d^{\text{los}}$, with measured values ranging from the noise limit up to a maximum, $(p_d^{\text{los}})_{\text{max}}$, which depends on the hydrogen column density, $N_H$: for $N_H \le 10^{21} \text{ cm}^{-2}$, $(p_d^{\text{los}})_{\text{max}} \approx 20\%$; for $10^{21} \text{ cm}^{-2} \le N_H \le 10^{22} \text{ cm}^{-2}$, $(p_d^{\text{los}})_{\text{max}}$ decreases steadily from $\approx 20\%$ to $\approx 12\%$; then $(p_d^{\text{los}})_{\text{max}}$ drops more steeply, down to $\le 5\%$ for $N_H \ge 2 \times 10^{22} \text{ cm}^{-2}$. Toward nearby dense cores,

with $N_H \ge 10^{22} \text{ cm}^{-2}$, $p_d^{\text{los}}$ also systematically decreases with increasing $N_H$.<fnref n="7" /> The observed scatter in $p_d^{\text{los}}$ arises from a combination of factors, including (1) spatial variations in the intrinsic polarization fraction and/or the efficiency of grain alignment, (2) spatial variations in the magnetic field inclination to the PoS, and (3) depolarization due to fluctuations in the PoS magnetic field orientation both along the LoS and across the telescope beam. The decrease of $(p_d^{\text{los}})_{\text{max}}$ with increasing $N_H$ can be attributed to either a gradual loss of dust grain alignment in dense regions shielded from the ambient UV radiation field (assuming radiative torques are largely responsible for dust grain alignment; see, e.g., Draine & Weingartner (1997); Hoang & Lazarian (2008)) or a gradual rise in the LoS + beam depolarization due to fluctuations in the PoS magnetic field orientation over long path lengths. In any case, the LoS-averaged $(p_d^{\text{los}})_{\text{max}}$ provides a lower limit to the local $(p_d^{\text{loc}})_{\text{max}}$.

In a companion paper, Planck Collaboration (2015a) analyzed the 353 GHz polarized dust emission from ten nearby molecular clouds located away from the Galactic plane and displaying a variety of physical conditions and polarization properties. They detected the largest $p_d^{\text{los}}$ in the most diffuse regions, and they confirmed the decrease of $(p_d^{\text{los}})_{\text{max}}$ with increasing $N_H$ above a certain threshold, the value of which actually depends on the considered cloud and is generally $\approx (1-3) \times 10^{21} \text{ cm}^{-2}$. The 2D distribution of $(p_d^{\text{los}})$, $N_H$) in each cloud falls below the upper envelope, $(p_d^{\text{los}})_{\text{max}}$ vs. $N_H$, of the all-sky distribution obtained by Planck Collaboration (2015b). Planck Collaboration (2015a) further compared their observed polarization maps to synthetic maps obtained with anisotropic, turbulent MHD simulations, in which they assumed a uniform polarization fraction parameter, $p_0 = 20\%$, corresponding to $(p_d^{\text{loc}})_{\text{max}} = 23\%$ (see their Eq. (8), where their intrinsic polarization fraction, $p_i$, corresponds to our $(p_d^{\text{loc}})_{\text{max}}$). The assumption of uniform $(p_d^{\text{loc}})_{\text{max}}$ is probably reasonable insofar as the dense cores appearing in the simulations are only weakly shielded from the ambient UV radiation. The simulations were able to reproduce the main statistical trends of the observed polarization maps (excluding dense cores), with the only source of depolarization arising from the magnetic field (both its inclination to the PoS and the spatial fluctuations in its PoS orientation), i.e., without invoking any loss of dust grain alignment. They also showed that a good estimate of $(p_d^{\text{loc}})_{\text{max}}$ could be inferred from the measured $(p_d^{\text{los}})_{\text{max}}$ in diffuse regions where B lies in the PoS and is ordered on large scales.

The statistical analysis of Planck Collaboration (2015b) was later refined by Planck Collaboration (2020a), who relied on the third public release of Planck data and debiased the estimates of the polarization fraction with the help of the modified asymptotic (MAS) estimator. The new results were found to be consistent with those of Planck Collaboration (2015b). In particular, the maximum value of the measured LoS-averaged polarization fraction was found to be $(p_d^{\text{los}})_{\text{max}} = 22.0^{+3.5}_{-1.4}\%$.

In dense molecular cores, $p_d^{\text{los}}$ appears to be much lower than in the large-scale diffuse medium. A number of studies reported observations of a drop in the polarization fraction toward dense cores (see, e.g., Ward-Thompson et al. (2000); Alves et al. (2014); and more references in the review by Pattle & Fissel (2019)). These so-called polarization holes could result

<footnote n="7">⁷ The derived values of $(p_d^{\text{los}})_{\text{max}}$ depend on the exact criterion used to define $(p_d^{\text{los}})_{\text{max}}$ and, more sensitively, on the angular resolution of the map. The values of $(p_d^{\text{los}})_{\text{max}}$ quoted here refer to the 99.99% percentile of $(p_d^{\text{los}})$ in the considered $N_H$ bin, at a resolution of 1° (see upper solid red curve in Fig. 2 of Planck Collaboration (2015a)).</footnote>
from a loss of grain alignment efficiency or from the tangling of magnetic field lines, within the beam or along the LoS. In a statistical analysis of Planck cold clumps (Planck Collaboration et al. 2016c), Ristorcelli et al. (in prep.) found a significant decrease of $p_d^{\text{los}}$ toward the clump centers, and they argued that this decrease is due to both grain misalignment and the magnetic field complex geometry. Typical values near clump centers are $p_d^{\text{los}} \sim 2-4\%$.

the LoS-averaged polarization fraction, $p_d^{\text{los,obs}}$, reduced by a depolarization factor, $\cos(2(\psi_d^{\text{los,obs}} - \bar{\psi}_{d,i}))$, due to fluctuations in the LoS-averaged polarization angle, $\psi_d^{\text{los,obs}}$, across the map. Eq. (C.5), for its part, shows that the best-fit polarization angle of cloud $C_i$, $\bar{\psi}_{d,i}$, is a weighted average over all the pixels of the LoS-averaged polarization angle, $\psi_d^{\text{los,obs}}$.

## Appendix C: Expressions of the best-fit polarization parameters

The procedure used to derive the polarization parameters of each cloud $C_i$ is to minimize the $\chi_r^2$ associated with the Stokes parameters for linear polarization (Eq. (40)), with $Q_d^{\text{mod}}$ and $U_d^{\text{mod}}$ given by Eqs. (38) and (39), respectively). Here, for simplicity, we assume $\sigma_Q = \sigma_U = \sigma$. Minimizing Eq. (40) with respect to $\bar{p}_{d,j}$, then with respect to $\bar{\psi}_{d,j}$, successively yields

$$ \sum_{n_{\text{pix}}} \frac{1}{\sigma^2} p_d^{\text{los,obs}} I_d^{\text{obs}} I_{d,j} \cos(2(\psi_d^{\text{los,obs}} - \bar{\psi}_{d,j})) = \sum_{n_{\text{pix}}} \frac{1}{\sigma^2} \sum_{i=1}^{n_{\text{cl}}} \bar{p}_{d,i} I_{d,i} I_{d,j} \cos(2(\bar{\psi}_{d,i} - \bar{\psi}_{d,j})) \quad (C.1) $$

and

$$ \sum_{n_{\text{pix}}} \frac{1}{\sigma^2} p_d^{\text{los,obs}} I_d^{\text{obs}} I_{d,j} \sin(2(\psi_d^{\text{los,obs}} - \bar{\psi}_{d,j})) = \sum_{n_{\text{pix}}} \frac{1}{\sigma^2} \sum_{i=1}^{n_{\text{cl}}} \bar{p}_{d,i} I_{d,i} I_{d,j} \sin(2(\bar{\psi}_{d,i} - \bar{\psi}_{d,j})). \quad (C.2) $$

These equations become physically more transparent in the case of a single cloud $C_i$, where they reduce to

$$ \bar{p}_{d,i} = \frac{\sum_{n_{\text{pix}}} p_d^{\text{los,obs}} \frac{I_{d,i}^2}{\sigma^2} \cos(2(\psi_d^{\text{los,obs}} - \bar{\psi}_{d,i}))}{\sum_{n_{\text{pix}}} \frac{I_{d,i}^2}{\sigma^2}} \quad (C.3) $$

and

$$ \sum_{n_{\text{pix}}} p_d^{\text{los,obs}} \frac{I_{d,i}^2}{\sigma^2} \sin(2(\psi_d^{\text{los,obs}} - \bar{\psi}_{d,i})) = 0, \quad (C.4) $$

when use is made of the identity $I_d^{\text{obs}} = I_{d,i}$. If $\psi_d^{\text{los,obs}}$ does not fluctuate too much across the map, such that $(\psi_d^{\text{los,obs}} - \bar{\psi}_{d,i})$ remains small, Eq. (C.4) can further be approximated by

$$ \bar{\psi}_{d,i} = \frac{\sum_{n_{\text{pix}}} \psi_d^{\text{los,obs}} \frac{p_d^{\text{los,obs}} I_{d,i}^2}{\sigma^2}}{\sum_{n_{\text{pix}}} \frac{p_d^{\text{los,obs}} I_{d,i}^2}{\sigma^2}}. \quad (C.5) $$

Equation (C.3) shows that the best-fit polarization fraction of cloud $C_i$, $\bar{p}_{d,i}$, is a weighted average over all the pixels of