© V. Pelgrims et al. 2021

1 Introduction

From a cosmological perspective, the characterization of the polarized diffuse Galactic emission from synchrotron and thermal dust is of prime importance. This emission dominates the signal in the frequency range of interest for the observation of the cosmic microwave background (CMB) polarization anisotropies (e.g., see Planck Collaboration X 2016).

In particular, possible contamination from the polarized diffuse Galactic emission has been shown to be one of the major limitations for the detection of primordial CMB polarization B-modes related to the inflationary era in the early universe (BICEP2/Keck & Planck Collaborations 2015; BICEP2 & Keck Array Collaborations 2016). Therefore, providing an accurate characterization and modeling of this polarized Galactic emission is of great importance as it would strengthen the level of confidence with regard to our understanding of its physical scope and would allow for accurate testing of the results obtained from elaborated component-separation techniques (e.g., Planck Collaboration IX 2016; Planck Collaboration X 2016) which are used to extract the cosmological CMB signal.

At microwave frequencies that are typically relevant for CMB experiments, the polarized sky is dominated by synchrotron emission below ~80 GHz and by thermal dust emission above that frequency. Both emission components result from a line-of-sight integration of local emission and these offer the possibility to infer the properties of the Galactic magnetic field (GMF), an important constituent and actor in the ecosystem of our Galaxy. The diffuse polarized Galactic synchrotron emission is produced by relativistic electrons that spiral along the GMF lines (see Ginzburg & Syrovatskiĭ 1966, for a review). Equivalently, the polarized thermal dust emission is produced by rotating aspherical dust grains that are totally or partially aligned with the GMF lines (Davis & Greenstein 1951; King & Harwit 1973; Onaka 1996; Lazarian et al. 1996; Onaka 2000; Efroimsky 2002; Jordan & Weingartner 2009; Vaillancourt et al. 2013; Andersson et al. 2015; Hoang 2017).

Over the last two decades, the polarized diffuse Galactic emission has been measured up to a relatively high level of accuracy and high angular resolution by the WMAP¹ and Planck² satellite experiments (Page et al. 2007; Bennett et al. 2013; Planck Collaboration X 2016), which have observed the sky in polarization in a large range of frequency going from 23 to 353 GHz. The wealth of information present in these full-sky observations is extremely valuable in modeling the various components of the Galaxy. In light of these polarization data, there have been attempts to constrain the geometry of the GMF at the largest Galactic scales and the relative amplitude of their main components (regular, turbulent,...) and of the Galactic matter content (Page et al. 2007; Ruiz-Granados et al. 2010; Jaffe et al. 2010, 2013; Fauvet et al. 2011, 2012; Jansson & Farrar 2012a,b).

Because full-sky polarization data first came at low frequencies, these investigations were driven by the study of the full-sky synchrotron emission. Page et al. (2007) used the three-year full-sky maps from the WMAP satellite at 22 GHz (the K-band) and fitted a parametric model using the polarization position angles of the emission. Ruiz-Granados et al. (2010) used the five-year WMAP polarization data at the same frequency and searched for the best fits of several parametric models on a grid-based exploration of the parameter space, thus providing the first systematic comparison of GMF models. Sun et al. (2008), Sun & Reich (2010), Jansson & Farrar (2012a,b) and Jaffe et al. (2010, 2013) built more sophisticated GMF models and used the same WMAP data to constrain them, complementing (for the most part) the synchrotron data with Faraday rotation or dispersion measures on Galactic or extragalactic sources. Recently, Planck Collaboration Int. XLII (2016) used the synchrotron data from the Planck satellite to update GMF models that were previously constrained from WMAP data and rotation measure data. The latter study has shown the limitations of the models at reproducing the current data sets. See also (Steininger et al. 2018) for a recent work.

Similar studies have been carried out using the polarized thermal dust emission. Page et al. (2007) showed that the 94 GHz band of the WMAP satellite measured the thermal dust emission and used it to constrain the GMF. Fauvet et al. (2011) used the 353-GHz data from the ARCHEOPS balloon experiment (Benoît & ARCHEOPS Collaboration 2004) in addition to the WMAP satellite 22-GHz channel data (tracing the polarized diffuse Galactic synchrotron emission) to constrain GMF models. Jaffe et al. (2013) used the full-sky WMAP 94-GHz polarization maps and showed that the diffuse Galactic emission observed at this frequency is not compatible with GMF configuration that fits best the polarized synchrotron emission as traced by the WMAP low frequency data. More recently, Planck Collaboration Int. XLII (2016) showed that the full-sky 353-GHz data from Planck is in conflict with predictions from the reconstructed large-scale GMF obtained from the synchrotron emission data. Other studies have considered the thermal dust emission only from the polar caps to adjust models of the magnetic field in the neighborhood of the Sun (Planck Collaboration Int. XLIV 2016; Alves et al. 2018; Pelgrims et al. 2020).

From these different works it has been established that the GMF is made of a regular (or coherent) component of an amplitude of a few micro Gauss and to which one or two turbulent components, random and ordered-random, are added (Jansson & Farrar 2012a,b; Jaffe et al. 2013; Planck Collaboration Int. XLII 2016). However, GMF models turn to be complex and highly dimensional (large number of parameters). As a consequence, they are highly degenerated and difficult to constrain, particularly as a result of the fact that according to some authors, the turbulent components (random or ordered-random) have amplitudes that might exceed that of the regular part and also because the density distributions of the Galactic matter content are poorly known.

In principle, a combined study of the diffuse synchrotron and thermal dust emission should allow us to constrain the GMF characteristics in three dimensions and to mitigate the loss of information introduced by the integration along the lines of sight that induces degeneracies. However, an additional complexity occurs. Three-dimensional (3D) density distributions of relativistic electron and of dust grains differ. Therefore, the GMF is not sampled in the exact same way along the lines of sight by the two species of the interstellar matter and the respective emission components may often reflect properties of the magnetized ISM at different locations. This characteristic, which could be used to our advantage to model the GMF, may likely be source of confusion and of apparent discrepancies if it is not properly accounted for.

To face this complexity, one approach would be to constrain models of the GMF and of the matter density distributions through a joint modeling of both emission components. Such an approach, followed by the IMAGINE Consortium (Boulanger et al. 2018), requires us to simultaneously constrain several models, with each residing in potentially high-dimensional space. Performing such an ambitious fit through, for instance, a Bayesian method should ideally be the aim in this instance.

In this paper we motivate and explore a somewhat different approach to infer the main characteristics of the large-scale GMF. In Sect. 2, we review the modeling of the polarized diffuse Galactic emission from the 3D models of the distribution of matter in the Galaxy and of the GMF. In Sect. 3, we discuss the main methodology used as well as the data and simulations used. Section 4 discusses the fitting procedure in intensity and polarization for the dust thermal and its application in simulating the data. In Sect. 5, we discuss the quality of our reconstruction of the regular part of the large-scale GMF, along with possible limitations and biases. Section 7 evaluates the impact of those limitations when independently tackling Galactic synchrotron emission. We apply our methodology to the full-sky polarization 353 GHz Planck maps in Sect. 8. Finally, we summarize our findings and present our conclusions in Sect. 9.

Fig. 1 Schematic view of the polarization direction of the Galactic synchrotron and dust thermal emission as a function of the Galactic magnetic field (GMF). For completeness, we also show the direction of the starlight polarization.

2 Modeling & implementation

2.1 Thermal dust polarized emission

The thermal dust emission dominates the polarized signal above 80 GHz and is the most significant foreground of the CMB at those frequencies (e.g., Planck Collaboration XXX 2014). If the GMF is coherent in a given region of the sky, dielectric aspherical dust grains tend to align with their major axis perpendicular to the field lines in the region and rotate with their spin axis parallel to the field lines (e.g., Martin 2007; Fauvet et al. 2011). As a consequence, the thermal dust emission arising at sub-millimeter and millimeter wavelengths is expected to be linearly polarized perpendicular to the GMF lines, as sketched in Fig. 1. Tracing the same dust, the optical polarized emission of stars is expected to be parallel to the GMF lines as this is due to anisotropic absorption in the plane of polarization.

To model the diffuse thermal emission from optically thin Galactic dust, we adopted a parameterization close to that of Fauvet et al. (2011), but which follows the physically motivated modeling of the emission introduced by Lee & Draine (1985) and reviewed by Pelgrims et al. (2020). We assume that the dust emission comes from a single population of dust grains heated at the same temperature from the interstellar radiation field, which implies a constant dust emissivity and a spatially constant intrinsic degree of polarization. We make the additional assumption that the degree of misalignment of the dust grains with respect to the GMF lines is spatially uniform.

Specifically, we model the intensity and the linear polarization Stokes parameters as $$\begin{array}{ll}{I}_{\text{d}}(n)=\hfill & {ϵ}_{\nu }^{\text{d}}{\displaystyle {\int }_0^{+\infty }\text{d}r{n}_{\text{d}}(r,n)\left\{1+p^{\text{d}}{f}_{\text{ma}}\left(\frac{2}{3}-{\sin }^2\alpha (r,n)\right)\right\}}\hfill \\ Q_{\text{d}}(n)=\hfill & {ϵ}_{\nu }^{\text{d}}{p}^{\text{d}}{f}_{\text{ma}}{\displaystyle {\int }_0^{+\infty }\text{d}r{n}_{\text{d}}(r,n){\sin }^2\alpha (r,n)\cos [2\gamma (r,n)]}\hfill \\ U_{\text{d}}(n)=\hfill & {ϵ}_{\nu }^{\text{d}}{p}^{\text{d}}{f}_{\text{ma}}{\displaystyle {\int }_0^{+\infty }\text{d}r{n}_{\text{d}}(r,n){\sin }^2\alpha (r,n)\sin [2\gamma (r,n)]},\hfill \end{array}$$(1)

where r is the radial distance from the observer along the line-of-sight at sky position, n. The different terms in the equation are:

• ${ϵ}_{\nu }^{\text{d}}$, the dust emissivity at observational frequency ν, which is linked to the dust temperature through a grey-body’s law (e.g., Planck Collaboration Int. XX 2015, Appendix B).

• p^d, the so-called intrinsic degree of linear polarization of the dust that depends on the properties of the dust grains. It represents the maximum value of the degree of linear polarization of the radiation emitted by an hypothetical ensemble of perfectly aligned dust grains from a small volume; it only depends on the geometry of the dust grains.

• f_ma, the misalignment term which generally should depend on the dust population. It characterizes the average tendency of the dust grains to align with the magnetic field line.

• n_d (r, n), the 3D Galactic dust grain density.

• α(r, n), the inclination angle of the GMF line with the line of sight at (r, n).

• γ(r, n), the so-called local polarization angle.

The local polarization angle is defined in the plane orthogonal to the line of sight as the angle between the polarization vector direction and the local meridian. Expressed in terms of the vector component of the ambient GMF, this angle is given as: $$\gamma (r,n)=\frac{1}{2}\arctan \left(\frac{-2B_{\theta }(r,n)B_{\varphi }(r,n)}{B_{\varphi }{(r,n)}^2-B_{\theta }{(r,n)}^2}\right),$$(2)

with B_θ and B_ϕ the local transverse components of the magnetic field in the local spherical coordinate basis (e_r, e_θ, e_ϕ) with e_θ pointing towards the South pole. Thus, Eq. (2) gives the polarization position angle of the polarization vectors stemming from the small space volume in the HEALPix (or COSMO) convention (Górski et al. 2005), which differs from the more commonly used IAU one. This angle is defined in the range $\left[0,180\right]$ degrees. We highlight that none of the values of I_d, Q_d and U_d depend on the amplitude of the magnetic field, but only on its geometrical structure through the angles α and γ.

Finally, in the following, we assume that the second term in the parenthesis of the top equation of Eq. (1) is negligible compared to the first term: $$I_{\text{d}}(n)={ϵ}_{\nu }^{\text{d}}{\displaystyle {\int }_0^{+\infty }\text{d}}r{n}_{\text{d}}(r,n).$$(3)

This simplifying assumption has been made in all previous works in the field (Fauvet et al. 2011; Jaffe et al. 2013; Planck Collaboration Int. XLII 2016). Relying on simulations, we found that in pixel space, the relative difference is at most of 10 percent (see also Sect. 3).

2.2 Synchrotron polarized emission

The diffuse Galactic synchrotron emission is produced by relativistic electrons that spiral around the GMF lines (Ginzburg & Syrovatskiĭ 1966; Rybicki & Lightman 1979). This emission is to be polarized perpendicularly to the GMF lines as sketched in Fig. 1. For the contribution of synchrotron emission to CMB frequencies (≳20 GHz) we follow the modeling of the emission presented by Page et al. (2007) and Fauvet et al. (2011) that relies on Rybicki & Lightman (1979). We adopt the notation of Fauvet et al. (2011). Assuming a power law for the relativistic electron energy the linear polarization Stokes parameters for the diffuse Galactic synchrotron emission is given as: $$\begin{array}{ll}{I}_{\text{s}}(n)=\hfill & {ϵ}_{\nu }^{\text{s}}{\displaystyle {\int }_0^{+\infty }\text{d}}r{n}_{\text{e}}(r,n){\left(B_{\perp }{(r,n)}^2\right)}^{(s+1)/4}\hfill \\ Q_{\text{s}}(n)=\hfill & {ϵ}_{\nu }^{\text{s}}{p}_{\text{s}}{\displaystyle {\int }_0^{+\infty }\text{d}}r{n}_{\text{e}}(r,n){\left(B_{\perp }{(r,n)}^2\right)}^{(s+1)/4}\cos [2\gamma (r,n)]\hfill \\ U_{\text{s}}(n)=\hfill & {ϵ}_{\nu }^{\text{s}}{p}_{\text{s}}{\displaystyle {\int }_0^{+\infty }\text{d}}r{n}_{\text{e}}(r,n){\left(B_{\perp }{(r,n)}^2\right)}^{(s+1)/4}\sin [2\gamma (r,n)],\hfill \end{array}$$(4)

where ${ϵ}_{\nu }^{\text{s}}$ is the synchrotron emissivity, n_e(r, n) is the local density of relativistic electrons, p_s is the intrinsic synchrotron polarization fraction which is related to the relativistic electron energy spectral index (s) as follows: $$p_{\text{s}}=\frac{s+1}{s+7/3}.$$(5)

