Issue 
A&A
Volume 652, August 2021



Article Number  A130  
Number of page(s)  31  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201833962  
Published online  23 August 2021 
Galactic magnetic field reconstruction using the polarized diffuse Galactic emission: formalism and application to Planck data
^{1}
Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSCIN2P3,
38000
Grenoble,
France
email: pelgrims@physics.uoc.gr
^{2}
Institute of Astrophysics, Foundation for Research and TechnologyHellas,
Vasilika Vouton,
70013
Heraklion,
Greece
^{3}
Department of Physics, and Institute for Theoretical and Computational Physics, University of Crete,
70013
Heraklion,
Greece
^{4}
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge,
MA
02139,
USA
Received:
26
July
2018
Accepted:
11
May
2021
The polarized Galactic synchrotron and thermal dust emission constitutes a major tool in the study of the Galactic magnetic field (GMF) and in constraining its strength and geometry for the regular and turbulent components. In this paper, we review the modeling of these two components of the polarized Galactic emission and present our strategy for optimally exploiting the currently existing data sets. We investigate a Markov chain Monte Carlo (MCMC) method to constrain the model parameter space through maximumlikelihood analysis, focusing mainly on dust polarized emission. Relying on simulations, we demonstrate that our methodology can be used to constrain the regular GMF geometry. Fitting for the reduced Stokes parameters, this reconstruction is only marginally dependent of the accuracy of the reconstruction of the Galactic dust grain density distribution. However, the reconstruction degrades, apart from the pitch angle, when including a turbulent component on the order of the regular one as suggested by current observational constraints. Finally, we applied this methodology to a set of Planck polarization maps at 353 GHz to obtain the first MCMC based constrains on the largescale regularcomponent of the GMF from the polarized diffuse Galactic thermal dust emission. By testing various models of the dust density distribution and of the GMF geometry, we prove that it is possible to infer the largescale geometrical properties of the GMF. We obtain coherent threedimensional views of the GMF, from which we infer a mean pitch angle of 27 degrees with 14% scatter, which is in agreement with results obtained in the literature from synchrotron emission.
Key words: polarization / dust, extinction / submillimeter: ISM / ISM: magnetic fields / methods: statistical / cosmic background radiation
© V. Pelgrims et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 Bmodes 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 componentseparation 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 lineofsight 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^{1} and Planck^{2} 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 fullsky 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; RuizGranados et al. 2010; Jaffe et al. 2010, 2013; Fauvet et al. 2011, 2012; Jansson & Farrar 2012a,b).
Because fullsky polarization data first came at low frequencies, these investigations were driven by the study of the fullsky synchrotron emission. Page et al. (2007) used the threeyear fullsky maps from the WMAP satellite at 22 GHz (the Kband) and fitted a parametric model using the polarization position angles of the emission. RuizGranados et al. (2010) used the fiveyear WMAP polarization data at the same frequency and searched for the best fits of several parametric models on a gridbased 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 353GHz data from the ARCHEOPS balloon experiment (Benoît & ARCHEOPS Collaboration 2004) in addition to the WMAP satellite 22GHz channel data (tracing the polarized diffuse Galactic synchrotron emission) to constrain GMF models. Jaffe et al. (2013) used the fullsky WMAP 94GHz 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 fullsky 353GHz data from Planck is in conflict with predictions from the reconstructed largescale 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 orderedrandom, 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 orderedrandom) 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. Threedimensional (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 highdimensional 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 largescale 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 largescale 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 fullsky 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 submillimeter 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 (1)
where r is the radial distance from the observer along the lineofsight at sky position, n. The different terms in the equation are:

, the dust emissivity at observational frequency ν, which is linked to the dust temperature through a greybody’s law (e.g., Planck Collaboration Int. XX 2015, Appendix B).
p^{d}, the socalled 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 socalled 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: (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 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: (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: (4)
where 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: (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 planeofsky component of the magnetic field (B_{⊥}). We note that , 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 equalarea 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 selfconsistent suite a Python codes that follows the above implementations. It is publicly available^{3} where we describe the architecture that is optimized for quick and userfriendly 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 largescale 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 largescale 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 nontrivial mixing of the matter density distribution and the geometrical structure of the GMF. However, we also notice that the thermal dust intensity is to firstorder 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 intensityweighted 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 mismodeling 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 freefree) 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 singlefrequency intensity and polarization maps at 353 GHz from the second release (PR2)^{4} that are available on the Planck Legacy Archive^{5}. 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^{6} of about 4.94′ (Planck Collaboration Int. XIX 2015) and are given in a HEALPix^{7} grid tessellation corresponding to N_{side} = 2048 (Górski et al. 2005). At the instrument resolution, the 353GHz polarization Stokes Q and U maps are noisedominated (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 fullsky dust extinction map derived in Planck Collaboration XIX (2014), specifically, the map of the dust optical depth (τ_{353}). 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 highresolution maps that simulate the Planck 353GHz linear polarization Stokes parameter maps. These sets, hereafter called S1, S2, and S2turb are defined as follows.
First, S1 is built by combining the exponentialdisk (ED) dust density distribution model with the LSA GMF model, with no turbulent component. Second, S2 is built by combining the fourspiralarms dust density model (ARM4) with the LSA GMF model; no turbulent component. And finally, S2turb 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 regularmodel 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^{8}. 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 353GHz are given by: (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 Highresolution simulations of the I, Q, U Stokes parameters (from left to right) obtained for S1 (top) and S2 (bottom). Units are supposedly K_{CMB}. Intensity map is given in coloredlogarithmic scale and polarization maps in linear scale. Fullsky maps are in Galactic coordinates, centered on the Galactic center. Galactic longitude decrease towards the righthand side of the plots. 
Free parameters of the parametric models for the two dust density distribution models (ED and ARM4) and for the GMF model (LSA).
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 largescale 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 bestfit 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 fittedmap 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 modelparameter 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 subdominant 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 MCMCbased exploration of the parameter spaces are thus fulfilled.
4.2 Likelihood definition
The bestfit parameters are obtained by maximizing a Gaussian likelihood function , defined as (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: .
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 signaldominated 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 signaldominated case close to the Galacticequator, while at high Galactic latitudes we are mainly in a noisedominated 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 perpixel covariance matrix maps released by the Planck collaboration at N_{side} = 2048. We neglect the offdiagonal terms^{9} and propagate the uncertainties at low resolution according to (8)
where 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: (9)
where D = {I, Q, U}^{10}. 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 lowresolution 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 χ^{2} 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 χ^{2} 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: (10)
These pseudo uncertainties are then propagated to the reduced Stokes parameters, q and u as: (11)
where x = X∕I and X = {Q, U}. And so, we write the covariance matrix for {q, u} as .
Reduced χ^{2} values for the bestfit in intensity and in polarization for the three study cases.
4.4 Implementation of the MCMC experiment
In order to explore the parameter space, determine the posterior distributions of the parameters, and find their bestfit values, we used the emcee MCMC Python software written by ForemanMackey et al. (2013), who implemented the AffineInvariant 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 χ^{2} hypersurface 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 χ^{2} and waited for the convergence criteria to be fulfilled to reach wellsampled 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 noninformative priors, adopting tophat 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 largescale 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 bestfit model of the dust density distribution (n_{d}). Then, we use the bestfit 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 lineofsight 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 bestfit parameter values that we obtained are reported in Table 1, both for the intensity and the polarization fits. The obtained values of thereduced χ^{2} 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 N_{side} = 64, see Sect. 4.5. The vertical and horizontal lightblue lines mark the input parameter values. Angles are given in degrees and the scale height (z_{0}) in kpc. 
5.1 Reconstruction of the dust density
The first steps in our fitting procedure is to provide bestfit 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 χ^{2} = 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 spacesampling resolution has larger effects.
For cases B and C, the intensity fit results are shown in Fig. 5. The bestfit 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 χ^{2} of 1.9 × 10^{5}. 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 χ^{2} do not account for mismodeling. 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 χ^{2} = 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 lineofsight 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 N_{side} = 64, middle: bestfit 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 bestfit models for the intensity map, we use the bestfit 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 bestfit 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 lineofsight, the elemental emitting volumes following Eq. (1) to produce the Q and U maps. Thenwe normalize these simulated maps by the bestfit 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 bestfit 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 bestfit (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 χ^{2} values of Table 2. The bestfit values of the GMF parameters are reported in Table 1 along with their marginalized uncertainties at the 1σ level. It is seen that the bestfit 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 bestfit 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 bestfit 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 Xshape (see Appendix A.2).
The reconstructed geometrical structures of the GMF corresponding to the cases A, B, and C bestfit 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 N_{side} = 64; second row: bestfit maps obtained while assuming the dust distribution follows an ED model (case B); third row: bestfit 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 outofplane 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 be any of the angles defined above and measured for the input GMF model and the same angle but for the reconstructed GMF models. We define the signed difference angle 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 for the angles between the field lines of a bestfit 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 bestfit 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 (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 bestfit 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): (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 twoslopes powerlaw 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 signalonly (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 S2turb (see Sect. 3.3) sets of simulated maps to estimate the uncertainties induced in the reconstruction of the largescale 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 S2turb 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 ψ_{0} 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 nonconsideration 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 S2turb 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 largescale 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 loglikelihood 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 modelparameter space that is induced by the nonconsideration 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 N_{side} = 64 for the cases with A_{turb} = 0.0, A_{turb} = 0.3 and A_{turb} = 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 (ρ_{0}, z_{0}) = (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 bestfit 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 MCMCfitting procedure described in Sect. 4 to the Planck polarization data at 353GHz (see Sect. 3.2) to infer the largescale structure of the GMF. We first rely on three dust density distribution models to fit for the intensity map. Then, we use those bestfit 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 reducedStokes 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 S2turb simulations. Angles are given in degree and the scale heights in kpc. From left to right: columns correspond to A_{turb} = 0.0, A_{turb} = 0.3 and A_{turb} = 0.9, respectively. From top to bottom for cases I to III defined as: case I: n_{d} ≡ ARM4 and gmf ≡ LSA (both models underlying the data); case II: n_{d} ≡ ED and gmf ≡ LSA; and case III: n_{d} ≡ ED and gmf ≡ ASS. The vertical and horizontal lightblue 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 bestfit 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 S2turb and the bestfit 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: n_{d} ≡ ARM4 and gmf ≡ LSA (both models underlying the data); case II: n_{d} ≡ ED and gmf ≡ LSA; and case III: n_{d} ≡ ED and gmf ≡ ASS. From left to right, the columns correspond to A_{turb} = 0.00, A_{turb} = 0.30, and A_{turb} = 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 Bestfit models
The bestfit parameters are reported in Table 3. In Fig. 14, we compare the data and the bestfit intensity maps obtained for the three dust density distribution models fitted to I_{353}. We also show χmaps, which give the statistical significance of the residuals. We observe that the bestfit 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 bestfits 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 smallscale feature turns out to be significant and the values of the reduced χ^{2} are ≫ 1 (see Table 5). This was to be expected from our simulationbased study where we show that such high values of reduced χ^{2} are reached as soon as we do not have the right parametric model at hand (see Sect. 5).
To verify the robustness of our bestfit 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 (τ_{353}) map presented inPlanck Collaboration XIX (2014). That observable, deduced from the multifrequency analysis, is known to trace the integrated dust density distribution at least as well as the intensity map does. The bestfit maps agree globally with those obtained from I_{353} fits. This, in turn, also justifies the assumption around Eq. (3).
In conclusion, none of the bestfit 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 largescale, 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 bestfit 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 ((I_{in} − Î)∕(I_{in} + Î)). 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: 353GHz map from Planck downgraded at N_{side} = 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 bestfits are shown in the first column and the statistical significance of the residual, perpixel, are shown in the second column. 
Free parameters of the dust density distribution models are given for the three explored modelings (ED, ARM4, and ARM4⊕ED).
Fig. 15 Histograms of the signed significance of the residuals of the bestfits 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 freeparameters 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 Bestfit models
In Fig. 16, we show the maps corresponding to the bestfit 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 bestfit 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 bestfits 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 χ^{2} 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_{353} 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 χ^{2} of the bestfit models are one order of magnitude better than those obtained for the fits of the intensity map.
Furthermore, the obtained bestfits 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 χ^{2} 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 bestfit 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 bestfit models, we also constrain the GMF parameters using the n_{d} bestfit models obtained from the fits on the dust opticaldepth map. The fitted (q, u) maps are obtained simply by replacement of the intensity map by the τ_{353} map, with the corresponding uncertainties. The bestfit 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_{353} 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 bestfit 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 mapmaking 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).
Free parameters of the GMF models are given for the three explored modelings (ASS, LSA, BSS, and QSS).
Fig. 16 Polarization maps. First row: 353GHz maps of the reduced Stokes parameters from Planck downgraded at N_{side} = 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 bestfit of the ED n_{d} model. The obtained bestfits are shown in the first and third columns and the statistical significance of their residuals, perpixel, 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 n_{d} model: ED, ARM4, and ARM4⊕ED from left to right. 
Reduced χ^{2} values of the bestfit models (i) of the dust density distribution on the intensity map and (ii) of the GMF models on the reducedjointedStokes maps (q, u) for each of thebestfit of the n_{d} model.
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 bestfit model. Solid, dashed, dashdotted, and dotted lines correspond to the ASS, LSA, BSS, and QSS GMF models. Blue, green, and red correspond to the bestfits obtained while assuming the n_{d} 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 bestfit 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 largescale 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 χ^{2} 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 Faradaydepth allsky 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 ψ_{1} parameter (see Eq. (A.11)) to span through negative values. An interesting feature of our bestfit 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 bestfit GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy. From topleft to bottomright: ASS, LSA, BSS, and QSS GMF models. The GMF reconstructed models corresponds to the bestfits obtained while assuming that n_{d} is our bestfit 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 counterclockwise). 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 n_{d} model and with the QSS GMF model. 
8.3 Bestfit maps: comparison and characteristics
Based on inspection of residual maps and on χ^{2} values, none of the twelve bestfit 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 bestfit 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 bestfit 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 reducedStokes 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 bestfits 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 crosscuts of the Galaxy obtained with the three bestfit 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 bestfit 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 bestfit from the least evolved modeling (i.e. n_{d} ≡ ED + GMF ≡ ASS) as a reference to which each bestfit 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 bestfit 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 bestfits 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 outofplane components among the ASS bestfits 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 lineofsight 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 bestfit 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 modelparameter 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 outofplane 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 bestfits exhibits the expected Xshape 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 bestfit ASS GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy for the three bestfit n_{d} models obtained from the adjustment of the I_{353} map. From left to right: with n_{d} ≡ 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 n_{d} = ED and the ASS GMF model. From top to bottom: pitch angle, tilt angle, inclination angles, and position angle. From left to right: n_{d} = ED, ARMϕ and ARM4⊕ED. The colors correspond to different GMF bestfit model. 
8.5 Alternative comparison of polarization maps
The degree of linear polarization () and the polarization position angle (ψ = 1∕2 arctan2(−u, q))^{11} 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^{12} 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 bestfits GMF reconstructions obtained with the n_{d} ≡ ARM4 model fitted to I_{353}. 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 lineofsight 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 largescale 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 ScoCen 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 largescale 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 NorthSouth 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 magnetohydrodynamic 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 bestfits of the GMF models, while assuming the bestfit model from n_{d} ≡ 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 largescale regular GMF; and (iii) the generalization of the largescale 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 largescale 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 decoherence 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 largescale geometrical features of the GMF by constraining models of the regular components of the GMF relying on a dedicated MCMC method and a maximumlikelihood 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 stepbystep 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 largescale geometrical features of the GMF.
Assuming a nonturbulent 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 mismodeling of the density distribution or of the mismodeling 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 maximumlikelihood uncertainties on the bestfit 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 bestfit 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 largescale 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 largescale 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 largescale 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 bellshape profile of the Galaxy. The second, more complex model draws spiralarm 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: (A.1)
This model has two free parameters, ρ_{0}, z_{0}, to be fitted to the data. The global amplitude is absorbed in the maximumprofile 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: (A.2)
where the function encodes the logarithmic spiral pattern as: (A.3)
with , ϕ_{00} 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 .
This parametric model has nine free parameters to be fitted to the data ρ_{0}, ρ_{c}, z_{0}, ϕ_{0}, ϕ_{00}, 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 AxiSymmetric Spiral, the Logarithmic Spiral Arm, the Bisymmetric Spiral, and the QuadriSymmetric Spiral models. Despite their simplicity, these four mathematical toy models are expected to capture, at least to some degree, the largescale geometrical features of the regular part of the magnetic field of our Galaxy with two to four freeparameters. 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 outoftheplane component, intended to produce the socalled Xshape of the GMF
 (iii)
the field lines have some degree of symmetry with respect to the zaxis 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. RuizGranados et al. (2010) for detailed discussion.
The four models take simple analytic expressions in the cylindrical coordinate system center on the Galactic center: (A.7)
where we introduced the orthonormal basis (e_{ρ}, e_{ϕ}, e_{z}) with e_{z} = e_{ρ} × e_{ϕ}, where the polar angle ϕ increases counterclockwise 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 nonprimordial 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: (A.8)
where B_{0} 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 nonzero z component. Different modeling of the field strength dependence exist. We set B_{0} 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: (A.9)
For the ASS model, we adopted the usual choice which is to fix the height scale at z_{0} = 1 kpc. However, it is to be noted that this choice impact the value of χ_{0} as these two parameters are not strictly independent.
According to the above implementation, the ASS model has two freeparameters (the pitch angle (p) and the tilt angle at large Galactic height (χ_{0})) 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 threeyear 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 axisymmetric 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 (A.10)
forces the magnetic field lines to follow a spiral pattern with varying pitch angle; ψ_{0} is the pitch angle at Sun radius (R_{⊙} = 8.0 kpc); and ψ_{1} the amplitude of the logarithmic radial modulation of the pitch angle. The tilt angle χ(z) is taken as in Eq. (A.9) with z_{0} = 1 kpc. Having fixed B_{0} to a constant, the LSA model has three free parameters to be fixed by the data: (ψ_{0}, ψ_{1} and χ_{0}). By construction, there is no reason to constrain the ψ_{1} parameter to only positive or negative values. Fixing ψ_{1} = 0 reduces the LSA model to the ASS model.
A.2.3 Bi Symmetric Spiral model: BSS
The class of bisymmetric 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 isoamplitude of the field have a constant and same pitch angle. The cylindrical components of this field model read (A.12)
where β = 1∕tan(p), with p the pitch angle, χ(z) the tilt angle defined as in Eq. (A.9), B_{0} the global amplitude of the field strength, which we assume to be a constant, and ρ_{0} 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 . 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 RuizGranados 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_{0}, 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 fourparameter 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 freeparameters, with field reversal and with outofplane component. This model has twice the number of field reversals than the BSS does. Such regular GMF models produce more efficient lineofsight 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 twoslopes powerlaw 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 powerlaw power spectrum with two spectral indices, as suggested by Han et al. (2004); Han (2017). The composite spatial magneticenergy 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 outerscale 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 highresolution and lowresolution simulations span the ranges of kvalues of and kpc^{−1}, 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 lowresolution simulations are normalized in the overlapping range of k values (i.e., kpc^{−1}). Inside the heliocentric radius of 2 kpc, the two turbulent components (high and low resolution) are considered in order to preserve continuity of largescale 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 lowresolution box but that are included in the highresolution one (namely, for ). 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 rescale 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 N_{side} (shown in the xaxis) and the high resolution map at N_{side} = 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 shadedblue (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 dashedred line represents the median of the distribution of the N_{pix} measurements of the relative difference, the greenshaded region contains 68% of the pixels and the blueshaded 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 highresolution maps at the sky coordinates of the centers of the lowresolutionmap 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 angularresolution 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 bestfit 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 lineofsight 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 bestfit parameters and the input ones, in sigma units,: . Here, p_{0} is the inputparameter value, is the bestfit parameter value corresponding to the minimum χ^{2} and 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 bestfit parameter values as compared to the input parameter values as a function of the N_{side} 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 bestfit 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 bestfit values may oscillate around the input across the N_{side} values. The value of the reduced χ^{2} 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 bestfit parameters for dust density distribution and the GMF, these models need to be taken with caution. Furthermore, the bestfit 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 lineofsight). 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 bestfit 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 bestfit 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 353GHz map from Planck downgraded at N_{side} = 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, perpixel, are shown in the second column. The north Galactic pole is on the righthand side panel and the vertical solid line is for Galactic longitude zero. Galactic latitude decrease counterclockwise 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 bestfits than the one presented in the core of the paper. Inspection of the maps of the significance of the residuals teaches us that our bestfit 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 353GHz maps of the reduced Stokes parameters from Planck downgraded at N_{side} = 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 bestfit of the ED n_{d} model. The obtained best fits are shown in the first and third columns and the statistical significance of their residuals, perpixel, 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 bestfit presented in the core of the paper (top row) and for the one obtained from the second local minimum (bottom row). We show the bestfit ASS GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy for the three bestfit n_{d} models obtained from the adjustment of the I_{353} map. From left to right: n_{d} ≡ ED, ARM4, and ARM4⊕ED, respectively. 
We illustrate here the existence of local minima in the hyper surface of the χ^{2} that we minimize to determine the bestfit 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 uninformative 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 bestfit model for this second solution depends more strongly on the underlying dust density distribution than the bestfit 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
 Alves, M. I. R., Boulanger, F., Ferrière, K., & Montier, L. 2018, A&A, 611, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Andersson, B.G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501 [Google Scholar]
 Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20 [Google Scholar]
 Benoît, A., & ARCHEOPS Collaboration 2004, Adv. Space Res., 33, 1790 [NASA ADS] [CrossRef] [Google Scholar]
 BICEP2/Keck & Planck Collaborations 2015, Phys. Rev. Lett., 114, 101301 [Google Scholar]
 BICEP2 & Keck Array Collaborations 2016, Phys. Rev. Lett., 116, 031302 [Google Scholar]
 Boulanger, F., Ensslin, T., Fletcher, A., et al. 2018, J. Cosmol. Astropart. Phys., 2018, 049 [Google Scholar]
 Chamandy, L., Shukurov, A., & Taylor, A. R. 2016, ApJ, 833, 43 [Google Scholar]
 Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Das, K. K., Zucker, C., Speagle, J. S., et al. 2020, MNRAS, 498, 5863 [Google Scholar]
 Davis, Jr. L., & Greenstein, J. L. 1951, ApJ, 114, 206 [Google Scholar]
 Efroimsky, M. 2002, ApJ, 575, 886 [NASA ADS] [CrossRef] [Google Scholar]
 Fauvet, L., MacíasPérez, J. F., Aumont, J., et al. 2011, A&A, 526, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fauvet, L., MacíasPérez, J. F., Jaffe, T. R., et al. 2012, A&A, 540, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ferrière, K., & Terral, P. 2014, A&A, 561, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ForemanMackey, D., Conley, A., Meierjurgen Farr, W., et al. 2013, Astrophysics Source Code Library [record ascl:1303.002] [Google Scholar]
 Gelman, A., & Rubin, D. B. 1992, Stat. Sci., 7, 457 [Google Scholar]
 Ginzburg, V. L., & Syrovatskiĭ S. I. 1966, Sov. Phys. Uspekhi, 8, 674 [Google Scholar]
 Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Computat. Sci., 5, 65 [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Han, J. L. 2017, ARA&A, 55, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Han, J. L., & Qiao, G. J. 1994, A&A, 288, 759 [NASA ADS] [Google Scholar]
 Han, J. L., Ferriere, K., & Manchester, R. N. 2004, ApJ, 610, 820 [NASA ADS] [CrossRef] [Google Scholar]
 Han, J. L., Manchester, R. N., Lyne, A. G., Qiao, G. J., & van Straten, W. 2006, ApJ, 642, 868 [NASA ADS] [CrossRef] [Google Scholar]
 Haverkorn, M., Gaensler, B. M., Brown, J. C., et al. 2006, ApJ, 637, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Haverkorn, M., Brown, J. C., Gaensler, B. M., & McClureGriffiths, N. M. 2008, ApJ, 680, 362 [NASA ADS] [CrossRef] [Google Scholar]
 Hildebrand, R. H., Kirby, L., Dotson, J. L., Houde, M., & Vaillancourt, J. E. 2009, ApJ, 696, 567 [NASA ADS] [CrossRef] [Google Scholar]
 Hoang, T. 2017, ArXiv eprints [arXiv:1704.01721] [Google Scholar]
 Jaffe, T. R., Leahy, J. P., Banday, A. J., et al. 2010, MNRAS, 401, 1013 [NASA ADS] [CrossRef] [Google Scholar]
 Jaffe, T. R., Ferrière, K. M., Banday, A. J., et al. 2013, MNRAS, 431, 683 [NASA ADS] [CrossRef] [Google Scholar]
 Jansson, R., & Farrar, G. R. 2012a, ApJ, 757, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Jansson, R., & Farrar, G. R. 2012b, ApJ, 761, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Jordan, M. E., & Weingartner, J. C. 2009, MNRAS, 400, 536 [NASA ADS] [CrossRef] [Google Scholar]
 King, L. W., & Harwit, M. O. 1973, IAU Symp., 52, 197 [Google Scholar]
 Lazarian, A., Efroimsky, M., & Ozik, J. 1996, ApJ, 472, 240 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, H. M., & Draine, B. T. 1985, ApJ, 290, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, P. G. 2007, EAS Pub. Ser., 23, 165 [CrossRef] [EDP Sciences] [Google Scholar]
 Onaka, T. 1996, ASP Conf. Ser., 97, 72 [Google Scholar]
 Onaka, T. 2000, ApJ, 533, 298 [NASA ADS] [CrossRef] [Google Scholar]
 Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJS, 170, 335 [Google Scholar]
 Pelgrims, V., Ferrière, K., Boulanger, F., Lallement, R., & Montier, L. 2020, A&A, 636, A17 [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2020, A&A, 641, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2014, A&A, 571, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2016, A&A, 594, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2016, A&A, 594, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2016, A&A, 594, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XIX. 2015, A&A, 576, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XX. 2015, A&A, 576, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXI. 2015, A&A, 576, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLII. 2016, A&A, 596, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLIV. 2016, A&A, 596, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poezd, A., Shukurov, A., & Sokoloff, D. 1993, MNRAS, 264, 285 [NASA ADS] [CrossRef] [Google Scholar]
 RuizGranados, B., RubiñoMartín, J. A., & Battaner, E. 2010, A&A, 522, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (Hoboken: Wiley) [Google Scholar]
 Steininger, T., Enßlin, T. A., Greiner, M., et al. 2018, ArXiv eprints [arXiv:1801.04341] [Google Scholar]
 Sun, X.H., & Reich, W. 2010, Res. Astron. Astrophys., 10, 1287 [NASA ADS] [CrossRef] [Google Scholar]
 Sun, X. H., Reich, W., Waelkens, A., & Enßlin, T. A. 2008, A&A, 477, 573 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vaillancourt, J., Andersson, B. G., & Lazarian, A. 2013, in Proceedings of The Life Cycle of Dust in the Universe: Observations, Theory, and Laboratory Experiments (LCDU2013). 1822 November, 2013. Taipei, Taiwan. Editors: Anja Andersen (University of Copenhagen, Denmark), Maarten Baes (Universiteit Gent, Belgium), Haley Gomez (Cardiff University, Ciska Kemper (Academia Sinica, Darach Watson (University of Copenhagen, Denmark). Online at http://pos.sissa.it/cgibin/reader/conf.cgi?confid=207, 3 [Google Scholar]
 Vallee, J. P. 1991, ApJ, 366, 450 [NASA ADS] [CrossRef] [Google Scholar]
 Vansyngel, F., Boulanger, F., Ghosh, T., et al. 2017, A&A, 603, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Waelkens, A., Jaffe, T., Reinecke, M., Kitaura, F. S., & Enßlin, T. A. 2009, A&A, 495, 697 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, J., Jaffe, T. R., Enßlin, T. A., et al. 2020, ApJS, 247, 18 [CrossRef] [Google Scholar]
 Zucker, C., Battersby, C., & Goodman, A. 2018, ApJ, 864, 153 [NASA ADS] [CrossRef] [Google Scholar]
These maps correspond to the dataset selected for the RADIOFOREGROUNDS project. We note that more recent sets of polarization maps have been delivered since the work and analyses presented in this paper were conducted (see, e.g., the third release of Planck data (PR3) Planck Collaboration III 2020).
Noise debiasing procedure is performed in the same way as in Appendix B.2 of Planck Collaboration Int. XIX (2015).
All Tables
Free parameters of the parametric models for the two dust density distribution models (ED and ARM4) and for the GMF model (LSA).
Reduced χ^{2} values for the bestfit in intensity and in polarization for the three study cases.
Free parameters of the dust density distribution models are given for the three explored modelings (ED, ARM4, and ARM4⊕ED).
Free parameters of the GMF models are given for the three explored modelings (ASS, LSA, BSS, and QSS).
Reduced χ^{2} values of the bestfit models (i) of the dust density distribution on the intensity map and (ii) of the GMF models on the reducedjointedStokes maps (q, u) for each of thebestfit of the n_{d} model.
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 Highresolution simulations of the I, Q, U Stokes parameters (from left to right) obtained for S1 (top) and S2 (bottom). Units are supposedly K_{CMB}. Intensity map is given in coloredlogarithmic scale and polarization maps in linear scale. Fullsky maps are in Galactic coordinates, centered on the Galactic center. Galactic longitude decrease towards the righthand 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 N_{side} = 64, see Sect. 4.5. The vertical and horizontal lightblue lines mark the input parameter values. Angles are given in degrees and the scale height (z_{0}) in kpc. 

In the text 
Fig. 4 From left to right, top: data from the S1 simulation downgraded to N_{side} = 64, middle: bestfit 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 N_{side} = 64; second row: bestfit maps obtained while assuming the dust distribution follows an ED model (case B); third row: bestfit 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 bestfit maps. 

In the text 
Fig. 9 For the S1 simulation, map of the reduced Stokes q downgraded at N_{side} = 64 for the cases with A_{turb} = 0.0, A_{turb} = 0.3 and A_{turb} = 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 S2turb simulations. Angles are given in degree and the scale heights in kpc. From left to right: columns correspond to A_{turb} = 0.0, A_{turb} = 0.3 and A_{turb} = 0.9, respectively. From top to bottom for cases I to III defined as: case I: n_{d} ≡ ARM4 and gmf ≡ LSA (both models underlying the data); case II: n_{d} ≡ ED and gmf ≡ LSA; and case III: n_{d} ≡ ED and gmf ≡ ASS. The vertical and horizontal lightblue lines mark input parameter values. 

In the text 
Fig. 11 Cumulative distributions the difference angles between the regular GMF model in simulations S2turb and the bestfit 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: n_{d} ≡ ARM4 and gmf ≡ LSA (both models underlying the data); case II: n_{d} ≡ ED and gmf ≡ LSA; and case III: n_{d} ≡ ED and gmf ≡ ASS. From left to right, the columns correspond to A_{turb} = 0.00, A_{turb} = 0.30, and A_{turb} = 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 bestfit 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 ((I_{in} − Î)∕(I_{in} + Î)). 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: 353GHz map from Planck downgraded at N_{side} = 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 bestfits are shown in the first column and the statistical significance of the residual, perpixel, are shown in the second column. 

In the text 
Fig. 15 Histograms of the signed significance of the residuals of the bestfits 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: 353GHz maps of the reduced Stokes parameters from Planck downgraded at N_{side} = 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 bestfit of the ED n_{d} model. The obtained bestfits are shown in the first and third columns and the statistical significance of their residuals, perpixel, 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 n_{d} 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 bestfit model. Solid, dashed, dashdotted, and dotted lines correspond to the ASS, LSA, BSS, and QSS GMF models. Blue, green, and red correspond to the bestfits obtained while assuming the n_{d} 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 bestfit GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy. From topleft to bottomright: ASS, LSA, BSS, and QSS GMF models. The GMF reconstructed models corresponds to the bestfits obtained while assuming that n_{d} is our bestfit 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 counterclockwise). 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 n_{d} model and with the QSS GMF model. 

In the text 
Fig. 21 Field line geometrical structures of bestfit ASS GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy for the three bestfit n_{d} models obtained from the adjustment of the I_{353} map. From left to right: with n_{d} ≡ 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 n_{d} = ED and the ASS GMF model. From top to bottom: pitch angle, tilt angle, inclination angles, and position angle. From left to right: n_{d} = ED, ARMϕ and ARM4⊕ED. The colors correspond to different GMF bestfit 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 bestfits of the GMF models, while assuming the bestfit model from n_{d} ≡ 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 N_{side} (shown in the xaxis) and the high resolution map at N_{side} = 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 shadedblue (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 bestfit parameter values as compared to the input parameter values as a function of the N_{side} 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 353GHz map from Planck downgraded at N_{side} = 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, perpixel, are shown in the second column. The north Galactic pole is on the righthand side panel and the vertical solid line is for Galactic longitude zero. Galactic latitude decrease counterclockwise 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 353GHz maps of the reduced Stokes parameters from Planck downgraded at N_{side} = 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 bestfit of the ED n_{d} model. The obtained best fits are shown in the first and third columns and the statistical significance of their residuals, perpixel, 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 bestfit presented in the core of the paper (top row) and for the one obtained from the second local minimum (bottom row). We show the bestfit ASS GMF models in the (x, y, z = 0) and (x, y = 0, z) planes of the Galaxy for the three bestfit n_{d} models obtained from the adjustment of the I_{353} map. From left to right: n_{d} ≡ ED, ARM4, and ARM4⊕ED, respectively. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.