© The Authors 2026

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 Fara-day 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.

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_⊥¹, 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_⊥ (or possibly B) in 3D. This type of tomographic mapping of B_⊥ 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°, 22.3°).

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, λ², times the rotation measure, $\text{RM\hspace{0.17em}}\propto {\displaystyle {\int }_0^L{n}_{\text{e}}}{B}_{\left|\right|}$, 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, $\text{DM\hspace{0.17em}=\hspace{0.17em}}{\displaystyle {\int }_0^L{n}_{\text{e}}}ds.$. The ratio of RM to DM directly yields the n_e-weighted average value of B_∥ between the pulsar and the observer. By combining the average values of B_∥ toward all pulsars with measured RM and DM (roughly 1 500 at the present time) with their measured distances, we can map out the large-scale 3D distribution of B_∥ (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 polar-ized radiation from a background radio source, we can exploit the λ²-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 λ² 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_⊥. 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_∥. 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_⊥, 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 (¹²CO and ¹³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^′ × 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 (¹²CO and ¹³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_{\text{d}}=\underset{0}{\overset{\infty }{{\displaystyle \int }}}{ℰ}_{\text{d}}(r)dr$$(1) and the two Stokes parameters for linear polarization, $$Q_{\text{d}}=\underset{0}{\overset{\infty }{{\displaystyle \int }}}{p}_{\text{d}}^{\text{loc}}(r){ℰ}_{\text{d}}(r)\cos \left(2{\psi }_{\text{d}}^{\text{loc}}(r)\right)dr$$(2) and $$U_{\text{d}}=\underset{0}{\overset{\infty }{{\displaystyle \int }}}{p}_{\text{d}}^{\text{loc}}(r){ℰ}_{\text{d}}(r)\sin \left(2{\psi }_{\text{d}}^{\text{loc}}(r)\right)dr,$$(3) where r denotes the LoS distance from the observer, subscript d stands for dust, E_d(r) is the local emissivity of the dust thermal emission, $p_{\text{d}}^{\text{loc}}(r)$ the local polarization fraction, and ${\psi }_{\text{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, $${\text{d}}_{}\equiv Q_{\text{d}}+i{U}_{\text{d}}=\underset{0}{\overset{\infty }{{\displaystyle \int }}}{p}_{\text{d}}^{loc}(r){ℰ}_{\text{d}}(r)e^{2i{\psi }_{\text{d}}^{loc}(r)}dr.$$(4)

The norm of the complex polarized intensity is generally referred to as the (real) polarized intensity, $$P_{\text{d}}=\left|{\rho }_{\text{d}}\right|=\sqrt{Q_{\text{d}}^2+U_{\text{d}}^2}.$$(5)

The local polarization fraction can be written as $$p_{\text{d}}^{\text{loc}}+{(}_{p_{\text{d}}^{\text{loc}}}{\text{cos2}}^{}{y}_{\text{B}},$$(6) where γ_B is the inclination angle of the local magnetic field to the PoS, ${(}_{p_{\text{d}}^{\text{loc}}}=(p_{\mathrm{i}\text{ntr}}R)$ is the theoretical maximum polarization fraction that can be achieved locally, p_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² (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, ψ_B, via $${\psi }_{\text{d}}^{loc}={\psi }_B\pm {90}^{\circ }.$$(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_{\text{d}}^{\text{los}}=\frac{P_{\text{d}}}{I_{\text{d}}}=\frac{\sqrt{Q_{\text{d}}^2+U_{\text{d}}^2}}{I_{\text{d}}},$$(8) and the LoS-averaged polarization angle by $${\psi }_{\text{d}}^{\text{los}}=\frac{1}{2}\arctan \left(\frac{U_{\text{d}}}{Q_{\text{d}}}\right),$$(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_{\text{d}}=p_{\text{d}}^{\text{los}}{I}_{\text{d}}\cos \left(2{\psi }_{\text{d}}^{\text{los}}\right)$$(10) and $$U_{\text{d}}=p_{\text{d}}^{\text{los}}{I}_{\text{d}}\sin \left(2{\psi }_{\text{d}}^{\text{los}}\right),$$(11) which can be rewritten in complex form as $${\text{d}}_{}\equiv Q_{\text{d}}+i{U}_{\text{d}}=p_{\text{d}}^{\text{los}}{I}_{\text{d}}{e}^{2i{\psi }_{\text{d}}^{\text{los}}}.$$(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_{\text{d}}^{\text{los}}=\frac{\left|{{\displaystyle \int }}_0^{\infty }{p}_{\text{d}}^{\text{loc}}(r){ℰ}_{\text{d}}(r)e^{2i{\psi }_{\text{d}}^{\text{loc}}(r)}dr\right|}{{{\displaystyle \int }}_0^{\infty }{ℰ}_{\text{d}}(r)dr}$$(13) and $$e^{2i{\psi }_{\text{d}}^{\text{los}}}=\frac{{{\displaystyle \int }}_0^{\infty }{p}_{\text{d}}^{loc}(r){ℰ}_{\text{d}}(r)e^{2i{\psi }_{\text{d}}^{\text{loc}}(r)}dr}{\left|{{\displaystyle \int }}_0^{\infty }{p}_{\text{d}}^{\text{loc}}(r){ℰ}_{\text{d}}(r)e^{2i{\psi }_{\text{d}}^{\text{loc}}(r)}dr\right|}.$$(14)

The physical meaning of Eqs. (13) and (14) is pretty straightforward: the LoS-averaged polarization fraction, $p_{\text{d}}^{\text{loc}}$ , is a dust-weighted LoS average of the local polarization fraction, $p_{\text{d}}^{\text{loc}}(r)$, reduced by a LoS depolarization factor, $e^{\text{2}i{\psi }_{\text{d}}^{\text{loc}}(r)}$, due to fluctua tions in the local polarization angle; the LoS-averaged polarization orientation, defined by the LoS-averaged polarization angle, ${\psi }_{\text{d}}^{\text{loc}}$, is a dust-weighted LoS average of the local polarization orientation, defined by the local polarization angle, ${\psi }_{\text{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_{\text{d}}^{\text{loc}}$ 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_{\text{d}}^{\text{loc}}(r)=p_{\text{d}}^{\text{loc}}$ and ${\psi }_{\text{d}}^{\text{loc}}(r)={\psi }_{\text{d}}^{\text{loc}}$. 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³, of the brightness temperature, T_b, of the Hi 21 cm, ¹²CO (1–0) 2.6 mm, and ¹³CO (1–0) 2.7 mm emission lines, based on the notion that the Hi line traces atomic gas, while the ¹²CO and ¹³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 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_{\text{b}}^{{12}_{\text{CO}}}}$ for CO (see paragraph preceding Eq. (A.20)). Once the values of T are set, we can convert the spectral cubes of $T_{\text{b}}^{{12}_{\text{CO}}}$, $T_{\text{b}}^{{13}_{\text{CO}}}$, and $T_{\text{b}\star }^{\text{HI}}$ into spectral cubes of ${(}_{T_{\text{b}\star }^{{12}_{\text{CO}}}}$, and ${(}_{T_{\text{b}\star }^{{12}_{\text{CO}}}}$, respectively, where ${(}_{T_{\text{b}\star }^{{12}_{\text{CO}}}}$ and ${(}_{T_{\text{b}\star }^{{12}_{\text{CO}}}}$ are two complementary estimates of $T_{\text{b}\star }^{{12}_{\text{CO}}}$. We can then combine the ${(}_{T_{\text{b}\star }^{{12}_{\text{CO}}}}$ and ${(}_{T_{\text{b}\star }^{{12}_{\text{CO}}}}$ spectral cubes into the spectral cube of a best-estimate $T_{\text{b}\star }^{{12}_{\text{CO}}}$ (Eq. (A.29)). From now on, superscript 12 in ¹²CO is dropped for notational simplicity.

The spectral cubes of $T_{\text{b}\star }^{\text{HI}}$ and $T_{\text{b}\star }^{\text{CO}}$ generally have different angular and spectral resolutions, so they first need to be brought to a common angular resolution, δθ, a common spectral resolution, δv, and a common (l, b, v) grid with, say, n₁ × n₂ pixels. The maps of I_d, Q_d, and U_d (see Sect. 2.1) also need to be brought to the angular resolution δθ and the n₁ × n₂ pixel grid.

We assume that our $T_{\text{b}\star }^{\text{HI}}$ and $T_{\text{b}\star }^{\text{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_{\text{H}}=N_{\text{H}}^{{\text{H}}_I}+N_{\text{H}}^{\text{CO}}.$$(15)

The intensity of the dust emission can be decomposed in a similar way, $$I_{\text{d}}=I_{\text{d}}^{{\text{HI}}_{}}+I_{\text{d}}^{\text{CO}}.$$(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_{\text{dA}}^{}$, is simply proportional to the hydrogen column density probed by that tracer, $N_{\text{HA}}^{}$ (e.g., Hildebrand 1983). If we denote the conversion factor from hydrogen column density to dust intensity by $X_{I_{\text{d}}/N_{\text{H}}}^{\text{A}}$, we can then write $$I_{\text{dA}}^{}=X_{I_{\text{d}}/N_{\text{H}}}^{\text{A}}{N}_{\text{HA}}^{}.$$(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_{\text{HA}}^{}$, is proportional to the opacity-corrected brightness temperature, $T_{\text{b}\star }^{\text{A}}(\upsilon )$, integrated over the RV (see Eq. (A.7) in Appendix A), $$N_{\text{HA}}^{}=X_{N_{\text{H}}/T_{\text{b}}}^{\text{A}}{{\displaystyle \int }}^{\text{}}{T}_{\text{b}\star }^{\text{A}}(v)dv,$$(18) where $X_{N_{\text{H}}/T_{\text{b}}}^{\text{A}}$ is the conversion factor from velocity-integrated opacity-corrected brightness temperature to hydrogen column density. Combining Eqs. (17) and (18) leads to $$I_{\text{dA}}^{}=X_{I_{\text{d}}/T_{\text{b}}}^{\text{A}}{{\displaystyle \int }}^{\text{}}{T}_{\text{b}\star }^{\text{A}}(v)dv$$(19) with $X_{I_{\text{d}}/T_{\text{b}}}^{\text{A}}=X_{I_{\text{d}}/N_{\text{H}}}^{\text{A}}{X}_{N_{\text{H}}/T_{\text{b}}}^{\text{A}}$. If we now insert Eq. (19) into Eq. (16), we finally obtain for the dust intensity $$I_{\text{d}}=X_{I_{\text{d}}/T_{\text{b}}}^{HI}{{\displaystyle \int }}^{\text{}}{T}_{\text{b}\star }^{HI}dv+X_{I_{\text{d}}/T_{\text{b}}}^{CO}{{\displaystyle \int }}^{\text{}}{T}_{\text{b}\star }^{CO}dv$$(20)

2.2.2 Conversion factors

The observables in Eq. (20) are the dust intensity, I_d, and the opacity-corrected brightness temperatures, $T_{\text{b}\star }^{\text{HI}}$ and $T_{\text{b}\star }^{\text{CO}}$. The conversion factors, $X_{I_{\text{d}}/T_{\text{b}}}^{\text{HI}}$ and $X_{I_{\text{d}}/T_{\text{b}}}^{\text{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_Id/NH and X_NH/Tb (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_{\text{d}}/T_{\text{b}}}^{\text{HI}}$ and $X_{I_{\text{d}}/T_{\text{b}}}^{\text{CO}}$, as free parameters, and we derive their best-fit values in the considered region by minimizing the reduced χ², defined by $${\chi }_{\text{r}}^2=\frac{1}{\left(n_{dat}-n_{par}\right)}{\displaystyle \sum }_{n_{pix}}\frac{{\left[I_{\text{d}}^{obs}-\left(I_{\text{d}}^{HI}+I_{\text{d}}^{CO}\right)\right]}^2}{{\sigma }_I^2},$$(21) where $I_{\text{d}}^{\text{obs}}$ is the observed dust intensity, $I_{\text{d}}^{\text{HI}}$ and $I_{\text{d}}^{\text{CO}}$ are the dust intensities associated with Hi and CO, respectively (Eq. (19) with A = Hi and CO), σ_I is the total observational uncertainty, n_dat = n_pix = n₁ × n₂ is the total number of data points, n_par = 2 is the number of free parameters, and the sum runs over the n_pix pixels. The total observational uncertainty is the quadratic sum of the measurement errors in $I_{\text{d}}^{\text{obs}},I_{\text{d}}^{\text{HI}}$, and $I_{\text{d}}^{\text{CO}}$, $${\sigma }_I^2={\sigma }^2\left(I_{\text{d}}^{\text{obs}}\right)+{\sigma }^2\left(I_{\text{d}}^{{\text{HI}}_{}}\right)+{\sigma }^2\left(I_{\text{d}}^{\text{CO}}\right),$$(22) where $${\sigma }^2\left(I_{\text{dA}}^{}\right)={\left(X_{I_{\text{d}}/T_{\text{b}}}^{\text{A}}\right)}^2\left({\displaystyle \sum }_v{\sigma }^2\left(T_{\text{b}\star }^{\text{A}}\right)\right)\delta v^2$$(23) for A = Hi and CO, ${\sigma }^2\left(T_{\text{b}\star }^{\text{HI}}\right)$ and ${\sigma }^2(T_{\text{b}\star }^{\text{CO}})$ are given by Eqs. (A.13) and (A.30), respectively, δv is the spectral resolution, and the sum runs over all the velocity bins.

Minimization of ${\chi }_{\text{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_{\text{d}}/T_{\text{b}}}^{\text{HI}})$ and $\sigma (X_{I_{\text{d}}/T_{\text{b}}}^{\text{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_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_cl clouds, $$I_{\text{d}}^{\mathrm{mod}}={\displaystyle \sum }_{i=1}^{n_{cl}}{I}_{\text{d},i},$$(24) and similarly for the two Stokes parameters for linear polarization, $$Q_{\text{d}}^{\mathrm{mod}}={\displaystyle \sum }_{i=1}^{n_{cl}}{Q}_{\text{d},i},$$(25) and $$U_{\text{d}}^{\mathrm{mod}}={\displaystyle \sum }_{i=1}^{n_{cl}}{U}_{\text{d},i},$$(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_{\text{b}\star }^{\text{HI}}$ and $T_{\text{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_{\text{b}\star }^{\text{HI}}$ and $T_{\text{b}\star }^{\text{CO}}$ spectral cubes separately. This algorithm is designed to decompose a spectral data cube, say, a cube of T_b⋆(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⋆(l, b, v) is approximated by a modeled ${\tilde{T}}_{\text{b}\star }(l,b,v)$, equal to the sum of n_G Gaussian components, $${\tilde{T}}_{\text{b}\star }(l,b,v)={\displaystyle \sum }_{j=1}^{n_{\text{G}}}{\tilde{T}}_{\text{b}\star ,j}(l,b,v),$$(27) where each component j is described by an amplitude Aj(l, b), a mean RV, ${\overline{v}}_j(l,b)$, and a standard deviation, σj(l, b), $${\tilde{T}}_{\text{b}\star ,j}(l,b,v)=A_j(l,b)\exp \left[-\frac{{\left(v-{\overline{v}}_j(l,b)\right)}^2}{2{\sigma }_j^2(l,b)}\right].$$(28)

The best-fit values of the 3n_G Gaussian parameters, Aj, ${\overline{v}}_j$, and σj, are derived through minimization of a cost function that includes the standard χ² 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, λ_A, ${\lambda }_{\overline{v}}$, and λ_σ. Since our ultimate purpose is to construct smooth, coherent clouds, we need to impose strong constraints of spatial coherence on A_j, ${\overline{v}}_j$, and σ_j, which is done by choosing large values for λ_A, ${\lambda }_{\overline{v}}$, and λ_σ. After verifying that the exact values are not critical, we adopt λ_A = λ_σ = 1000 and ${\lambda }_{\overline{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}}_{\text{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_{\text{G}}^{\text{tot}}=n_{\text{G}}^{\text{HI}}+n_{\text{G}}^{\text{CO}}$ Gaussian components, which we consider two by two. For every pair of components j and j^′ (j ≠ j^′), we compute a velocity correlation coefficient at each pixel, (l, b), $$C_{j{j}^{\prime }}(l,b)=\frac{{{\displaystyle \sum }}_v{\tilde{T}}_{\text{b}\star ,j}{\tilde{T}}_{\text{b}\star ,j^{\prime }}}{\sqrt{{\displaystyle \sum _v{\tilde{T}}_{\text{b}\star ,j}^2}}\sqrt{{\displaystyle \sum _v{\tilde{T}}_{\text{b}\star ,j^{\prime }}^2}}},$$(29) where the sums run over all the velocity bins. Clearly, C_jj′ depends only on the velocity profiles (mean velocities, ${\overline{v}}_j$ and ${\overline{v}}_{j^{\prime }}$ , and standard deviations, σ_j and σ_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_{j{j}^{\prime }}\rangle =\frac{{\displaystyle \sum _{n_{pix}}\left({\tilde{I}}_{\text{d},j}+{\tilde{I}}_{\text{d},j^{\prime }}\right)C_{j{j}^{\prime }}}}{{\displaystyle \sum _{n_{pix}}\left({\tilde{I}}_{\text{d},j}+{\tilde{I}}_{\text{d},j^{\prime }}\right)}},$$(30) where ${\tilde{I}}_{\text{d},j}$ is the dust intensity of component j (defined by Eq. (33) below) and the sum runs over the n_pix pixels.

We consider that components j and j^′ are correlated, and thus parts of a same cloud, if ⟨𝒞_jj’⟩ lies above a certain threshold, 𝒞_th $$\langle C_{j{i}^{\prime }}\rangle \ge C_{\text{th}}.$$(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 𝒞_th. For our application in Sect. 3, we tested all the integer values of 𝒞_th between 0 and 100%. Based on the results (summarized in Table 3), we decided to present a detailed analysis in the case 𝒞_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 𝒞_th as well as their sensitivity to the exact value of 𝒞_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 ⟨𝒞_jk⟩ ≥ 𝒞_th and ⟨𝒞_kl⟩ ≥ 𝒞_th, they are grouped into one cloud enclosing components j,k, and l, even if ⟨𝒞_jl⟩ < 𝒞_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_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_{\text{d},i}={\displaystyle \sum }_j{\tilde{I}}_{\text{d},j},$$(32) where the dust intensity of component j, Ĩ_d,j, is related to its opacity-corrected brightness temperature, ${\tilde{T}}_{\text{b}\star ,j}$ (Eq. (28)), through an equation similar to Eq. (19): $${\tilde{I}}_{\text{d},j}=X_{I_{\text{d}}/T_{\text{b}}}^{\text{A}}{{\displaystyle \int }}^{\text{}}{\tilde{T}}_{\text{b}\star ,j}(v)dv$$(33)

Superscript A in the conversion factor $X_{I_{\text{d}}/T_{\text{b}}}^{\text{A}}$ refers to the gas tracer ( Hi or CO) associated with component j. As a reminder, the best-fit value of $X_{I_{\text{d}}/T_{\text{b}}}^{\text{A}}$ and the attendant uncertainty, $\sigma (X_{I_{\text{d}}/T_{\text{b}}}^{\text{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_{\text{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_{\text{d}}^{\mathrm{mod}}$ and $U_{\text{d}}^{\mathrm{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_{\text{d},i}=p_{\text{d},i}{I}_{\text{d},i}\cos \left(2{\psi }_{\text{d},i}\right)$$(34) and $$U_{\text{d},i}=p_{\text{d},i}{I}_{\text{d},i}\sin \left(2{\psi }_{\text{d},i}\right),$$(35) with $$p_{\text{d},i}=\frac{\left|{{\displaystyle \int }}_{L_i}{p}_{\text{d}}^{\text{loc}}(r){ℰ}_{\text{d}}(r)e^{2i{\psi }_{\text{d}}^{\text{loc}}(r)}dr\right|}{{{\displaystyle \int }}_{L_i}{ℰ}_{\text{d}}(r)dr}$$(36) and $$e^{2i{\psi }_{\text{d},i}}=\frac{{{\displaystyle \int }}_{L_i}{p}_{\text{d}}^{\text{loc}}(r){ℰ}_{\text{d}}(r)e^{2i{\psi }_{\text{d}}^{\text{loc}}(r)}dr}{\left|{{\displaystyle \int }}_{L_i}{p}_{\text{d}}^{\text{loc}}(r){ℰ}_{\text{d}}(r)e^{2i{\psi }_{\text{d}}^{\text{loc}}(r)}dr\right|}.$$(37)

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

To proceed, we consider that p_d,i and ψ_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 ψ_d,i at the cloud scale, which we denote with an overscore⁴. Thus, the pair of equations we work here with is obtained by inserting Eqs. (34) and (35), with p_d,i and ψ_d,i overscored, into Eqs. (25) and (26): $$Q_{\text{d}}^{\mathrm{mod}}={\displaystyle \sum }_{i=1}^{n_{\text{cl}}}{\overline{p}}_{\text{d},i}{I}_{\text{d},i}\cos \left(2{\overline{\psi }}_{\text{d},i}\right)$$(38) and $$U_{\text{d}}^{\mathrm{mod}}={\displaystyle \sum }_{i=1}^{n_{\text{cl}}}{\overline{p}}_{\text{d},i}{I}_{\text{d},i}\sin \left(2{\overline{\psi }}_{\text{d},i}\right).$$(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, ${\overline{p}}_{\text{d},i}$ and ${\overline{\psi }}_{\text{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 ${\overline{p}}_{\text{d},i}$ and ${\overline{\psi }}_{\text{d},i}$ are the values that minimize $${\chi }_{\text{r}}^2=\frac{1}{\left(n_{\text{dat}}-n_{\text{par}}\right)}{\displaystyle \sum }_{n_{pix}}\left[\frac{{\left(Q_{\text{d}}^{\text{obs}}-Q_{\text{d}}^{\mathrm{mod}}\right)}^2}{{\sigma }_Q^2}+\frac{{\left(U_{\text{d}}^{\text{obs}}-U_{\text{d}}^{\mathrm{mod}}\right)}^2}{{\sigma }_U^2}\right],$$(40) where $Q_{\text{d}}^{\text{obs}}$ and $U_{\text{d}}^{\text{obs}}$ are the observed Stokes parameters, $Q_{\text{d}}^{\mathrm{mod}}$ and $U_{\text{d}}^{\mathrm{mod}}$ are the modeled Stokes parameters given by Eqs. (38) and (39), respectively, σ_Q and σ_U are the total “observational uncertainties” associated with Q_d and U_d, respectively, n_dat = 2 n_pix is the total number of data points, n_par = 2 n_cl is the number of free parameters, and the sum runs over the n_pix pixels. The total “observational uncertainties” are the quadratic sums of the measurement errors in the observed Stokes parameters, $(Q_{\text{d}}^{\text{obs}})$ and $(U_{\text{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_{\text{d}}^{\mathrm{mod}}$ and $U_{\text{d}}^{\mathrm{mod}}$ are both equal to the error in the modeled intensity, $I_{\text{d}}^{\text{mod}}$, times the observed LoS-averaged polarization fraction, $p_{\text{d}}^{\text{los,obs}}$. We may further approximate the error in $I_{\text{d}}^{\text{mod}}$ by the residual $\text{Δ(}{I}_{\text{d}})=I_{\text{d}}^{\text{obs}}-I_{\text{d}}^{\mathrm{mod}}$. Altogether, we have $${\sigma }_Q^2={\sigma }^2\left(Q_{\text{d}}^{\text{obs}}\right)+{\left(p_{\text{d}}^{\text{los,obs}}\right)}^2{\text{Δ}}^2\left(I_{\text{d}}\right)$$(41) and $${\sigma }_U^2={\sigma }^2\left(U_{\text{d}}^{\text{obs}}\right)+{\left(p_{\text{d}}^{\text{los,obs}}\right)}^2{\text{Δ}}^2\left(I_{\text{d}}\right).$$(42)

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

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

The inclination angle of ${\overline{B}}_i$ to the PoS, ${\overline{\gamma }}_{B,i}$, can in principle be inferred from the polarization fraction, ${\overline{p}}_{\text{d},i}$, using Eq. (6) applied to cloud i, $${\cos }^2{\overline{\gamma }}_{B,i}=\frac{{\overline{p}}_{\text{d},i}}{{\left({\overline{p}}_{\text{d},i}\right)}_{\max }}.$$(44) together with an adopted value of the maximum polarization fraction averaged at the cloud scale, ${\left({\overline{p}}_{\text{d},i}\right)}_{\max }$. Following Planck Collaboration (2015a), we could, for instance, take ${\left({\overline{p}}_{\text{d},i}\right)}_{\max }=23\%$ (see Appendix B). Because of the uncertainty in ${\left({\overline{p}}_{\text{d},i}\right)}_{\max }$, the derived value of ${\overline{\gamma }}_{B,i}$ is much less reliable than the derived value of ${\overline{\psi }}_{B,i}$. Moreover, the existence of two opposite-signed solutions for ${\overline{\gamma }}_{B,i}$ implies that the orientation of ${\overline{B}}_i$ can only be determined with a mirror ambiguity with respect to the PoS. Zee-man observations, which are sensitive to B_∥, would be needed to set the sign of ${\overline{\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°30^′, −3°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, δθ, and (when relevant) their spectral resolution, δv, and their spectral extent, ∆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, ¹²CO 2.6 mm, and ¹³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 δ ≥ −5° (Winkel et al. 2016), and the ¹²CO and ¹³CO cubes from the Milky Way Imaging Scroll Painting (MWISP) survey of the Galactic plane, with 10° ≤ l ≤ +250° and |b| ≲ 5.2° (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 ¹²CO (middle panel) and ¹³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 ¹²CO peak does not stand out as prominently as the ¹³CO peak, which is most likely because ¹²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_{\text{b}*}^{\text{HI}}$ using Eq. (A.5) with T = 80 K. For CO, we derive two complementary estimates of $T_{\text{b}*}^{12\text{CO}}$ assuming $T={\left(T_{\text{b}}^{12\text{CO}}\right)}_{\max }:{\left(T_{\text{b}*}^{12\text{CO}}\right)}_1$ based on the ¹²CO cube (Eq. (A.18)) and ${\left(T_{\text{b}*}^{12\text{CO}}\right)}_2$ based on the rescaled ¹³CO cube (Eq. (A.19)); we then combine ${\left(T_{\text{b}*}^{12\text{CO}}\right)}_1$ and ${\left(T_{\text{b}*}^{12\text{CO}}\right)}_2$ to obtain a best-estimate $T_{\text{b}*}^{12\text{CO}}$ (Eq. (A.29)), from now on referred to as $T_{\text{b}*}^{\text{CO}}$. We find that $T_{\text{b}*}^{\text{CO}}$ is nearly equal to the estimate with the smaller uncertainty, which is ${\left(T_{\text{b}*}^{12\text{CO}}\right)}_1$ away from the CO peak and ${\left(T_{\text{b}*}^{12\text{CO}}\right)}_2$ toward the CO peak. This dichotomy results from the huge jump in the opacity correction to $T_{\text{b}}^{\text{CO}}$ toward the CO peak, which renders ${\left(T_{\text{b}*}^{12\text{CO}}\right)}_1$ extremely uncertain.

The maps of the velocity-integrated $T_{\text{b}*}^{\text{HI}}$ and $T_{\text{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 ¹²CO and ¹³CO maps (not shown) are quite similar, except toward the CO peak, where only $T_{\text{b}*}^{\text{13CO}}$ is reliable. Altogether, the combined opacity-corrected CO map (right panel in the second row) is very close to the opacity-corrected ¹²CO map, with only the CO peak taken from the rescaled opacity-corrected ¹³CO map.

Before jointly exploiting the $T_{\text{b}*}^{\text{HI}}$ and $T_{\text{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, δθ = 10.8^′, together with a pixel size at the (rounded) Nyquist limit, δθ_pix = 4^′. Accordingly, the (l, b) grid needed to cover our 104^′ × 104^′ G139 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 spectral resolution δv = 0.3 km s⁻¹, which is much finer than the spectral resolution of the Hi cube and roughly twice the spectral resolution of the CO cube. This spectral oversampling of $T_{\text{b}*}^{\text{HI}}$ makes it easier to extract meaningful Gaussian kinematic components with ROHSA. The re-gridded maps of I_d, Q_d, and U_d and those of the velocity-integrated $T_{\text{b}*}^{\text{HI}}$ and $T_{\text{b}*}^{\text{CO}}$ are displayed in the second row of Fig. 1 and the third row of Fig. 2, respectively.

Fig. 1 Maps of the intensity, Id (left), and of the two Stokes parameters for linear polarization, Qd (middle) and Ud (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 Id for Id(Planck Collaboration 2016b, 2020b) and 0.015 Pd for Qd and Ud (Planck Collaboration 2020b,c).

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_{I_{\text{d}}/T_{\text{b}}}^{\text{A}}$ (see Eq. (19)). As explained in Sect. 2.2, the best-fit values of $X_{I_{\text{d}}/T_{\text{b}}}^{\text{HI}}$ and $X_{I_{\text{d}}/T_{\text{b}}}^{\text{CO}}$, together with their uncertainties, are obtained by minimizing ${\chi }_{\text{r}}^2$ in Eq. (21) through bootstrap sampling + MCMC simulations.

We find that the best fit has χ_r = 1.63. This small value of χ_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, σ_I (Eq. (22)). The latter is dominated by $\sigma \left(I_{\text{d}}^{\text{HI}}\right)$, which generally exceeds $\sigma \left(I_{\text{d}}^{\text{obs}}\right)$ and $\sigma \left(I_{\text{d}}^{\text{CO}}\right)$ by a factor ∼2–10, except toward the CO peak, where $\sigma \left(I_{\text{d}}^{\text{CO}}\right)~2\sigma \left(I_{\text{d}}^{\text{HI}}\right)~3.5\sigma \left(I_{\text{d}}^{\text{obs}}\right)$. The dominant contribution to $\sigma \left(I_{\text{d}}^{\text{HI}}\right)$, 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 \left(I_{\text{d}}^{\text{CO}}\right)$ 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}^{\text{HI}}$ and $X_{I_d/T_b}^{\text{CO}}$. Their best-fit values and standard deviations, written above the corresponding 1D histograms (upper-left and lower-right panels), are $X_{I_d/T_b}^{\text{HI}}=\left(0.00081\pm 0.000006\right)\left({\text{MJy\hspace{0.17em}sr}}^{-1}\right){\left({\text{K\hspace{0.17em}km\hspace{0.17em}s}}^{-1}\right)}^{-1}$ and $X_{I_d/T_b}^{\text{CO}}=\left(0.0685\pm 0.0019\right)\left({\text{MJy\hspace{0.17em}sr}}^{-1}\right){\left({\text{K\hspace{0.17em}km\hspace{0.17em}s}}^{-1}\right)}^{-1}$, respectively. The shape of the joint distribution (lower-left panel) indicates that $X_{I_d/T_b}^{\text{HI}}$ and $X_{I_d/T_b}^{\text{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 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_{\text{d}}^{\text{HI}}$ and $I_{\text{d}}^{\text{CO}}$ as well as the best-fit map of the reconstructed dust intensity, $I_{\text{d}}^{\text{HI}}+I_{\text{d}}^{\text{CO}}$. The latter is displayed in the top-middle panel of Fig. 4, where it can be compared to the observational Planck map, $I_{\text{d}}^{\text{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, $\left[I_{\text{d}}^{\text{obs}}-\left(I_{\text{d}}^{\text{HI}}+I_{\text{d}}^{\text{CO}}\right)\right]$, 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 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_{\text{d}}/T_{\text{b}}}^{\text{HI}}$ and $X_{I_{\text{d}}/T_{\text{b}}}^{\text{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_{\text{d}}/N_{\text{b}}}^{\text{HI}}$, $X_{N_{\text{H}}/T_{\text{b}}}^{\text{HI}}$ , $X_{I_{\text{d}}/N_{\text{H}}}^{\text{12CO}}$, and $X_{N_{\text{H}}/T_{\text{b}}}^{\text{12CO}}$ defined in Sect. 2.2.

Regarding H_I, the value of the conversion factor from velocity-integrated brightness temperature to hydrogen column density is generally taken as $X_{N_{\text{H}}/T_{\text{b}}}^{\text{HI}}=0.0182\left({10}^{20}{\text{cm}}^{-2}\right){\left({\text{K\hspace{0.17em}km\hspace{0.17em}s}}^{-1}\right)}^{-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_{\text{d}}/N_{\text{H}}}^{\text{HI}}=\left(0.0484\pm 0.0010\right)\left({\text{MJy\hspace{0.17em}sr}}^{-1}\right){\left({10}^{20}{\text{cm}}^{-2}\right)}^{-1}$ in their reference region defined by |b| > 20° and N_HI < 1.2×10²¹ cm⁻², while Planck Collaboration (2011a) derived values in the range [0.021, 0.071] (MJy sr⁻¹) (10²⁰ cm⁻²)⁻¹ for local (low-velocity) clouds toward 14 high-latitude fields covering ≃825 deg² on the sky. Combining the above values of $X_{N_{\text{H}}/T_{\text{b}}}^{\text{HI}}$ and $X_{I_{\text{d}}/N_{\text{H}}}^{\text{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_{\text{d}}/T_{\text{b}}}^{\text{HI}}=X_{I_{\text{d}}/N_{\text{H}}}^{\text{HI}}{X}_{N_{\text{H}}/T_{\text{b}}}^{\text{HI}}=\left(0.00088\pm 0.000018\right)\left({\text{MJy\hspace{0.17em}sr}}^{-1}\right){\left({\text{K\hspace{0.17em}km\hspace{0.17em}s}}^{-1}\right)}^{-1}$ and [0.00038, 0.00129] (MJy sr⁻¹) (K km s⁻¹)⁻¹, respectively, where the uncertainties include only uncertainties in $X_{I_{\text{d}}/N_{\text{H}}}^{\text{HI}}$, not uncertainties in $X_{N_{\text{H}}/T_{\text{b}}}^{\text{HI}}$ associated with, for instance, opacity saturation. Both estimates are consistent. Our best-fit value, $X_{I_{\text{d}}/T_{\text{b}}}^{\text{HI}}=\left(0.00081\pm 0.000006\right)\left({\text{MJy\hspace{0.17em}sr}}^{-1}\right){\left({\text{K\hspace{0.17em}km\hspace{0.17em}s}}^{-1}\right)}^{-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 (τ₃₅₃/N )^CO H for six nearby molecular clouds between longitude l = 139° and 191° and between latitude b = −3° and −56°. Relying on 353 GHz dust emission data and on [0.4, 100] GeV γ-ray emission data, they obtained values of X_CO in the range [0.86, 1.73] (10²⁰ cm⁻²) (K km s⁻¹)⁻¹ and [0.36, 1.14] (10²⁰ cm⁻²) (K km s⁻¹)⁻¹, respectively. They further obtained values of (τ₃₅₃/N_H)^CO in the range [12.6, 44] (10²⁷ cm⁻²)⁻¹. Noting that $X_{N_{\text{H}}/T_{\text{b}}}^{\text{CO}}=2X_{\text{CO}}$ and I_d = (4.76 × 10⁴ MJy sr⁻¹) τ₃₅₃ for a dust temperature T_d = 19.7 K (Planck Collaboration 2014), we find that the corresponding ranges of $X_{I_{\text{d}}/T_{\text{b}}}^{\text{CO}}=X_{I_{\text{d}}/N_{\text{H}}}^{\text{CO}}{X}_{N_{\text{H}}/T_{\text{b}}}^{\text{CO}}$ are [0 10, 0 73] (MJy sr⁻ ) (K km s⁻¹)⁻¹ and [0.043, 0.48] (MJy sr⁻¹) (K km s⁻¹)⁻¹, 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_{\text{d}}/T_{\text{b}}}^{\text{CO}}=\left(0.0685\pm 0.0019\right)\left({\text{MJy\hspace{0.17em}sr}}^{-1}\right){\left({\text{K\hspace{0.17em}km\hspace{0.17em}s}}^{-1}\right)}^{-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_{\text{d}}/T_{\text{b}}}^{\text{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.

Fig. 2 Top row: observational maps of the velocity-integrated brightness temperatures, Tb, of the Hi 21 cm (left), 12CO 2.6 mm (middle), and 13CO 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, Tb⋆, of the Hi (left) and 12CO (right) lines, where $T_{\text{b}*}^{12\text{CO}}$ is derived based on the combined 12CO and 13CO data. Third row: same maps resampled to the common 26×26 pixel grid (pixel size = 4’). 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−1] for Hi and ≃[0, 2.5 MJy sr−1] for CO.

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_{\text{d}}/T_{\text{b}}}^{\text{HI}}$ and $X_{I_{\text{d}}/T_{\text{b}}}^{\text{CO}}$ [in (MJy sr−1) (K km s−1)−1], for the two gas tracers, Hi and CO, in the G139 region.

Fig. 4 Maps of the intensity, Id, 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.

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×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_{\text{b}*}^{\text{HI}}$ and $T_{\text{b}*}^{\text{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 ¹²CO and ¹³CO cubes.

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

The CO decomposition (bottom panel of Fig. 5) leads to 7 components peaking between ≃ −32 and 0 km s⁻¹. The four dominant components (CO [3] - CO [6]) are clustered around v ∼ −10 km s⁻¹, with peak velocities between ≃−16 and 0 km s⁻¹; 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 ≃ 0 km s⁻¹ as well as two very close components (CO [1] and CO [2]) centered at v ≃ −32 km s⁻¹. The latter probably form a single physical entity, which was artificially split by the sharp transition between our two estimates of $T_{\text{b}*}^{\text{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 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 ≃ −32 km s⁻¹.

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 ¹²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 ≃ [−5, +10] km s⁻¹, the Camelopardalis (Cam) OB1 association layer, with v ≃ [−20, −5] km s⁻¹, and the Perseus arm, with v ≃ [−60, −30] km s⁻¹. 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 ≈ [150, 300] pc, [800, 900] pc, and [2, 3] kpc, respectively. Later, Montillaud et al. (2015) identified three kinematic components, with RVs v ≈ −34 km s⁻¹, −16 km s⁻¹, and −9 km s⁻¹, in the same ¹²CO data from Dame et al. (2001) as those studied by Straižys & Laugalys (2007). The v ≈ − 34 km s⁻¹ component clearly peaks at the position of the bright core in the dust intensity map (top-left panel of Fig. 1), the v ≈ −16 km s⁻¹ component has a more extended and more diffuse emission, and the v ≈ −9 km s⁻¹ component is fainter. Referring to the work of Straižys & Laugalys (2007), they concluded that the v ≈ −34 km s⁻¹ component must be part of the Perseus arm, while the v ≈ −16 km s⁻¹ and v ≈ −9 km s⁻¹ components must belong to the Cam OB1 layer. In addition, their identification of the v ≈ −34 km s⁻¹ component with the bright core led them to place the latter at a distance d = 2.5 ± 0.5 kpc. To make the link with our own study, the v ≈ −34 km s⁻¹, −16 km s⁻¹, and −9 km s⁻¹ components of Montillaud et al. (2015) presumably correspond, in our decomposition, to the CO[1] - CO[2] pair (centered at v ≃ −32 km s⁻¹) and to the two peaks formed by the clustered CO[3] - CO[6] components (at v ≃ −16 and −7 km s⁻¹).

Fig. 5 Spectra of the opacity-corrected brightness temperatures, $T_{\text{b}*}^{\text{HI}}$ and $T_{\text{b}*}^{\text{CO}}$, 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_{\text{b}*}^{\text{HI}}$ and $T_{\text{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).

Fig. 6 Reconstructed maps of the opacity-corrected brightness temperatures, $T_{\text{b}*}^{\text{HI}}$ (top) and $T_{\text{b}*}^{\text{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).

3.3.2 Cloud reconstruction

Moving on to the second step of the procedure described in Sect. 2.3, we examine the $n_{\text{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 𝒞_jj′ (Eq. (29)) at each of the 26 × 26 pixels of our grid, retain the weighted average value of 𝒞_jj′ over all the pixels, ⟨𝒞_jj′⟩ (Eq. (30)), and compare ⟨𝒞_jj′⟩ to the velocity-coherence threshold, 𝒞_th, entering Eq. 31). Here, for the purpose of illustration, we adopt 𝒞_th = 50%. The results obtained for this value of 𝒞_th are presented in Sects. 3.3–3.4, while an overview of the results obtained for all the integer values of 𝒞_th between 0 and 100% is provided in Sect. 3.5. When ⟨𝒞_jj′ ⟩ ≥ 𝒞_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 𝒞_jj′ displayed in the upper-right half and the derived values of ⟨𝒞_jj′⟩ indicated in the lower-left half.

It emerges from Fig. 7 that 83 pairs satisfy the condition ⟨ 𝒞_jj′ ⟩ ≥ 𝒞_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 ⟨𝒞_jj′⟩ increases), as opposed to light blue for the other pairs. Amongst the pairs with ⟨𝒞_jj′⟩ ≥ 𝒞_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 ℂ₁, ℂ ₂, ℂ ₃, ℂ ₄, ℂ ₅, ℂ ₆, and ℂ ₇, 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.

ℂ₁, ℂ ₂, and ℂ₃ are three purely atomic clouds, containing 1, 2, and 1 Hi components, respectively. ℂ₁ is faint, with an enhancement in the southwest corner; ℂ₂ is much brighter, especially in the northwest part; ℂ₃ is faint, knotty, and mostly confined to the northeast corner. ℂ₁ and ℂ₂ both cover a rather broad velocity range, which could perhaps indicate that they are quite extended along the LoS.

ℂ₄ and ℂ₅ are two mixed atomic-molecular clouds. ℂ₄ encloses an oblique, elongated CO structure (2 CO components), partially surrounded by an extended Hi envelope (2 Hi components). ℂ₅ 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.

ℂ₆ and ℂ₇ are two purely molecular clouds, with two very close CO components for ℂ₆ and a single CO component for ℂ₇. ℂ₆ 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). ℂ₇ is much fainter, confined to the northeast corner, and also localized in velocity.

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

Fig. 7 Correlation coefficient, 𝒞jj′ (Eq. (29)), expressed in percent, for the 325 possible pairs of Gaussian kinematic components j and j′ (j ≠ j′) extracted with ROHSA 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 𝒞jj′ are displayed in the upper-right half of the figure. The weighted average values of 𝒞jj′ over the 26 × 26 pixels, ⟨𝒞jj′ ⟩ (Eq. (30)), are indicated in the lower-left half, with red or light-blue shading according to whether ⟨𝒞jj′ ⟩ ≥ 50% or ⟨𝒞jj′ ⟩ < 50% and with a level of shading increasing with increasing ⟨𝒞jj′⟩. Pairs for which ⟨𝒞jj′⟩ ≥ 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, ℂ1, …, ℂ 7, as indicated along the diagonal.

Fig. 8 Reconstructed spectra (left) and maps (right) of the 353 GHz dust intensities, Id,i, of the seven clouds, ℂ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).

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, ℂ₁, …, ℂ₇. For each cloud ℂ_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, ${\overline{p}}_{\text{d},i}$, and polarization angle, ${\overline{\psi }}_{\text{d},i}$. Following the procedure described in Sect. 2.4, we derive the best-fit values of ${\overline{p}}_{\text{d},i}$ and ${\overline{\psi }}_{\text{d},i}$, together with their uncertainties, $\sigma \left({\overline{p}}_{\text{d},i}\right)$ and $\sigma \left({\overline{\psi }}_{\text{d},i}\right)$, by minimizing ${\chi }_{\text{r}}^2$ in Eq. (40) through MCMC simulations. As a prior, we require ${\overline{p}}_{\text{d},i}$, as suggested by Planck Collaboration (2015a) (see Appendix B).

The best fit has χ_r = 3.17. This fairly large value indicates that our model composed of seven clouds with uniform polarization parameters (${\overline{p}}_{\text{d},i}$ and ${\overline{\psi }}_{\text{d},i}$) does not reproduce the Planck polarization maps within the total “observational uncertainties”, σ_Q and σ_U, defined below Eq. (40). These uncertainties are dominated by “decomposition errors” in $Q_{\text{d}}^{\mathrm{mod}}$ and $U_{\text{d}}^{\mathrm{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_{\text{d}}^{\text{obs}}$ and $U_{\text{d}}^{\text{obs}}$ (first term) by a factor ∼2–3. The main reason for the fairly large value of χ_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, ${\overline{p}}_{\text{d},i}$ and ${\overline{\psi }}_{\text{d},i}$, together with their standard deviations, are listed in Table 2. It appears that the polarization fractions vary widely, from ${\overline{p}}_{\text{d},4}\simeq 2\%$ to ${\overline{p}}_{\text{d},1}$, ${\overline{p}}_{\text{d},7}$ hitting our imposed upper limit of 23%. Clearly, ℂ₇ 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 ℂ₁ 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 $\left|{\overline{\psi }}_{\text{d},1}\right|$, $\left|{\overline{\psi }}_{\text{d},5}\right|\underset{˜}{<}5°$ to $\left|{\overline{\psi }}_{\text{d},3}\right|\simeq 90°$. 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 ${\overline{p}}_{\text{d},i}$ and ${\overline{\psi }}_{\text{d},i}$ (third and fourth columns of Table 2) to compute the Stokes parameters Q_d,i and U_d,i of each cloud ℂ_i (Eqs. (34) and (35)) as well as the resulting Stokes parameters, $Q_{\text{d}}^{\mathrm{mod}}$ and $U_{\text{d}}^{\mathrm{mod}}$ (Eqs. (25) and (26)), and the associated polarized intensity, $P_{\text{d}}^{\mathrm{mod}}$ (Eq. (5)), LoS-averaged polarization fraction, $p_{\text{d}}^{\text{los},\mathrm{mod}}$ (Eq. (8)), and LoS-averaged polarization angle, ${\psi }_{\text{d}}^{\text{los},\mathrm{mod}}$ (Eq. (9)). The best-fit maps of $Q_{\text{d}}^{\mathrm{mod}}$, $U_{\text{d}}^{\mathrm{mod}}$, $P_{\text{d}}^{\mathrm{mod}}$, $p_{\text{d}}^{los,\mathrm{mod}}$ and ${\psi }_{\text{d}}^{los,\mathrm{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_{\text{d}}^{los,\text{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_{\text{d}}^{\text{obs}}$ (left panel in the third row of Fig. 9) is globally more uniform than the map of $I_{\text{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_{\text{d}}^{\text{obs}}$ (left panel of Fig. 4) and $p_{\text{d}}^{\text{los,obs}}$ (left panel in the fourth row of Fig. 9). The map of $P_{\text{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_{\text{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_{\text{d}}^{\text{mod}}$ and I_d,5 (fifth panel in the right column of Fig. 8) suggest that the bright and strongly polarized $\left({\overline{p}}_{\text{d},5}\simeq 18\%\right)$ cloud ℂ₅ provides a dominant contribution to $P_{\text{d}}^{\text{mod}}$.

The map of $p_{\text{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 ℂ₆, which is very bright and has a polarization angle very different from those of the other bright clouds, ℂ₅ and ℂ₂, detected along its LoS (see fourth column of Table 2). As a result, ℂ₆ contributes constructively to $I_{\text{d}}^{\text{mod}}$, leading to a bright spot in the map of $I_{\text{d}}^{\text{mod}}$ (Fig. 4), and contributes destructively to $P_{\text{d}}^{\text{mod}}$, leading to a dip in the map of $P_{\text{d}}^{\text{mod}}$ and a low-polarization spot in the map of $p_{\text{d}}^{\text{los,mod}}$ (Fig. 9).

The map of ${\psi }_{\text{d}}^{\text{los,obs}}$ (left panel in the fifth row of Fig. 9) is fairly homogeneous, with $\left|{\psi }_{\text{d}}^{\text{los,obs}}\right|\underset{˜}{<}10°$ everywhere. The four large-${\psi }_{\text{d}}^{\text{los,obs}}$ spots (orange regions) near the middle are associated with low $p_{\text{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_{\text{d}}^{\text{los,obs}}$ is not associated with a particularly large ${\psi }_{\text{d}}^{\text{los,obs}}$ (it actually sits between two large-${\psi }_{\text{d}}^{\text{los,obs}}$ spots) is probably because the putative depolarizing cloud has a polarization angle relatively close to 90^°5. The model map of ${\psi }_{\text{d}}^{\text{los,mod}}$ (middle panel) is even more homo geneous than its observational counterpart, except toward ℂ₆, where ${\psi }_{\text{d}}^{\text{los,mod}}$ reaches ≃16°. The homogeneous, low-$\left|{\psi }_{\text{d}}^{\text{los,mod}}\right|$ background can be mostly attributed to the widespread, bright, and strongly polarized cloud ℂ₅, which dominates the polarized emission and has a small polarization angle $\left({\overline{\psi }}_{\text{d},5}\simeq -3°\right)$. The jump in ${\psi }_{\text{d}}^{\text{los,mod}}$ toward ℂ₆ can be explained by this very bright and moderately polarized cloud having a much larger polarization angle $\left({\overline{\psi }}_{\text{d},6}\simeq 66°\right)$. The discrepancies between ${\psi }_{\text{d}}^{\text{los,obs}}$ and ${\psi }_{\text{d}}^{\text{los,mod}}$ are mostly due to our assumption that each cloud has a uniform polarization angle. In particular, the behavior of ${\psi }_{\text{d}}^{\text{los,obs}}$in the vicinity of ℂ₆ suggests that the true polarization angle of ℂ₆ has values closer to 90° in the bright center of the cloud (hence a small impact on ${\psi }_{\text{d}}^{\text{los,obs}}$) and values closer to 45° on either side (hence a more significant increase of ${\psi }_{\text{d}}^{\text{los,obs}}$).