The angle γ found in the expression of Q_s and U_s is the same as in Eq. (2). It marks the orientation perpendicular to the plane-of-sky component of the magnetic field (B_⊥). We note that $B_{\perp }{(r,n)}^2=B{(r,n)}^2{\sin }^2[\alpha (r,n)]$, where the angle α(r, n) is the same as in Eq. (1); namely, it is the inclination angle of the GMF vector with respect to the line of sight.

2.3 Implementation

To simulate the Stokes parameters I, Q, and U of the polarized diffuse thermal dust and synchrotron Galactic emission, we follow Eqs. (1) and (4). We chose to sample the Galactic space according to a spherical coordinate system centered on the Sun. The radial distance to the observer is linearly sampled and we adopt the equal-area HEALPix tessellation (Górski et al. 2005) to uniformly sample the angular coordinates as this is the most commonly used format when dealing with CMB data. We thus consider as many spherical shells as radial bins. Then, the matter density distribution and the GMF models are evaluated at each 3D location along the lines of sight.

In this work, we consider parametric models both for the matter density distribution and for the regular and stochastic (turbulent) GMF components. We assume that these quantities are constant within the elemental volume and we integrate them according tothe midpoint rule. Our choice of the sampling of the Galactic space and the overall implementation are thus similar to the one adopted in the hammurabi code implementation (Waelkens et al. 2009), except that we do not consider refining the angular gridas the radial distance increases. At this stage, we do not consider this simplification as a substantial issue.

Within the framework of the RADIOFOREGROUNDS project, we have developed gpempy, a self-consistent suite a Python codes that follows the above implementations. It is publicly available³ where we describe the architecture that is optimized for quick and user-friendly simulations of the polarized sky. These codes are not optimized for an MCMC exploration of parameter spaces of models.

2.4 Parametric models

In this paper, we consider gentle parametric models for the matter density distribution and for the large-scale regular GMF. The models are reviewed in details in Appendix A.

3 Methodology, data and simulations

3.1 Methodology

To efficiently constrain, from the existing data, the large-scale GMF features relying on parametric models, it would be best to search for a combination of observables that allows us to tackle the lowest levels of complexity at a time. Based on the modeling of the emission presented in Eqs. (1) and (4), we observe that the Stokes parameters of the linear polarization result from a non-trivial mixing of the matter density distribution and the geometrical structure of the GMF. However, we also notice that the thermal dust intensity is to first-order independent from the GMF (see Eq. (3)). Furthermore, the dust polarization emission is only affected by the geometry of the GMF and not by its intensity. These two facts open up the possibility to constrain the dust density distribution separately from the GMF and to considerably reduce the number of free parameter in the fitting procedure. Ideally armed with a good model for the dust density distribution, we can use it to constrain GMF models using polarization data. Here, we consider using the intensity normalized (or reduced) Stokes parameters q_d = Q_d∕I_d and u_d = U_d∕I_d, rather than the Stokes Q_d and U_d. This choice is motivated by the fact that the reduced Stokes parameters q_d and u_d may be regarded as the intensity-weighted mean of the GMF geometry. Therefore, we argue that the integrated GMF geometry dominates these quantities more than the dust density. For the same reason, we expect that possible biases or mis-modeling of the dust densitydistribution have less of an effect on the final reconstruction of the GMF geometry (see also our Sect. 5.2).

In the case of the synchrotron emission, it is more difficult to find an approach that separates the matter density distribution from the GMF. Furthermore, unlike the case of thermal dust emission, it is risky to consider thereduced Stokes parameters at low frequency as other Galactic emission components (e.g., from anomalous microwave emission and free-free) are also important in intensity at CMB frequencies. Thus, a careful separation of these Galactic components is needed to safely compute the reduced synchrotron Stokes parameters q_s = Q_s∕I_s and u_s = U_s∕I_s. In contrast, assuming that a good model of the geometrical structure of the GMF can be obtained from thermal dust emission, this model couldthen be used to constrain the relativistic electron density distribution and of the GMF strength through a fit of the synchrotron data.

3.2 Data

In this paper, we mainly concentrate on the diffuse polarized thermal dust emission. We use the Planck single-frequency intensity and polarization maps at 353 GHz from the second release (PR2)⁴ that are available on the Planck Legacy Archive⁵. We refer the reader to Planck Collaboration Int. XIX (2015); Planck Collaboration Int. XXI (2015) for details and discussions regarding these data. The Planck HFI 353 GHz maps have a native resolution⁶ of about 4.94′ (Planck Collaboration Int. XIX 2015) and are given in a HEALPix⁷ grid tessellation corresponding to N_side = 2048 (Górski et al. 2005). At the instrument resolution, the 353-GHz polarization Stokes Q and U maps are noise-dominated (Planck Collaboration Int. XIX 2015). This is particularly true at high Galactic latitudes. The dispersion arising from CMB polarization anisotropies is much lower than the instrumental noise for Q and U (Planck Collaboration VI 2014). Thus, its impact on our analysis is expected to be negligible. The cosmic infrared background (CIB) is considered to be unpolarized (Planck Collaboration XXX 2014). At this frequency, subdominant contributions are expected in the intensity map either from the CMB, Galactic, and extragalactic point sources, the CIB and the zodiacal light (e.g., Planck Collaboration Int. XXI 2015). We note that we carry out our work at low resolution so that contributions from point sources and the CIB are not expected to be significant. The CMB temperature anisotropies and zodiacal light contributions are expected to be subdominant with respect to the thermal dust emission in intensity. To consolidate the results of our analysis for the intensity part, we also make use of the full-sky dust extinction map derived in Planck Collaboration XIX (2014), specifically, the map of the dust optical depth (τ₃₅₃). This quantity is related to the dust column density integrated along the line of sight. It has been shown to be a good tracer of the dust column density at high Galactic latitudes. Deviation occurs in denser molecular region and, thus, towards regions of the Galactic plane (Planck Collaboration XIX 2014).

3.3 Simulated maps

For the simulated data used in the following, we consider two dust density distribution models: the regular Logarithmic Spiral Arm (LSA) GMF model and the (random) turbulent GMF model, both presented in Appendix A. Using these parametric models and realistic Planck noise maps, we produced different sets of high-resolution maps that simulate the Planck 353-GHz linear polarization Stokes parameter maps. These sets, hereafter called S1, S2, and S2-turb are defined as follows.

First, S1 is built by combining the exponential-disk (ED) dust density distribution model with the LSA GMF model, with no turbulent component. Second, S2 is built by combining the four-spiral-arms dust density model (ARM4) with the LSA GMF model; no turbulent component. And finally, S2-turb is the same as S2 but including different degree of turbulence in the GMF (A_turb = 0, 0.3, 0.9) (see Sect. 6.2).

The maps are shown in Fig. 2 for N_side = 2048. We give the input values of the regular-model parameters in Table 1. An illustration of the input GMF is given in Fig. A.1.

To construct these simulated maps, we adopted a sampling of the Galactic space defined by N_side = 1024 and a radial step of 0.2 kpc⁸. We then upgraded the maps at N_side = 2048 to directly calibrate the Galactic thermal dust emission of our simulated maps to the measured one at 353 GHz by Planck through a linear fit. To add realistic Planck noise, we consider the one derived from the difference between the Stokes parameter maps obtained from the data corresponding to the first and the second half mission of the Planck individual pointings.

In summary, the simulated Planck maps, M, of the polarized thermal dust emission at 353-GHz are given by: $$M_X=A_{\text{cal}}{S}_X^{\text{dust}}+N_X^{\text{Planck}},$$(6)

where X refers either to I_d, Q_d, U_d, N is the noise map, S is computed using gpempy to estimate Eqs. (1) and (3), and A_cal is the global calibration factor that is computed by a linear fit to the Planck data independently for the intensity and the polarization data. The same A_cal is applied to both Q_d and U_d to preserve the polarization position angle. In the remainder of this paper, we drop the subscript _d as we are only dealing with the thermal dust emission, apart from Sect. 7.

Fig. 2 High-resolution simulations of the I, Q, U Stokes parameters (from left to right) obtained for S1 (top) and S2 (bottom). Units are supposedly KCMB. Intensity map is given in colored-logarithmic scale and polarization maps in linear scale. Full-sky maps are in Galactic coordinates, centered on the Galactic center. Galactic longitude decrease towards the right-hand side of the plots.

4 Fitting procedure in intensity and polarization for the Galactic thermal dust emission

In this section, we describe the fitting procedure that we use to infer large-scale GMF features using the dust Galactic thermal dust emission. The fit to the (simulated) data is performed in two steps. First, we fit for the intensity in order to recover the dust density distribution. Second, we use the best-fit model of the dust density distribution and fit for the polarization maps (q, u) to recover the GMF. We rely on MCMC technique to recover both models of the dust density distribution and of the GMF.

4.1 Choice of fitted-map resolution

Although we have optimized our codes, producing simulated maps of the diffuse Galactic emission at the instrument resolution is extremely time consuming. To effectively adjust models to data and efficiently explore the different parameter spaces, a large number of simulations is needed. Inevitably we thus have to work at a lower resolution than the native one of the data. In Appendix B we explore in detail the consequences of this requirement both at the map level and in term of the model-parameter posterior distribution reconstruction. We demonstrate that the necessity to work at lower resolution induces a bias in the parameter space. However, this bias is currently, and for many practical cases, not the dominant source of uncertainty in the kind of modeling tackled in this study. Nonetheless, it is best to ensure that this source of uncertainty is sub-dominant by working at the highest resolution made possible by the computing facilities.

Inside the fitting procedure, we chose to implement the model maps at N_side = 64 and with a radial bin of 0.2 kpc. For the models we consider in this paper, the computation time of a full set of Stokes parameters maps is always on the order of few tenths of a second. The requirement to run a MCMC-based exploration of the parameter spaces are thus fulfilled.

4.2 Likelihood definition

The best-fit parameters are obtained by maximizing a Gaussian likelihood function $\mathcal{L}$, defined as $$\ln \mathcal{L}=-\frac{1}{2}\left(D-\alpha M\right)C_{\text{D}}^{-1}{\left(D-\alpha M\right)}^T\text{}\text{},$$(7)

where D and M represent either the intensity, I, or the concatenated reduced polarization parameters (q, u) for all pixels in the data (or simulated data) and the model maps, respectively; C_D is the covariance matrix associated to the uncertainties in the data and α is an overall normalization factor. It is estimated for each model, and independently for I and (q, u), at each MCMC step via a simple linear fit: $\alpha =({\sum }_i(D_i{M}_i/{\sigma }_i^2))/({\sum }_i(M_i^2/{\sigma }_i^2))$.

4.3 Estimation of the covariance matrix

The choice of the exact form for the noise covariance matrix is not trivial when we think about the properties of the Planck data at 353 GHz. In the case of the intensity we are in a highly signal-dominated case over most of the sky and there is significant intrinsic dispersion in the signal with respect to the low degree of complexity allowed by the models currently used in the literature. With regard to polarization, the situation is slightly more complex as we are mainly in a signal-dominated case close to the Galacticequator, while at high Galactic latitudes we are mainly in a noise-dominated regime. To account for these peculiarities, we developed a hybrid approach considering both statistical instrumental uncertainties and intrinsic signal dispersion.

In the case of the instrumental uncertainties in the Stokes parameters, {I, Q, U}, we consider the block diagonal per-pixel covariance matrix maps released by the Planck collaboration at N_side = 2048. We neglect the off-diagonal terms⁹ and propagate the uncertainties at low resolution according to $${\sigma }_{D_{64};\text{stat},i}^2=\frac{1}{N_{[2048,64]}^2}{\sum }_{j\in i}{\sigma }_{D_{2048};\text{stat},j}^2,$$(8)

where ${\sigma }_{D_{N_{\text{side}}};\text{stat},i}^2$ is the noise variance for a pixel i in the data maps at the resolution given by N_side; and, N_[2048,64] is the number of pixels in a N_side = 2048 map corresponding to a given pixel in a N_side = 64 map. Furthermore, a first estimate of the signal intrinsic dispersion can be obtained using the normalized dispersion of the data at high resolution (N_side = 2048) corresponding to pixel, i, of the low resolution (N_side = 64) map, as: $${\sigma }_{D_{64};\text{disp},i}^2=\frac{1}{N_{[2048,64]}^2}{\sum }_{j\in i}{\left(D_{2048;j}-D_{64;i}\right)}^2,$$(9)

where D = {I, Q, U}¹⁰. Let us note that in the synchrotron literature (e.g. Jansson & Farrar 2012a,b), the signal dispersion is often taken as simply the dispersion of the data in low-resolution pixels. The two estimates of the variance due to signal dispersion differ by a factor of N_[2048,64] = 1024. We adopt the normalized dispersion, also referred to as the standard error on the mean, in order to obtained a reduced χ² on the order of one when the data are fitted with the appropriate models (see Table 2) and when the uncertainties are dominated by instrumental noise. The drawback of this choice, however, is that the reduced χ² values rapidly increases as soon as the models do not fully correspond to the fitted data.

Finally, we can define “pseudo” uncertainties in the low resolution {I, Q, U} maps, whichwould account for the intrinsic signal dispersion and the noise, as: $${\sigma }_{D_{64}}=\max \left({\sigma }_{\text{stat}},{\sigma }_{\text{disp}}\right).$$(10)

These pseudo uncertainties are then propagated to the reduced Stokes parameters, q and u as: $${\sigma }_x^2=\frac{1}{I^2}{\sigma }_X^2+\frac{X^2}{I^4}{\sigma }_I^2,$$(11)

where x = X∕I and X = {Q, U}. And so, we write the covariance matrix for {q, u} as $C_{x{x}^{\prime }}=\text{diag}\left\{{\sigma }_x^2\right\}$.

4.4 Implementation of the MCMC experiment

In order to explore the parameter space, determine the posterior distributions of the parameters, and find their best-fit values, we used the emcee MCMC Python software written by Foreman-Mackey et al. (2013), who implemented the Affine-Invariant sampler proposed by Goodman & Weare (2010). We run the MCMC code for several Markov chains until the convergence criteria proposed by Gelman & Rubin (1992) is fulfilled for all the model parameters with a threshold value of 1.03. We tested for the convergence every 100 MCMC steps. Due to the complex nature of the problem at hand, local minima can be encountered by some of the chains. The reason is that the χ² hyper-surface may exhibit multiple minima that might be sharp in the case of thermal dust emission (mainly for dust intensity). To remedy this problem, we stopped the MCMC after several thousand steps and relaunched the experiment from the location of the minimum χ² and waited for the convergence criteria to be fulfilled to reach well-sampled posterior distributions.

