Galactic magnetic field reconstruction using the polarized diffuse Galactic emission: Formalism and application to $Planck$ data
Vincent Pelgrims Juan Francisco Macías-Pérez Florian Ruppin
 


   

   

   Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 38000 Grenoble, France
   Institute of Astrophysics, Foundation for Research and Technology-Hellas, GR-71110 Heraklion, Greece
   Department of Physics, and Institute for Theoretical and Computational Physics, University of Crete, GR-70013 Heraklion, Greece
   Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
   

              pelgrims@physics.uoc.gr ; macias@lpsc.in2p3.fr
             

   
        GMF Reconstruction

  

   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 maximum-likelihood
   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 large-scale regular-component 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 large-scale geometrical
        properties of the GMF.
        We obtain coherent three-dimensional (3D) 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.
   

   

   Galactic magnetic field reconstruction using the polarized diffuse Galactic emission: Formalism and application to Planck data
    V. Pelgrims
          1,2,3
          J.F. Macías-Pérez
          1
          F. Ruppin
          1,4
          
    Received July 26, 2018; accepted May 11, 2021
==============================================================================================================================







§ 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 ).

In particular, possible contamination from the polarized diffuse
Galactic emission has been shown to be one of the major limitations for
the detection of primordial CMB polarization B-modes related to the
inflationary era in the early universe
<cit.>.
Therefore, providing an accurate characterization and modeling of this
polarized Galactic emission is of great importance as it would strengthen the level of confidence
with regard to our understanding of its physical scope and would allow for accurate testing of 
the results obtained from elaborated component-separation techniques 
<cit.> which are
used to extract the cosmological CMB signal. 


