© ESO, 2010

1. Introduction

Spiral galaxies exhibit large-scale magnetic fields. The Milky Way is not an exception, but obtaining its spatial distribution is extremely difficult. Most methods for observing magnetic fields are based either on radio observations of the synchrotron emission (Wolleben et al. 2006; Reich 2006; Testori et al. 2008, and references therein), or on the Faraday rotation (hereafter FR) of pulsars (e.g. Weisberg et al. 2004; Han et al. 2006; Noutsos et al. 2008) and extragalactic radio sources (hereafter EGRS) (e.g. Gaensler et al. 2001; Brown et al. 2007; Haverkorn et al. 2008a; Carretti et al. 2008). In addition, several radio lines show Zeeman splitting, which can be used also to directly constrain the strength of the magnetic field (e.g. Fish et al. 2003; Han & Zhang 2007).

Despite the large numbers of pulsars and EGRS for which the rotation measure (hereafter RM) has been recently determined, there is no consensus about the large-scale pattern of the Galactic magnetic field (hereafter GMF). It is probably more complex than previously expected, as pointed out in Men et al. (2008). Recent results by Sun et al. (2008) show an axisymmetric disk distribution with reversals inside the solar circle as the best model to describe the GMF when all-sky maps at 1.4 GHz from DRAO and Villa Elisa, the 22 GHz map from WMAP satellite, and the Effelsberg RM survey of EGRS are combined. Results derived from Brown et al. (2007) by using RM of EGRS suggest an axisymmetric pattern of the disk magnetic field. Vallée (2008) proposed the inclusion of a ring model to describe the field. The interest on the structure of this large scale GMF is justified by a number of reasons. First of all, this field might be important when considering the dynamics of the galaxy at those large scales (Nelson 1988; Battaner et al. 1992; Battaner & Florido 1995; Kutschera & Jalocha 2004; Battaner & Florido 2007). A good characterization of the large-scale GMF pattern would allow a detailed correction of the galactic contribution for a better understanding of cosmological magnetic fields (see Battaner & Florido 2009, for a recent review), which would be potentially observed with upcoming CMB missions such as PLANCK (The Planck Collaboration 2006).

In addition, the GMF modifies the trajectory of high energy cosmic rays, therefore its knowledge is crucial for understanding their distribution in energy and direction, as it is obtained in experiments like Auger (Bluemer et al. 2008), Milagro (Abdo et al. 2009) and others. The global anisotropies found by Milagro have been interpreted by Battaner et al. (2009) as produced by Galactic magnetic fields. Moreover, GMF could be important to explain the well known knee in the spectrum of cosmic rays energies around 10⁶ GeV (Masip & Mastromatteo 2008), because sub-knee-energy cosmic protons trapped by GMF should have a much higher optical depth for interactions with WIMPs.

The extraction of structure of the Galactic magnetic field from the measurements of the polarized intensity is extremely difficult. The Galactic magnetic field probably cannot be described by a single model because there would need to be at least three components, the thin disk, the thick disk (Beuermann et al. 1985) and the halo, each one defined by its own model and parameters. Taking into account that each component is largely unknown, an analysis of a multi-component magnetic field renders an objective estimation of the best configuration a task that is too complicated.

The thin disk is characterized by the highest field strengths, and we are embedded in it. However, magnetic fields in the thin disk are predominantly random, and therefore they barely contribute to the observed net polarized emission, even if a z-component of the regular fields exists (Han et al. 1999). On the other hand, local spurs (Berkhuijsen et al. 1971) like the North Polar Spur, highly distort the main field configuration. Even at this high frequency Faraday depolarization cannot be neglected in particular regions. Observing at high galactic latitudes the thin disk field is contaminated by the thick disk and the halo fields, which could even become dominant. Models for the thick disk are scarce in the literature, but those proposed for the thin disk could be tested, even if different parameters characterized both disks. Our results for the disk inferred from the polarized synchrotron emission will give a complementary insight into the field structure. The halo structure remains unknown and has a very different structure, consisting in a double torus in two hemispheres with opposite directions. After some pioneer detections (Simard-Normandin & Kronberg 1980; Han & Qiao 1994), it has been modeled by Han et al. (1997), Harari et al. (1999), Tinyakov & Tkachev (2002), Prouza & Šmída (2003), and others. To illustrate the current uncertainty on this component, we can compare the maximum value of the magnetic strength, 1 μG following Prouza & Šmída (2003), and 10 μG following Sun et al. (2008), although these last authors propose a maximum strength of only 2 μG when the thermal electron scale height is increased by a factor of 2. The contribution of the halo field to the polarized emission is therefore difficult to estimate. Most models take only RM from extragalactic sources into account to estimate the halo structure.

We here carry out a systematic comparison of a number of possible GMF models, exploring which one is providing a better fit to the large-scale polarization map at 22 GHz.

Our analysis is based on the 5-year WMAP data release, and extends the work by Page et al. (2007) by considering a detailed comparison not only with the polarization angle, but with the polarized intensity (i.e. the Stokes Q and U parameters). Although the polarization angle can be used to described some properties of the large scale pattern of the GMF, it is not sensitive to some parameters (e.g. the field strength) and it also contains some intrinsic degeneracy with respect to the direction of the field lines. Because of this reason, our main results will be obtained with the analysis of the Stokes Q and U maps, although the independent analysis based on the position angle will be done in some cases for comparison.

We used eight disk models and one halo model. Different masks enabled us in an indirect way to estimate the different contributions of the galactic components at different galactic latitudes, but the consideration of a multicomponent field has not been fully undertaken. Below, we will not make an explicit distictintion between thin and thick disk.

For obtaining the emission we also need a model of the distribution and spectrum of cosmic rays, which is another important source of uncertainties. Here we have assumed that the cosmic ray structure follows that of the gas, as they are produced by supernovae and cannot travel far away from the birth place. This assumption is rather usual in the literature but also rather questionable.

Here we focus on the study of the large-scale field. In principle, we could divide the galactic magnetic field into three components: a random component for scales lower than about 100 pc (see Haverkorn et al. 2008a), for which some works have published the turbulence spectrum (Han et al. 2004; Han 2009); a main “spiral” field for scales typical of spiral arms, and a “galactic scale” main field. Some authors consider that 1 kpc is a length defining the large scale (e.g. Han 2008a), therefore taking spiral arms as a large scale phenomenon. Indeed, the field direction is opposite in arms and in inter-arms regions (e.g. Beck et al. 1996; Han et al. 2006). Here, however, we are interested in scales of the galaxy itself, which is why the spiral arms are considered as wavy perturbations. It is clear that we need to investigate these three components, but this separation is important because the tools, mainly the interpretation in terms of the generation mechanism, could be completely different. The random field should be produced by turbulence in a magnetized medium, supernova explosions and other local mechanisms. Spiral arms produce characteristic motions that enhance magnetic fields in a high conductivity medium. Fields at the galactic scale would be interpreted in terms of galaxy formation and/or dynamo effects.

The polarized synchrotron emission is a valuable tool for investigating the overall field pattern at large scales and specially at high galactic latitudes. The total emission is much affected by random fields and by non-polarized galactic emission as free-free (see e.g. Miville-Deschênes et al. 2008). Faraday rotation of pulsars is a powerful technique but it is very much affected by enhancements and tangling of the field by the passage of spiral waves. Faraday rotation of extragalactic sources would inform about the larger scale fields, but here the poor knowledge of the intrinsic FR is a heavy problem. When using all-sky observations such as provided by WMAP, and PLANCK in the future, the observed emission is integrated along a path traversing the galaxy, therefore the wavy effect of spiral arms is at least in part smoothed.

This could explain why RM of pulsars and WMAP polarization have given somewhat inconsistent results. There are other all-sky surveys at lower frequencies (see Reich 2006), but Faraday depolarization and dust emission are very strong at latitudes below 30deg. The detailed study of the WMAP polarization measurements in terms of the Galactic magnetic fields is more a complementary than an additional technique. It is relatively free of Faraday depolarization, of dust contamination, of random fields and of spiral waves perturbations.