We used 250 Markov walkers for the intensity fits and 100 for the polarization fits. The Markov chains are initialized according to uniform distributions. For the exploration of the parameter space, we considered non-informative priors, adopting top-hat distributions. The explored ranges of values of the free parameters are given in Table 1 for the different fitted models. The results of different adjustments are presented in the next sections for the simulated and true data sets. In the more evolved cases, it was necessary to run the MCMC for several thousands of steps. This implied the generation of millions of magnetized dusty Milky Ways.

We optimized the codes such that a MCMC realization of a map is always a fraction of a second at N_side = 64. We are limited by the efficiency of the basic Python functions to compute the functional forms of the models, such as the hyperbolic cosine. The computing time required for such a converged MCMC fit to be reached is at the day scale running on twelve CPU cores.

4.5 Validation of the MCMC algorithm

We validated our MCMC implementation on simple simulations of the thermal dust emission. We used the same 3D dust density distribution and regular GMF models that those discussed in Sect. 3.3 (specifically S1 and S2). However, for this test, the I, Q, and U maps are simulated and fitted at the same HEALPix N_side = 64 resolution. The results of the fitting procedure are summarized in Fig. 3 which shows the marginalized 1D and 2D posterior probability distributions for the parameters of the regular LSA GMF model. It turns out that our MCMC reconstruction algorithm performs correctly with the input model parameters being well within the 95% C.L. contours. Similar to better results are obtained while fitting for the dust density distribution.

5 Reconstruction accuracy for the dust density and regular GMF component

In this section, we use simulations S1 and S2 to demonstrate that we can retrieve the input GMF model from the “observed” thermal dust emission using the approach described above. In other words, we show that the first step of our global attempt to infer the large-scale GMF features from the Galactic polarized diffuse emission is achievable. Thus, a first fit is performed on the intensity map in order to obtain the best-fit model of the dust density distribution (n_d). Then, we use the best-fit n_d model as an input to constrain the GMF model by fitting the (q, u) polarization maps.

We consider the following cases. First, Case A, we fit the S1 simulations using the ED model for the dust density and the LSA model for the GMF. Second, Case B, we fit the S2 simulations using the ED model for the dust density and the LSA model for the GMF. And finally, Case C, we fit the S2 simulations using the ARM4 model for the dust density and the LSA model for the GMF.

Case A allows us to validate the methodology in a relatively simple model and is thus appropriated to diagnose possible issues. Case B is of interest because it helps us to evaluate the impact of the poor knowledge of the dust density distribution on the GMF reconstruction. This situation is likely to occur when tackling real data given the convoluted nature of the real data set. Case C helps us to evaluate the possible effect of the loss of information due to line-of-sight integration (i) on the modeling of the intensity map and (ii) to evaluate the effect of the propagation of this source of uncertainty on the reconstruction of the 3D GMF.

The best-fit parameter values that we obtained are reported in Table 1, both for the intensity and the polarization fits. The obtained values of thereduced χ² for each case are reported in Table 2.

Fig. 3 Corner plot showing the marginalized 1D and 2D posterior probability distributions for the parameters of the regular LSA GMF model in the case where the simulated data and the fit are performed at Nside = 64, see Sect. 4.5. The vertical and horizontal light-blue lines mark the input parameter values. Angles are given in degrees and the scale height (z0) in kpc.

5.1 Reconstruction of the dust density

The first steps in our fitting procedure is to provide best-fit parameter values for the dust density distribution model from a fit on the intensity map using the MCMC procedure for intensity described in Sect. 4.

The main results for the intensity fit for case A are shown in the first column of Fig. 4. We find that we obtained a good fit to the data with and overall reduced χ² = 1.02  and residuals consistent with noise at high Galactic latitudes. Nevertheless, we also see that the significance of the residuals increases (first towards positive values and suddenly towards negative values) while getting closer to the Galactic equator. This behavior actually takes its origin in the resolution bias discussed in Appendix B. It is indeed close to the Galactic plane that the functional forms of the dust density are allowed to vary quickly. Therefore, it is from those locations (then projected on the sky) that the loss in space-sampling resolution has larger effects.

For cases B and C, the intensity fit results are shown in Fig. 5. The best-fit model and residuals are shown in rows 2 (3) and 4 (5) for case B (C), respectively. For case B, the fitting to the intensity data is very poor with a reduced χ² of 1.9 × 10⁵. This is as expected because of the very small error bars, the large number of data points and principally because the density distribution model cannot reproduce the complexity of the intensity signal injected in the S2 simulation. The uncertainties used to computed the reduced χ² do not account for mis-modeling. The significance of the residuals are large both in and outside the Galactic equatorial region.

For case C, we obtain a good fit to the data with reduced χ² = 1.33. The residuals are consistent with noise at high Galactic latitudes. In the Galactic plane similar comments as for case A hold. The significance of the residuals even shows some structures and it is slightly larger than for case A due to the somewhat larger complexity of the underlying model.

Finally, from the results for case A and C we conclude that, despite the line-of-sight integration, and assuming we have the right model at hand, we are able to retrieve the parameters of a complex dust density distribution with reasonable confidence. In contrast, as shown in case B, if the applied dust density model does not reflect reality, large residuals can be found. Although we cannot directly extrapolate the results of case B, it allows us to anticipate on the kind of issues we would face with real data sets, as discussed in Sect. 8.

Fig. 4 From left to right, top: data from the S1 simulation downgraded to Nside = 64, middle: best-fit maps obtained while assuming the dust distribution follows an ED model and the GMF the LSA model (Case A in the text), and bottom: significance of the residuals. The color scale of the last row ranges from −5 to 5.

5.2 Reconstruction of the GMF

Having found the best-fit models for the intensity map, we use the best-fit values of the parameters to populate the Galaxy with dust densitydistribution and then to constrain the underlying GMF model using MCMC fit on the polarization maps.

In practice, we use the best-fit parameter values obtained above to evaluate the dust density at all locations of the sampled space. Then, using a MCMC sampling, we vary the four parameters of the GMF model (given in Eq. (A.10)) and integrate, along the line-of-sight, the elemental emitting volumes following Eq. (1) to produce the Q and U maps. Thenwe normalize these simulated maps by the best-fit intensity map to obtain the q and u maps. Finally, we proceed to a comparison of the simulated data and the computed model by maximizing the likelihood function described in Sect. 4.2. The best-fit maps are presented in the second and third columns of Fig. 4, second row, for case A, and on Fig. 5, second and third rows, for cases B and C, respectively. The maps of the significance of the residuals are also shown in the last rows of the figures.

From the comparison of the input and best-fit (q, u) maps of case A and case C, it appears that we are able to provide good fit to the data sets. This isconfirmed by the maps of the significance of the residuals, and by the reduced χ² values of Table 2. The best-fit values of the GMF parameters are reported in Table 1 along with their marginalized uncertainties at the 1σ level. It is seen that the best-fit parameter values and the input ones are very close in terms of relative differences even if the latter do not stand within the 1σ confidence interval. Again, the observed shift is most likely due to the resolution bias which turns out to be significant in this case since it is the only important source of uncertainty (see Appendix B). Unlike what we obtained for the intensity fit, we do not observe any particular trend or feature in the residual maps near the Galactic equator. At the map level, the fitting issue from the resolution bias is likely dimmed by the fact that we consider the reduced Stokes parameters.

For case B, we note that the residuals are significantly larger than for the other two cases. We observe a lack in the signal amplitude of the best-fit polarization maps with respect to the input ones. Nevertheless, we see that we are able to capture the correct sky location of the extrema of the reduced Stokes parameters. An inspection of Table 1 shows us that the best-fit values of the GMF parameter are close to the input ones. While the agreement is globally worse here than in case A and C, we see that the agreement is better for the parameters of the spiral shape of the GMF than those of the X-shape (see Appendix A.2).

The reconstructed geometrical structures of the GMF corresponding to the cases A, B, and C best-fit parameters are shown in Fig. 6, and can be compared to the input one shown in Fig. A.1. These figures qualitatively illustrate the fact that we are indeed able to reconstruct (i.e., constrain) the geometry of the GMF from the observation of the polarized thermal dust emission in the three cases discussed above.

Fig. 5 Maps corresponding to the modified Stokes I, q, and u from left to right. Top: data from the S2 simulation downgraded to Nside = 64; second row: best-fit maps obtained while assuming the dust distribution follows an ED model (case B); third row: best-fit maps obtained while assuming the dust distribution follows an ARM4 model (case C). The two last rows show the significance of the residuals for the case B (4th row) and the case C (5th row). For the case B, the color scales range from −500 and 500 and from −15 to 15 for the intensity and the polarization, respectively. For the case C, color scales range from −5 to 5, as in Fig. 4.

5.3 Quality of the reconstructed GMF

In this subsection, we quantitatively compare the GMF geometrical structures that we reconstructed by fitting the polarization maps for case A, B, and C (see Fig. 6) to the input model shown in Fig. A.1. At each location of the Galactic space considered when computing the polarization observables (see Sect. 2), we consider two sets of angles that determine the orientation of the GMF lines. First, in the cylindrical coordinate system centered on the Galactic center (with the z axis perpendicular to the Galactic equator), wedefine the angles p and χ, and second in the spherical coordinate system centered on the Sun, the angles α and γ. In the first scheme, p, the pitch angle, is the angle that makes the GMF line with e_ϕ, and χ, the tilt angle, is the angle that makes the field line with planes parallel to the Galactic plane (z constant); χ is the complementary to the angle made by theGMF line with e_z and characterizes the out-of-plane GMF component. In the second scheme, α is the inclination angle that makes the field line relative to the line of sight and γ the position angle in the plane orthogonal to it. Let ${\widehat{\xi }}_{\text{r}}$ be any of the angles defined above and measured for the input GMF model and $\widehat{\xi }$ the same angle but for the reconstructed GMF models. We define the signed difference angle ${\Delta }_2(\widehat{\xi },{\widehat{\xi }}_{\text{r}})$ as the difference of the two angles that range between −90 and 90 degrees (−90 and 90 corresponding to the same configuration).

To quantify the accuracy of the reconstructed GMF structures, we compute at each location the ${\Delta }_2(\widehat{\xi },{\widehat{\xi }}_{\text{r}})$ for the angles between the field lines of a best-fit GMF model and the field lines of the input GMF model. We generate histograms of the angle differences. An hypothetical perfect reconstruction would lead to a single bin centered in zero. The histograms for the three reconstruction cases are presented in Fig. 7 for the four angles discussed above. Cases A, B, and C are shown on the same plots for comparison.

For case A and case C, the difference between the pitch and the tilt angles measured at all locations of the sampled Galaxy never exceeds one degree. The two reconstructions are very good. For case B, the departure is slightly larger. The distributions for the whole sampled space peak at about two degree and are centered on zero, and can go up to four degrees. We consider the reconstruction of the GMF geometry in case B as being fair given the large difference between the modeled dust density distribution and the recovered one. Such an achievement is possible because we fit the GMF parametric model to the data using reduced Stokes parameters.

We can also compare the best-fit models to the input one via the differences of the polarization position angles in the maps. For this, we use the acute angle, defined between 0 and 90 degree as $${\Delta }_{{\widehat{\psi }}_{\text{r}},\widehat{\psi }}={90}^{°}-|{90}^{°}-|\widehat{\psi }-{\widehat{\psi }}_{\text{r}}\left|\right|,$$(12)

where the consecutive absolute values take into account the axial nature of the polarization vectors, that is, that only their orientation are relevant. In Fig. 8 we show the difference of polarization position angles as defined before. We observe that the reconstruction is excellent for case A and very good for case C, with a distribution peaked at zero and about 95% of the pixels below 0.1 and 4 degrees, respectively. For case B, the results are worse with a distribution also peaked at zero but with a larger scatter (95% of the pixels below 20 degrees).

Fig. 6 Reconstructed GMF structures in the XY plane (top) and XZ plane (bottom) in case A, case B, and case C, respectively, shown from top to bottom. See text for details. Similar plot for the input GMF model is shown in Fig. A.1.

6 Impact of turbulence in the reconstruction of the regular GMF

In the context of a first study, we consider in this paper only the reconstruction of the regular GMF and do not explicitly account for a stochastic component in the likelihood function. Nevertheless, in this section, we study the consequences of ignoring the turbulent component of the GMF in our analysis.

Fig. 7 Histograms of the angle differences (in degree) between the best reconstructions of the GMF and the input one computed at every location of the sampled Galactic space (left). From top to bottom: differences of the pitch angles, the tilt angles, the inclination angles, and the position angles (see text). Case C overlaps case A in tilt and position angle differences. Right: same as for left column but we present the cumulative distributions of the absolute values of the angle differences.

Fig. 8 Histogram of the differences of the polarization position angles between the data and the best fits. The polarization position angles are deduced from the q and u maps both from the downgraded simulations (S1 and S2 but without the Planck noise) and from the best-fit maps.

6.1 Modeling of the turbulent component

We consider an overall turbulent contribution so that the total 3D GMF can be written as (Chandrasekhar & Fermi 1953; Hildebrand et al. 2009): $$B^{\text{tot}}(r)=B^{\text{reg}}(r)+A_{\text{turb}}{B}^{\text{turb}}(r),$$(13)

where the regular (B^reg) and turbulent (B^turb) GMF components are normalized so that A_turb represents the relative amplitude between them (e.g. A_turb = 1 means that contributions to the GMF of the regular and turbulent components are equal). In this paper B^reg is assumed to be one of the regular models depicted in Appendix A.2 and B^turb is assumed to be a Gaussian realization of a random isotropic 3D vector field that follows a two-slopes power-law spectrum as described in detailed in Appendix A.3.

Using polarization data only, Fauvet et al. (2011) obtained an upper limit for the turbulent relative amplitude, A_turb ≲ 0.25, by combining synchrotron and dust polarization data. In more recent studies (e.g., Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XLIV 2016 and references therein) which include also data from pulsars, it is commonly admitted that the turbulent component might be of the same order as the regular one, A_turb ≲ 1.0. Accounting for these results we limit our investigations to the comparison of results for three different cases: A_turb = 0 (regular component only as in Sect. 5.2), 0.3 (comparable to the upper limit of A_turb from Fauvet et al. 2011), and 0.9 (comparable regular and turbulent components as indicated by the more recent studies). In all the cases, the turbulent component will not be considered in the fits described below.