At microwave frequencies that are typically relevant for CMB experiments, the polarized sky is dominated by synchrotron emission below  ∼ 80 GHz and by thermal dust emission above that frequency.
Both emission components result from a line-of-sight integration of local emission and these offer the possibility to infer the properties of the Galactic magnetic field (GMF), an important constituent and actor in the ecosystem of our Galaxy.
The diffuse polarized Galactic synchrotron emission is produced by relativistic electrons that spiral along the GMF  lines <cit.>.
Equivalently, the polarized thermal dust emission is produced by rotating aspherical dust grains that are totally or partially aligned with the GMF lines (; ; ; ; ; ; ; ; ; ).



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[<https://map.gsfc.nasa.gov/>]
and Planck[<http://www.esa.int.Planck>] satellite experiments (; ; ), which have observed the sky in polarization in a large range of frequency going from 23 to 353 GHz.
The wealth of information present in these full-sky observations is extremely valuable in modeling the various components of the Galaxy.
In light of these polarization data, there have been attempts to constrain the geometry of the GMF at the largest Galactic scales
and the relative amplitude of their main components (regular, turbulent,...) and of the Galactic matter content
(; ; ; ; ; ; ; ).

Because full-sky polarization data first came  at low frequencies, these
investigations were driven by the study of the full-sky synchrotron
emission.
<cit.> used the three-year full-sky maps from the WMAP
satellite at 22 GHz (the K-band) and fitted a parametric model using the
polarization position angles of the emission. 
<cit.> used the five-year WMAP polarization data at the
same frequency and searched for the best fits of several parametric models
on a grid-based exploration of the parameter space,
thus providing the first systematic comparison of GMF models.
<cit.>, <cit.>, <cit.>,
<cit.>, <cit.>, and <cit.>
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, <cit.> 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 <cit.> for a recent work.

Similar studies have been carried out
using the  polarized thermal dust emission. <cit.> showed that the 94 GHz band of the WMAP satellite measured the thermal dust emission and used it to constrain the GMF. <cit.> used the 353-GHz data from the
ARCHEOPS balloon experiment <cit.> in addition to the WMAP satellite 22-GHz channel data (tracing the polarized diffuse Galactic synchrotron emission) to constrain GMF models.
<cit.> used the full-sky WMAP 94-GHz polarization maps and showed that the diffuse Galactic emission observed at this frequency is not compatible with GMF configuration that fits best the polarized synchrotron emission as traced by the WMAP low frequency data. More recently, <cit.> showed that the full-sky 353-GHz data from Planck is in conflict with predictions from the reconstructed large-scale GMF obtained from the synchrotron emission data.
Other studies have considered the thermal dust emission only from the polar caps to adjust models of the magnetic field in the neighborhood of the Sun (; ; ).



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



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

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

In this paper we motivate and explore a somewhat different approach to infer the main characteristics of the large-scale GMF.
In Sect. <ref>, 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. <ref>, we discuss the main methodology used as well as the data and simulations used. 
Section <ref> discusses the fitting procedure in intensity and polarization for the dust thermal and its application in simulating the data.
In Sect. <ref>, we discuss the quality of our reconstruction of the regular part of the large-scale GMF, along with possible limitations and biases. 
Sect. <ref> evaluates the impact of those limitations when independently tackling Galactic synchrotron emission.
We apply our methodology to the full-sky polarization 353 GHz Planck maps in Sect. <ref>. 
Finally, we summarize our findings and present our conclusions in Sect. <ref>. 




§ MODELING & IMPLEMENTATION





 §.§ 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 <cit.>.
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 <cit.>.
As a consequence, the thermal dust emission arising at sub-millimeter and millimeter wavelengths is expected to be linearly polarized perpendicular to the GMF lines, as sketched in Fig. <ref>.
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 <cit.>, but which follows the physically motivated modeling of the emission introduced by <cit.> and reviewed by <cit.>.
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


    I_d(𝐧) =     ϵ_ν^d ∫_0^+∞dr   n_d(r,𝐧)  {
    1 + p^d  f_ma(2/3 - sin^2 α(r,𝐧))
    }
    
    Q_d(𝐧) =     ϵ_ν^d  p^d  f_ma∫_0^+∞dr   n_d(r,𝐧)  sin^2 α(r,𝐧)  cos[2 γ(r,𝐧)]
    
    U_d(𝐧) =     ϵ_ν^d  p^d  f_ma∫_0^+∞dr   n_d(r,𝐧)  sin^2 α(r,𝐧)  sin[2 γ(r,𝐧)]  ,


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


  * ϵ_ν^d, the dust emissivity at observational frequency ν, which is linked to the dust temperature through a grey-body's law (e.g.,, Appendix B). 

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

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

  * n_d(r, 𝐧), the 3D Galactic dust grain density.

  * α(r, 𝐧), the inclination angle of the GMF line with the line of sight at (r, 𝐧).

  * γ(r, 𝐧), the so-called local polarization angle.

The local polarization angle is defined in the plane orthogonal to the line of sight as the angle between the polarization vector direction and the local meridian. Expressed in terms of the vector component of the ambient GMF, this angle is given as:
 

    γ(r,  𝐧) = 1/2arctan( -2   B_θ(r,  𝐧)   B_ϕ(r,  𝐧)/B_ϕ(r,  𝐧)^2 - B_θ(r,  𝐧)^2)   ,

with B_θ and B_ϕ the local transverse components of the magnetic field in the local spherical coordinate basis (𝐞_r, 𝐞_θ, 𝐞_ϕ) with 𝐞_θ pointing towards the South pole. Thus, Eq. <ref>  gives the polarization position angle of the polarization vectors stemming from the small space volume in the HEALPix (or COSMO) convention <cit.>, which differs from the more commonly used IAU one. This angle is defined in the range [ 0 ,   180 [ 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 Eqs. <ref> is negligible compared to the first term:

    I_d(𝐧) =  ϵ_ν^d ∫_0^+∞ dr   n_d(r,𝐧)  .

This simplifying assumption has been made in all previous works in the field (<cit.>). Relying on simulations, we found that in pixel space, the relative difference is at most of 10 percent (see also Sect. <ref>).




 §.§ Synchrotron polarized emission


The diffuse Galactic synchrotron emission is produced by relativistic electrons that spiral around the GMF lines <cit.>.
This emission is to be polarized perpendicularly to the  GMF lines as sketched in Fig. <ref>. For the contribution of synchrotron emission to CMB frequencies (≳ 20 GHz) we follow the modeling of the emission presented by <cit.> and <cit.> that relies on <cit.>. We adopt the notation of <cit.>.

Assuming a power law for the relativistic electron energy the linear polarization Stokes parameters for the diffuse Galactic synchrotron emission is given as:


    I_s(𝐧) =     ϵ_ν^s ∫_0^+∞ dr  n_e(r,𝐧) (𝐁_⊥(r,𝐧) ^2 )^(s+1)/4
    
    Q_s(𝐧) =     ϵ_ν^s  p_s ∫_0^+∞ dr  n_e(r,𝐧) (𝐁_⊥(r,𝐧) ^2 )^(s+1)/4 cos[2γ(r,𝐧)] 
    
    U_s(𝐧) =     ϵ_ν^s  p_s ∫_0^+∞ dr  n_e(r,𝐧) (𝐁_⊥(r,𝐧) ^2 )^(s+1)/4 sin[2γ(r,𝐧)],


where ϵ_ν^s is the synchrotron emissivity, n_e(r,𝐧) is the local density of relativistic electrons, p_s is the intrinsic synchrotron polarization fraction which is related to the relativistic electron energy spectral index (s) as follows:

    p_s = s+1/s+7/3 .


The angle γ found in the expression of Q_s and
U_s is the same as in Eq. <ref>.
It marks the orientation perpendicular to the plane-of-sky component
of the magnetic field (𝐁_⊥).
We note that 𝐁_⊥(r,𝐧) ^2 = 𝐁(r,𝐧)^2 sin^2[α(r,𝐧)], where the angle α(r,𝐧) is the same as in Eq. <ref>; namely, it is the inclination angle of the GMF vector with respect to the line of sight.







 §.§ Implementation


To simulate the Stokes parameters I, Q, and U of the polarized diffuse thermal dust and synchrotron Galactic emission, we follow Eqs. <ref> and <ref>.
We chose to sample the Galactic space according to a spherical coordinate system centered on the Sun. The radial distance to the observer is linearly sampled and we adopt the equal-area HEALPix tessellation <cit.> 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 to the midpoint rule.
Our choice of the sampling of the Galactic space and the overall implementation are thus similar to the one adopted in the  code implementation <cit.>, except that we do not consider refining the angular grid as 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 , a self-consistent suite a Python codes that follows the above implementations.
It is publicly available[<http://www.radioforegrounds.eu/pages/software/gmf-reconstruction.php>]
where we describe the architecture that is optimized for quick and user-friendly simulations of the polarized sky.
These codes are not optimized for an MCMC exploration of parameter spaces of models.



 §.§ Parametric models

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





§ METHODOLOGY, DATA & SIMULATIONS





 §.§ Methodology


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



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 the reduced Stokes parameters at low frequency as other Galactic emission components (e.g., from anomalous microwave emission and free-free) are also important in intensity at CMB frequencies.
Thus, a careful separation of these Galactic components is needed to safely compute the reduced synchrotron Stokes parameters
q_s = Q_s/I_s and u_s = U_s/I_s.
In contrast, assuming that a
good model of the geometrical structure of the GMF can be obtained from thermal dust emission, this model could then be used to constrain
the relativistic electron density distribution and of the GMF strength through a fit of the synchrotron data.



 §.§ Data


In this paper, we mainly concentrate on the diffuse polarized thermal dust emission.
We use the Planck single-frequency intensity and polarization maps at 353 GHz from the second release (PR2)[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) ).] that are available on the Planck Legacy Archive[http://pla.esac.esa.int/pla/#maps].
We refer the reader to <cit.> for details and discussions regarding these data.
The Planck HFI 353 GHz maps have a native resolution[https://wiki.cosmos.esa.int/planckpla2015/index.php/HFI_ 
 performance_summary] of about 4.94' <cit.> and are given in a HEALPix[See http://healpix.sourceforge.net]
grid tessellation corresponding to N_side = 2048 <cit.>.
At the instrument resolution, the 353-GHz polarization Stokes Q and U maps are noise-dominated <cit.>.
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 <cit.>. Thus, its impact on our analysis is expected to be negligible. The cosmic infrared background (CIB) is considered to be unpolarized <cit.>.
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., ). We note that we  carry out our work at low resolution so that contributions from point sources and the CIB are not expected to be significant. The CMB temperature anisotropies and zodiacal light contributions are expected to be subdominant with respect to the thermal dust emission in intensity. To consolidate the results of our analysis for the intensity part, we also make use of the full-sky dust extinction map derived in <cit.>, 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  <cit.>.




 §.§ 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 <ref>.
Using these parametric models and realistic Planck noise maps, we produced different sets of high-resolution maps that simulate the Planck 353-GHz linear polarization Stokes parameter maps. These sets, hereafter called ,  and  are defined as follows.

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

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


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.[Using the  software this realization takes up to 300 gigabyte of the RAM before the line-of-sight integration performed.]
We then upgraded the maps at N_side = 2048 to directly calibrate the Galactic thermal dust emission of our simulated maps to the measured one at 353 GHz by Planck through a linear fit.
To add realistic Planck noise, we consider the one derived from the difference between the Stokes parameter maps obtained from the data corresponding to the first and the second half mission of the Planck individual pointings. 


In summary, the simulated Planck maps, M, of the polarized thermal dust emission at 353-GHz are given by:

    M_X =  A_cal S_X^dust + N_X^Planck ,

where X refers either to I_d, Q_d, U_d, N is the noise map, S is computed using  to estimate Eqs. <ref> and <ref>, 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. <ref>.




§ FITTING PROCEDURE IN INTENSITY AND POLARIZATION FOR THE GALACTIC THERMAL DUST EMISSION


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




 §.§ Choice of fitted-map resolution


Although we have optimized our codes, producing simulated maps of the diffuse Galactic emission at the instrument resolution is extremely time consuming.
To effectively adjust models to data and efficiently explore the different parameter spaces, a large number of simulations is needed.
Inevitably we thus have to work at a lower resolution than the native one of the data.

In Appendix <ref> we explore in detail the consequences of this requirement both at the map level and in term of the model-parameter posterior distribution reconstruction. We demonstrate that the necessity to work at lower resolution induces a bias in the parameter space. However, this bias is currently, and for many practical cases, not the dominant source of uncertainty in the kind of modeling tackled in this study. Nonetheless, it is best to ensure that this source of uncertainty is sub-dominant by working at the highest resolution made possible by the computing facilities.

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



 §.§ Likelihood definition


The best-fit parameters are obtained by maximizing a Gaussian likelihood function ℒ, defined as

    lnℒ  = -1/2 ( D - α M )   C^-1_D ( D - α M )^T
    
    ,

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:
α = (∑_i (D_i   M_i/σ^2_i))/(∑_i (M_i^2/σ^2_i)).




 §.§ Estimation of the covariance matrix


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



In the case of the instrumental 
uncertainties in the Stokes parameters, {I,Q,U }, we consider the block diagonal per-pixel covariance matrix maps released by the Planck collaboration at N_side = 2048. We neglect the off-diagonal terms[To cross-check this working simplification, we evaluated the full noise covariance matrix for the reduced Stokes parameters q and u, relying on MC simulations. We find that they distribute around zero and that for 90 percent of the pixels we have |C_qu|/ C_qq≤ 16.4% and |C_qu|/ C_uu≤ 16.3%.]
and propagate the uncertainties at low resolution according to

    σ^2_D_64;stat,i = 1/N_[2048,64]^2∑_j ∈ iσ^2_D_2048;stat,j ,

where σ^2_D_N_side;stat,i 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:

    σ^2_D_64;disp,i = 1/N_[2048,64]^2∑_j ∈ i( D_2048;j  -  D_64;i)^2
    ,

where D={I,Q,U}[Notice that, we do not consider the mixing of Q and U due to parallel transport when averaging over pixels.].
Let us note that in the synchrotron literature (e.g. ), the signal dispersion is often taken as simply the dispersion of the data in low-resolution pixels.
The two estimates of the variance due to signal dispersion differ by a factor of N_[2048,64] = 1024. We adopt the normalized dispersion, also referred to as the standard error on the mean, in order to obtained a reduced χ^2 on the order of one when the data are fitted with the appropriate models (see Table <ref>) 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, which
would account for the intrinsic signal dispersion and the noise, as:

    σ_D_64 = max(σ_stat, σ_disp)  .


These pseudo uncertainties are then propagated to the reduced Stokes parameters, q and u as:

    σ^2_x = 1/I^2σ^2_X + X^2/I^4σ^2_I
    ,

where x = X/I and X={Q,U}. And so, we write the covariance matrix for { q,  u } as C_xx' = diag{σ^2_x}.






 §.§ Implementation of the MCMC experiment

In order to explore the parameter space, determine the posterior distributions of the parameters, and find their best-fit values,
we used the  MCMC Python software written by <cit.>, who implemented the Affine-Invariant sampler proposed by <cit.>.
We run the MCMC code for several Markov chains until the convergence criteria proposed by <cit.> 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 hyper-surface may exhibit multiple minima that might be sharp in the case of thermal dust emission (mainly for dust intensity).
To remedy this problem, we stopped the MCMC after several thousand steps and relaunched the experiment from the location of the minimum χ^2 and waited for the convergence criteria to be fulfilled to reach well-sampled posterior distributions.

We used 250 Markov walkers for the intensity fits and 100 for the polarization fits. The Markov chains are initialized according to uniform distributions.
For the exploration of the parameter space, we considered non-informative priors, adopting top-hat distributions. The explored ranges of values of the free parameters are given in Table <ref> 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.









 §.§ 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. <ref> (specifically  and ).
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. <ref> 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.




§ RECONSTRUCTION ACCURACY FOR THE DUST DENSITY AND REGULAR GMF COMPONENT


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

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


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


The best-fit parameter values that we obtained are reported in Table <ref>, both for the intensity and the polarization fits.
The obtained values of the reduced χ^2 for each case are reported in Table <ref>.





 §.§ Reconstruction of the dust density

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



The main results for the intensity fit for case A are shown in the first
column of Fig. <ref>.
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 <ref>.
It is indeed close to the Galactic plane that the functional forms of the dust density are allowed to vary quickly. Therefore, it is from those locations (then projected on the sky) that the loss in space-sampling resolution has larger effects.



For cases B and C, the intensity fit results are shown in Fig. <ref>.
The best-fit model and residuals are shown in rows 2 (3)  and 4 (5) for case B (C), respectively.
For case B, the fitting to the intensity data is very poor with a reduced χ^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  simulation.
The uncertainties used to computed the reduced χ^2 do not account for mis-modeling.
The significance of the residuals are large both in and
outside the Galactic equatorial region. 

For case C, we obtain a good fit to the data with reduced χ^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 line-of-sight integration, and assuming we have the right model at hand, we are able to retrieve the parameters of a complex dust density distribution with reasonable confidence. In contrast, as shown in case B, if the applied dust density model does not reflect reality, large residuals can be found. Although we cannot directly extrapolate the results of case B, it allows us to anticipate on the kind of issues we would face with real data sets, as discussed in Sect. <ref>.



 §.§ Reconstruction of the GMF



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

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



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



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



The reconstructed geometrical structures of the GMF corresponding to the cases A, B, and C best-fit parameters
are shown in Fig. <ref>, and can be compared to the input one shown in Fig. <ref>. 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.











 §.§ 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. <ref>) to the input
model shown in Fig. <ref>.
At each location of the Galactic space considered when computing
the polarization observables (see Sect. <ref>), 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), we
define 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 𝐞_ϕ, 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 the GMF
line with 𝐞_z and characterizes the out-of-plane GMF
component.
In the second scheme, α is the inclination angle that makes the
field line relative to the line of sight and γ the position angle in
the plane orthogonal to it.
Let ξ̂_r 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 Δ_2(ξ̂, ξ̂_r)
as the difference of the two angles that range between -90 and 90
degrees (-90 and 90 corresponding to the same configuration).




To quantify the accuracy of the reconstructed GMF structures, we compute at each location the Δ_2(ξ̂, ξ̂_r) for the angles between the field lines of a best-fit GMF model and the field lines of the input GMF model. We generate histograms of the angle differences.
An hypothetical perfect reconstruction would lead to a single bin centered in zero. The histograms for the three reconstruction cases are presented in Fig. <ref> for the four angles discussed above. Cases A, B, and C are shown on the same plots for comparison.


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



We can also compare the best-fit models to the input one via the differences of the polarization position angles in the maps. For this, we use the acute angle, defined between 0 and 90 degree as

    Δ_ψ̂_r, ψ̂ = 90^∘ - | 90^∘ - |ψ̂ - ψ̂_r | |,

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. <ref> 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).