To complement our study, for some of the models we also investigate if a modification of the radial variation of the strength of the GMF is producing an impact on the quality of the fit. This radial variation is also largely unknown and may be important for the production of the all-sky polarization maps.

The paper is organized as follows. In Sect. 2 we describe the set of GMF models that have been investigated in this analysis. In Sect. 3, we describe the numerical method used to produce the template maps, as well as the the model selection method used to determine our preferred model and the corresponding set of parameters. We present the results in Sect. 4 and discuss them in the following section.

2. Models of Galactic magnetic field

We have focused our analyses on the constraints derived from the polarized intensity map at 22 GHz. At these frequencies, the physical process which dominates the polarized intensity is the synchrotron radiation emitted by the population of relativistic cosmic ray (hereafter CR) electrons with energies between 400 MeV and 25 GeV (Strong et al. 2007), as they interact with the GMF.

Therefore, in order to obtain a prediction of the polarization pattern at this frequency, we first require a description of the distribution of the relativistic CR electrons in our Galaxy. For the purposes of this paper, in which we are interested only in the large-scale pattern of the GMF, it is sufficient to consider a simplified description of the CR electron population. Here, we use the spatial distribution of relativistic electrons as in Drimmel & Spergel (2001), i.e. ${\mathit{N}}_{\mathrm{e}}\mathrm{=}{\mathit{N}}_{\mathrm{0}}\exp \left(\frac{\mathrm{-}\mathit{r}}{\mathrm{5} \mathrm{kpc}}\right){\mathrm{sech}}^{\mathrm{2}}\left(\frac{\mathit{z}}{\mathrm{1} \mathrm{kpc}}\right)\mathit{,}$(1)where N₀ ≈ 3.2 × 10^-4   cm^-3 is derived from the value for the CR electron density on Earth (Sun et al. 2008), and r and z are the radial and the vertical coordinates in cylindrical galactocentric coordinates, respectively. This model essentially assumes the same spatial density distribution for the CR electrons as for the interstellar gas. At first order, this is what we would expect, because the higher the interstellar gas density, the higher the star-formation rate, and therefore the higher the supernova production, which gives a higher relativistic electron density. These cosmic electrons lose energy (via synchrotron itself) on short distances of less than 1 kpc.

The value for N₀ in Eq. (1) is very uncertain. This value is usually obtained by assuming that the CR electron spectrum can be described by a power-law with constant spectral index p = 3. However, observations in the last few years as well as numerical simulations (see Strong et al. 2007, and references therein) suggest that this assumption is not appropriate for the entire spectrum. Moreover, different observations show variations of the order of 50% (or even higher) for this quantity. Therefore, we expect that this uncertainty might introduce a bias in the recovered amplitude of the field strength for the different models, and will be taken into account as explained below.

In the following sub-sections we present the set of GMF models we investigate in this work, all of them taken from the literature. Most of these models have been proposed/constrained using the analysis of Faraday rotation of pulsars and EGRSs. Thus, it is also interesting to explore if these models are also able to describe the large-scale polarization pattern seen in WMAP 22 GHz maps. In total we investigate eight models describing the disk and halo fields, and one model conceived to describe the halo field. The disk models are: 1) axisymmetric (ASS), 2) axisymmetric with radial dependence of the strength (ASS(r)), 3) bisymmetric positive (BSS₊), 4) bisymmetric negative (BSS₋), 5) bisymmetric positive with radial dependence of the strength (BSS₊(r)), 6) bisymmetric negative with radial dependence of the strength (BSS₋(r)), 7) concentric circular rings (CCR) and 8) logarithmic spiral arms (LSA). Moreover we shall consider for the halo component the bi-toroidal (BT) model. Throughout this section, all coordinates (r,z) refer to cylindrical galactocentric coordinates.

2.1. Axisymmetric model

The axisymmetric model (see e.g. Vallee 1991; Poezd et al. 1993) is one of the simplest descriptions of the GMF. It is compatible with a non-primordial origin of the galactic magnetism, based on the dynamo theory. There are several possible models of the family of the axisymmetric distribution. The components for this model of GMF are given by $\begin{array}{ccc}{\mathit{B}}_{\mathit{r}}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\sin \mathrm{\left(}\mathit{p}\mathrm{\right)}\cos \mathrm{\left(}\mathit{\chi }\mathrm{\right(}\mathit{z}\mathrm{\left)}\mathrm{\right)}\\ {\mathit{B}}_{\mathit{\phi }}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\cos \mathrm{\left(}\mathit{p}\mathrm{\right)}\cos \mathrm{\left(}\mathit{\chi }\mathrm{\right(}\mathit{z}\mathrm{\left)}\mathrm{\right)}\\ {\mathit{B}}_{\mathit{z}}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\sin \mathrm{\left(}\mathit{\chi }\mathrm{\right(}\mathit{z}\mathrm{\left)}\mathrm{\right)}\mathit{,}\end{array}$where p is the pitch angle¹, which is considered to be constant; B₀(r) is the field strength (which in principle might be a function of the radial distance), and χ(z) corresponds to a “tilt angle”. For our study, we have adopted the following functional dependence $\mathit{\chi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{\chi }}_{\mathrm{0}}\tanh \left(\frac{\mathit{z}}{{\mathit{z}}_{\mathrm{0}}}\right)\mathit{,}$(3)where we use a value of z₀ = 1   kpc for the characteristic scale of variation in the vertical direction.

For the computations in this paper, we shall consider two cases for the radial dependence of B₀(r). The first case (hereafter ASS) corresponds to B₀(r) = B₀ (constant value). As a second case (hereafter ASS(r)), we consider a radial variation of the type ${\mathit{B}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{1}}}{\mathrm{1}\mathrm{+}\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{1}}}}\mathit{,}$(4)where r₁ represents the characteristic scale at which B₀(r) decreases to half its value at the galactic centre. This radial variation is based on a possible extension of the model by Poezd et al. (1993). The expression has appropriate asymptotic behaviors, in the sense that we obtain a finite value when r is close to the galactic center (r → 0), and asymptotically tends to  ∝ 1/r when r → ∞, as suggested in Battaner & Florido (2007). We must note that in Eq. (4), B₁ and r₁ are not independent, provided that we fix the value of the GMF strength in the solar neighborhood. For example, using R₀ = 8 kpc for the Sun galactocentric distance, and B_⊙ = 3 μG for the magnetic field strength in the solar neighborhood, we can re-write Eq. (4) in terms of a single free-parameter as ${\mathit{B}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}{\mathit{r}}_{\mathrm{1}}\mathrm{+}\mathrm{24}}{{\mathit{r}}_{\mathrm{1}}\mathrm{+}\mathit{r}}\mathit{,}$(5)where r is given in kpc and B₀(r) in μG.

Summarizing, a particular ASS model is fully described once these three parameters are given:  [B₀,p,χ₀] , while in the ASS(r) model, one would in principle require four parameters. However, for the ASS(r) family of models, we will use the constraint given in Eq. (5), which in practice means that we will only have three free parameters to represent a certain model:  [r₁,p,χ₀] .

The typical range of variation of the p values found in the literature for the ASS model of the disk (see Vallee 1991) is shown in Table 1. In general, different pitch angle and field strength values are derived for the different spiral arms of the Galaxy.

2.2. Bisymmetric model