In practice, assuming a model for the regular part of the GMF and a realization of our turbulent model, the regular and turbulent components of the GMF are combined at every sampled positions in the Galaxy according to Eq. (13). Then, the thermal dust emission I, Q and U Stokes parameters are computed by integrating Eq. (1) considering a thermal dust density distribution model and the total GMF. For illustration, we show in Fig. 9 the simulated q maps (downgraded at N_side = 64) of the thermal dust emission as observed by Planck at 353 GHz when considering the turbulent component of the GMF, and, for the S1 simulation regular GMF component. We represent the signal-only (top) and signal plus noise maps (bottom) for A_turb = 0.0 (top), 0.3 (middle), and 0.9 (bottom).

6.2 Quality of the reconstruction of the regular GMF in the presence of turbulence

We rely on the S2-turb (see Sect. 3.3) sets of simulated maps to estimate the uncertainties induced in the reconstruction of the large-scale regular component of the GMF by ignoring the turbulent component. We apply our MCMC procedure, which does not account for the turbulent component, on these three sets of polarization maps for three different cases:

Case I: n_d ≡ ARM4 and GMF ≡ LSA,

Case II: n_d ≡ ED and GMF ≡ LSA,

Case III: n_d ≡ ED and GMF ≡ ASS,

where ASS is a more constrained version of the LSA model for which the opening angle of the spiral is fixed to a constant (see Appendix A.2).

In Fig. 10, we show the 1D and 2D marginalized posterior distribution obtained while fitting the S2-turb simulations following cases I, II, and III. From left to right, we show the results for the three turbulent contributions A_turb = 0.0, 0.3, 0.9.

For case I (first row) without turbulence (A_turb = 0.0, thus corresponding the case C of Sect. 5), the input parameters for the regular GMF model are recovered well inside the 95% CL of the posterior distribution. However, this is not the case when turbulence is included and we observe a bias in the recovered regular GMF model parameters with respect to the input ones. Significant biases are also observed for case II and case III both for the cases with and without turbulence. Interestingly, though, we see that differences in the ψ₀ angle valuesare only of few degrees for all cases. Inspection of Fig. 10 also reveals that the biases in recovered parameter values induced by the non-consideration of the turbulence is roughly of the same order of the one induced by the mismodeling of the density distribution or of the mismodeling of the regular part of the GMF. In this paper, we concentrate in the latter effect and do not attempt to fit for the turbulence component in the likelihood.

In Fig. 11, we present the cumulative histograms of the angle differences (in degrees) measured at every location of the Galactic sampled space between the input regular GMF model in the S2-turb simulations and the reconstructed ones obtained from the Cases I, II, and III (top to bottom). From left to right we show results for the three turbulent contributions A_turb = 0.0, 0.3, 0.9. Uncertainties at the 68% C.L. are overplotted as shaded regions and are computed by sampling the posterior MCMC probability distributions. We observe that the direction of the regular GMF component can be reconstructed to a good precision, better than 10 degrees, in 70–80% of the space even in the presence of a moderate amount of turbulence and for the three cases considered. However, when the fractional amount of turbulence is high, A_turb = 0.9, the reconstruction degrades significantly. Nevertheless, for cases II and III we observe a surprisingly good reconstruction when comparing to case I. This may point to the conclusion that using an incorrect model for the dust density allows to better reconstruct the 3D geometry of the large-scale regular part of the GMF. Such a result requires further investigation and might be an unexpected advantage of reconstructing the regular GMF using the reduced Stokes parameters. However, the fact that the reconstruction in case III is better than in case II is expected. Indeed, the LSA model is a generalization of the ASS model and the LSA model realization in the simulations is not very different from an ASS model. Therefore, the regular GMF model that is adjusted in case III leaves less degrees of freedom for the field lines but constrains the field to be close to the simulated one. We verify this interpretation by fitting another regular GMF model that is intrinsically different from the LSA. In that case we obtained a poor reconstruction of the 3D geometry of the GMF. Nevertheless, we observe that reliable constraints can be obtained for the pitch angle even for the case of comparable regular and turbulent components, which is favored by current studies.

As expected the effect of not considering the turbulent part in the modeling induces a bias in the parameter space. This bias appears to be very significant due to the large number of data points, the relatively small error values and, more importantly, that uncertainties from mismodeling or incomplete modeling are not included in our log-likelihood definition and so, not included in our posterior distributions. This stresses the need for introducing the turbulent component in the modeling. Furthermore, inspection of both the marginalized posterior distributions and the cumulative histograms reveals that the bias in the model-parameter space that is induced by the non-consideration of the turbulence is roughly of the same order of the one induced by the mismodeling of the density distribution or of the mismodeling of the regular part. This shows that improvements in the modelings are mandatory both on the regular and turbulent contributions to the GMF.

Fig. 9 For the S1 simulation, map of the reduced Stokes q downgraded at Nside = 64 for the cases with Aturb = 0.0, Aturb = 0.3 and Aturb = 0.9, from left to right. Top panels results from regular and turbulent field; Planck noise is added in the bottom panels.

7 Proof of concept for future synchrotron analysis

As discussed in Sect. 2.2 for the case of the synchrotron emission, the matter density distribution and the GMF are deeply intertwined. Therefore, they have to be constrained simultaneously along with the synchrotron emission spectral index. Here, relying on simulations, we investigate in a toy model if we can safely use the constraints obtained on the geometry of the GMF from the Galactic thermal dust emission analysis to help the fit of the Galactic synchrotron emission maps.

We assume a relatively simple density distribution for the Galactic relativistic electrons. As the density distribution of the relativistic electron is expected to be smoother than the dust density distribution, for instance, due to the diffusion mechanism, we adopt the same exponential disk model as the one used to produce the S1 thermal dust simulations. We adopted the most used values for the two scale lengths of the model (ρ₀, z₀) = (8.0, 1.0) kpc in the case of synchrotron emission. We then produce, directly at N_side = 64, the intensity and polarization maps of the synchrotron emission at a reference frequency according to Eq. (4). We consider two cases in terms of the GMF model parameters. In the first one, we use the input parameter of the GMF given in Table 1, and in the second, the best-fit parameters from the fit of the thermal dust emission discussed in Sect. 5.2. We use the GMF reconstruction from case B, where the incorrect dust density distribution model was used, which it is more likely to correspond to what we would obtain in the case of real data sets and produces the worst reconstruction of the GMF geometrical structure of Sect. 5.2.

Using gpempy, we computed the intensity and polarization maps corresponding to the Galactic synchrotron emission in the two cases presented above. We show them in Fig. 12. In that figure we also show the position angles of the polarization vectors deduced from the Stokes Q and U parameters. The maps obtained from the two GMF look very similar. In pixel space, the relative difference of the intensity exceeds the ten percent threshold for only five percent of the sky. The relative difference of the polarization degree (not shown) exceeds the ten percent threshold for less than one percent of the sky. The agreement in the polarization position angle is also remarkable as also inferred from Fig. 13 where we show the histogram of the acute angles (Eq. (12)) between the polarization vectors corresponding to the two realizations of the synchrotron sky. It shows that the position angles agree at the 3 degree level almost for 95% of the lines of sight.

From these results, we conclude that the GMF model reconstructed from the Galactic thermal dust emission analysis, even when using an oversimplified model for the dust density distribution, can be used at lower frequencies as an input, that is, a prior, to model the Galactic synchrotron emission and obtain constraints on the population of relativistic electrons. This constitutes a proof of concept of the overall methodology proposed in Sect. 3. Full validation and application of this methodology will be presented elsewhere.

8 Reconstruction of the GMF structure from the Planck 353 GHz data

In this section, we apply the MCMC-fitting procedure described in Sect. 4 to the Planck polarization data at 353-GHz (see Sect. 3.2) to infer the large-scale structure of the GMF. We first rely on three dust density distribution models to fit for the intensity map. Then, we use those best-fit models to constrain four different models of the regular component of the GMF through a fit to the reduced Stokes parameter maps. Both the dust density and GMF distribution models are detailed in Appendix A and encode different degree of complexity whereas sharing geometrical features. In this first attempt, the contribution of the turbulent component of the GMF is not accounted for, as discussed above, and local structures of the ISM are not included.

The fits, as discussed below, are performed at the resolution of N_side = 64. We thus downgrade the I, Q and U 353 GHz Planck maps from N_side = 2048 to N_side = 64 and compute the reduced-Stokes parameters q and u and their uncertainties (see Eq. (11)). These data are shown in the first row of Fig. 16. Maps at N_side = 32 are also used in the case of the intensity.

Fig. 10 Corner plots of the 1D and 2D marginalized posterior distributions of the regular GMF parameter for the three reconstruction cases discussed in Sect. 6.2 applied to the S2-turb simulations. Angles are given in degree and the scale heights in kpc. From left to right: columns correspond to Aturb = 0.0, Aturb = 0.3 and Aturb = 0.9, respectively. From top to bottom for cases I to III defined as: case I: nd ≡ ARM4 and gmf ≡ LSA (both models underlying the data); case II: nd ≡ ED and gmf ≡ LSA; and case III: nd ≡ ED and gmf ≡ ASS. The vertical and horizontal light-blue lines mark input parameter values.

8.1 Dust density reconstruction

8.1.1 Fitting procedure

For the fit of the intensity map, we first fit the data at low resolution, N_side = 32 to explore the full parameter space more quickly and determine regions of the parameter space favored by the data. At this resolution, the Markov chains are initialized according to a uniform distribution for each of the parameters. Then, when a converged solution is reached at that low resolution, we start a new exploration of the parameter space by comparing simulations to data at N_side = 64. In this case, the Markov chains are initialized according to Gaussian centered on the best-fit parameter values obtained at N_side = 32 and with a width of a few percent (typically 10%).

The initially explored ranges of values of the free parameters are given in Table 3 for the different fitted models. We consider the two dust density distribution models used for the simulation data (see Sect. 5), ED and ARM4, plus and an extra one formed from the combination of these two: ARM4⊕ED (see Appendix A).

Fig. 11 Cumulative distributions the difference angles between the regular GMF model in simulations S2-turb and the best-fit reconstructed models using the MCMC likelihood analysis. From top to bottom for the pitch angles, the tilt angles, the inclination angles, and the position angles. The three reconstruction cases corresponds to different models assumed to model the data: case I: nd ≡ ARM4 and gmf ≡ LSA (both models underlying the data); case II: nd ≡ ED and gmf ≡ LSA; and case III: nd ≡ ED and gmf ≡ ASS. From left to right, the columns correspond to Aturb = 0.00, Aturb = 0.30, and Aturb = 0.90. The shaded area corresponds to the 68% CL extracted from the posterior of the probability distribution as sampled from the MCMC analysis.

8.1.2 Best-fit models

The best-fit parameters are reported in Table 3. In Fig. 14, we compare the data and the best-fit intensity maps obtained for the three dust density distribution models fitted to I₃₅₃. We also show χ-maps, which give the statistical significance of the residuals. We observe that the best-fit models cannot account for all the complexity and the richness contained in the data. This is not surprising since our models do not include patterns such as clumps, filaments, or bones that can be directly spotted by eye and that require dedicated extra modeling (see e.g. Zucker et al. 2018). Accounting for such features would significantly increase the number of parameters and will not be easily tractable within a MCMC approach.

In Fig. 15, we show the histograms of the χ_i corresponding to the three best-fits to the intensity map. The distributions are nearly Gaussian with very large spread. This is again due to the fact that the uncertainties in the data (even including the intrinsic dispersion) cannot overcome the residuals to such simplistic parametric models as discussed before. As a result, any small-scale feature turns out to be significant and the values of the reduced χ² are ≫ 1 (see Table 5). This was to be expected from our simulation-based study where we show that such high values of reduced χ² are reached as soon as we do not have the right parametric model at hand (see Sect. 5).

To verify the robustness of our best-fit models of the dust density distribution (n_d) with respectto contamination from other emission components (e.g., point sources, CIB), we also fit for the dust optical depth (τ₃₅₃) map presented inPlanck Collaboration XIX (2014). That observable, deduced from the multi-frequency analysis, is known to trace the integrated dust density distribution at least as well as the intensity map does. The best-fit maps agree globally with those obtained from I₃₅₃ fits. This, in turn, also justifies the assumption around Eq. (3).

In conclusion, none of the best-fit models is satisfying on its own. Nevertheless, as discussed in Sect. 4, it is still possible to constrain the GMF geometry even in the case of any mismodeling of the Galactic dust grain density. Therefore, in the following, we use the recovered dust density distribution models to constrain the geometry of the large-scale, regular, GMF models. In order to propagate the uncertainties due to the assumed n_d model, we find relevant to fit the GMF models with all the three n_d models and to take the dispersion on the parameters of the reconstructed GMF as a conservative estimate of the uncertainties on those parameters.

Fig. 12 Maps of the synchrotron emission. Columns correspond to the maps of I, q, u, and polarization position angles (in the IAU convention). The first row correspond to synchrotron emission when the original GMF model is adopted. For the second row, the adopted GMF corresponds to the best-fit obtained from the analysis of dust maps, case B in Sect. 5. In the third row, we present the relative difference of the intensity maps ((Iin − Î)∕(Iin + Î)). For the polarization angles, we display the acute angles between the polarization vectors.

Fig. 13 Acute angle histogram comparing the polarization vector orientations from the synchrotron emission assuming the original GMF model and the worst one from dust analysis (case B). Angles are in degree.

Fig. 14 Intensity maps. First row: 353-GHz map from Planck downgraded at Nside = 64 and the corresponding map of uncertainties that we use to compute the χ2. Rows two to four correspond to dust density distribution models labeled ED, ARM4, and ARM4⊕ED, respectively.The obtained best-fits are shown in the first column and the statistical significance of the residual, per-pixel, are shown in the second column.

Fig. 15 Histograms of the signed significance of the residuals of the best-fits obtained by adjusting the three dust density distribution models (ED, ARM4, and ARM4⊕ED) to the thermal dust intensity map.

8.2 3D regular GMF reconstruction

8.2.1 Fitting procedure

Similar to what we did to fit the intensity map, we adopted flat priors on the free-parameters of the different models of the GMF, which are presented in Table 4. Details on the parametrization of these GMF models are given in Appendix A.2. We initialized the several hundred walkers according to uniform distributions and let the system evolve until convergence is reached. The parameter space for each class of model is also given in Table 4. Unlike the case of the intensity fit, we work only at the N_side = 64 because our implementation is fast enough and there are only a few free parameters of the GMF models that it can make vary.

8.2.2 Best-fit models