§ 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.



 §.§ Modeling of the turbulent component

We consider an overall turbulent contribution so that the total 3D GMF can be written as <cit.>:

    𝐁^tot(𝐫) =  𝐁^reg(𝐫) + A_turb𝐁^turb(𝐫),

where the regular (𝐁^reg) and turbulent (𝐁^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 𝐁^reg is assumed to be one of the regular models depicted in Appendix <ref> and 𝐁^turb is assumed to be a Gaussian realization of a random isotropic 3D vector field that follows a two-slopes power-law spectrum as described in detailed in Appendix <ref>.

Using polarization data only, <cit.> 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., ;  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. <ref>), 0.3 (comparable to the upper limit of A_turb from <cit.>), 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. <ref>. Then, the thermal dust emission I, Q and U Stokes parameters are computed by integrating Eq. <ref> considering a thermal dust density distribution model and the total GMF.
For illustration, we show in Fig. <ref> 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  simulation regular GMF component. We represent the signal-only (top) and signal plus noise maps (bottom) for A_turb = 0.0 (top), 0.3 (middle), and 0.9 (bottom). 








 §.§ Quality of the reconstruction of the regular GMF in the presence of turbulence


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


  * Case I:  n_d≡ ARM4 and GMF≡ LSA,

  * Case II:  n_d≡ ED and GMF≡ LSA,

  * Case III:  n_d≡ ED and GMF≡ ASS,

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