This model is compatible with a primordial origin of the cosmic magnetism. It could explain the reversals of the magnetic field derived from observations of RM of pulsars (e.g. Han & Qiao 1994; Han 2001; Han et al. 2006). The three components for this model are given by $\begin{array}{ccc}{\mathit{B}}_{\mathit{r}}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\cos \left(\mathit{\phi }\mathrm{\pm }\mathit{\beta }\ln \left(\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{0}}}\right)\right)\sin \mathrm{\left(}\mathit{p}\mathrm{\right)}\cos \mathrm{\left(}\mathit{\chi }\mathrm{\right(}\mathit{z}\mathrm{\left)}\mathrm{\right)}\\ {\mathit{B}}_{\mathit{\phi }}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\cos \left(\mathit{\phi }\mathrm{\pm }\mathit{\beta }\ln \left(\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{0}}}\right)\right)\cos \mathrm{\left(}\mathit{p}\mathrm{\right)}\cos \mathrm{\left(}\mathit{\chi }\mathrm{\right(}\mathit{z}\mathrm{\left)}\mathrm{\right)}\\ {\mathit{B}}_{\mathit{z}}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\sin \mathrm{\left(}\mathit{\chi }\mathrm{\right(}\mathit{z}\mathrm{\left)}\mathrm{\right)}\mathit{,}\end{array}$where B₀(r) is the field strength; β = 1/tan(p), with p the pitch angle; R₀ is the distance Sun-Galactic center (≈8 kpc); and finally χ(z) is the tilt angle, which we assume is also given by Eq. (3). Note that in Eq. (6) we are considering two possible families of bisymmetric models, which correspond to the positive (negative) sign inside the big parenthesis for B_r and B_φ. In this way, we are considering the two conventions that are found in the literature. “Positive” bisymmetric models (hereafter BSS₊) are defined with the same convention to recover the model proposed by Han & Qiao (1994), while “negative” bisymmetric models (hereafter BSS₋) use the same sign convention of Jansson et al. (2009). Note that for BSS₋ the spiral has the opposite sign for the magnetic direction.

As in the previous case, we consider two sub-families of models. The first one, noted as BSS_±, corresponds to constant strength of the magnetic field. In this case, a model is completely defined by giving three parameters:  [B₀,p,χ₀] . The second family, noted as BSS_±(r), includes a radial variation of the field strength according to Eq. (4). For this case, we again fix the field strength to the value in the Solar neighborhood (Eq. (5)), so a model is specified by giving three parameters:  [r₁,p,χ₀] .

Typical parameter values for the disk model found in the literature, after conversion to our convention for the sign of the pitch angle, are p = (8.2 ± 0.5)°, B₀ = 1.4 μG and B_global = (1.8 ± 0.3) μG (Han & Qiao 1994); p = 14° and B₀ = 1 μG (Simard-Normandin & Kronberg 1980); and p = 11° and B₀ = (2.1 ± 0.3) μG (Han et al. 2006).

2.3. Concentric circular ring (CCR) model

This model was proposed by Rand & Kulkarni (1989) to fit the RM of their pulsars catalog, with the aim of providing a way to take into account reversals of the magnetic field at different radii. In their original expressions, they do not consider a vertical component for the magnetic field (i.e. B_z = 0). Here, we shall extend the equations for this model presented in Indrani & Deshpande (1999) to account for a vertical dependence in the following way $\begin{array}{ccc}{\mathit{B}}_{\mathit{r}}& \mathrm{=}& \mathrm{0}\\ {\mathit{B}}_{\mathit{\phi }}& \mathrm{=}& \frac{{\mathit{B}}_{\mathrm{0}}}{\sin \mathrm{(}\mathit{\pi }{\mathit{D}}_{\mathit{r}}\mathit{/}\mathit{w}\mathrm{)}}\sin (\frac{\mathit{\pi }\mathrm{(}\mathit{r}\mathrm{-}{\mathit{R}}_{\mathrm{0}}\mathrm{+}{\mathit{D}}_{\mathit{r}}\mathrm{)}}{\mathit{w}})\cos \mathrm{\left(}\mathit{\chi }\mathrm{\right(}\mathit{z}\mathrm{\left)}\mathrm{\right)}\\ {\mathit{B}}_{\mathit{z}}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\sin \mathrm{\left(}\mathit{\chi }\mathrm{\right(}\mathit{z}\mathrm{\left)}\mathrm{\right)}\mathit{,}\end{array}$where w is the space between reversals; D_r is the distance to the first reversal; B₀ is the field strength, and finally χ(z) is given by Eq. (3). All distances (w and D_r) are given in kiloparsecs. Note that we introduce an additional factor sin(πD_r/w) in the definition of B_φ, so B₀ still preserves the meaning of the magnetic field strength in the solar neighborhood.

In our analyses we did not consider any radial dependence of the field strength, so for the CCR model the parameter space is defined by  [D_r,w,B₀,χ₀] .

The best-fit values obtained by Rand & Kulkarni (1989) for the disk component were D_r = (0.6 ± 0.08) kpc , w = (3.1 ± 0.08) kpc and B₀ = (1.6 ± 0.2) μG.

2.4. Logarithmic spiral arms (LSA) model