In Fig. 16, we show the maps corresponding to the best-fit of each GMF and the maps corresponding to the significance of the residuals in the case where the n_d model is thesimplest (ED). In that figure, the different degree of complexity of the pattern appearing in the maps between the different best-fit models is simply due to the GMF model. It is quite remarkable how the simple modeling adopted here are capable of roughly reproducing the largest patterns observed in the data. However, an inspection of the residuals maps calls for further developments in the polarization map modeling, especially with regard to tackling low and intermediate Galactic latitudes. The agreement between the models and the data are indeed visually found to be better at high Galactic latitudes. This is further discussed in the next section.

The parameter values of all the best-fits of the different GMF models considered are reported in Table 4 for each of the n_d models that we investigated in Sect. 8.1. In Table 5, we also present the values of the reduced χ² that quantify thegoodness of our fits to the GMF models. In Fig. 17, we show the histograms of the χ_i corresponding to the best adjustments of the four investigated GMF models and for the three n_d models fitted to the I₃₅₃ map. These distributions are nearly Gaussian. Again, their widths are large as the uncertainties in the data can not entirely cope with the residual to such simplistic models. We note, however, that the values of the reduced χ² of the best-fit models are one order of magnitude better than those obtained for the fits of the intensity map.

Furthermore, the obtained best-fits of the GMF are generally not compatible within the uncertainties when moving from one dust densitydistribution to another. The main reason for this is that the data are complex and signal dominated in some regions and that our models are far from being able to account for all the complexity and the richness contained in the maps. Indeed, despite the conservative determination of the errors that we adopt (see Eq. (10)), the posterior distributions of the model parameters are too narrow, as they do not reflect the uncertainties that come from mismodeling. The use of a MCMC technique, however, is fully justified by the fact that the multidimensional χ² surfaces, overall, have convoluted geometries and present a large number of local minima that would compromise the performance of other minimization technique. For now, we consider that the scatter of the parameter values obtained with the different n_d models is more representative of the uncertainties on those best-fit parameters. This again stresses the need for more complex and realistic models of the dust density distribution and of the GMF in order to optimally exploit the currently available data.

To test the robustness of our best-fit models, we also constrain the GMF parameters using the n_d best-fit models obtained from the fits on the dust optical-depth map. The fitted (q, u) maps are obtained simply by replacement of the intensity map by the τ₃₅₃ map, with the corresponding uncertainties. The best-fit polarization maps and the geometrical structure of the GMF are found to be in good agreement with those obtained when the modeling of n_d is from the I₃₅₃ fits. This test makes us confident that the GMF constraints obtained through the fit of the reduced Stokes parameter maps are robust against possible residuals from point sources or other components of the diffuse Galactic emission.

In order to test further the stability of our GMF reconstruction against systematic effects, such as the leakage of the intensity towards the Q and U polarization channels, we apply the Global Generalized Fit correction to the data, as suggested by Planck Collaboration VIII (2016). We proceed to the fit of the GMF model for the least evolved and the most evolved modelings (ED + ASS and ARM4⊕ED + QSS) with and without applying the quoted corrections. The best-fit parameter values that we obtain are in agreement at the ten percent level with the previous values. These two tests make us confident that our GMF reconstruction with the Planck data is not strongly biased by instrumental and map-making systematics and shows that those systematics in Planck polarization maps have a significant contribution to the modeling uncertainty budget whereas not accounted for in the covariance matrix. This conclusion was already reached in Alves et al. (2018) and Pelgrims et al. (2020).

Fig. 16 Polarization maps. First row: 353-GHz maps of the reduced Stokes parameters from Planck downgraded at Nside = 64 and the corresponding map of uncertainties that we use to compute the χ2 (q, σq, u, σu). Rows two to five: GMF models labeled ASS, LSA, BSS and QSS using the best-fit of the ED nd model. The obtained best-fits are shown in the first and third columns and the statistical significance of their residuals, per-pixel, are shown in the second and fourth columns.

Fig. 17 Histograms of the significance of the residuals. x stands for the concatenation of the q and u parameters. Each GMF model is represented by a different color and each panel corresponds to a different underlying nd model: ED, ARM4, and ARM4⊕ED from left to right.

Fig. 18 Functional forms of the pitch (top) and tilt (bottom) angles as a function of the Galactic cylindrical coordinate; ρ and z, respectively.Each line corresponds to a best-fit model. Solid, dashed, dash-dotted, and dotted lines correspond to the ASS, LSA, BSS, and QSS GMF models. Blue, green, and red correspond to the best-fits obtained while assuming the nd model to be ED, ARM4 or ARM4⊕ED, respectively.The pitch angle is constant for the ASS, BSS, and QSS models.

8.2.3 Constraints on the GMF geometry

The radial and height evolution of the pitch and the tilt angles that correspond to the best-fit parameters for all the models considered are shown on Fig. 18. For a given GMF model (same line style but different color in Fig. 18), the agreement among the different choice of n_d model is remarkable. The values of the pitch angles and of the tilt angles agree fairly well through the different modelings. The relative difference between the parameters characterizing those angles range from few percent to about few tens of a percent when comparing very different modelings, such as that of the model having the lowest degree of complexity (ED + ASS) and the one with the largest one (ARM4⊕ED + QSS). This level ofuncertainty is of the same order as the one attributable to the Planck residual systematics. As inferred from Fig. 18, the agreement on the pitch angle is the strongest for a Galactic radius of about 8 kpc, with a mean pitch angle of about 27 degrees and a mean scatter among the twelve computed models of about five degrees, corresponding to a mean relative difference of 14 percent between the reconstructions. However, we also observe that in the case of the LSA model, for which the pitch angle can vary with radius, the evolution from the inner part of the Galaxy to the outskirt is very large. Finally, we also find an inversion of the sign of the tilt angle for the QSS model with respect to the other ones.

In Fig. 19 we show the direction of the large-scale GMF lines in horizontal and vertical planes of the Galaxy, namely the (x, y, z = 0) and the (x, y = 0, z) planes. Stream plots of the field lines are shown for the best fits obtained for each GMF models assuming the n_d to follow the ED model. Very similar plots are obtained with the other n_d models, as we discuss later. In terms of the visuals and given that the functional forms of the models may be different, the similarity among the GMF reconstructions is remarkable. We note that for the LSA model, a second minimum of the χ² is found with a distant solution in the parameter space but that is similar to the solution obtained by Steininger et al. (2018) in fitting Faraday-depth all-sky data and synchrotron data. We present this solution in Appendix D. The reason why our global minimum is not reached by these authors is likely due to the fact that they do not allow the ψ₁ parameter (see Eq. (A.11)) to span through negative values. An interesting feature of our best-fit LSA model is that the field lines appear to wind up on themselves at a cylindrical radius of about 18 kpc.

Fig. 19 Fieldline geometrical structures of best-fit GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy. From top-left to bottom-right: ASS, LSA, BSS, and QSS GMF models. The GMF reconstructed models corresponds to the best-fits obtained while assuming that nd is our best-fit ED model. The black dot at the center of the figure is the Galactic center and the black star is at the Sun position. The color scales offer information about the relative strength of GMF (i.e., the norm of the vector field at each location) and of its direction (red going clockwise, blue going counter-clockwise). In the case of ASS and LSA models, the (constant) norm has been fixed to 2.1.

Fig. 20 Corner plot showing the correlations between the significance of the residuals in I, q, and u. This figure corresponds to the modeling with the ARM4⊕ED nd model and with the QSS GMF model.

8.3 Best-fit maps: comparison and characteristics

Based on inspection of residual maps and on χ² values, none of the twelve best-fit models seems to provide a significantly better fit than the others. Each combination of dust density and GMF models results on different possible features on the best-fit maps. On the maps, the differences between models are more evident at low and intermediate Galactic latitudes, where the 3D models (n_d and GMF) differ more strongly and where the uncertainties are very small. At the contrary, at high Galactic latitudes, at about |b_gal|≥ 60°, every best-fit maps seem to converge towards the same solution. This is well illustrated in Figs. C.1 and C.2, where we present orthographic projections of the maps presented in Figs. 14 and 16, respectively.

In Fig. 20, we show, a corner plot illustrating the possible correlations between the significance of the residuals in I, q and u. We show them for the modeling corresponding to n_d ≡ ARM4⊕ED and GMF ≡ QSS; similar figures and results are obtained for all the other cases. We observe that the residuals on the polarization maps are not strongly correlated to the residuals on the intensity map. This reinforces our view that the reduced-Stokes parameter maps are more sensitive to the GMF than to the dust density distribution. Furthermore, we find that the polarization residuals are centered on zero and not correlated. While this is not shown on the figure, we observe that the scatter of the significance of the residuals (in I, q, and u) decreases with increasing Galactic latitude. This is observed for all of the models and we understand it as being due to the small uncertainties on q and u at low |b_gal|.

8.4 Comparison of reconstructed GMF

The fact that the best-fits GMF are very consistent with respect to the choice of the dust density distribution model can be observed, for example, in Fig. 21. In this figure we show the geometrical structure of the GMFin the horizontal and vertical cross-cuts of the Galaxy obtained with the three best-fit models of n_d when inferring theGMF geometry constraining the ASS model. The similarity of the recovered field geometrical structures is striking, especially in the (x, y, z = 0) plane. Equivalent figures are obtained for the other GMF models.

In order to compare quantitatively the GMF geometrical structures that we reconstructed when fitting the polarization maps, we proceed as in Sect. 5.3. For each best-fit model, we compute the pitch (p), tilt (χ), inclination (α) and position (γ) angles at each position of the Galactic space used for our MC realizations. We then compute the signed differences of each of the four angles considering the best-fit from the least evolved modeling (i.e. n_d ≡ ED + GMF ≡ ASS) as a reference to which each best-fit GMF model is compared with. The histograms of the signed differences are presented in Fig. 22.

In terms of the pitch and tilt angles, inspection of the figure leads to the following main results:

• (i)

  The pitch angles of all the best-fit GMF obtained for the ASS, BSS, and QSS agree within less than five degrees (with the exception of ARM4⊕ED + BSS, which is about nine degrees away).

• (ii)

  The histogram of the pitch angles of the LSA model (allowed to vary with Galactic radius) peaks at value near the values of the other GMF models. At the Sun radius, the pitch angle of this model is between six to eight degrees away from ourreference model.

• (iii)

  The tilt angles of all the models are centered around the same value. The scatter is however larger than that for the pitch angle when comparing BSS and QSS to ASS. This was expected due to the modulation of the field amplitude in planes parallel to the Galactic disk (z constant) that produces a change of the tilt angle of the GMF lines.

• (iv)

  Comparison of the best-fits of the same GMF model but for different n_d reveals a sizable coherence. We have checked that this conclusion, already informed around Fig. 18, equally holds for the other GMF than for the ASS, which is further discussed below.

As an example, the good agreement observed in the spirals of Fig. 21 for the ASS model is quantified by the blue histograms in the middle and right panels of the first row of Fig. 22, except that, here, we are not limited to (x, y, z = 0) or (x, y = 0, z) planes but account for the whole sampled space. The histograms show that through the whole sampled space the differences of the measured pitch angle are below one degree. The blue histograms in the middle and right panels of the second row quantify the small differences in the out-of-plane components among the ASS best-fits observed in lower panels of Fig. 21. The measured tilt angle differences (compared to the reference GMF) are less than about two degrees throughout the whole space.

In terms of the inclination (α) and position (γ) angles, which are the angles directly linked to the observable through line-of-sight integration (see Eq. (1)), the distributions of the relative angles shown in Fig. 22 (last two rows) also show an overall agreement between the model reconstructions. However, the wings of the distributions are more important for the angles α and γ than those of p and χ.

These results illustrate the non trivial mapping between the two set of angles, p and χ on one side and α and γ on the other side. Despite their similarities, the best-fit GMF models look different from the observer’s view, whereas when combined with the dust density and integrated along the lines of sight, all sets of polarization maps are equivalently a good (bad) fit to the data. This also illustrates the degree of degeneracy of the inclination and position angle of the GMF lines along the lines of sight.

In the above, it is the comparison of model-parameter values (and corresponding geometrical features) that is relevant, not the absolutevalues of the parameters since they are expected to be biased, despite the fact that we ignore the turbulent component of the GMF in the fit (see Sect. 6). However, we have shown that reliable constraints can be obtained for the pitchangle and therefore we limit our final conclusions to this parameter. We find with our analysis that the spiral pattern of the GMF has a mean pitch angle (at the Sun radius for the LSA model) of about 27 degrees, with a 14 percent scatter, and takes 17.5 and 35.4 degree as the two extreme values. Our pitch angle values are rather large but are still compatible with what can be expected from dynamo theory (e.g., Chamandy et al. 2016). We also find that, when allowed, the pitch angle vary dramatically from the inner part of the Galaxy to the outside region. A value of 27 degree for the pitch angle was already reported in Page et al. (2007). The out-of-plane pattern of the GMF is poorly constrained in our analysis. This is likely because it depends more on the particular GMF model under study and, at the same time, because its influence on maps is more noticeable at low and intermediate Galactic latitudes where our fits are not satisfying. Let us emphasize, though, that none of our best-fits exhibits the expected X-shape of the GMF lines that has been found in radio observations of external galaxies (see e.g., Jansson & Farrar 2012a for a discussion).

Fig. 21 Field line geometrical structures of best-fit ASS GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy for the three best-fit nd models obtained from the adjustment of the I353 map. From left to right: with nd ≡ ED, ARM4 and ARM4⊕ED, respectively.The color scale indicates the strength of GMF, i.e. the norm of the vector field at each location, which is fixed to 2.1 in the case of the ASS model.

Fig. 22 Histograms of the relative (local) angles as compared to our reference model defined as being the best fit obtained with nd = ED and the ASS GMF model. From top to bottom: pitch angle, tilt angle, inclination angles, and position angle. From left to right: nd = ED, ARMϕ and ARM4⊕ED. The colors correspond to different GMF best-fit model.

8.5 Alternative comparison of polarization maps

The degree of linear polarization ($p_{\text{lin}}=\sqrt{q^2+u^2}$) and the polarization position angle (ψ = 1∕2 arctan2(−u, q))¹¹ are commonly used to study and characterize the magnetized interstellar medium (e.g., Planck Collaboration Int. XIX 2015). In Fig. 23we show the noise debiased¹² p_lin and ψ maps from the Planck data at 353 GHz (top row) and those deduced from the polarization maps (q and u) corresponding to our best-fits GMF reconstructions obtained with the n_d ≡ ARM4 model fitted to I₃₅₃. These maps are alternative ways at looking at the q and u maps presented in Fig. 16 (but for another n_d model, which does not change the maps substantially). The joint consideration of all of the maps (q, u, p_lin and ψ) can be of interest, for instance, in diagnosing the limitation of our modeling and proposing further developments.