In Fig. <ref>, we show the 1D and
2D marginalized posterior distribution obtained while fitting the  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. <ref>), 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 values are only of few degrees for all cases.
Inspection of Fig. <ref> also reveals that the biases in recovered parameter values induced by the non-consideration of the turbulence is roughly of the same order of the one induced by the mismodeling of the density distribution or of the mismodeling of the regular part of the GMF. In this paper, we concentrate in the latter effect and do not attempt to fit for the turbulence component in the likelihood.


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



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










§ PROOF OF CONCEPT FOR FUTURE SYNCHROTRON ANALYSIS



As discussed in Sect. <ref> 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  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 Eqs. <ref>.
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 <ref>, and in the second, the best-fit parameters from the fit of the thermal dust emission discussed in Sect. <ref>.
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. <ref>.



Using , we computed the intensity and polarization maps corresponding to the Galactic synchrotron emission in the two cases presented above.
We show them in Fig. <ref>. 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. <ref> where we show the histogram of the acute angles (Eq. <ref>) 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. <ref>. Full validation and application of this methodology will be presented elsewhere.








§ RECONSTRUCTION OF THE GMF STRUCTURE FROM THE PLANCK 353 GHZ DATA



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

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



 §.§ Dust density reconstruction






  §.§.§ Fitting procedure

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



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



  §.§.§ Best-fit models