This model was used by Page et al. (2007) to describe the distribution of the large-scale pattern of the polarization angle in the 22 GHz WMAP 3-year data. The equations which describe the field are $\begin{array}{ccc}{\mathit{B}}_{\mathit{r}}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\sin \mathit{\psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\cos \mathit{\chi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\\ {\mathit{B}}_{\mathit{\phi }}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\cos \mathit{\psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\cos \mathit{\chi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\\ {\mathit{B}}_{\mathit{z}}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\sin \mathit{\chi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{,}\end{array}$where

$\mathit{\psi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\psi }}_{\mathrm{0}}\mathrm{+}{\mathit{\psi }}_{\mathrm{1}}\ln \left(\frac{\mathit{r}}{\mathrm{8} \mathrm{kpc}}\right)$and

$\mathit{\chi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{\chi }}_{\mathrm{0}}\tanh \left(\frac{\mathit{z}}{\mathrm{1} \mathrm{kpc}}\right)\mathrm{·}$Note that the LSA model is essentially an axisymmetric model in which the pitch angle is not constant, because magnetic field lines are logarithmic spirals. According to our definition of pitch angle for the ASS model, ψ(r) would play the role of pitch angle, in which we would have a constant part, ψ₀, and a characteristic amplitude for the logarithmic dependence of the arms, ψ₁. Following Page et al. (2007), B₀ is assumed constant with a value of 3 μG.

Thus, the parameter space is defined by  [ψ₀,ψ₁,χ₀] . The proposed values for the different parameters are  [ψ₀,ψ₁,χ₀]  =  [27°,0.9°,25°]  (see Page et al. 2007, and also the erratum available at the LAMBDA web site²).

2.5. Bi-toroidal (BT) model (halo model)

Some authors have found hints of a halo or thick disk field component with opposite directions in both hemispheres. For example, Han & Wielebinski (2002) and Prouza & Šmída (2003) detected it with maximum strengths at a large height of 1.5 kpc at both sides above and below the plane, with maximum strengths of approximately 1 μG. Sun et al. (2008) give a more complete description of this double torus, with the height of maximum strength again at 1.5 kpc, but its maximum is much stronger (around 10 μG).

Following this scenario, we propose a possible configuration for the components of the magnetic field which do contain a different sign in both hemispheres, and it is given by $\begin{array}{ccc}{\mathit{B}}_{\mathit{r}}& \mathrm{=}& \mathrm{0}\\ {\mathit{B}}_{\mathit{\phi }}& \mathrm{=}& {\mathit{B}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\arctan \left(\frac{\mathit{z}}{{\mathit{\sigma }}_{\mathrm{1}}}\right)\exp \left(\frac{\mathrm{-}{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}\right)\\ {\mathit{B}}_{\mathit{z}}& \mathrm{=}& \mathrm{constant}\mathit{,}\end{array}$where σ₁ and σ₂ are two constants (measured in kpc) which encode the characteristic scales of variation of the field with the vertical distance, and do take into account in a simplified way the change in the sign. For our computations, we fix the B_z value to be 0.2 μG (Han & Qiao 1994), and we only consider the case of a radial variation of B₀(r) given by Eq. (5). Thus, the parameter space for this model is specified by  [r₁,σ₁,σ₂] .

Fig. 1 Observed Stokes Q (top) and U (bottom) maps at 22 GHz from WMAP5 data. The maps are shown in a Mollweide projection of HEALPix in Galactic coordinates. The center of the map corresponds to (l,b) = (0°,0°), and the graticule increases with Δl = Δb = 45°. The maps are degraded to nside = 16 to enhance the large scale pattern. Units are mK (thermodynamic).

3. Methodology

In order to obtain constraints on the different parameters for each one of the models described above, we performed a systematic comparison of the predicted polarized intensity due to synchrotron emission with the observed map at 22 GHz by WMAP satellite. Here, we describe the relevant details of the dataset, the numerical procedure to compute the synchrotron maps for a certain GMF model, and the model selection criterion that we adopted for our analyses.

3.1. Description of the K-band WMAP5 data

The analysis of this paper is based on a comparison with the K-band (equivalent to 22 GHz) polarization map obtained by the WMAP satellite after five years of operation (Hinshaw et al. 2009). This map is publicly available in the LAMBDA website³, and it is given in HEALPix⁴ format (Górski et al. 2005).

Figure 1 shows the all-sky Stokes Q and U maps at 22 GHz, degrading the resolution to enhance the large-scale pattern, using an nside = 16 pixelization (which corresponds to a pixel size of 3.66°). For all the computations in this paper, we will use these degraded maps as the input data. As shown by the WMAP team (Page et al. 2007), at this frequency (22 GHz) the large-scale pattern observed in polarization is completely dominated by the galactic contribution, and the CMB component is sub-dominant. Thus, we can safely neglect the contribution of the CMB to the polarization map in our analyses.

From these two observables (Stokes Q and U parameters), one can obtain the map of the direction of the polarization angle (hereafter PA) as ${\mathit{\gamma }}_{\mathrm{obs}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\arctan \left(\frac{\mathit{U}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}}{\mathit{Q}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}}\right)\mathrm{+}\frac{\mathit{\pi }}{\mathrm{2}}\mathit{,}$(10)where is the direction of the line of sight. Note that we include in our definition the π/2 factor, so Eq. (10) represents the angle of the magnetic field, and takes values in the  [ 0,π ]  region. Note that the polarization convention adopted here is the one described in HEALPix, labelled as COSMO, which differs from the usual IAU convention in a minus sign for Stokes U parameter. In addition, the WMAP definition of Stokes parameters includes an additional 1/2 factor with respect to the definition used by Chandrasekhar (1960), which will be the one adopted here. All these quantities (U, Q and PA) are defined in a galactic (heliocentric) coordinate system. This has to be taken into account when comparing the observed maps with the models. Figure 2 shows the observed direction of the PA at the same 3.66° (nside = 16) resolution.

Fig. 2 Observed direction of the polarization angle (PA) at 22 GHz, obtained from the two maps shown in Fig. 1. The PA is defined here as the local direction of the magnetic field (see Eq. (10)). Units are degrees.

3.1.1. Noise maps

In order to perform the model selection, we also need to estimate the noise maps associated to the data. It is important to note that those noise maps should account for the covariances introduced by all components which are present in the observed data but are not included in the theoretical model. In particular, it should account for the instrumental noise component and for the random component of the GMF, which are not included in our theoretical model. The first one could be easily estimated based on the information provided by the WMAP team about instrument sensitivity and the overall integration time that the satellite has spent on each particular pixel. However, this would not account for the second part of the covariance. Thus, in order to model all the different contributions to the covariance, we follow a different procedure, which makes use of the fact that the original WMAP-K band maps have a much better angular resolution than the pixel size that is finally adopted in our analyses.

Noise maps for Q and U.

In this paper we probed three different methods to characterize the noise maps associated to Q and U. In addition to the pure instrumental noise, all these three methods attempt to estimate the contribution to the total covariance matrix of the residual astrophysical components which are not included in our model (for example, the random component of the magnetic field, or the variance introduced by point-to-point variations of the mean level of the galactic emission). All methods produce a noise map at nside = 16 resolution, which provides a pixel size of  ~3.66°. The three methods to obtain the σ_Q map (or equivalently, the σ_U), are

• Method 1. This corresponds to the same procedure asdescribed in Jansson etal. (2008). We start from the observed WMAPStokes-Q map at full resolution (nside = 512), degrading it to  ~0.92° pixel resolution (nside = 64). For a given pixel i within our nside = 16 pixelization scheme, we obtain the associated noise σ_Q(i) by computing the square root of the variance of the  ~0.92° pixels inside a radius of 2° from the center of our pixel i.

• Method 2. We start from the observed Q-map at nside = 512, and we convolve it with a Gaussian of FWHM = 1°. For a given pixel i within our nside = 16 pixelization scheme, the associated noise σ_Q(i) is computed as the square root of the variance in each pixel of the smoothed map within a radius of 2° from the center of our pixel i.

• Method 3. We start from the observed WMAP Stokes-Q map degraded at nside = 16. For each pixel i, the noise σ_Q(i) is computed by obtaining the square root of the variance inside a circle of radius r ~ 7.4° (i.e. twice the pixel radius). This method produces an unbiased estimate of the noise map in the case of uncorrelated noise, provided that the scale of variation of the noise map is larger than 7.4°.

Each one of these three methods produces different results for the noise map in certain areas, specially those dominated by galactic features. In turn, these differences imply significant changes in the goodness-of-fit statistic (up to a factor of 2−3 in some cases).

In general, we find that the larger noise amplitudes are obtained with method 1, while the smaller amplitudes are obtained with method 2.

Thus, and as a conservative approach, we decided to present the computations of this paper using the method 3.

The implications of this uncertainty in the determination of the noise maps are discussed in Sect. 5.

Noise maps for the Stokes Q and U parameters, obtained with method 3, are shown in Fig. 3. These maps have been used to obtain the final results presented in Tables 3 and 4.

Noise maps for PA.

Once these two noise maps (σ_Q and σ_U) have been obtained, one can in principle obtain the noise map associated to the PA (σ_γ) from them. However, there is an important point here. Equation (10), which defines the PA, is not linear in Q and U. This implies that the variance map for the PA will depend on the particular model which is used to compute the average value; or in other words, the noise map for PA will depend not only on σ_Q and σ_U, but also on Q and U themselves. Therefore, when doing the model selection, the noise map for PA will also be a function of the model.

To illustrate this, we compute here the noise map for PA, using as reference model the observed WMAP-K map. We use here a Monte Carlo method, drawing N_sim realizations of pairs of (Q,U) maps, with mean equal to the observed (Q,U) maps, and variance given by σ_Q and σ_U, as computed in previous step. We have checked that using N_sim = 5000 noise realizations is enough to get convergence on the σ_γ map at the per cent level. Figure 3 also shows the noise map for the PA, obtained for this particular case.

Fig. 3 Noise maps for the WMAP5 datasets presented in Figs. 1 and 2. They show the Stokes Q (top), Stokes U (middle) and position angle (PA, bottom) maps. Units for Stokes Q and U maps are mK. Units for the PA map are degrees. The noise map for PA has been computed as fluctuations around the observed WMAP5 K-band map.

3.2. Producing the synchrotron polarized emission for a given magnetic field model

To obtain the predicted polarized synchrotron emission for any of the Galactic magnetic field models described above, we have developed a code which only includes the relevant physics, and it is optimized in terms of computational time with the aim of efficiently performing all the computations. The code works directly within the HEALPix pixelization scheme, and obtains the predictions for the synchrotron emission directly at the resolution level that we have chosen (i.e. nside = 16). As discussed above, this resolution is enough as long as we are interested in the large-scale pattern of the Galactic emission.

In general, assuming that the cosmic ray spectrum is a power-law distribution of spectral index p, we can predict the Stokes Q and U parameters that characterize the polarization of the synchrotron emission at a certain frequency by computing the emissivity (energy per unit time per unit volume per frequency per solid angle) in the two orthogonal directions, parallel and perpendicular to the projection of the magnetic field on the plane of the sky. Following Rybicki & Lightman (1986), we have $\begin{array}{ccc}& & {\mathit{ϵ}}_{\mathrm{\perp }}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\mathit{N}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\frac{\sqrt{\mathrm{3}}{\mathit{e}}^{\mathrm{3}}}{\mathrm{8}\mathit{\pi m}{\mathit{c}}^{\mathrm{2}}}{\left(\frac{\mathrm{4}\mathit{\pi mc}}{\mathrm{3}\mathit{e}}\right)}^{\frac{\mathrm{1}\mathrm{-}\mathit{p}}{\mathrm{2}}}{\mathit{\nu }}^{\frac{\mathrm{1}\mathrm{-}\mathit{p}}{\mathrm{2}}}\mathrm{(}\mathit{B}\sin \mathit{\alpha }{\mathrm{)}}^{\frac{\mathit{p}\mathrm{+}\mathrm{1}}{\mathrm{2}}} \\ & & \end{array}$(11) $\begin{array}{ccc}& & {\mathit{ϵ}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\mathit{N}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\frac{\sqrt{\mathrm{3}}{\mathit{e}}^{\mathrm{3}}}{\mathrm{8}\mathit{\pi m}{\mathit{c}}^{\mathrm{2}}}{\left(\frac{\mathrm{4}\mathit{\pi mc}}{\mathrm{3}\mathit{e}}\right)}^{\frac{\mathrm{1}\mathrm{-}\mathit{p}}{\mathrm{2}}}{\mathit{\nu }}^{\frac{\mathrm{1}\mathrm{-}\mathit{p}}{\mathrm{2}}}\mathrm{(}\mathit{B}\sin \mathit{\alpha }{\mathrm{)}}^{\frac{\mathit{p}\mathrm{+}\mathrm{1}}{\mathrm{2}}}\\ & & \end{array}$(12)where B is the magnetic field, ν is the frequency, p is the spectral index, and e and m are the electron charge and mass, respectively. The function N(r,z) represents the electron number density at the corresponding position (r,z) in the Galaxy, and it is obtained from Eq. (1). From these two equations, the polarized intensity at a given frequency is obtained by integrating the emissivity along the line of sight ${\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z,}\mathit{n̂}\mathrm{\right)}\mathrm{\equiv }\mathrm{(}{\mathit{Q}}_{\mathit{\nu }}\mathrm{-}\mathit{i}{\mathit{U}}_{\mathit{\nu }}\mathrm{)}\mathrm{=}\mathrm{\int }\mathrm{\left[}{\mathit{ϵ}}_{\mathrm{\perp }}\mathrm{\right(}\mathit{\nu ,z,}\mathit{n̂}\mathrm{)}\mathrm{-}{\mathit{ϵ}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{(}\mathit{\nu ,z,}\mathit{n̂}\mathrm{\left)}\mathrm{\right]}{\mathrm{e}}^{\mathrm{-}\mathit{i}\mathrm{2}\mathit{\gamma }\mathrm{\left(}\mathit{z,}\mathit{n̂}\mathrm{\right)}}\mathrm{d}\mathit{z}\mathit{,}$(13)where we have set the coordinate system in a way that the z-axis represents the line-of-sight direction, and the other two directions are contained in the plane on the sky, with the y-axis pointing east and x-axis pointing south (i.e. HEALPix coordinate convention, as explained in Sect. 3.1). With these definitions, the Stokes Q and U parameters are given by $\begin{array}{ccc}{\mathit{Q}}_{\mathit{\nu }}& \mathrm{=}\mathrm{\int }\mathrm{\left[}{\mathit{ϵ}}_{\mathrm{\perp }}\mathrm{\right(}\mathit{\nu ,z,}\mathit{n̂}\mathrm{)}\mathrm{-}{\mathit{ϵ}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{(}\mathit{\nu ,z,}\mathit{n̂}\mathrm{\left)}\mathrm{\right]}\cos \mathrm{\left(}\mathrm{2}\mathit{\gamma }\mathrm{\right)}\mathrm{d}\mathit{z}& \\ {\mathit{U}}_{\mathit{\nu }}& \mathrm{=}\mathrm{\int }\mathrm{\left[}{\mathit{ϵ}}_{\mathrm{\perp }}\mathrm{\right(}\mathit{\nu ,z,}\mathit{n̂}\mathrm{)}\mathrm{-}{\mathit{ϵ}}_{\mathrm{\left|}\mathrm{\right|}}\mathrm{(}\mathit{\nu ,z,}\mathit{n̂}\mathrm{\left)}\mathrm{\right]}\sin \mathrm{\left(}\mathrm{2}\mathit{\gamma }\mathrm{\right)}\mathrm{d}\mathit{z}\mathit{,}& \end{array}$By assuming a spectral index p = 3, we obtain the simulated Q and U components along the line of sight (z-axis) by numerical integration $\begin{array}{ccc}{\mathit{Q}}_{\mathit{\nu }}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}& \mathrm{=}{\mathit{K}}_{\mathit{Q}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathrm{\int }}_{\mathrm{LOS}}{\mathit{N}}_{\mathrm{e}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathrm{[}{\mathit{B}}_{\mathit{x}}^{\mathrm{2}}\mathrm{-}{\mathit{B}}_{\mathit{y}}^{\mathrm{2}}\mathrm{]}\mathrm{d}\mathit{z}& \\ {\mathit{U}}_{\mathit{\nu }}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}& \mathrm{=}\mathrm{-}{\mathit{K}}_{\mathit{U}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathrm{\int }}_{\mathrm{LOS}}{\mathit{N}}_{\mathrm{e}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathrm{2}{\mathit{B}}_{\mathit{x}}{\mathit{B}}_{\mathit{y}}\mathrm{d}\mathit{z}\mathit{,}& \end{array}$where we explicitly introduce a minus sign in the equation for the U component, in order to convert from the IAU convention for the polarization to the HEALPix convention, which is used in the WMAP maps. The K_U(ν) and K_Q(ν) constants also include the conversion factors between brightness and temperature. At 22 GHz, we can safely use the Rayleigh-Jeans approximation. Substituting the numerical values, we obtain K_Q(ν) = 1.41 × 10¹¹ mK cm³ (μG)^-2 kpc^-1, and K_U(ν) = 1.25 × 10¹¹ mK cm³ (μG)^-2 kpc^-1. The simulated PA map is derived from these two Eq. (15) as $\mathit{\gamma }\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathrm{=}\mathrm{0.5}\arctan \left(\frac{\mathrm{-}{\mathit{K}}_{\mathit{U}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathrm{\int }}_{\mathrm{LOS}}{\mathit{N}}_{\mathrm{e}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathrm{2}{\mathit{B}}_{\mathit{x}}{\mathit{B}}_{\mathit{y}}\mathrm{d}\mathit{z}}{{\mathit{K}}_{\mathit{Q}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathrm{\int }}_{\mathrm{LOS}}{\mathit{N}}_{\mathrm{e}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathrm{[}{\mathit{B}}_{\mathit{x}}^{\mathrm{2}}\mathrm{-}{\mathit{B}}_{\mathit{y}}^{\mathrm{2}}\mathrm{]}\mathrm{d}\mathit{z}}\right)\mathrm{+}\frac{\mathit{\pi }}{\mathrm{2}}\mathit{,}$(16)where B_x and B_y represent in our coordinate system the two components of the magnetic field which are perpendicular to the line of sight.

Finally, when predicting the expected synchrotron emission for a particular model, we have included an additional restriction on the line-of-sight integration by excluding those points whose galactocentric radial coordinate r_G is smaller than 3 kpc or larger than 20 kpc. The first restriction excludes the inner region of the Galaxy, where large deviations from the regular pattern are expected (La Rosa et al. 2006), while the second one introduces a radial cut-off. In any case, we have checked that the results are robust against changes in these numbers.

3.3. Exploration of the parameter space

As described in Sect. 2, for each one of the families of GMF models we have a set of parameters which define each particular model. Given that in all cases the dimension of the parameter space is small (there are, at the most, four parameters describing a particular model), we decided to carry out the exploration of the parameter space using a grid-based approach. For higher dimensions, a Monte Carlo method would be more appropriate.

For each one of the different GMF models we constructed three different grids of models, which we label as “literature”, “blind” and “non-blind”. The first one is centered around the average values which are found in the literature for each one of the different parameters. The second one spans the maximum range which is reasonably expected for each particular parameter. Finally, the third one is built a-posteriori, once the model selection has been performed on the previous grid, by centering the new grid around the best-fit parameters for each case. Table 1 summarizes all the relevant parameters for each one of these three grids. In total we computed more than one million models (290 000 for the blind grid, 970 000 for the non-blind, and 51 000 for the literature) for all the different GMF models described in Sect. 2. Each one of those models corresponds to a set of three maps (Q, U and PA) of the expected synchrotron polarized emission of the sky at 22 GHz.

To give an indication, the average execution time in a standard desktop computer for the computation of a particular model requires  ≲ 4 s of CPU time. Thus, the total CPU time for the construction of all the grids is around 1500 CPU hours.

3.4. Model selection and parameter estimation for each GMF model

Once we explored the parameter space with these three grids, we derived the best-fit parameters for each one of the GMF models using a Bayesian approach. To this end, we have to both compute the likelihood function (ℒ), and to provide an expression for the priors. Once we obtained the best-fit parameters for each GMF model, different models will be compared in terms of the reduced-χ² statistic.

3.4.1. Likelihood function

In general, we may assume that the likelihood function is defined as a multivariate Gaussian when it is written in terms of the observables, i.e. $\ln \mathrm{ℒ}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\chi }}^{\mathrm{2}}\mathit{.}$(17)If we assume that the correlations between the different pixels are negligible, we have ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}}\frac{\mathrm{(}{\mathit{x}}_{\mathit{i}}\mathrm{-}{\mathit{k}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}\mathit{,}$(18)where x_i represents the observational data, k_i the simulated data and , the associated noise covariance. In our case, we performed two different evaluations of the likelihood:

• the first one corresponds to a direct comparison of Stokes Q and U parameters. In this case, we have i = 1,...,2N_pix, and x_i = Q_i for i = 1,N_pix; and x_i = U_i − Npix for i = N_pix + 1,...,2N_pix. This case will be noted as ;

• the second case corresponds to the comparison of PA, so now we have i = 1,...,N_pix, and x_i = γ_i. Note that in this case, σ_i will depend on k_i, and thus the PA noise map needs to be computed for each particular model. This case will be noted as , and will be used for comparison with the results of the previous case.

Once these two functions ( and ) are evaluated in all the data-points of the different grids, the posterior distributions are obtained, and then marginalized⁵ over all the unwanted parameters. At the end, we end up with marginal probability distribution functions for each one of the parameters. From these, confidence intervals are derived as the 0.16, 0.5 and 0.84 points of the cumulative probability distribution function. Thus, our parameter estimate is the median of the marginalized posterior probability distribution function, and the confidence interval encompasses 68 per cent of the probability.

3.4.2. Priors

For the analyses in this paper, we have not introduced prior information on the parameter values describing any of the models. This is equivalent to say that we have implicitly adopted a top-hat prior in all the parameters, where the top-hat function is defined by the parameter ranges presented in Table 1. Thus, in all cases, the evaluation of the posterior would reduce to the computation of the likelihood function (ℒ).

However, for the case of (Q,U) analysis, we slightly modified the standard analysis in the following way. As discussed in Sect. 2, the amplitude of the CR electron spectrum in the solar neighborhood is highly uncertain. This will in turn imply a large uncertainty in the recovered strength for the magnetic field, and moreover, it might produce a bias on the recovered GMF amplitude. To account for this additional degree of uncertainty (at least at first order), we have introduced an additional parameter ϵ, which multiplies to both the predicted Q and U maps for a given GMF model. Note that this parameter has no impact on PA. If the CR electron density were entirely correct, we would have ϵ = 1. If there is an uncertainty in this parameter due to the modeling of the CR distribution, we account for it by introducing a Gaussian prior for this additional parameter $\mathrm{-}\mathrm{2}\ln {\mathrm{ℒ}}_{\mathit{QU}}\mathrm{=}\frac{\mathrm{(}\mathit{ϵ}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{ϵ}}^{\mathrm{2}}}\mathrm{+}{\mathit{\chi }}_{\mathit{Q,U}}^{\mathrm{2}}\mathit{,}$(19)and we marginalize over it. Given the existing uncertainty on the CR electron density, we have chosen a conservative value of σ_ϵ = 0.8. The marginalization over this additional parameter, ϵ, can be done analytically, yielding $\ln {\mathrm{ℒ}}_{\mathit{QU}}\mathrm{=}\frac{\mathit{A}\mathrm{+}\mathrm{2}\mathit{C}}{\mathrm{4}\mathit{B}}\mathrm{+}\frac{\ln \mathit{B}}{\mathrm{4}}\mathit{,}$(20)where $\mathit{A}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i}}\frac{{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{ϵ}}^{\mathrm{2}}}$(21)$\mathit{B}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i}}\frac{{\mathit{k}}_{\mathit{i}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{ϵ}}^{\mathrm{2}}}$(22)$\mathit{C}\mathrm{=}\sum _{\mathit{i}}\frac{{\mathit{x}}_{\mathit{i}}{\mathit{k}}_{\mathit{i}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{ϵ}}^{\mathrm{2}}}\mathrm{·}$(23)Note that this scheme is completely equivalent to the marginalization over calibration uncertainties which is also adopted in CMB experiments (see e.g. Bridle et al. 2002).

3.5. Masks

Nearby structures in the Galaxy (e.g. supernova remnants) might distort the regular pattern of the GMF as seen in the synchrotron polarized emission, thus biasing the determination of some parameters of the GMF model. To probe the robustness of the estimates as a function of the region which is used for the fit, we have considered different masks to constrain the global, halo and disk components of the GMF. In all cases, we exclude the galactic center region defined as the region in which the line of sight crosses a region with radius 3 kpc in cylindrical coordinates, that is not taken into account in our integration scheme. To constrain the global component we considered two masks, listed in decreasing order of covered sky fraction:

• Mask 1. It corresponds to the Galactic center region.

• Mask 2. It combines mask 1 with the well-known local regions with a strong polarized intensity that could be interpreted as part of the polarized intensity produced by the regular GMF. Here we considered the four “loops” described in Berkhuijsen et al. (1971), and we added the polarization mask from WMAP team, which includes a better masking of the North Galactic spur region and several other small objects (e.g. LMC).

To constrain the halo component we exclude the following regions:

• Mask 3. It excludes the disk defined as the emission contains in |b|  < 10°.

• Mask 4. It excludes the emission of the disk (mask 3) and the “loops” defined in mask 2.

For completeness, we also considered two masks to constrain the field pattern in the disk. They are

• Mask 5. It excludes the halo region defined as the emissionobtained for |b|  > 10°.

• Mask 6. It excludes the halo and the “loops” described in Berkhuijsen et al. (1971).

In Fig. 4 we show all regions which have been described in this subsection. Figure 5 shows the polarization mask used by the WMAP team⁶. Table 2 presents the detailed information about the sky coverage of each one of these masks. It also indicates the number of available pixels for the analysis (i.e. number of terms in the summation in Eq. (18). This quantity is relevant in order to compute the reduced-χ² for the best-fit models. As a reference, note that in this pixelization (nside = 16), a whole-sky map contains 3072 pixels. Therefore, we would have in this case N_pix(QU) = 6144.

Note that even if we use these masks to eliminate the effects of the random local spurs, other undetected random local features could be a source of errors in our large-scale models and these errors are difficult to quantify. But our spur masks probably eliminate the major contribution of local features.

Fig. 4 Regions used for the definition of the six masks adopted for the analyses (see text for details).

Fig. 5 Mask for the galactic polarized emission used by the WMAP team.

4. Results and discussion

For each one of the masks described in Table 2 and for each one of the GMF models presented in Sect. 2, we evaluate the posterior distribution in each one of the three grids, and marginalizing over the relevant parameters, we obtain the corresponding confidence regions. Our default analysis uses the Q − U maps, and the noise maps obtained with method 3 in Sect. 3.1.1. The results are summarized in Tables 3 and 4.

In order to evaluate which model better reproduces the data, we used as a goodness-of-fit the reduced-χ² statistic, which is obtained as the ratio of the the minimum value for χ²(≡  −2lnℒ_QU) to the number of degrees of freedom (hereafter d.o.f.). The number of dof is obtained as N_pix − M, being N_pix the number of terms in Eq. (18), and M the number of parameters for the considered GMF model. The last two columns in each table present the values for the minimum χ², and the reduced-χ².

4.1. The magnetic field in the Galactic halo

The results for the halo field are summarized in Table 3.

Our reference computation corresponds to the case labelled as mask 4, in which we exclude the emission of the disk, the “loops” (including the polarization mask provided by the WMAP team), and the galactic center region (which is not accounted for in our analysis). For this reference case, the model with the minimum reduced-χ² is the ASS(r), although the other two axisymmetric models (LSA and ASS) provide approximately the same goodness-of-fit. Thus, the large-scale 22 GHz polarized synchrotron emission seems to be more compatible with some type of axisymmetry, a conclusion also reached by Page et al. (2007), and also compatible with the results shown by Sun et al. (2008). For illustration, Fig. 6 shows the field pattern of the best-fit ASS(r) model at z = 4 kpc; Fig. 7 shows the marginalized one-dimensional posteriors distributions for the parameters of this model; and Fig. 8 shows the predicted Stokes Q and U and PA maps for the same best-fit model. We now discuss each one of the relevant parameters separately.

Radial scale.

For this ASS(r) model, the derived constraint on r₁ is  <2.5 kpc (95% confidence level). This parameter essentially controls the distance at which the magnetic field is no longer constant and begins to decrease proportional to r^-1 (see Eq. (5)). The obtained value is indeed very small compared with the radial scale of the electron density or any other scale distance in our galaxy, suggesting that the data require an important variation of the field in the inner part of the halo, probably due to the presence of stronger magnetic fields at the galactic center (see e.g. Roy et al. 2008). Indeed, in the literature, the halo model proposed in Prouza & Šmída (2003) requires also a small radial scale of r₁ = 4 kpc, where the radial dependence is  ∝  in that case.

Field strength.

For the ASS(r), we adopted a fixed value for the magnetic field strength of B₀ = 3 μG at the solar neighborhood. The original LSA model proposed by Page et al. (2007) also assumed this fixed value. However, it was considered as a free parameter for the ASS model, and was found to agree with that value ( μG).

Pitch angle.

The pitch angle was considered in the ASS(r) and ASS models as a constant free parameter. In both cases, we consistently obtain a value between 24° and 26° respectively. For the LSA model, it is considered as a radial function with a logarithmic dependence, and in this case the derived pitch angle at the solar neighborhood is again  ≈26°, a value which is consistent with that given by Page et al. (2007) of p = 27°. As discussed below, these values are higher than the typical pitch angles obtained for the disk field.

Tilt angle.

This parameter controls the vertical structure of the field. The derived tilt angle in all the axysimmetric models is of the order of 30°, which implies a vertical component close to 1 μG at z = 1 kpc. This vertical field component could be identified as the poloidal component corresponding to the dipole field responsible of the mG vertical component at the very center (Han & Qiao 1994). It could correspond as well to the vertical component of the cluster field diffused into the disk by turbulent magnetic diffusion (e.g. Battaner & Florido 2000), in this case to the Local Group field.

Other models.

Table 3 shows a good consistency between the three axisymmetric models. The rest of the families considered here provide slightly worse reduced-χ² values, although none can be clearly rejected. In general, the two families of bisymmetric models provide a slightly poorer goodness-of-fit. The derived radial scales, field strengths and pitch angles are in general similar to those found for the axisymmetric models. However, the tilt values (χ₀) are considerably lower, and in some cases negligible. This could be due to a compensation produced by the inherent reversals of the fields for a given direction in a bisymmetric configuration.

The two remaining models (CCR and BT) produce the poorer fits (χ² = 2.1), but again they cannot be rejected. With regard to these values, the CCR points to the existence of a reversal at r ~ 3 kpc from the galactic center. The bi-toroidal model has been suggested as a possible explanation for the halo double torus that it is interpreted by Han et al. (1997, 1999) and Han (2009) as a consequence of an α−Ω effect. Another possibility is that the aforementioned vertical field diffused from the galaxy cluster could be twisted by differential rotation in the vertical direction, producing toroidal fields above and below the plane with opposite directions in both hemispheres (e.g. Battaner & Florido 2000).

Fig. 6 Large-scale pattern of the ASS(r) model at z = 4 kpc. This model provides the best fit for the halo field.

Fig. 7 One dimensional (marginalized) posterior distribution functions for the parameters of ASS(r) halo model (top: r1, middle: p and bottom: χ0) when mask 4 is considered in the Stokes Q and U analysis.

Effect of the loops and the disk emission on the determination of the halo field.

As described in Sect. 3.5, nearby structures in the Galaxy might introduce biases on the recovered parameters of a given GMF model. Moreover, the emission of the disk could contaminate the halo field as discussed in Sect. 1. If the disk emission is not excluded by the mask, we could constrain the global component of the GMF and quantify the impact of the disk emission on the parameters describing halo field. The corresponding results are also shown in Table 3.

The influence of disk emission can be seen by comparing the results from masks 2 and 4. In both cases, the best-fits are obtained for axisymmetric models but when disk emission is not masked out, the reduced-χ² becomes slightly poorer. In general, the pitch and tilt angles and the field strength remain unchanged, with the important exception of the radial scale factor for ASS(r) and BSS_±(r) models. The inclusion of the disk emission in the analysis drastically increases this radial scale, probably because as shown below, the disk does not require strong radial variations.

Finally, we can also evaluate the impact of the loops on the fit by comparing the results of masks 3 and 4, or 1 and 2. The basic conclusion in this case is that including the loop regions in the analysis do not significantly bias the results (even for the radial scale parameter), but the quality of the fits get worse in all cases.

4.2. The magnetic field in the disk

For completeness, we also used two masks (5 and 6) to study the magnetic field in the Galactic plane by masking the halo emission. Table 4 summarizes the constraints on the different parameters for this case. The reference mask now is number 6, which also excludes the contribution of loops. However, the available sky area for the fit in this case is very small (7.8% in total, which corresponds to approximately 40% of the total area of the disk). Because of this limited area, the conclusions on the magnetic field parameters might be uncertain. Nevertheless, we consider that it is still important to compare these results with the numbers obtained with other methods.

Focusing on mask 6 alone, the three axisymmetric models (LSA, ASS and ASS(r)) provide practically identical values of the goodness-of-fit, which is slightly better than for the other cases. For these three models, the derived pitch angle values are lower than for the halo. In the solar neighborhood the pitch angle of the spiral arm is  ~18° for the stars and  ~13° for all gaseous components⁷ (see Vallee 1995; Vallée 2002). Therefore, our best-fit suggests that the magnetic arms follow the gas structure (p ~ 14°−15°). However, we note that Jansson et al. (2009) found a value of p ≈ 35° for the ASS+RING model proposed by Sun et al. (2008), which is not compatible with the one obtained here.

The two families of bisymmetric models provide good results, and again, the pitch angle values are low. Moreover, the derived constraints (between 7°−10°) are fully compatible with the results obtained by other authors. In the literature, values also range from 7.2° to 11° (see Han & Qiao 1994 and Han 2001 for values obtained with Faraday rotation of pulsars; and Heiles (1996), for values obtained with polarized starlight). For comparison purposes, Fig. 9 shows the pattern of our best-fit BSS₊ model in the galactic disk (z = 0), which is similar in shape to the results obtained by other authors (see e.g. Fig. 5 in Han & Qiao 1994).

We note that the values derived for the tilt angle in all models are on the order of  ≲ 19°, which again implies a z-dependence of the field strength within the disk which is compatible with the values observed by Han & Qiao (1994) and up to 0.4   μG for the thin disk.

Finally, note that we do not expect changes in these results if we include a more refined treatment of the turbulent magnetic field in the analysis, since the polarized synchrotron emission comes from the regular pattern of the GMF, which is located in the inter-arms regions (Beck 2007). Indeed, we find that the results obtained with the other two methods for the noise determination (see Sect. 3.1.1) are fully consistent with those presented here for all masks and all parameters. We only found a small dependence of the constrained value for χ₀ for axisymmetric and bisymmetric models with the three noise maps, which is on the order of  ≲ 15 per cent.

Fig. 8 Best-fit for the halo field by using the Stokes Q and U parameters and PA excluding the disk, the galactic center and the “loops” (ASS model with radial dependence).

Fig. 9 Large-scale pattern of the best fit BSS+ model for the disk emission in the galactic disk (z = 0).

4.3. Comparison of the result with the PA analysis

For comparison purposes, we also performed the analyses using the PA information alone. As discussed elsewhere, this analysis provides a limited amount of information because the PA is not sensitive to some parameters, like a constant field strength.

In general, the results using the PA provide compatible results to those of the Stokes Q and U analysis, although there is a larger number of unconstrained parameters and a poorer goodness-of-fit. To illustrate this, the ASS(r) analysis with mask 4 gives r₁ < 45.9 kpc, p ~ 20.0°, and χ₀ ~ 21.0° (χ² ~ 2.3); while for the ASS case gives p ~ 20.5°, χ₀ ~ 19.0° and no constraint on B₀ (χ² ~ 2.3). The pitch and tilt angles are compatible with those derived from the Stokes Q and U analysis. The best fit is given by the LSA with ψ₀ > 22.3°, ψ₁ > 8.6° and χ₀ < 18.1° (χ² ~ 2.3).

5. Conclusions

We have constrained the regular Galactic magnetic field by using the polarized emission at 22 GHz. To this aim we considered nine models of the Galactic magnetic field, each defined by three or four free parameters and for six different masks to interpret the polarized maps from WMAP5. The combination of models, free parameters, and masks produced a huge number of simulated maps to be compared with the observational ones in Stokes Q and U parameters, which in turn provide valuable constraints to determine the three dimensional configuration of the magnetic field of our Galaxy.

The family of GMF models that better describes the halo emission is the axisymmetric one, although any of the other considered models can be rejected based on their goodness-of-fit. The magnetic spiral arms have a pitch angle of p ≈ 24° and a tilt angle of χ₀ ≈ 30°, implying a strong vertical field of  ~1 μG at z = 1 kpc. When a radial variation is fitted, the models generally require a fast variation in the inner part of the galaxy (r₁ < 2.5 kpc).

We stress that an accurate determination of the covariance matrix, which accounts for both the noise and the residual astrophysical components, is very complicated. We explored in detail three different methods and found that they may lead to differences of a factor of  ~2−3 in the goodness-of-fit values, while the values of the best-fit parameters for each model do not vary significantly. In practice, this means that rejecting a model based on the goodness-of-fit could be inappropriate, although the relative differences of the χ² statistic between models can be used for a comparison and for selecting preferred models.

We have also tried to constrain the disk field parameters using the polarized synchrotron emission, even though the fitted region is very small (~8% of the sky). Here, the obtained results are remarkably consistent with those obtained with other methods. In this case, all the considered models give a very similar goodness-of-fit, with a (very small) preference for an axisymmetric model. The data do not require a radial dependence of the strength, and the pitch angle values are much lower than in the halo case. Indeed, the constrained pitch angles are compatible with those found by using other observational methods like Faraday rotation of pulsars (Han et al. 2006) However, we note that this conclusion has not been reached by recent results where polarized synchrotron emission is used (see e.g. Jansson et al. 2009). The tilt angle in the disk is χ₀ ~ 17°, which implies a vertical field structure with B_z ~ 0.1 μG at z = 200 pc. This value is compatible with that found by Han & Qiao (1994), based on the rotation measure of pulsars.

We remark that there are still some important uncertainties in the modeling of the synchrotron emission. In particular, we probably need a better knowledge of the distribution of cosmic rays, which is today a clear source of indeterminations in this type of analyses. The detailed modeling of the halo field could have an influence on the cosmic ray trajectories and could be crucial for the direct detection of the primordial magnetic fields.

Finally, we expect that the higher sensitivity and angular resolution in the polarized channels of the PLANCK telescope (The Planck Collaboration 2006; Tauber et al. 2010), the low-frequency channels of the QUIJOTE-CMB experiment (Rubiño-Martín et al. 2010) and experiments like LOFAR and SKA (Beck 2009b), will provide a greater improvement of our knowledge of the Galactic magnetic field.

Acknowledgments

We acknowledge Dr. Marco Tucci and Dr. Juan Betancort-Rijo for useful discussions and Dr. Ángel de Vicente for computing support. We acknowledge the MPIfR “magnetic group” and specially Prof. Richard Wielebinski for help and useful discussions. We acknowledge the use 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. To reduce the computational time needed for the production of the polarized synchrotron maps in this paper, we made use of the Condor workload management system (http://www.cs.wisc.edu/condor/) installed at the IAC. This project has been partially supported by Spanish MEC Grant ESP 2004-06870-C02. This work has been partially funded by project AYA2007-68058-C03-01 of the Spanish Ministry of Science and Innovation (MICINN). JAR-M is a Ramón y Cajal fellow of the MICINN.

References

All Tables

All Figures

	Fig. 1 Observed Stokes Q (top) and U (bottom) maps at 22 GHz from WMAP5 data. The maps are shown in a Mollweide projection of HEALPix in Galactic coordinates. The center of the map corresponds to (l,b) = (0°,0°), and the graticule increases with Δl = Δb = 45°. The maps are degraded to nside = 16 to enhance the large scale pattern. Units are mK (thermodynamic).
In the text

	Fig. 2 Observed direction of the polarization angle (PA) at 22 GHz, obtained from the two maps shown in Fig. 1. The PA is defined here as the local direction of the magnetic field (see Eq. (10)). Units are degrees.
In the text

	Fig. 3 Noise maps for the WMAP5 datasets presented in Figs. 1 and 2. They show the Stokes Q (top), Stokes U (middle) and position angle (PA, bottom) maps. Units for Stokes Q and U maps are mK. Units for the PA map are degrees. The noise map for PA has been computed as fluctuations around the observed WMAP5 K-band map.
In the text

	Fig. 4 Regions used for the definition of the six masks adopted for the analyses (see text for details).
In the text

	Fig. 5 Mask for the galactic polarized emission used by the WMAP team.
In the text

	Fig. 6 Large-scale pattern of the ASS(r) model at z = 4 kpc. This model provides the best fit for the halo field.
In the text

	Fig. 7 One dimensional (marginalized) posterior distribution functions for the parameters of ASS(r) halo model (top: r1, middle: p and bottom: χ0) when mask 4 is considered in the Stokes Q and U analysis.
In the text

	Fig. 8 Best-fit for the halo field by using the Stokes Q and U parameters and PA excluding the disk, the galactic center and the “loops” (ASS model with radial dependence).
In the text

	Fig. 9 Large-scale pattern of the best fit BSS+ model for the disk emission in the galactic disk (z = 0).
In the text