Inspection of the left column of Fig. 23 clearly illustrates qualitatively the limitation of the models at reproducing the data. Despite the noise dominated map of the data, it is nevertheless interesting to note how the modeling can qualitatively reproduce low p_lin values in some region of the Galactic equator where we know the intensity to be high. These low degrees of linear polarization results from line-of-sight depolarization induced by the varying GMF orientation. This can be inferred from the comparison, at low |b_gal|, of the BSS and QSS models (two last rows) to the ASS model (second row), for instance. Let us stress that this depolarization is obtained using the regular GMF only.

Furthermore, for lines of sight with p_lin ≲ 5%, the relatively simple models of the large-scale regular GMF that we consider are able to reproduce qualitatively the low and high values of the degree of linear polarization for low and high values of the intensity as suggested by the data (Planck Collaboration Int. XIX 2015). However, these models fail at reproducing very high degree of linear polarization (p_lin > 5%), where large difference are observed between model maps and data. These difference are well localized in the sky towards regions of nearby known structures (spurs, arms, and clouds) mainly attributable to the nearby structures such as the Sco-Cen association (e.g., Das et al. 2020). We used the gpempy software to investigate this further and we reached the conclusion that, in general, p_lin is not boosted by the simple addition of a clump of dust towards such a maximum to solve the problem. The addition of a clump actually produces the opposite. To reproduce the observed p_lin in these objects, it is likely that the effective degree of polarization of the dust population in these structures will have to be enhanced. This can be achieved by having locally either (i) a more efficient alignment of the dust grains along the field line or (ii) a higher intrinsic degree of linear polarization of the grain population, or (iii) a GMF particularly well aligned and coherent with respect to the line of sight (see Eq. (1)).

As opposed to p_lin, it is remarkable how the studied GMF models are able to capture qualitatively the broad tendency in the polarization position angle data (right column of Fig. 23). Indeed, the agreement between the models and the data is fair in terms of ψ, namely, in terms of the projected (and weighted) orientation of the GMF lines integrated along the lines of sight. Beside the fair resemblance, the modeled ψ maps show more coherence and smooth spatial variations than observed in the data. Also, inherently to the large-scale regular GMF considered in this study, the modeled maps have a mirror symmetry about the Galactic equator that is accompanied by a change of sign. In the data, that symmetry is disturbed on small scale and twisted and skewed on global scale. This is well illustrated in the upper right quadrant of the maps and comparing the ‘lines’ for which ψ = 0° = 180°, respectively. The inclusion of magnetic fields attached to local structures, such as depicted by Alves et al. (2018) and Pelgrims et al. (2020) in the case of the Local Bubble, addresses the latter issue in breaking the North-South symmetry. The modeling of (relatively) small scale perturbations would require the modeling of magnetized structures well localized in the Galaxy and will necessitate the modeling of the turbulent component, perhaps in the form of parameterized magneto-hydrodynamic turbulence, such as that depicted in Wang et al. (2020). Any such investigation is beyond the scope of this paper and will be addressed somewhere else.

Fig. 23 Degree of linear polarization (left) and polarization position angle (right), as deduced from the Planck data at 353 GHz (first row) and as deduced from our best-fits of the GMF models, while assuming the best-fit model from nd ≡ ARM4. Rows two tofive: ASS, LSA, BSS, and QSS models for the GMF, respectively.

8.6 Prospects

Improvements in the modelings of the regular part of the GMF with respect to this paper will require (i) the implementation of local magnetic field structure (e.g., Alves et al. 2018; Pelgrims et al. 2020); (ii) the interconnection of the latter with the large-scale regular GMF; and (iii) the generalization of the large-scale regular GMF models. As a straightforward generalized model, we would propose a model that contains both the radial variation of the pitch angle, such as encoded in the LSA model, and the strength field modulation, such as in the BSS or QSS model. Independently, both features appear indeed to be tentatively favored by dust data, as compared to the ASS model. The consideration of the other observables such as synchrotron emission and Faraday rotation measurement data promise to help significantly at constraining and modeling further the large-scale regular part of GMF and to solve possible ambiguities. In this respect,it might be worth undertaking physical models of the GMF such as those derived by Ferrière & Terral (2014).

A second stage of improvement in the GMF modeling will be to include the turbulent part of the GMF in the reconstruction. This step is of interest for CMB foreground studies as it allows for the statistical characterization of the polarized Galactic emission from thermal dust (e.g., see Vansyngel et al. 2017). However, in order to avoid overestimating the amplitude of the turbulent component compared to the regular one, we would stress the importance to progress as far as possible in the modeling of the regular component of the field without the introduction of random components. Indeed, given its effect at the map level (also discussed in Steininger et al. 2018), a random (turbulent) component can effectively help reconcile the model with the data. Therefore, we have to be cautious to not overestimate the turbulent component to compensate for a poor modeling of the regular part. The turbulence is certainly present in the Galaxy and must leave its imprint in the data, but perhaps, with a lower relative amplitude compared to the regular component than previously reported (e.g., Jansson & Farrar 2012a,b; Planck Collaboration Int. XLII 2016), especially at the physical scales probed by the Galactic diffuse emission considered at low resolution.

At the map level, the inclusion of a turbulent component also results in the two following effects: (i) an angular spatial de-coherence of the polarization position angles and (ii) an overall lowering of the degree of linear polarization, due to line of sight depolarization. The latter would call for an overall increase of the global amplitude of the modeled polarization maps and therefore would tend to enlarge the range of values span by Q and U. Such an increase might help at recovering the missing amplitude of the field strength, that can be visually inferred from comparing the data with the modeled polarization maps (see Figs. 16 or 23). We found that this apparent missing amplitude decrease the larger the pixels (alternatively the smaller the N_side) with which the fits are performed. This points to the fact that turbulence in the GMF is not the dominant factor when considering sufficiently large pixels in the sky. We found that the factor by which we should multiply the q and u maps for a better visual fit is, averaged on the 12 GMF reconstruction, of about 1.08 at N_side = 32 and 1.13 at N_side = 64. This observational fact opens new possible ways of tackling the modeling of the GMF and its different components. We have not investigated this possibility further and we leave this issue for a future work.

9 Summary and conclusions

In this paper, we present a methodology to infer the large-scale geometrical features of the GMF by constraining models of the regular components of the GMF relying on a dedicated MCMC method and a maximum-likelihood analysis of polarization data of the diffuse Galactic emission at microwave and millimeter frequencies, and more specifically, from the thermal dust emission.

Based on the current understanding of the modeling of the Galactic synchrotron and the thermal dust emission mechanisms in intensity and polarization, we built a step-by-step methodology to optimally divide and to simplify the fitting procedure. First, assuming the Galactic thermal dust intensity to be independent of the GMF, we can obtain a fair representation of the dust density distribution across the Galaxy. Second, we use the recovered thermal dust density distribution to fit the Galactic thermal dust polarization and so, to constrain the geometrical structure of the regular GMF component, assuming only this component is present in the data. Finally, the latter could be used as an input to constrain simultaneously the GMF amplitude and the Galactic relativistic electron density distribution by fitting synchrotron maps in intensity and polarization.

Relying on two sets of realistic simulations of the thermal polarized Galactic thermal dust emission, we demonstrate that it is possible to provide constraints on models of the regular component of the GMF and thus that it is possible to infer large-scale geometrical features of the GMF.

Assuming a non-turbulent contribution to the GMF in the simulations, we were able to provide excellent to fair reconstruction of the regular GMF geometrical structure, even in the case where we adopt an overly simplistic model (different and simpler parameterization with respect to the input one) for the Galactic dust density distribution. The latter case is likely to be representative of what happens when dealing with real data sets given the convoluted nature of the Galactic thermal dust emission in intensity, in particular. These relatively good results are possible thanks to the fitting of the reduced Stokes parameters (i.e., normalized to the intensity) that reduces the variations in the thermal dust polarization data induced by variations of the dust density along the lines of sight. In the case of simulations including a turbulent component of the GMF, we found that it is possible to satisfactorily recover the regular GMF component in the case of moderate turbulent contribution. When the turbulence component is on the same order of the regular one, as suggested by observation of the synchrotron sky, the reconstruction of the regular GMF degrades significantly, requiring the inclusion of the turbulent component in the fitting procedure. Nevertheless, reliable constraints can be obtained for the pitch angle of the regular GMF.

As a general result, we found that the uncertainties, in terms of model parameters and GMF geometry, arising from the lack of a consideration of the turbulence component are roughly of the same order of the ones induced by the mis-modeling of the density distribution or of the mis-modeling of the regular part of the GMF. This shows that improvements in the modelings are mandatory both for the deterministic (regular) and stochastic (turbulent) parts of the models. We also find that the maximum-likelihood uncertainties on the best-fit parameters are underestimated as they do not account for any mismodeling effects. Then, we argue that the scatter resulting from different modelings may be more representative of the true uncertainties.

Additionally, in the case of simulations without turbulence, we show that it appears possible to tackle the reconstruction of the distribution of the Galactic relativistic electrons by fitting the Galactic synchrotron emission using the best-fit of the GMF geometry obtained from the dust analysis as a prior.

Finally, we applied our methodology to a set of Planck maps of the polarized thermal dust emission at 353 GHz. We have been able toinfer large-scale geometrical features of the GMF by constraining different models for the dust density distributions and for the regular component of the GMF. We find that we are able to infer large-scale geometrical properties of the GMF and obtained coherent 3D views of the GMF, which, in terms of absolute parameter values, are limited by turbulence and the mismodeling of the dust density The most robust geometrical features that we infer is that at the Sun radius, the spiral of theGMF has a mean pitch angle of about 27 degrees with a 14 percent scatter depending on the precise modeling.

Despite some caveats, our analysis demonstrates that the Galactic thermal dust polarized emission is a powerful probe of the large-scale GMF, which is complementary to the more conventional ones (Galactic synchrotron data, Faraday rotation measures, and Faraday dispersion measures) and that it promises to play an important role in the overall efforts undertaken to reconstruct the actual 3D structure of the GMF (see e.g., Boulanger et al. 2018).

Acknowledgements

Special thanks to Céline Combet and David Maurin for the numerous discussions and encouragements during the realization of this work. This work has been funded by the European Union’s Horizon 2020 research and innovation program under grant agreement number 687312. V.P. also acknowledges support from the European Research Council under the European Union’s Horizon 2020 research and innovation program, under grant agreement No 771282. We acknowledge the use of data from the Planck/ESA mission, downloaded from the Planck Legacy Archive, and of the Legacy Archive for Microwave Background Data Analysis (LAMBDA). Support for LAMBDA is provided by the NASA Office of Space Science. Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package.

Appendix A 3D parametric models

In this paper, we consider a particular set of models for the dust density distribution and for the regular and random components of the GMF. These parametric models are described below.

A.1 Dust density distribution model

For the dust density distribution, we adopt two parametric models and their superposition. The first, simple, model reproduces the bell-shape profile of the Galaxy. The second, more complex model draws spiral-arm pattern similar as the model studied in Jaffe et al. (2013). The two models are parameterized in cylindrical coordinate (ρ, ϕ, z) centered on the Galactic center.

A.1.1 Exponential Disk (ED)

According to this model, the dust density distribution takes the parametric form: $$n_{\text{d}}(\rho ,z)=A_0\frac{\exp (-\rho /{\rho }_0)}{{(\cosh (z/z_0))}^2}.$$(A.1)

This model has two free parameters, ρ₀, z₀, to be fitted to the data. The global amplitude is absorbed in the maximum-profile likelihood computation and is thus irrelevant.

A.1.2.4 Spiral Arms (ARM4)

In this case, the dust density distribution is obtained by summing the contributions of four logarithmic spiral arms. In this implementation, and up to a relative amplitude, the arms are assumed to be identical with a rotation of ninety degrees. Our implementation is given as: $$n_{\text{d}}(\rho ,\varphi ,z)=A_0\frac{\exp (-{(\rho -{\rho }_c)}^2/(2{\rho }_0^2))}{{\left(\cosh (z/z_0)\right)}^2}\mathcal{S}(\rho ,\varphi )$$(A.2)

where the function $\mathcal{S}$ encodes the logarithmic spiral pattern as: $$\mathcal{S}(\rho ,\varphi )={\sum }_i\left(A_i\exp \left(\frac{-{(\varphi -{\varphi }_{s,i})}^2}{2{\varphi }_0^2}\right)\right)$$(A.3)

and $$\begin{array}{ll}{\varphi }_{s,i}=\hfill & \frac{1}{\tan (p)}\log (\rho /{\rho }_{0,i}),\hfill \\ {\rho }_{0,i}=\hfill & \exp \left({\varphi }_{0,i}\tan (p)\right),\hfill \\ {\varphi }_{0,i}=\hfill & {\varphi }_{00}+i\pi /2,\hfill \end{array}$$

with $i=\left\{1,2,3,4\right\}$, ϕ₀₀ the angular starting point of the fourth arm and p the pitch angle of the logarithmic spirals. Particular attention has to be made regarding the 2π ambiguity of polar angle during the evaluation of the function $\mathcal{S}$.

This parametric model has nine free parameters to be fitted to the data ρ₀, ρ_c, z₀, ϕ₀, ϕ₀₀, p, and three relative amplitudes for the spiral arms. As for the ED model, the global amplitude is irrelevant.

A.1.3.4 Spiral Arms plus 1 Exponential Disk (ARM4⊕ED)

We further consider a dust density model that results from the superposition of the ED and ARM4 models. This model has twelve free parameters to be fitted to the data. The radial and height scale lengths for the disk part and for the arm part are allowed to take different values.

Fig. A.1 GMF structures in the XY plane (top) and XZ plane (bottom) for the input LSA model used as an input to build the S1and S2 realistic simulations.

A.2 Regular Galactic Magnetic Field model

Throughout this paper, we make use of four parametric models for the regular component of the GMF. These are the Axi-Symmetric Spiral, the Logarithmic Spiral Arm, the Bi-symmetric Spiral, and the Quadri-Symmetric Spiral models. Despite their simplicity, these four mathematical toy models are expected to capture, at least to some degree, the large-scale geometrical features of the regular part of the magnetic field of our Galaxy with two to four free-parameters. The four GMF parametric models share the following main features:

• (i)

  field lines have spiral geometries in planes parallel to the Galactic disk (at z constant)

• (ii)

  field lines have out-of-the-plane component, intended to produce the so-called X-shape of the GMF

• (iii)

  the field lines have some degree of symmetry with respect to the z-axis of cylindrical Galactic coordinate system

The first three models have already been discussed in the literature. We briefly present the models below but refer the reader to e.g. Ruiz-Granados et al. (2010) for detailed discussion.

The four models take simple analytic expressions in the cylindrical coordinate system center on the Galactic center: $$B=B_{\rho }{e}_{\rho }+B_{\varphi }{e}_{\varphi }+B_z{e}_z,$$(A.7)

where we introduced the orthonormal basis (e_ρ, e_ϕ, e_z) with e_z = e_ρ × e_ϕ, where the polar angle ϕ increases counter-clockwise in the (x, y, z = 0) plane of the Galaxy. The four parametric models of the regular GMF component differ by their functional forms of the cylindrical components of the field.

A.2.1 Axi Symmetric Spiral model: ASS

The ASS model (see e.g., Vallee 1991; Poezd et al. 1993) is one of the simplest descriptions of the Galactic magnetic field. It is compatible with a non-primordial origin of the Galactic magnetism, based on the dynamo theory. The field lines follow logarithmic spirals with constant pitch angle. The cylindrical components of this model are: $$\begin{array}{ll}{B}_{\rho }=\hfill & B_0\sin (p)\cos (\chi (z))\hfill \\ B_{\varphi }=\hfill & B_0\cos (p)\cos (\chi (z))\hfill \\ B_z=\hfill & B_0\sin (\chi (z)),\hfill \end{array}$$(A.8)

where B₀ is the field strength that might, in general, depend on the distance to the Galactic center, p is the (constant) pitch angle and χ(z) is a “tilt angle” that allows the field lines to have a non-zero z component. Different modeling of the field strength dependence exist. We set B₀ to 2.1 μG throughout space because the polarized thermal dust emission is not sensitive to the field strength.

The pitch angle is defined as the angle between the azimuthal direction, e_ϕ, and the magnetic field. p = 0° corresponds to circle and p = 90° produces (cylindrical) radial lines.

The tilt angle is defined as the angle between the field line and the plane parallel to the Galactic plane. This angle usually takes the parametric form: $$\begin{array}{ll}\chi (z)=\hfill & {\chi }_0\tan h\left(\frac{z}{z_0}\right).\hfill \end{array}$$(A.9)

For the ASS model, we adopted the usual choice which is to fix the height scale at z₀ = 1 kpc. However, it is to be noted that this choice impact the value of χ₀ as these two parameters are not strictly independent.

According to the above implementation, the ASS model has two free-parameters (the pitch angle (p) and the tilt angle at large Galactic height (χ₀)) that can be adjusted to fit the polarization maps to the data.

A.2.2 LSA model

This parametric model, originally named the Logarithmic Spiral Arm, was introduced in Page et al. (2007) to describe the three-year WMAP data at 22 GHz for the synchrotron emission and at 94 GHz for the dust emission. Basically, this model implements an extension of the axi-symmetric class of models where the pitch angle defining the spiral structure depends on the cylindrical radius (ρ). Consequently, and despite the name of the model, the spirals are not logarithmic. As for the ASS model, there is no arm structure in the field amplitude.

The parametric form of the field components read $$\begin{array}{ll}{B}_{\rho }=\hfill & B_0\sin (\psi (\rho ))\cos (\chi (z))\hfill \\ B_{\varphi }=\hfill & B_0\cos (\psi (\rho ))\cos (\chi (z))\hfill \\ B_z=\hfill & B_0\sin (\chi (z)),\hfill \end{array}$$(A.10)

where the function $$\psi (\rho )={\psi }_0+{\psi }_1\ln \left(\frac{\rho }{R_{\odot }}\right)$$(A.11)

forces the magnetic field lines to follow a spiral pattern with varying pitch angle; ψ₀ is the pitch angle at Sun radius (R_⊙ = 8.0 kpc); and ψ₁ the amplitude of the logarithmic radial modulation of the pitch angle. The tilt angle χ(z) is taken as in Eq. (A.9) with z₀ = 1 kpc. Having fixed B₀ to a constant, the LSA model has three free parameters to be fixed by the data: (ψ₀, ψ₁ and χ₀). By construction, there is no reason to constrain the ψ₁ parameter to only positive or negative values. Fixing ψ₁ = 0 reduces the LSA model to the ASS model.

A.2.3 Bi Symmetric Spiral model: BSS

The class of bi-symmetric spiral models produces GMF that could be compatible with a primordial origin. This model includes magnetic field line reversion as suggested from pulsar rotation measurements (e.g., Han & Qiao 1994; Han et al. 2006). This model was also referred to as the Modified Logarithmic Spiral model in Fauvet et al. (2011). The field lines and the spirals drawn by the surfaces of iso-amplitude of the field have a constant and same pitch angle. The cylindrical components of this field model read $$\begin{array}{ll}{B}_{\rho }=\hfill & B_0\cos \left(\varphi \pm \beta \ln \left(\frac{\rho }{{\rho }_0}\right)\right)\sin (p)\cos (\chi (z))\hfill \\ B_{\varphi }=\hfill & B_0\cos \left(\varphi \pm \beta \ln \left(\frac{\rho }{{\rho }_0}\right)\right)\cos (p)\cos (\chi (z))\hfill \\ B_z=\hfill & B_0\sin (\chi (z)),\hfill \end{array}$$(A.12)

where β = 1∕tan(p), with p the pitch angle, χ(z) the tilt angle defined as in Eq. (A.9), B₀ the global amplitude of the field strength, which we assume to be a constant, and ρ₀ a radial scaling factor dictating the position of the spiral arms in the (x, y, z = 0) plane. In this model, the amplitude of the GMF is shaped in spirals by the term $\cos \left(\varphi \pm \beta \ln \left(\rho /{\rho }_0\right)\right)$. Because this modulation appears only in the components B_ρ and B_ϕ, it is not a simple scaling of the GMF strength. It produces changes in field line directions and thus produces the reversal of the field. The ± in the parentheses was introduced by Ruiz-Granados et al. (2010) to account for the different conventions met in the literature. Given our working scheme, we adopt the “−” sign so that the spiral pattern of the GMF amplitude and of the GMF lines coincide, which is reasonable to postulate.

Unlike for the ASS and LSA models, we let the parameter z₀, involved in the tilt angle modeling, to be free. The BSS model has therefore four free parameters that can be fitted by the observations.

A.2.4 Quadri Symmetric Spiral model: QSS

This is a class of GMF models is a slight modification of the BSS in which we multiply the argument responsible for the amplitude modulation by a factor of 2. This small change results in a class of GMF models having four spiral arms located ninety degrees away from one another. This four-parameter model therefore shares the four spiral arm features of the model presented in Jaffe et al. (2010, 2013) but with a much lower number of free-parameters, with field reversal and with out-of-plane component. This model has twice the number of field reversals than the BSS does. Such regular GMF models produce more efficient line-of-sight depolarization of the dust emission than the other models, a feature that might be called by the data.

A.3 Turbulent Galactic Magnetic Field model

We considered modeling the 3D turbulent part of the GMF at each position in the Galaxy (B^turb(r)), as derived from an isotropic random component that follows a two-slopes power-law magnetic energy spectrum.

To simulate the turbulent part of the GMF, each component of the turbulent GMF vector is obtained from a 3D Gaussian realization that follows a power-law power spectrum with two spectral indices, as suggested by Han et al. (2004); Han (2017). The composite spatial magnetic-energy spectrum takes the form P(k) = C k^α with α =−0.37 for k < k_out and α =−5∕3 for k ≥ k_out with the turnover at k_out = 10 corresponding to the outer-scale of the turbulence estimated at 0.1 kpc (Haverkorn et al. 2006, 2008).

Following Fauvet et al. (2011); Waelkens et al. (2009); Planck Collaboration Int. XLII (2016), the turbulent Gaussian random field, is computed in Fourier space leading to a 3D cube that contains the sampled Galaxy. In order to compute the turbulent component with the highest possible resolution in reasonable computing time, we perform both a low and a high spatial resolution simulation of the turbulent GMF in a cubic Cartesian grid of 1024 cells on a side. The physical resolution of the grids are 0.39 and3.9 kpc, respectively. In the Fourier space the high-resolution and low-resolution simulations span the ranges of k-values of $\left[0.25,256\right]$ and $\left[0.025,25.6\right]$ kpc⁻¹, respectively.The vector components of the turbulent magnetic field are then extracted at each location of our spherical sampling of the Galactic space (with N_side = 1024 as in Sect. 3.3), assuming a nearest neighbor scheme to map the two space samplings. The high- and low-resolution simulations are normalized in the overlapping range of k values (i.e., $k\in \left[0.25,25.6\right]$ kpc⁻¹). Inside the heliocentric radius of 2 kpc, the two turbulent components (high and low resolution) are considered in order to preserve continuity of large-scale fluctuations in the neighborhood of the Sun. Outside the heliocentric radius of 2 kpc, only the low resolution box is used, we add random white noise with an rms that corresponds to the spatial frequencies not sampled in the low-resolution box but that are included in the high-resolution one (namely, for $k\in \left[25.6,256\right]$). This procedure ensures an identical rms of the turbulent GMF vectors throughout the sampled space. Unlike Planck Collaboration Int. XLII (2016); Steininger et al. (2018), we do not re-scale the turbulent GMF component as a function of position in the Galaxy. This is because the regular GMF model of Page et al. (2007), the only one to which the turbulent component is added in this study, has a constant amplitude.

Appendix B Resolution bias

To efficiently explore the model parameter spaces while running MCMC method, large number of simulations are necessary. It is therefore often necessary to work at resolution lower than that of the original data to keep the computing time reasonable. Here, we explore possible consequences of this choice in the final results.

Fig. B.1 Distribution of the relative difference per pixel between low resolution maps produced at several HEALPix Nside (shown in the x-axis) and the high resolution map at Nside = 1024 downgraded to equivalent resolution for the signal term of the S2 simulated maps. From left to right: distribution for Stokes parameter I, Q and U. The shaded-blue (-green) regions contains the 95% (68%) of the total number of pixels.

First, we have tested the impact of the choice of the angular resolution of the HEALPix tessellation on the simulated observables for the models presented in Sect. 3.3. This comparison study is illustrated in Fig. B.1, where we represent the distribution of the relative differences of the models as computed for N_side = 8, 16, 32, 64, 128, 256, 512 with respect to the ones computed at N_side = 1024 and downgraded to the respective lower N_side values. For each pixel, the relative differences is computed as 2(S_HR − S_LR)∕(S_HR + S_LR), where the subscripts HR and LR stand for high resolution and low resolution, respectively; and S represents any of the linear polarization Stokes parameters.

In Fig. B.1, the dashed-red line represents the median of the distribution of the N_pix measurements of the relative difference, the green-shaded region contains 68% of the pixels and the blue-shaded region 95%. In that figure we see that, at the pixel level, a scatter arises due to the resolution at which the simulations are computed. This is a sampling issue. The lower the map resolution, the coarser the sampling of the 3D space. At the map level this leads to features that are either missed or grossly represented. As a consequence, other values of the model parameters may lead to a better fit to the observed downgraded features, thus producing an overall shift (or bias) in the parameter space.

To remedy this issue, Wang et al. (2020) suggested to interpolate high-resolution maps at the sky coordinates of the centers of the low-resolution-map pixels. Interpolation can indeed reduce the magnitude of this sampling issue but will not make it disappear in general. Moreover, interpolation has its drawbacks. The resulting maps depend on the interpolation scheme adopted and that the propagation of uncertainties is more complex and in any case not analytic. In this work, we do not consider this possibility any further. We note that, as shown in Fig. B.1, the high and low angular-resolution simulations agree on pixel values at the percent level starting from N_side = 64 and the agreement increases as N_side increases.

Second, we investigate the bias introduced in the best-fit parameters when working at resolution lower than that of the native data. For this purpose we consider the case of fitting the simulation S1 because the underlying dust density distribution model is the simplest and the loss of information due to line-of-sight integration is expected to be negligible, if any. That is, for this simple case the source of bias in recovered parameter values mainly comes from the effect of the simulation resolution. To demonstrate the significance of the resolution bias, we performed the MCMC fits of the intensity map and of the polarization maps at the N_side values of 8, 16, and 32, in addition to the fits in N_side = 64 presented in the core of the paper. For the two parameters of the dust density distribution and for the four parameters of the LSA GMF model, we then evaluate the difference between the best-fit parameters and the input ones, in sigma units,: $|p_0-\widehat{p}|/{\sigma }_{\widehat{p}}$. Here, p₀ is the input-parameter value, $\widehat{p}$ is the best-fit parameter value corresponding to the minimum χ² and ${\sigma }_{\widehat{p}}$ is computed as being the standard deviation of the marginalized distribution of the considered parameter as sampled by the MCMC algorithm.

Fig. B.2 Significance of the resolution bias of the best-fit parameter values as compared to the input parameter values as a function of the Nside value used to generate the models in the MCMC algorithm. for simulation S1. Top panel: two parameters of the ED dust densitydistribution. Bottom panel: behavior for the four GMF parameters.

The results of this analysis are shown in Fig. B.2. It is seen that the best-fit parameters converge toward the input values as the resolution of the model used in the MCMC analysis increases. The convergence is however not trivial and depends upon the considered parameter and also on its input values. To this respect, we note that for some parameters, the best-fit values may oscillate around the input across the N_side values. The value of the reduced χ² obtained at the different resolutions go from 984 to 1.02 for the intensity fits and from 24.2 to 1.935 for the fits of the (q, u) maps.

We conclude that due to the resolution bias the values of the best-fit parameters for dust density distribution and the GMF, these models need to be taken with caution. Furthermore, the best-fit parameters obtained at a given resolution should not be used at other resolutions. Formally speaking, it should be possible to account for such a bias in the definition of the likelihood function. However, we have observed that the bias depends strongly both on the parametric form of the models and on the parameter values themselves. Therefore, accounting for it is cumbersomeand goes beyond the scope of this paper. Overcoming the resolution bias will be mandatory once the models to be fitted to the real data sets will be sufficiently evolved to account for all the complexity and the richness contained in them. Indeed, the resolution bias can be significant only if it turns out to be the leading source of uncertainties. However, this can happen in the case signal dominated data, as it is already the case for the 353 GHz Planck data, and when the model accurately represents the data, what is clearly not the case as of today (see also Sect. 8). Therefore, for now, we notice that this bias exists and that the larger the N_side value of the working resolution, the smaller the bias is.