The best-fit parameters are reported in Table <ref>.
In Fig. <ref>, we compare the data and the best-fit 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 best-fit models cannot account for all the complexity and the richness contained in the data. This is not surprising since our models do not include patterns such as clumps, filaments, or bones that can be directly spotted by eye and that require dedicated extra modeling (see e.g. ). Accounting for such features would significantly increase the number of parameters and will not be easily tractable within a MCMC approach.

In Fig. <ref>, we show the histograms of the χ_i corresponding to the three best-fits to the intensity map. The distributions are nearly Gaussian with very large spread. This is again due  to the fact that the uncertainties in the data (even including the intrinsic dispersion) cannot overcome the residuals to such simplistic parametric models  as discussed before.
As a result, any small-scale feature turns out to be significant and the values of the reduced χ^2 are ≫1 (see Table <ref>). 
This was to be expected from our simulation-based 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. <ref>).







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



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





 §.§ 3D regular GMF reconstruction











  §.§.§ Fitting procedure

Similar to what we did to fit the intensity map, we adopted flat priors on the free-parameters of the different models of the GMF, which are presented in Table <ref>. Details on the parametrization of these GMF models are given in Appendix <ref>. 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 <ref>.
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.




  §.§.§ Best-fit models

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



The parameter values of all the best-fits of the different GMF models considered are reported in Table <ref> for each of the n_d models that we investigated in Sect. <ref>.

In Table <ref>, we also present the values of the reduced χ^2 that quantify the goodness of our fits to the GMF models.
In Fig. <ref>, 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 best-fit models are one order of magnitude better than those obtained for the fits of the intensity map.



Furthermore, the obtained best-fits of the GMF are generally not compatible within the uncertainties when moving from one dust density distribution 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. <ref>),
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 best-fit parameters. This again stresses the need for more complex and realistic models of the dust density distribution and of the GMF in order to optimally exploit the currently available data.