Third, we tested the impact of the adopted value for the radial sampling (sampling along the line-of-sight). As far as we do not pretend to model the very nearby structures, we found that the radial step can be as large as few hundreds of parsecs. The difference inpixel space is more pronounced due to the adopted angular sampling than due to the radial sampling. However, a fine radialsampling would be essential as soon as nearby and physically small features are introduced in the modeling.

In summary, it turns out to be important to work at the highest resolution allowed by the computing facilities. As a result, in this paper we choose to simulate the maps at N_side = 64 and with a radial bin of 0.2 kpc. For the models considered the computation time of a full set of Stokes parameters maps is always of the order of few tenths of a second. A greater computing time would limit significantly our ability to run the MCMC algorithm.

Appendix C Orthographic view of the Planck data and of the best-fit models

In Fig. C.1, we provide an alternative view of the intensity data and obtained best fits than the one presented inthe core of the paper. The agreement between the best-fit models at high Galactic latitudes is salient in these orthographic views. We note that here we maintain the same color scale as in Fig. 14.

Fig. C.1 Orthographic view of the intensity maps. First row show the 353-GHz map from Planck downgraded at Nside = 64 and the corresponding map of uncertainties that we use to compute the χ2. Rows two to four correspond to dust density distribution models labeled ED, ARM4 and ARM4⊕ED, respectively.The obtained best fits are shown in the first column and the statistical significance of the residual, per-pixel, are shown in the second column. The north Galactic pole is on the right-hand side panel and the vertical solid line is for Galactic longitude zero. Galactic latitude decrease counter-clockwise on the right panel and clockwise in the left panel.

In Fig. C.2 we provide an alternative view of the fitted polarization data and of the obtained best-fits than the one presented in the core of the paper. Inspection of the maps of the significance of the residuals teaches us that our best-fit models agree well at high Galactic latitudes. The same color scale as in Fig. C.2 is adopted.

Fig. C.2 Orthographic view of the polarization maps. First row shows the 353-GHz maps of the reduced Stokes parameters from Planck downgraded at Nside = 64 and the corresponding map of uncertainties that we use to compute the χ2 (q, σq, u, σu). Rows 2–5: GMF models labeled ASS, LSA, BSS, and QSS using the best-fit of the ED nd model. The obtained best fits are shown in the first and third columns and the statistical significance of their residuals, per-pixel, are shown in the second and fourth columns. The same convention is used as for Fig. C.1.

Appendix D Second solution for the LSA GMF model

Fig. D.1 GMF geometrical structure for the best-fit presented in the core of the paper (top row) and for the one obtained from the second local minimum (bottom row). We show the best-fit ASS GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy for the three best-fit nd models obtained from the adjustment of the I353 map. From left to right: nd ≡ ED, ARM4, and ARM4⊕ED, respectively.

We illustrate here the existence of local minima in the hyper surface of the χ² that we minimize to determine the best-fit to the Planck data. For this purpose, we present the second solution that we obtained when fitting the LSA GMF model in the second row of Fig. D.1. We also present this solution because it illustrates the importance of considering large and as much un-informative prior as possible when fittingmodels to the existing dataset. Indeed, the LSA GMF model has recently been fitted to Faraday rotation measures and synchrotron data by Steininger et al. (2018), but while exploring a reduced parameter space compared tothe one we considered in this paper. It is interesting to notice that the second solution of the GMF, which we present here as a local minimum, appears to be consistent with their solution.

The best-fit model for this second solution depends more strongly on the underlying dust density distribution than the best-fit models that we discussed in the core of the paper and that we present (for the three dust density model) in the first row of Fig. D.1.

References

All Tables

All Figures

	Fig. 1 Schematic view of the polarization direction of the Galactic synchrotron and dust thermal emission as a function of the Galactic magnetic field (GMF). For completeness, we also show the direction of the starlight polarization.
In the text

	Fig. 2 High-resolution simulations of the I, Q, U Stokes parameters (from left to right) obtained for S1 (top) and S2 (bottom). Units are supposedly KCMB. Intensity map is given in colored-logarithmic scale and polarization maps in linear scale. Full-sky maps are in Galactic coordinates, centered on the Galactic center. Galactic longitude decrease towards the right-hand side of the plots.
In the text

	Fig. 3 Corner plot showing the marginalized 1D and 2D posterior probability distributions for the parameters of the regular LSA GMF model in the case where the simulated data and the fit are performed at Nside = 64, see Sect. 4.5. The vertical and horizontal light-blue lines mark the input parameter values. Angles are given in degrees and the scale height (z0) in kpc.
In the text

	Fig. 4 From left to right, top: data from the S1 simulation downgraded to Nside = 64, middle: best-fit maps obtained while assuming the dust distribution follows an ED model and the GMF the LSA model (Case A in the text), and bottom: significance of the residuals. The color scale of the last row ranges from −5 to 5.
In the text

Fig. 5 Maps corresponding to the modified Stokes I, q, and u from left to right. Top: data from the S2 simulation downgraded to Nside = 64; second row: best-fit maps obtained while assuming the dust distribution follows an ED model (case B); third row: best-fit maps obtained while assuming the dust distribution follows an ARM4 model (case C). The two last rows show the significance of the residuals for the case B (4th row) and the case C (5th row). For the case B, the color scales range from −500 and 500 and from −15 to 15 for the intensity and the polarization, respectively. For the case C, color scales range from −5 to 5, as in Fig. 4.

In the text

	Fig. 6 Reconstructed GMF structures in the XY plane (top) and XZ plane (bottom) in case A, case B, and case C, respectively, shown from top to bottom. See text for details. Similar plot for the input GMF model is shown in Fig. A.1.
In the text

Fig. 7 Histograms of the angle differences (in degree) between the best reconstructions of the GMF and the input one computed at every location of the sampled Galactic space (left). From top to bottom: differences of the pitch angles, the tilt angles, the inclination angles, and the position angles (see text). Case C overlaps case A in tilt and position angle differences. Right: same as for left column but we present the cumulative distributions of the absolute values of the angle differences.

In the text

	Fig. 8 Histogram of the differences of the polarization position angles between the data and the best fits. The polarization position angles are deduced from the q and u maps both from the downgraded simulations (S1 and S2 but without the Planck noise) and from the best-fit maps.
In the text

	Fig. 9 For the S1 simulation, map of the reduced Stokes q downgraded at Nside = 64 for the cases with Aturb = 0.0, Aturb = 0.3 and Aturb = 0.9, from left to right. Top panels results from regular and turbulent field; Planck noise is added in the bottom panels.
In the text

Fig. 10 Corner plots of the 1D and 2D marginalized posterior distributions of the regular GMF parameter for the three reconstruction cases discussed in Sect. 6.2 applied to the S2-turb simulations. Angles are given in degree and the scale heights in kpc. From left to right: columns correspond to Aturb = 0.0, Aturb = 0.3 and Aturb = 0.9, respectively. From top to bottom for cases I to III defined as: case I: nd ≡ ARM4 and gmf ≡ LSA (both models underlying the data); case II: nd ≡ ED and gmf ≡ LSA; and case III: nd ≡ ED and gmf ≡ ASS. The vertical and horizontal light-blue lines mark input parameter values.

In the text

Fig. 11 Cumulative distributions the difference angles between the regular GMF model in simulations S2-turb and the best-fit reconstructed models using the MCMC likelihood analysis. From top to bottom for the pitch angles, the tilt angles, the inclination angles, and the position angles. The three reconstruction cases corresponds to different models assumed to model the data: case I: nd ≡ ARM4 and gmf ≡ LSA (both models underlying the data); case II: nd ≡ ED and gmf ≡ LSA; and case III: nd ≡ ED and gmf ≡ ASS. From left to right, the columns correspond to Aturb = 0.00, Aturb = 0.30, and Aturb = 0.90. The shaded area corresponds to the 68% CL extracted from the posterior of the probability distribution as sampled from the MCMC analysis.

In the text

Fig. 12 Maps of the synchrotron emission. Columns correspond to the maps of I, q, u, and polarization position angles (in the IAU convention). The first row correspond to synchrotron emission when the original GMF model is adopted. For the second row, the adopted GMF corresponds to the best-fit obtained from the analysis of dust maps, case B in Sect. 5. In the third row, we present the relative difference of the intensity maps ((Iin − Î)∕(Iin + Î)). For the polarization angles, we display the acute angles between the polarization vectors.

In the text

	Fig. 13 Acute angle histogram comparing the polarization vector orientations from the synchrotron emission assuming the original GMF model and the worst one from dust analysis (case B). Angles are in degree.
In the text

Fig. 14 Intensity maps. First row: 353-GHz map from Planck downgraded at Nside = 64 and the corresponding map of uncertainties that we use to compute the χ2. Rows two to four correspond to dust density distribution models labeled ED, ARM4, and ARM4⊕ED, respectively.The obtained best-fits are shown in the first column and the statistical significance of the residual, per-pixel, are shown in the second column.

In the text

	Fig. 15 Histograms of the signed significance of the residuals of the best-fits obtained by adjusting the three dust density distribution models (ED, ARM4, and ARM4⊕ED) to the thermal dust intensity map.
In the text

Fig. 16 Polarization maps. First row: 353-GHz maps of the reduced Stokes parameters from Planck downgraded at Nside = 64 and the corresponding map of uncertainties that we use to compute the χ2 (q, σq, u, σu). Rows two to five: GMF models labeled ASS, LSA, BSS and QSS using the best-fit of the ED nd model. The obtained best-fits are shown in the first and third columns and the statistical significance of their residuals, per-pixel, are shown in the second and fourth columns.

In the text

	Fig. 17 Histograms of the significance of the residuals. x stands for the concatenation of the q and u parameters. Each GMF model is represented by a different color and each panel corresponds to a different underlying nd model: ED, ARM4, and ARM4⊕ED from left to right.
In the text

Fig. 18 Functional forms of the pitch (top) and tilt (bottom) angles as a function of the Galactic cylindrical coordinate; ρ and z, respectively.Each line corresponds to a best-fit model. Solid, dashed, dash-dotted, and dotted lines correspond to the ASS, LSA, BSS, and QSS GMF models. Blue, green, and red correspond to the best-fits obtained while assuming the nd model to be ED, ARM4 or ARM4⊕ED, respectively.The pitch angle is constant for the ASS, BSS, and QSS models.

In the text

Fig. 19 Fieldline geometrical structures of best-fit GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy. From top-left to bottom-right: ASS, LSA, BSS, and QSS GMF models. The GMF reconstructed models corresponds to the best-fits obtained while assuming that nd is our best-fit ED model. The black dot at the center of the figure is the Galactic center and the black star is at the Sun position. The color scales offer information about the relative strength of GMF (i.e., the norm of the vector field at each location) and of its direction (red going clockwise, blue going counter-clockwise). In the case of ASS and LSA models, the (constant) norm has been fixed to 2.1.

In the text

	Fig. 20 Corner plot showing the correlations between the significance of the residuals in I, q, and u. This figure corresponds to the modeling with the ARM4⊕ED nd model and with the QSS GMF model.
In the text

Fig. 21 Field line geometrical structures of best-fit ASS GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy for the three best-fit nd models obtained from the adjustment of the I353 map. From left to right: with nd ≡ ED, ARM4 and ARM4⊕ED, respectively.The color scale indicates the strength of GMF, i.e. the norm of the vector field at each location, which is fixed to 2.1 in the case of the ASS model.

In the text

	Fig. 22 Histograms of the relative (local) angles as compared to our reference model defined as being the best fit obtained with nd = ED and the ASS GMF model. From top to bottom: pitch angle, tilt angle, inclination angles, and position angle. From left to right: nd = ED, ARMϕ and ARM4⊕ED. The colors correspond to different GMF best-fit model.
In the text

	Fig. 23 Degree of linear polarization (left) and polarization position angle (right), as deduced from the Planck data at 353 GHz (first row) and as deduced from our best-fits of the GMF models, while assuming the best-fit model from nd ≡ ARM4. Rows two tofive: ASS, LSA, BSS, and QSS models for the GMF, respectively.
In the text

	Fig. A.1 GMF structures in the XY plane (top) and XZ plane (bottom) for the input LSA model used as an input to build the S1and S2 realistic simulations.
In the text

Fig. B.1 Distribution of the relative difference per pixel between low resolution maps produced at several HEALPix Nside (shown in the x-axis) and the high resolution map at Nside = 1024 downgraded to equivalent resolution for the signal term of the S2 simulated maps. From left to right: distribution for Stokes parameter I, Q and U. The shaded-blue (-green) regions contains the 95% (68%) of the total number of pixels.

In the text

	Fig. B.2 Significance of the resolution bias of the best-fit parameter values as compared to the input parameter values as a function of the Nside value used to generate the models in the MCMC algorithm. for simulation S1. Top panel: two parameters of the ED dust densitydistribution. Bottom panel: behavior for the four GMF parameters.
In the text

Fig. C.1 Orthographic view of the intensity maps. First row show the 353-GHz map from Planck downgraded at Nside = 64 and the corresponding map of uncertainties that we use to compute the χ2. Rows two to four correspond to dust density distribution models labeled ED, ARM4 and ARM4⊕ED, respectively.The obtained best fits are shown in the first column and the statistical significance of the residual, per-pixel, are shown in the second column. The north Galactic pole is on the right-hand side panel and the vertical solid line is for Galactic longitude zero. Galactic latitude decrease counter-clockwise on the right panel and clockwise in the left panel.

In the text

Fig. C.2 Orthographic view of the polarization maps. First row shows the 353-GHz maps of the reduced Stokes parameters from Planck downgraded at Nside = 64 and the corresponding map of uncertainties that we use to compute the χ2 (q, σq, u, σu). Rows 2–5: GMF models labeled ASS, LSA, BSS, and QSS using the best-fit of the ED nd model. The obtained best fits are shown in the first and third columns and the statistical significance of their residuals, per-pixel, are shown in the second and fourth columns. The same convention is used as for Fig. C.1.

In the text

Fig. D.1 GMF geometrical structure for the best-fit presented in the core of the paper (top row) and for the one obtained from the second local minimum (bottom row). We show the best-fit ASS GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy for the three best-fit nd models obtained from the adjustment of the I353 map. From left to right: nd ≡ ED, ARM4, and ARM4⊕ED, respectively.

In the text