To test the robustness of our best-fit models, we also constrain the GMF parameters using the n_d best-fit models obtained from the fits on the dust optical-depth map.
The fitted (q, u) maps are obtained simply by replacement of the intensity map by the τ_353 map, with the corresponding uncertainties. The best-fit polarization maps and the geometrical structure of the GMF are found to be in good agreement with those obtained when the modeling of n_d is from the I_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 <cit.>. We proceed to the fit of the GMF model for the least evolved and the most evolved modelings (ED + ASS and ARM4⊕ED + QSS) with and without applying the quoted corrections.
The best-fit parameter values that we obtain are in agreement at the ten percent level with the previous values. These two tests make us confident that our GMF reconstruction with the Planck data is not strongly biased by instrumental and map-making systematics and shows that those systematics in Planck polarization maps have a significant contribution to the modeling uncertainty budget whereas not accounted for in the covariance matrix. This conclusion was already reached in <cit.> and <cit.>.




  §.§.§ Constraints on the GMF geometry

The radial and height evolution of the pitch and the tilt angles that correspond to the best-fit parameters for all the  models considered are shown on Fig. <ref>.



For a given GMF model (same line style but different color in Fig. <ref>), 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 of uncertainty is of the same order as the one attributable to the Planck residual systematics.
As inferred from Fig. <ref>, 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. <ref> we show the direction of the large-scale GMF lines in horizontal and vertical planes of the Galaxy, namely the (x, y, z=0) and the (x, y=0, z) planes. Stream plots of the field lines are shown for the best fits obtained for each GMF models assuming the n_d to follow the ED model.
Very similar plots are obtained with the other n_d models, as we discuss later.
In terms of the visuals and given that the functional forms of the models may be different,
the similarity among the GMF reconstructions is remarkable.
We note that for the LSA model, a second minimum of the χ^2 is found with a distant solution in the parameter space but that is similar to the solution obtained by <cit.> in fitting Faraday-depth all-sky
data and synchrotron data.
We present this solution in Appendix <ref>. 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 Eqs. <ref>) to span through negative values.
An interesting feature of our best-fit LSA model is that the field lines appear to wind up on themselves at a cylindrical radius of about 18 kpc.




 §.§ Best-fit maps: Comparison and characteristics

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

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




 §.§ Comparison of reconstructed GMF


The fact that the best-fits GMF are very consistent with respect to the choice of the dust density distribution model can be observed, for example, in Fig. <ref>.
In this figure we show the geometrical structure of the GMF in the horizontal and vertical cross-cuts of the Galaxy obtained with the three best-fit models of n_d when inferring the GMF 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. <ref>.
For each best-fit model, we compute the pitch (p), tilt (χ), inclination (α) and position (γ) angles at each position of the Galactic space used for our MC realizations.
We then compute the signed differences of each of the four angles considering the best-fit from the least evolved modeling (i.e. n_d≡ ED + GMF≡ ASS) as a reference to which each best-fit GMF model is compared with.
The histograms of the signed differences are presented in Fig. <ref>.




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


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

  * 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 our reference model.

  * 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.

  * Comparison of the best-fits of the same GMF model but for different n_d reveals a sizable coherence. We have checked that this conclusion, already informed around Fig. <ref>, 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. <ref> for the ASS model is quantified by the blue histograms in the middle and right panels of the first row of Fig. <ref>, except that, here, we are not limited to (x, y, z=0)
or (x, y=0, z) planes but account for the whole sampled space.
The histograms show that through the whole sampled space the differences of the measured pitch angle are below one degree.
The blue histograms in the middle and right panels of the second row quantify the small differences in the out-of-plane components among the ASS best-fits observed in lower panels of Fig. <ref>. The measured tilt angle differences (compared to the reference GMF) are less than about two degrees throughout the whole space.




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

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



In the above, it is the comparison of model-parameter values (and corresponding geometrical features) that is relevant, not the absolute values 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. <ref>).
However, we have shown that reliable constraints can be obtained for the pitch angle 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 <cit.>. 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 <cit.>.
The out-of-plane pattern of the GMF is poorly constrained in our analysis. This is likely because it depends more on the particular GMF model under study and, at the same time, because its influence on maps is more noticeable at low and intermediate Galactic latitudes where our fits are not satisfying.
Let us emphasize, though, that none of our best-fits exhibits the expected X-shape of the GMF lines that has been found in radio observations of external galaxies (see e.g.,  for a discussion).




 §.§ Alternative comparison of polarization maps




The degree of linear polarization (p_lin = √(q^2 + u^2)) and the polarization position angle (ψ = 1/2 arctan2(-u,q))[
The polarization position angles are given in the IAU convention.]
are commonly used to study and characterize the magnetized interstellar medium <cit.>.
In Fig. <ref> we show the noise debiased[Noise debiasing procedure is performed in the same way as in Appendix B.2 of <cit.>.] p_lin and ψ maps from the Planck data at 353 GHz (top row) and those deduced from the polarization maps (q and u) corresponding to our best-fits GMF reconstructions obtained with the n_d≡ ARM4 model fitted to I_353.
These maps are alternative ways at looking at the q and u maps presented in Fig. <ref> (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. <ref> clearly illustrates qualitatively the limitation of the models at reproducing the data.
Despite the noise dominated map of the data, it is nevertheless interesting to note how the modeling can qualitatively reproduce low p_lin values in some region of the Galactic equator where we know the intensity to be high. These low degrees of linear polarization results from line-of-sight depolarization induced by the varying GMF orientation.
This can be inferred from the comparison, at low |b_gal|, of the BSS and QSS models (two last rows) to the ASS model (second row), for instance. Let us stress that this depolarization is obtained using the regular GMF only.


Furthermore, for lines of sight with p_lin≲ 5%, the relatively simple models of the large-scale regular GMF that we consider are able to reproduce qualitatively the low and high values of the degree of linear polarization for low and high values of the intensity as suggested by the data <cit.>.
However, these models fail at reproducing very high degree of linear polarization (p_lin > 5%), where large difference are observed between model maps and data. These difference are well localized in the sky towards regions of nearby known structures (spurs, arms, and clouds) mainly attributable to the nearby structures such as the Sco-Cen association (e.g., <cit.>).
We used the  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. <ref>).



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. <ref>). Indeed, the agreement between the models and the data is fair in terms of ψ, namely, in terms of the projected (and weighted) orientation of the GMF lines integrated along the lines of sight.
Beside the fair resemblance, the modeled ψ maps show more coherence and smooth spatial variations than observed in the data. Also, inherently to the large-scale regular GMF considered in this study, the modeled maps have a mirror symmetry about the Galactic equator that is accompanied
by a change of sign. In the data, that symmetry is disturbed on small scale and twisted and skewed on global scale. This is well illustrated in the upper right quadrant of the maps and comparing the `lines' for which ψ = 0^∘ = 180^∘, respectively.
The inclusion of magnetic fields attached to local structures, such as depicted by <cit.> and <cit.> in the case of the Local Bubble, addresses the latter issue in breaking the North-South symmetry.
The modeling of (relatively) small scale perturbations would require the modeling of magnetized structures well localized in the Galaxy and will necessitate the modeling of the turbulent component, perhaps in the form of parameterized magneto-hydrodynamic turbulence, such as that depicted in <cit.>.
Any such investigation is beyond the scope of this paper and will be addressed somewhere else.




 §.§ 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., ; );
(ii) the interconnection of the latter with the large-scale regular GMF; and (iii) the generalization of the large-scale regular GMF models.
As a straightforward generalized model, we would propose a model that contains both the radial variation of the pitch angle, such as encoded in the LSA model, and the strength field modulation, such as in the BSS or QSS model. Independently, both features appear indeed to be tentatively favored by dust data, as compared to the ASS model.
The consideration of the other observables such as synchrotron emission and Faraday rotation measurement data promise to help significantly at constraining and modeling further the large-scale regular part of GMF and to solve possible ambiguities. In this respect, it might be worth undertaking physical models of the GMF such as those derived by <cit.>.

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 ).

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 ), 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., ; ; ), especially at the physical scales probed by the Galactic diffuse emission considered at low resolution.

At the map level, the inclusion of a turbulent component also results in the two following effects:
(i) an angular spatial de-coherence of the polarization position angles and
(ii) an overall lowering of the degree of linear polarization, due to line of sight depolarization. The latter would call for an overall increase of the global amplitude of the modeled polarization maps and therefore would tend to enlarge the range of values span by Q and U. Such an increase might help at recovering the missing amplitude of the field strength, that can be visually inferred from comparing the data with the modeled polarization maps (see Fig. <ref> or Fig. <ref>).
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.




§ SUMMARY & CONCLUSIONS




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

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



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

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


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


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





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



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



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.








aa 
 





  


§ 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.



 §.§ Dust density distribution model



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



  §.§.§ Exponential Disk (ED)


According to this model, the dust density distribution takes the parametric form:

    n_d(ρ,  z) = A_0  exp(- ρ / ρ_0)/(cosh(z / z_0))^2 .

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




  §.§.§ 4 Spiral Arms (ARM4)

In this case, the dust density distribution is obtained by summing the contributions of four logarithmic spiral arms. In this implementation, and up to a relative amplitude, the arms are assumed to be identical with a rotation of ninety degrees.
Our implementation is given as:

    n_d(ρ, ϕ,  z) = A_0  exp(- (ρ - ρ_c)^2 / (2 ρ_0^2))/(cosh(z / z_0))^2 𝒮(ρ, ϕ)

where the function 𝒮 encodes the logarithmic spiral pattern as:

    𝒮(ρ, ϕ) = ∑_i ( A_i  exp( - (ϕ - ϕ_s,i )^2/2 ϕ_0^2) )

and

    ϕ_s,i =    1/tan(p)log (ρ / ρ_0,i)         , 
    ρ_0,i =    exp( ϕ_0,i tan(p) )         , 
    ϕ_0,i =    ϕ_00 + i  π/2             ,

with i = { 1, 2, 3, 4 }, ϕ_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.



  §.§.§ 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.




 §.§ Regular Galactic Magnetic Field model


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


  * field lines have spiral geometries in planes parallel to the Galactic disk (at z constant)
 
  * field lines have out-of-the-plane component, intended to produce the so-called X-shape of the GMF

  * the field lines have some degree of symmetry with respect to the z-axis of cylindrical Galactic coordinate system
The first three models have already been discussed in the literature. We briefly present the models below but refer the reader to e.g. <cit.> for detailed discussion.



The four models take simple analytic expressions in the cylindrical coordinate
system center on the Galactic center:

    𝐁 = B_ρ𝐞_ρ + B_ϕ𝐞_ϕ +B_z 𝐞_z,

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



  §.§.§ Axi Symmetric Spiral model: ASS

The ASS model (see e.g., ; ) is one of the simplest descriptions of the Galactic magnetic field.
It is compatible with a non-primordial origin of the Galactic magnetism, based on the dynamo theory. The field lines follow logarithmic spirals with constant pitch angle. The cylindrical components of this model are:

    B_ρ  =           B_0 sin(p) cos(χ(z))                
    
    B_ϕ  =           B_0 cos(p) cos(χ(z))                
    
    B_z             =           B_0 sin(χ(z)),

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 non-zero 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, 𝐞_ϕ, 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:

    χ(z) =     χ_0 tanh( z/z_0) .

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 free-parameters (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.





  §.§.§ LSA model

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

The parametric form of the field components read

    B_ρ =        B_0 sin( ψ(ρ)) cos(χ(z))     
    
    B_ϕ =        B_0 cos(ψ(ρ)) cos(χ(z))      
    
    B_z     =      B_0 sin(χ(z)),

where the function 

    ψ(ρ) = ψ_0 + ψ_1  ln( ρ/R_⊙)

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. <ref>  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.




  §.§.§ Bi Symmetric Spiral model: BSS

The class of bi-symmetric spiral models produces GMF that could be compatible with a primordial origin. This model includes magnetic field line reversion as suggested from pulsar rotation measurements <cit.>.
This model was also referred to as the Modified Logarithmic Spiral model in
<cit.>.
The field lines and the spirals drawn by the surfaces of iso-amplitude of
the field have a constant and same pitch angle.
The cylindrical components of this field model read

    B_ρ  =      B_0 cos(
                                    ϕ±βln(ρ/ρ_0)) sin(p) cos(χ(z))    
    
    B_ϕ =      B_0  cos(
                                    ϕ±βln(ρ/ρ_0)) cos(p) cos(χ(z))    
    
    B_z =      B_0 sin(χ(z)),

where β = 1/tan(p), with p the pitch angle, χ(z) the tilt angle defined as in Eq. <ref>, 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 cos(ϕ±βln(ρ/ρ_0)). 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 <cit.> 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.




  §.§.§ Quadri Symmetric Spiral model: QSS

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



 §.§ Turbulent Galactic Magnetic Field model


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


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

Following <cit.>, 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 and 3.9 kpc, respectively. In the Fourier space the high-resolution and low-resolution simulations span the ranges of k-values of [ 0.25,  256 ] and [0.025,  25.6] 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. <ref>), assuming a nearest neighbor scheme to map the two space samplings.

The high- and low-resolution simulations are normalized in the overlapping range of k values (i.e., k ∈[0.25,  25.6] 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 large-scale fluctuations in the neighborhood of the Sun. Outside the heliocentric radius of 2 kpc, only the low resolution box is used, we add random white noise with an rms that corresponds to the spatial frequencies not sampled in the low-resolution box but that are included in the high-resolution one (namely, for k ∈[25.6,  256 ]). This procedure ensures an identical rms of the turbulent GMF vectors throughout the sampled space. Unlike <cit.>, we do not re-scale the turbulent GMF component as a function of position in the Galaxy. This is because the regular GMF model of <cit.>, the only one to which the turbulent component is added in this study, has a constant amplitude.




§ 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. 

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


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


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

The results of this analysis are shown in Fig. <ref>. It is seen that the best-fit parameters converge toward the input values as the resolution of the model used in the MCMC analysis increases. The convergence is however not trivial and depends upon the considered parameter and also on its input values.
To this respect, we note that for some parameters, the best-fit values may oscillate around the input across the N_side values.
The value of the reduced χ^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 best-fit parameters for dust density distribution and the GMF, these models need to be taken with caution.
Furthermore, the best-fit parameters obtained at a given resolution should not be used at other resolutions.
Formally speaking, it should be possible to account for such a bias in the definition of the likelihood function. However, we have observed that the bias depends strongly both on the parametric form of the models and on the parameter values themselves.
Therefore, accounting for it is cumbersome and 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. <ref>).
Therefore, for now, we notice that this bias exists and that the larger the N_side value of the working resolution, the smaller the bias is.




Third, we tested the impact of the adopted value for the radial sampling (sampling along the line-of-sight). As far as we do not pretend to model the very nearby structures, we found that the radial step can be as large as few hundreds of parsecs.
The difference in pixel space is more pronounced due to the adopted angular sampling than due to the radial sampling. However, a fine radial sampling 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.







§ ORTHOGRAPHIC VIEW OF THE PLANCK DATA AND OF THE BEST-FIT MODELS





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

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




§ SECOND SOLUTION FOR THE LSA GMF MODEL




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



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


