# The Coherent Magnetic Field of the Milky Way

Michael Unger  Institute for Astroparticle Physics (IAP), Karlsruhe Institute
of Technology (KIT), Karlsruhe, Germany  Institutt for fysikk, Norwegian
University of Science and Technology (NTNU), Trondheim, Norway  Glennys R.
Farrar  Center for Cosmology and Particle Physics, Department of Physics, New
York University, New York, NY 10003, USA

###### Abstract

We present a suite of models of the coherent magnetic field of the Galaxy
(GMF) based on new divergence-free parametric functions describing the global
structure of the field. The model parameters are fit to the latest full-sky
Faraday rotation measures of extragalactic sources (RMs) and polarized
synchrotron intensity maps from WMAP and Planck (PI). We employ multiple
models for the density of thermal and cosmic-ray electrons in the Galaxy,
needed to predict the skymaps of RMs and PI for a given GMF model. The
robustness of the inferred properties of the GMF is gauged by studying many
combinations of parametric field models and electron density models.

We determine the pitch angle of the local magnetic field ($(11\pm
1)^{\circ}$), explore the evidence for a grand-design spiral coherent magnetic
field (inconclusive), determine the strength of the toroidal magnetic halo
below and above the disk (same within 10%), set constraints on the half-height
of the cosmic-ray diffusion volume ($\geq 2.9$ kpc), investigate the
compatibility of RM\- and PI-derived magnetic field strengths (compatible
under certain assumptions) and check if the toroidal halo field could be
created by the shear of the poloidal halo due to the differential rotation of
the Galaxy (possibly).

A set of eight models is identified to broadly encompass the present
uncertainties in the coherent GMF – spanning different functional forms, data
products and auxiliary input, and maximizing the differences in their
predictions. We present the corresponding skymaps of rates for axion-photon
conversion in the Galaxy, and deflections of ultra-high energy cosmic rays.
Based on comparing the predicted deflections of UHECRs traversing the Galaxy
in these models, we conclude that the deflection uncertainties are smaller
than the deflections themselves for most of the sky and that the rigidity
(${\mathcal{R}}=E/Z$) threshold for nuclear astronomy can be lowered by a
factor of two if GMF deflections are taken into account.

Galactic magnetic field, Galactic physics, Milky Way, Cosmic rays

††software: The results of this paper have been derived using the software
packages HEALPix (Górski et al., 2005), $\overline{\textrm{Off}}$
$\underline{\textrm{line}}$ (Argiro et al., 2007) and Root (Brun & Rademakers,
1997).

##  1 Introduction

Spiral galaxies are known to be permeated by large-scale magnetic fields, with
energy densities comparable to the turbulent and thermal energy densities of
the interstellar medium; see e.g. Beck (2016) for a recent review. A good
knowledge of the global structure of these fields is important for
understanding their origin, infering their effect on galactic dynamics,
estimating the properties of diffuse motion of low-energy Galactic cosmic rays
and studying the impact of magnetic deflections on the arrival directions of
extragalactic ultrahigh-energy cosmic rays. It is also important for new
physics studies, for instance axion-photon conversion in the GMF or the
interpretation of possible signatures of astrophysical dark matter
annihilation.

The determination of the large-scale structure of the magnetic field of our
Galaxy is particularly challenging since one must infer it from from the
vantage point of Earth, located inside the field. Previous attempts to model
the Galactic magnetic field (GMF) are summarized by Jaffe (2019). In this
paper, we focus on the coherent magnetic field of the Galaxy, leaving the
study of its turbulent component to the near future. Following Jansson &
Farrar (2012a) (hereafter JF12), we derive the GMF by fitting suitably general
parametric models of its structure to the two astrophysical data sets which
are the most constraining of the coherent magnetic fields: the rotation
measures (RMs) of extragalactic polarized radio sources and the polarized
intensity (PI) of the synchrotron emission of cosmic-ray electrons in the
Galaxy.

The relation of these two astrophysical observables to the magnetic field is
detailed in Sec. 2, followed by a description of the RM and PI data in Sec. 3.
The interpretation of this data relies on the knowledge of the three-
dimensional density of thermal electrons and cosmic-ray electrons in the
Galaxy. We discuss these auxiliary models in Sec. 4. The parametric models of
the GMF investigated in this paper are introduced in Sec. 5 and the model
optimization is described in Sec. 6.

The combination of different data sets, auxiliary models and parametric
functions yields an ensemble of GMF models that reflect the uncertainties and
degeneracies inherent in the inference of the global field structure from the
limited information provided by the RM and PI data. In Sec. 7 we narrow down
these model variants to a few benchmark models that encompass the largest
differences within the ensemble.

A typical application of the global magnetic field models is the inference of
the arrival direction of cosmic rays at the edge of the Galaxy and the
determination of the conversion probability of axions in the magnetic field of
the Galaxy. In Sec. 8 we briefly comment on the implcations of our study on
these topics.

We conclude this paper by giving a brief summary of our findings in Sec. 9 and
addressing the question: What is known and not known about the coherent
magnetic field of the Galaxy?

##  2 Observables

In this work, we derive the global structure of the GMF from two astrophysical
observations, Faraday rotation measures and polarized synchrotron intensity.
These two quantities provide complementary information on the parallel and
perpendicular components of the coherent magnetic field along the line of
sight from the observer at Earth through the Galaxy.

The rotation measure (RM) relates the Faraday-rotated polarization angle,
$\chi$, of an astrophysical source observed at a wavelength $\lambda$ to its
intrinsic polarization angle, $\chi_{0}$, after passing through the magnetized
plasma of the Galaxy via

$\chi(\lambda)=\chi_{0}+\text{RM}\,\lambda^{2}.$ (1)

RM is related to the parallel magnetic field via the line-of-sight integral
from the observer to the source (distance $\approx\infty$ for extragalactic
objects),

$\text{RM}=-C_{\text{RM}}\int_{0}^{\infty}n_{e}(\bm{x}(r))\,B_{\parallel}(\bm{x}(r))\,\text{d}r,$
(2)

where the position $\bm{x}$ at distance $r$ from an observer located at
$\bm{x}_{0}$ is $\bm{x}(r)=\bm{x}_{0}+\bm{u}_{r}\,r$111Throughout this work,
we assume an observer at
$\bm{x}_{0}=\bm{x}_{\odot}=(-r_{\odot},\,0,\,z_{\odot})$, in a right-handed
Galactic coordinate system. The vertical distance of the Sun from the Galactic
plane, $z_{\odot}$, is somewhere between $0.006\pm 0.001$ kpc (Joshi et al.,
2016) and $0.026\pm 0.003$ kpc (Majaess et al., 2009) and can thus be
neglected for the purpose of this work, i.e., we set $z_{\odot}\rightarrow 0$.
For the distance of Earth to the Galactic center we adopt the precise
measurement of the Gravity collaboration, $r_{\odot}=8.178\pm 0.025$ kpc
(Abuter et al., 2019)., $n_{e}$ denotes the density of the thermal electron
plasma of the warm ionized medium, $B_{\parallel}=\bm{B}\,\cdot\bm{u}_{r}$ is
the projection of the magnetic field onto the line-of-sight with unit vector
$\bm{u}_{r}$ pointing away from the observer and
$\smash[b]{C_{\text{RM}}=e^{3}/(8\,\varepsilon_{0}\,\pi^{2}\,m_{e}^{2}\,c^{3})\approx
0.8119\,(\nicefrac{{\text{rad}}}{{\text{m}^{2}}})\,(\nicefrac{{\text{cm}^{3}}}{{\text{pc}\,\mu\text{G}}})}$
(e.g. Bradt, 2008). Thus, RM is negative for a magnetic field oriented away
from the observer towards the source and RM is positive if the magnetic field
points from the source to the observer.

The interpretation of RM becomes more complicated if one considers the
possibility of a small-scale correlation or anti-correlation of magnetic
fields and thermal electrons. An anti-correlation could arise if the magnetic
field is in pressure equilibrium with the thermal electron plasma. In this
case, the local Faraday rotation would be systematically diminished and
interpreting the RM assuming no anti-correlation would underestimate the
magnitude of the integrated line-of-sight magnetic field. A positive
correlation could be caused by compression enhancing both the magnetic field
and electron density. These effects were studied by Beck et al. (2003) who
provide the following approximate relation between the rotation measure for
the uncorrelated case, $\text{RM}_{0}$, and the general case with a
correlation coefficient $\kappa$,

$\text{RM}=\text{RM}_{0}\left(1+\frac{2}{3}\,\kappa\,\frac{\langle
b^{2}\rangle}{B^{2}+\langle b^{2}\rangle}\right),$ (3)

where $B$ denotes the coherent field and $\langle b^{2}\rangle$ is the mean
squared field strength of the turbulent magnetic field. In one of our model
variants we will allow for a non-zero $\kappa$.

The polarized synchrotron intensity originates from cosmic-ray electrons and
positrons spiraling in the coherent magnetic field of the Galaxy. The received
intensity depends on the sum of the CR electron-density-weighted synchrotron
emission along the line of sight. Because the emission is strongly forward-
beamed along the helical electron trajectory, the observed emission from a
particular direction is dominated by the electron motion along the magnetic
field component perpendicular to the line of sight. At the magnetic field
strengths and frequencies relevant for this analysis (few $\mu$G and 30 GHz),
the typical cosmic-ray energy responsible for synchrotron emission is of the
order of tens of GeV (e.g. Longair, 2011).

Figure 1: Overview of the RM data used in this analysis. Measurements are
shown in Galactic coordinates using a cylindrical projection that visually
preserves the point density per solid angle. RMs published before 2012 (and
thus available for the JF12 analysis) are shown as gray and black circles.

The relationship between the synchrotron volume emissivity and the magnetic
field strength is particularly simple for the case of electrons with an energy
distribution following a power-law, $n_{\text{cre}}(E)=n_{0}\,E^{-p}$:

$j_{\nu}\varpropto
n_{0}\,\nu^{\frac{-(p-1)}{2}}\,B_{\perp}^{\frac{p+1}{2}}\overset{p=3}{=}n_{0}\,\nu^{-1}\,B_{\perp}^{2},$
(4)

where in the last step the approximation of $p\approx 3$ was used, which is
applicable for a canonical $E^{-2}$ spectrum at the source softened to
$E^{-3}$ due to energy losses from synchrotron cooling and inverse Compton
scattering.

However, the cosmic-ray electron spectrum is not expected to be an exact power
law (e.g. Orlando & Strong, 2013). To obtain the local volume emissivity at
the position $\bm{x}$ we therefore integrate the single-electron emissivity
over energy,

$j_{\nu}(\bm{x})=\int_{0}^{\infty}j(E,B_{\perp}(\bm{x}))\,n_{\rm
cre}(\bm{x},E)\,{\text{d}}E,$ (5)

where the emissivity $j(E)$ of an electron of energy $E$ is given in Eqs.
(8.56) and (8.57) of Longair (2011) and can be conveniently evaluated using
the synchrotron functions provided by the GSL library. The integral over the
emissivities along the line of sight given the three-dimensional distribution
of $n_{\text{cre}}(\bm{x},E)$ and $B_{\perp}(\bm{x})$ yields the Stokes
parameters $Q$ and $U$ of the polarized synchrotron intensity. $Q$ and $U$ add
up quadratically to the total polarized intensity via

$\text{PI}^{2}=Q^{2}+U^{2}$ (6)

and their ratio defines the observed polarization angle

$\psi_{\rm PA}=\frac{1}{2}\,\arctan(U/Q).$ (7)

$\psi$ is perpendicular to the line-of-sight average of the
$n_{\text{cre}}$-weighted magnetic field angle

$\langle\psi_{\rm mag}\rangle=\psi_{\rm PA}+\pi/2,$ (8)

i.e. the analysis of $Q$ and $U$ is sensitive to both the strength and the
direction of perpendicular magnetic along the line of sight.

Isotropic random magnetic fields do not contribute to the polarized intensity,
but an anisotropic random field results in a net-polarized intensity, because
$Q$ and $U$ define only the plane of polarization, and are thus invariant
under a 180∘ rotation of the magnetic field. Anisotropic random magnetic
fields can originate from magneto-hydrodynamic turbulence (Goldreich &
Sridhar, 1995) or from the stretching or compression of a flux-frozen
isotropic random magnetic field (Laing, 1980). The latter can lead to
basically one-dimensional fluctuations along a preferred direction, termed
“ordered random” by Jaffe et al. (2010) or “striated” by Jansson & Farrar
(2012a). A striated random field does not make a net contribution to the RM
since the contributions of opposite-sign regions cancel out. However, it does
contribute to the polarized synchrotron emission since that depends only on
the orientation and not the direction of the field. Here we model the possible
striation via a multiplicative factor to the coherent field,

$B^{\prime}=(1+\xi)B.$ (9)

The striation factor $\xi$ is related to the striation factor $\beta$ used in
JF12, for which a cosmic-ray electron energy spectrum with a spectral index of
$p=3$ was assumed, via $(1+\xi)=\sqrt{1+\beta}$.

##  3 Data Products

###  3.1 Rotation Measures

We assemble a sky map of extragalactic RMs based on individual measurements
from various catalogues and surveys as illustrated in Fig. 1. Most of the
currently known extragalactic RMs were derived by Taylor et al. (2009) from
the two-band polarization data of the NVSS survey (Condon et al., 1998)
leading to 37543 RMs with a declination of $>-40^{\circ}$. Since only two
frequencies are available in this data set, the derived RMs are susceptible to
$n\,\pi$ ambiguities in Eq. 1 for high values of $|\text{RM}|$. We replace the
values of 20 of these high-RM sources by the ones obtained by the broadband
follow-up observations of Ma et al. (2019), three are discarded based on these
re-observations.

For measurements predating the year 2014, we use the catalogue of 4553 high-
quality RMs curated by Xu & Han (2014). Notable surveys with
$N_{\text{RM}}>100$ collected in this catalogue are from Simard-Normandin et
al. (1981), Broten et al. (1988), Brown et al. (2003), Klein et al. (2003),
Broderick et al. (2007), Brown et al. (2007a), Heald et al. (2009), Feain et
al. (2009), Mao et al. (2010), Van Eck et al. (2011), Mao et al. (2012a) and
Mao et al. (2012b).

3220 RMs in the Southern equatorial hemisphere are taken from the S-PASS/ATCA
wide-band radio polarimetry survey (version 0.9) of Schnitzeler et al. (2019)
(applying the quality cuts given in Sec. 4 of that paper and requiring
emission at a single Faraday depth as well as model fits with
$\chi^{2}/\text{ndf}<10$). This survey is of particular significance, since it
complements the sky coverage of the NVSS survey. GMF model fits predating the
S-PASS/ATCA release, e.g. JF12, were constrained by only a very small amount
of RMs in the 18% southernmost equatorial sky.

Further 5999 RMs published in the years 2014-2022 are taken from the data
compilation of Van Eck et al. (2023) (v1.1.0 at Van Eck et al., 2022). Notable
surveys with $N_{\text{RM}}>100$ in this compilation are from Farnes et al.
(2014), Anderson et al. (2015), O’Sullivan et al. (2017), Kaczmarek et al.
(2017), Betti et al. (2019), Riseley et al. (2020), Ma et al. (2020) and Van
Eck et al. (2021). In addition, we include the recent 2461 extragalactic RMs
of O’Sullivan et al. (2023).

(a) Sky map of extragalactic rotation measures (44857 RMs averaged over
$N_{\text{side}}=32$ HEALPix pixels). The color scale is saturated at
$|\text{RM}|\geq\,150\text{rad}/\text{m}^{2}$)

(b) Masked region for the RM analysis. Selected astrophysical objects (red)
and HII regions (blue). The gray-scale background shows the sky map of
emission measures.

Figure 2: Sky map of rotation measures (left) and masked regions (right).

In this data set of 53773 RMs, we identify obvious multiple measurements of
the same extragalactic object if their coordinates correspond to the same
pixel id on a high-resolution HEALPix sky map with
$N_{\text{side}}=2048$222The HEALPix package (Hierarchical, Equal Area, and
iso-Latitude Pixelation of the sphere) is used to subdivide the sky into
equal-area pixels (Górski et al., 2005). The resolution, $N_{\text{side}}$, of
a HEALPix map is related to the number of pixels via
$N_{\text{pix}}=\smash[b]{12\,N_{\text{side}}^{2}}$. For a map with
$N_{\text{side}}=16$, the standard resolution used in this paper, the number
of pixels is $N_{\text{pix}}=3072$ with an angular width of approximately
$\smash[b]{\theta\approx\sqrt{4\pi/N_{\text{pix}}}\;180^{\circ}/\pi=\sqrt{3/\pi}\;60^{\circ}/N_{\text{side}}=3.66^{\circ}}$.
All duplicate RMs but the measurement with the best frequency coverage are
dropped, leaving 47054 RMs of unique extragalactic objects.

In the next step of data selection, we apply a two-pass algorithm to reject
outliers caused by either a large source-intrinsic RM or a wrong resolution of
the $n\pi$-ambiguity. The procedure is similar in spirit but different in
detail to that used for the data set of JF12. For each RM measurement we
accumulate surrounding measurements within an angular collection radius
starting at 1 degree. The collection radius is increased until a sample of at
least 10 independent surrounding RM measurements are found that differ by less
than three median absolute deviations from the median RM of that sample. In
the second pass, these sources are then used to calculate the mean and
standard deviation of the RM in the region surrounding each object. We keep
objects if their RM value is closer than three standard deviations from the
surrounding mean.

The remaining 44857 RMs are displayed in Fig. 2a. For this figure, we show the
average of RMs in an $N_{\text{side}}=32$ HEALPix map333This and subsequent
maps shown in this paper are Mollweide projections of the sky in Galactic
coordinates with the Galactic center at the origin. The longitude increases
towards the left from $0^{\circ}$ at the center towards $+180^{\circ}$ and
decreases towards the right from $+360^{\circ}$ at the center to
$+180^{\circ}$ (or, equivalently, from $0^{\circ}$ to $-180^{\circ}$). , in-
painting 1524 pixels (12%) without a measurement. At this somewhat higher
resolution than the one used for the analysis (see below), both the large and
small-scale features of the extragalactic rotation measures are visible. The
most salient large-scale features of the RM sky are known since the early
advent of RM catalogues, see e.g. the discussion in Andreasyan & Makarov
(1988): For the inner Galaxy (longitudes between -90 to 90∘) there is an anti-
symmetry in both longitude and latitude of the sign of the average rotation
measures, whereas in the outer Galaxy (longitudes from 90 to 270∘), the
average RMs have the same sign above and below the plane, as summarized by the
following cartoon of an RM sky map,

$\text{RM}(l,b)\approx\mathchoice{\leavevmode\hbox to25.23pt{\vbox
to22.08pt{\pgfpicture\makeatletter\hbox{\hskip 12.6138pt\lower-11.04227pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{12.4138pt}{0.0pt}\pgfsys@curveto{12.4138pt}{5.98808pt}{6.85603pt}{10.84227pt}{0.0pt}{10.84227pt}\pgfsys@curveto{-6.85603pt}{10.84227pt}{-12.4138pt}{5.98808pt}{-12.4138pt}{0.0pt}\pgfsys@curveto{-12.4138pt}{-5.98808pt}{-6.85603pt}{-10.84227pt}{0.0pt}{-10.84227pt}\pgfsys@curveto{6.85603pt}{-10.84227pt}{12.4138pt}{-5.98808pt}{12.4138pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-7.77788pt}{-3.33333pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\\!\\!\frac{-+-+}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to15.8pt{\vbox to16.43pt{\pgfpicture\makeatletter\hbox{\hskip
7.89973pt\lower-8.21384pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{7.69974pt}{0.0pt}\pgfsys@curveto{7.69974pt}{4.42596pt}{4.2525pt}{8.01384pt}{0.0pt}{8.01384pt}\pgfsys@curveto{-4.2525pt}{8.01384pt}{-7.69974pt}{4.42596pt}{-7.69974pt}{0.0pt}\pgfsys@curveto{-7.69974pt}{-4.42596pt}{-4.2525pt}{-8.01384pt}{0.0pt}{-8.01384pt}\pgfsys@curveto{4.2525pt}{-8.01384pt}{7.69974pt}{-4.42596pt}{7.69974pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.44453pt}{-2.33333pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\\!\\!\frac{-+-+}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to9.51pt{\vbox to12.66pt{\pgfpicture\makeatletter\hbox{\hskip
4.75702pt\lower-6.32823pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{4.55702pt}{0.0pt}\pgfsys@curveto{4.55702pt}{3.38457pt}{2.5168pt}{6.12823pt}{0.0pt}{6.12823pt}\pgfsys@curveto{-2.5168pt}{6.12823pt}{-4.55702pt}{3.38457pt}{-4.55702pt}{0.0pt}\pgfsys@curveto{-4.55702pt}{-3.38457pt}{-2.5168pt}{-6.12823pt}{0.0pt}{-6.12823pt}\pgfsys@curveto{2.5168pt}{-6.12823pt}{4.55702pt}{-3.38457pt}{4.55702pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.2223pt}{-1.66666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\\!\\!\frac{-+-+}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to9.51pt{\vbox to12.66pt{\pgfpicture\makeatletter\hbox{\hskip
4.75702pt\lower-6.32823pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{4.55702pt}{0.0pt}\pgfsys@curveto{4.55702pt}{3.38457pt}{2.5168pt}{6.12823pt}{0.0pt}{6.12823pt}\pgfsys@curveto{-2.5168pt}{6.12823pt}{-4.55702pt}{3.38457pt}{-4.55702pt}{0.0pt}\pgfsys@curveto{-4.55702pt}{-3.38457pt}{-2.5168pt}{-6.12823pt}{0.0pt}{-6.12823pt}\pgfsys@curveto{2.5168pt}{-6.12823pt}{4.55702pt}{-3.38457pt}{4.55702pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.2223pt}{-1.66666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\\!\\!\frac{-+-+}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\,.$ (10)

In addition, many small-scale features can be identified, some of which can be
attributed to foreground objects. We remove some of the most obvious regions,
indicated by circles in Fig. 2b. Firstly, we discard lines of sight if they
pass through magnetized objects with a large angular size, i.e. M31, the Small
or Large Magellanic Cloud, Centaurus A (core and lobes) and the Galactic
Center. The excluded regions around these objects are shown in red in Fig. 2b.
Secondly, we deselect RMs if the thermal electron density along the line of
sight is dominated by a single object. Such local overdensities in the diffuse
WIM are caused by the ultraviolet light of young massive stars ionizing the
ISM around them, creating an HII region. RMs in lines of sights passing
through HII regions are often dominated by the product
$n_{e}B_{\parallel}\Delta$ inside the region of thickness $\Delta$, see e.g.
Heiles et al. (1981) and Harvey-Smith et al. (2011). These RMs are thus not
representative of the large-scale GMF. Regions of locally enhanced thermal
electron density can e.g. be identified by the emission measure,
$\text{EM}=\int_{0}^{\infty}n_{e}\,{\text{d}}l$, which is displayed in Fig. 2b
based on the composite of H$\alpha$ surveys (VTSS, SHASSA and WHAM) from
Finkbeiner (2003) . More than 8000 Galactic HII regions are known in the
Galaxy (Anderson et al., 2014), but most of them are either far away, and thus
of small angular extent, and/or at low Galactic latitudes, where the overall
RM integral through the plane is large enough that it is not overshadowed by
the contribution of a single HII region. We therefore deselect only the lines
of sight that overlap with an HII region if it is either of large angular
extent ($r>10^{\circ}$) or at high latitudes, $|b|>5^{\circ}$. These regions
are displayed in blue in Fig. 2b if their size is $>0.5^{\circ}$.

The final 41686 RMs are then binned in 3072 angular pixels of an
$N_{\text{side}}=16$ HEALPix map with an angular diameter of $\sim
3.7^{\circ}$. For each pixel, we calculate the mean and unbiased sample
variance to be used in the model optimization (see Eq. 49 below). If the
number of RMs within one pixel is $<10$, then the variance (but not the mean)
is calculated using the 10 RMs with the closest angular distance to the center
direction of the pixel. This leads to 2838 pixels with RM data, where 46
pixels are excluded because no data was observed in this direction and another
188 pixels are excluded by the masks displayed in Fig. 2b. The final masked
and binned RM data are shown in Fig. 10 in Sec. 7.

###  3.2 Polarized Synchrotron Emission

The polarized synchrotron emission from the Galaxy is best observed at high
frequencies, where the depolarization of the signal due to Faraday rotation is
negligible. At Faraday depths typical for the Galaxy, the change in
polarization angle becomes negligible in the tens of GHz range, c.f. Eq. 1.
The received polarized intensity in units of Rayleigh-Jeans antenna
temperature, $T=c^{2}/(2\nu^{2}\,k_{\text{B}})\,\text{PI}$, decreases as
$\nu^{\beta_{\text{S}}}$ with frequency. Here the synchrotron spectral index
$\beta_{\text{S}}\approx-3$ for a cosmic-ray electron spectrum with index
$p=3$, see Eq. 4. Another source of polarized emission from the Galaxy
originates from thermal dust with a spectral index of
$\beta_{\text{D}}\approx+1.6$ (Ade et al., 2015). Empirically, the cross-over
between the two components is at around 100 GHz.

The “sweet spot” for the observation of synchrotron emission is therefore at
high enough frequencies such that Faraday depolarization is negligible, but
not too high, such that the signal is not dominated by polarized dust
emission. The Wmap and Planck satellites operated in this frequency regime,
with a frequency threshold of 20 and 30 GHz respectively.

From Wmap we use the final nine-year results (DR5) on the polarized
synchrotron emission at 22.5 GHz provided at a resolution of
$N_{\text{side}}=64$ (Bennett et al., 2013). Four variants of $Q$ and $U$ sky
maps are available, differing mainly by the constraints placed on
$\beta_{\text{S}}$ during analysis. These four variants of the derived
polarized synchrotron intensity are in good agreement outside of the Galactic
plane, which is anyway masked in our analysis (see below). We use the “base”
model, which is the most data-driven variant in which the synchrotron spectral
index was allowed to float freely in each sky pixel. We in-paint the pixels
that are flagged as having an erroneous component separation (2.7%) with the
$Q$ and $U$ averages of their eight surrounding neighbor pixels.

Figure 3: Stokes $Q$ (top panel) and $U$ (middle panel) parameters and
polarized intensity PI (bottom panel) from the averaged synchrotron maps of
Planck and Wmap at 30 GHz.

From Planck we use the third release (R3.0) of the polarized synchrotron
foreground at 30 GHz derived from component separation using the Commander
software (Akrami et al., 2020). These sky maps are provided at high resolution
of $N_{\text{side}}=2048$ and we average the $Q$ and $U$ values to obtain maps
with $N_{\text{side}}=64$ to match the resolution of Wmap. During this
averaging, spurious outliers in the high-resolution maps are identified and
discarded if the value is more than $5~{}\sigma$ away from the median of the
values within the low-resolution pixel. This procedure removes 0.18% of the
high-resolution data points.

We then combine the $Q$ and $U$ sky maps from Wmap and Planck into a lower
resolution map at $N_{\text{side}}=16$ by taking the simple arithmetic average
of all 32 measurements in each pixel (16 from each experiment). Wmap
intensities are extrapolated to the Planck frequency with a spatially constant
synchrotron spectral index of $\beta_{\text{S}}=-3.15$ based on the mean value
predicted by our simulations described in Sec. 4.2. Note that due to the
proximity of the Wmap and Planck frequencies, even a large difference between
the actual and assumed of $\Delta\beta_{\text{S}}=0.2$ would introduce a
variation in extrapolated intensity of only
$1-(22.5~{}\text{GHz}/30~{}\text{GHz})^{\pm\Delta\beta_{\text{S}}}=\pm 6\%$.
The variance of the 32 data points per pixel is used to calculate the weights
in the $\chi^{2}$ minimization as discussed in Sec. 6.

Figure 4: Polarized synchrotron mask used in this analysis. Masked regions are
shown in red, blue and green (see text). For comparison, the outline of the
Fermi and eROSITA bubbles are shown in violet and orange.

The Stokes parameters derived by the Wmap and Planck collaborations exhibit
systematic large-scale differences, see Fig. 22 in the appendix. These
differences could either be attributed to large-scale spatial variation of the
synchrotron spectral index or to residual calibration uncertainties in one or
both of the data sets. A combined analysis of Wmap and Planck data was
recently presented by the CosmoGlobe Collaboration (Watts et al., 2023a, b),
performing the component separation and the calibration of the data sets
simultaneously. We will use the CosmoGlobe results as an alternative to our
average and study the effect on the parameters of GMF models in Sec. 7.

Our combined maps of Stokes $Q$ and $U$ and polarized intensity PI of the
Galactic synchrotron emission at 30 GHz are shown in Fig. 3. As can be seen,
both $Q$ and $U$ exhibit large-scale features, that we will interpret in the
following as imprints of the large-scale features of GMF. However, care must
be taken not to include regions with strong local features that could bias our
fits. We therefore mask-out some regions of the sky in fitting to the GMF
parameters.

The elements of our mask when fitting the polarized intensity data are shown
in Fig. 4. Red indicates the map used in the original Wmap analysis. It
excludes regions of high polarized intensity along the Galactic plane and the
North Polar Spur, and pixels containing strong extragalactic sources such as
the radio lobes of Centaurus A. The blue region depicts the additional
polarized intensity mask introduced by Jansson & Farrar (2012a) mainly to
remove a presumably local, high-latitude polarized emission at
$90^{\circ}<l<180^{\circ}$. Finally, we also mask the green regions as an
attempt to remove further large circular arcs (called loops, spurs or
filaments) visible in PI. Here we exclude data in the direction of loops I to
IV, as defined by Berkhuijsen et al. (1971) and further filaments identified
by Vidal et al. (2015). The union of the above leads to our final PI mask. The
cumulative application of these three masks leaves 73.1% (Wmap), 63.4% (Wmap
+JF12) and 57.8% (Wmap +JF12+loops) of the sky for analysis.

It is worthwhile noting that the exact attribution of features in the
polarized radio sky to local or global phenomena in the Galaxy is still under
debate, see e.g. Lallement (2022). Some of the loops and filaments could be
caused by local supernova remnants expanding into the surrounding ambient
magnetic field (Spoelstra, 1973), whereas others might be related to large-
scale magnetized outflow from the Galactic center (Carretti et al., 2013),
related to the so-called “Fermi bubbles” observed in gamma rays (Ackermann et
al., 2014) and surrounded by “eROSITA bubbles” in X-rays (Predehl et al.,
2020). These are shown as violet and orange regions in Fig. 4, demonstrating
that our mask retains some directions which contain contributions from the
Fermi and eROSITA bubbles.

##  4 Auxiliary Models

###  4.1 Thermal Electrons

Figure 5: Left: Thermal electron density, $n_{e}$, for the Ymw16 and Ne2001
model. The top panel shows the mid-plane density at $z=0$ and the bottom panel
gives an edge-on view of the density at $y=0$ in the $x-z$ plane. Structures
visible at around the sun position ($\sim(-8.2,0,0)$ kpc) are due to the
modeled under- and over-densities of the local interstellar medium. Right:
Thermal electron density along the Galactic $x$-axis at $z=y=0$ (top) and
along the $z$-axis at $x=-6$ kpc and $y=0$.

The magnetized plasma responsible for the Faraday rotation of extragalactic
radio sources resides mostly in the warm ionized medium (WIM) of the Galaxy
(e.g. Ryden & Pogge, 2021). To interpret the rotation measures, a three-
dimensional model of the density, $n_{e}$, of these free thermal electrons is
needed. Here we use two models of $n_{e}$: Ne2001 of Cordes & Lazio (2002) and
Ymw16 of Yao et al. (2017). Both models were tuned to describe the dispersion
measure, DM, of Galactic pulsars, that is given by the line-of-sight integral
from Earth to the pulsar at a distance $d$,
$\text{DM}=\int_{0}^{d}n_{e}(\bm{x}(r))\,{\text{d}}r$. More than 3000 pulsars
with measured DMs are listed in the current version (1.70) of the ATNF pulsar
catalogue (Manchester et al., 2005), but only for a few the distance $d$ is
well known: Ne2001 was tuned to the DMs of 112 pulsars and Ymw16 used 189.
Given the large number of parameters of these models (e.g. the Ymw16 model has
82 fixed parameters and 32 fitted parameters) and the scarcity of data, they
rely to a large extent on astrophysical priors for the geometrical topology of
the thermal electron density.

Particularly important for the modeling of the large-scale structure of the
GMF, the vertical structure of the Galactic WIM is relatively well constrained
by the DMs of pulsars in high-latitude globular clusters Gaensler et al.
(2008). Since the fit of Ne2001 pre-dates most of this data, we replace the
original value of the exponential scale height of the diffuse WIM of
$h_{\text{WIM}}=0.95$ kpc with $h_{\text{WIM}}=1.3\pm 0.2$ kpc as derived by
Schnitzeler (2012) for this model. A larger value of $h_{\text{WIM}}=1.67\pm
0.05$ kpc was inferred by Yao et al. (2017) for their Ymw16 model.

An illustration of the thermal electron densities of Ymw16 and Ne2001 is shown
in Fig. 5. As can be seen, the two models differ substantially, especially
regarding the positions and widths of spiral arms, the density in the
molecular ring at $r\sim 5$ kpc and the Galactic Center and concerning the
scale-height of the thick disk of the WIM. Whereas Ymw16 can predict pulsar
DMs with a somewhat higher fidelity than Ne2001, the latter model remains a
viable alternative to describe the large-scale features of the WIM (e.g. Price
et al., 2021). By using both models in our GMF fits, we can study the
systematic effects arising from different assumptions on the thermal electron
density.

Figure 6: Left: Flux of cosmic-ray electrons and positrons at Earth. Black
points show the measurements. The calculated flux in the local interstellar
medium (LISM), for the $h_{D}=6$ kpc case, is shown as a blue band and the
contribution from secondary production
($\text{p}_{\text{CR}}+\text{p}_{\text{ISM}}$ interactions) is indicated by a
dashed line; the solar-modulated flux at Earth is shown as a red line. The
additional $x$-axis at the top of the figure shows the magnetic field strength
needed to obtain a synchrotron peak at $\nu_{c}=30$ GHz (i.e., the frequency
for which the Planck collaboration reports the polarized synchrotron
emission), for the electron energy given in the $x$-axis at the bottom of the
figure. Middle: Density of 10 GV cosmic-ray electrons at a galactocentric
radius $r=4$ kpc, as a function of Galactic height $z$ and Right: as a
function of radius $r$ at mid-plane ($z=0$). Different choices of the half-
height $h_{D}$ of the diffusion volume are shown as lines with different
colors. The dashed line is the result of Orlando & Strong (2013) for
$h_{D}=10$ kpc.

###  4.2 Cosmic Ray Electrons

Calculating Galactic synchrotron emission requires not only a model of the GMF
but also a model of the three-dimensional density distribution and energy
spectrum of cosmic-ray electrons in the Galaxy, $n_{\text{cre}}({\bf x})$.

In contrast to the purely phenomenological thermal electron models discussed
in the previous section, predictions of $n_{\text{cre}}$ are based on detailed
modeling of the production and propagation of electrons and positrons in the
Galaxy (e.g. Moskalenko & Strong, 1998), to obtain a steady-state solution of
the diffusion equation of cosmic rays in the Galaxy (Ginzburg & Syrovatskii,
1964). In its simplest version (sometimes referred to as “plain diffusion”),
the defining quantity of the model is the diffusion coefficient of charged
particles in the interstellar medium. Most calculations assume a homogeneous
and isotropic Galactic diffusion coefficient $D$ within a “diffusion volume”
approximated as a cylinder of half-height $h_{D}$.

The normalization and rigidity444The rigidity of a particle with charge $Z\,e$
and momentum $p$ (energy $E$) is ${\mathcal{R}}=p\,c/(Z\,e)\simeq E/(Z\,e)$
dependence of $D$ is determined from measurements of the fluxes of secondary
and primary cosmic-ray nuclei at Earth. The most precise estimates (in terms
of the uncertainty of both cosmic-ray flux and nuclear cross sections) are
derived from the ratio of the flux of secondary boron nuclei and (mostly)
primary carbon nuclei. However, secondary-to-primary ratios can only constrain
the ratio, $r_{D}=D/h_{D}$, (e.g. Maurin et al., 2001). Estimates of $r_{D}$
in a plain diffusion scenario range from 0.03 kpc/Myr (Génolini et al., 2019;
Yuan et al., 2017) to 0.10 kpc/Myr (Cummings et al., 2016) at a reference
rigidity of 10 GV.

The degeneracy between $h_{D}$ and $D$ can in principle be broken by data on
“cosmic clocks” (e.g. 10Be/9Be or Be$/$C), but due to the poor quality of
current data and uncertainties in the spallation cross sections, only mild
constraints on the halo height can be derived and the current estimates of
$h_{D}$ are in the range 2 to 10 kpc (e.g. Weinrich et al., 2020; Evoli et
al., 2020; Maurin et al., 2022). We therefore consider diffusion volumes
having $h_{D}=2,4,6,8$ and 10 kpc and derive the corresponding
$n_{\text{cre}}$ for each, to allow us to assess the uncertainty in the GMF
due to the present uncertainty in $h_{D}$. For each of the five different
values of the height of the diffusion volume, $h_{D}$, we solve the cosmic-ray
diffusion equation with the Dragon program (Evoli et al., 2008)555We choose
Dragon over the other widely used solver Galprop (Strong & Moskalenko, 1998)
because at the time of setting up our analysis chain only Dragon supported
spatially varying diffusion coefficients. In this paper we, however,
investigate only homogenous diffusion..

For each value of $h_{D}$ we obtain the diffusion coefficient at 10 GV from
$r_{D}=D/h_{D}=0.046$ kpc/Myr and we take the rigidity dependence of
$D({\mathcal{R}})$ to be the one derived for the “PD2” model of Cummings et
al. (2016). At high rigidities (${\mathcal{R}}\gtrsim 5~{}\text{GV}$) the
diffusion coefficient scales as $\smash[b]{D\varpropto{\mathcal{R}}^{-0.578\pm
0.073}}$, i.e., compatible with the power-law scaling
$\smash[b]{{\mathcal{R}}^{-\nicefrac{{1}}{{2}}}}$ typical for a turbulent
magneto-hydrodynamical cascade (Iroshnikov, 1964; Kraichnan, 1965). The
spatial distribution of Galactic cosmic-ray sources is taken to follow the
radial distribution of pulsars in the Galaxy (Lorimer et al., 2006), used as a
proxy for the distribution of supernova remnants, and we use the JF12 magnetic
field model to calculate the cooling of cosmic-ray electrons due to
synchrotron radiation. The maximum galactocentric radius of the diffusion
volume was set to $R_{\text{max}}=18$ kpc.

We constrain the solutions of the diffusion equations by the measurements of
cosmic-ray proton and lepton fluxes, in the local interstellar medium from the
Voyager I satellite (Cummings et al., 2016), and inside the heliosphere from
AMS-02 orbiting Earth on the International Space Station (Aguilar et al.,
2014, 2015, 2019). We first find the proton spectrum at the source; this then
fixes the contribution to the lepton flux at Earth coming from proton
interactions with cosmic rays and the ISM
($\text{p}_{\text{CR}}+\text{p}_{\text{ISM}}$), accounting for a small
adjustment due to contributions from nuclei. We then attribute the remaining
lepton flux to the injected lepton spectrum. Thus, in this simplified ansatz,
the well-known “positron anomaly” (Adriani et al., 2009), i.e., the excess of
positron flux beyond expectations from secondary production at energies
$\gtrsim 10$ GeV, is attributed to the injected lepton spectrum. Such an
ansatz is plausible if astrophysical sources of primary positrons have a
similar spatial distribution as the sources of primary electrons, e.g., if
pulsars are the sources of the “anomalous” positrons and electrons (Hooper et
al., 2009) and the bulk of electrons is accelerated in supernova remnants.
Then, if the interest is only in the sum of electrons and positrons as it is
in our case, the two sources can be lumped together into a single source
class. Due to the fast cooling time of electrons and positrons at high
energies, the high energy flux at Earth might be dominated by a local source
and thus not representative of the average lepton flux in the Galaxy (e.g. Di
Mauro et al., 2014; Joshi & Razzaque, 2017; Mertsch, 2018). However these
effects are expected to play a role only above 100 GeV and are thus not
important for the synchrotron frequencies of interest in this study. An
example of our fit of the cosmic-ray lepton flux is shown in the left panel of
Fig. 6.

As a further systematic check, we also interpret the synchrotron data using
the z10LMPDE cosmic-ray electron model, a Galprop calculation of Orlando &
Strong (2013) (OS13) used in (Adam et al., 2016) for GMF modeling (Jaffe,
2017). This model assumes a half-cylinder height of $h_{D}=10$ kpc, a maximum
radius of $R_{\text{max}}=20$ kpc and the source distribution from Strong &
Moskalenko (1998) with a radial cutoff at 15 kpc. Further differences to our
fiducial Galprop calculation include a different magnetic field (Sun & Reich,
2010) for electron cooling and a different value of $r_{D}=0.031$ at 10 GV. We
normalized this model to match the AMS flux of electrons and positrons at 35
GV.

The vertical and radial distribution of the cosmic-ray lepton density at 10 GV
for different values of $h_{D}$ and the OS13 model are shown in the middle and
right panel of Fig. 6. As expected, a larger half-height of the diffusion
volume leads to a larger cosmic-ray occupation in the halo, while the density
in the disk is approximately constant due to the normalization at Earth. The
different source distribution assumed for OS13 results is mainly responsible
for the very different radial distribution of cosmic ray leptons.

Of course, further variations beyond the scale height are possible and we
studied the effect on synchrotron maps of several other type of variations, as
described below. However only the variation of $h_{D}$ makes a significant
change in the PI predictions outside of the mask we use, so we only include
the six $n_{\text{cre}}$ models outlined above in our GMF model
fitting.666Note that in the framework of a tunable striation parameter, an
overall factor in the polarized synchrotron intensity does not affect the
derived coherent magnetic field, but only the striation factor.. A spiral
distribution of sources can affect the PI up to 25% in the Galactic plane, but
has negligible effects outside of the polarization mask used in this work.
Using the Planck-tune of the JF12 random field (Adam et al., 2016) for
synchrotron cooling instead of the original model (Jansson & Farrar, 2012b)
affects the intensity by up to 16% in the outer Galaxy (where a spiral arm
with a large random field is present in the original model to describe the
intensity from the “fan region”), but the differences are again negligible
outside the masked region used in this paper. The same holds for
$n_{\text{cre}}$ models using a three-dimensional model of the interstellar
radiation field from Porter et al. (2017). Moreover, in this study we did not
iteratively re-adjust the magnetic field used for the synchrotron energy
losses in the $n_{\text{cre}}$ calculations. But the resulting refitted
coherent fields are of the same order of magntitude as the ones of JF12 and
therefore the electron cooling is dominated by the sum of contributions from
the random magnetic field and inverse Compton scattering. Finally, it was
shown by Orlando (2018) that the cosmic-ray electron spectrum above a few GeV
is insensitive to the inclusion of re-acceleration and/or convection in the
diffusion equation, therefore we do not consider it here.

The six models shown in Fig. 6 will be used in this paper to study the
systematic effect of the three-dimensional model of $n_{\text{cre}}$ on the
derived structure of the GMF.

We note, that future systematic studies should include an investigation of the
effects of a more realistic particle propagation with an inhomogenous and/or
anisotropic diffusion coefficent (e.g. Di Bernardo et al., 2013; Giacinti et
al., 2018; Merten et al., 2017; AL-Zetoun & Achterberg, 2018), ideally with a
self-consistent description of relation between the magnetic field and the
diffusion coefficient (e.g. Kuhlen et al., 2022; Blasi, 2023).

##  5 Magnetic Field Models

Figure 7: Left: Fixed-arm spiral disk field of the JF12 model. Right: Fourier
spiral of the “neCL” model, cf. Sec. 7. The position of the sun is indicated
as a white star.

Inspired by previous models of the coherent Galactic magnetic field (e.g.
Stanev, 1997; Prouza & Šmída, 2003; Sun et al., 2008; Pshirkov et al., 2011;
Jansson & Farrar, 2012a) we describe its global structure as a superposition
of a spiral field close to the Galactic disk and a large-scale halo field. The
halo field is composed of a toroidal and a poloidal field as introduced by
Jansson & Farrar (2012a).

In the following, we give a brief overview of the different variants of
parametric descriptions investigated in this paper. The necessity of each
component to describe large-scale features of the RM and PI sky will be
discussed in Sec. 7.

###  5.1 Disk Field

####  5.1.1 Fourier Spiral

Observations of synchrotron polarization patterns in face-on external spiral
galaxies and the RMs of Galactic pulsars suggest that the coherent field close
to the disk of the Milky Way follows a spiral pattern which can be
approximated as a logarithmic spiral with pitch angle $\alpha$. In previous
GMF models the pitch angle was often fixed to a value close to the one of the
spiral matter density of the Milky Way and the magnetic field is assumed to be
organized in “arms” with a pre-defined geometry. The strength and direction
within these arms are adjusted to match the data. Within one arm, the field is
approximated to be azimuthally constant, leading to strong discontinuities at
the boundaries of the arms (see e.g. Brown et al., 2007b; Sun et al., 2008;
Jaffe et al., 2010; Pshirkov et al., 2011; Jansson & Farrar, 2012a; Han et
al., 2018).

Here we adopt a new description of the disk field that provides the
flexibility to fit the geometry of the magnetic arms and that avoids
discontinuities between them. For this purpose, we decompose the magnetic
field strength at reference radius $r_{0}$ and azimuthal angle $\phi_{0}$ into
$n$ modes of strength $B_{m}$ and phase $\phi_{m}$,

$B(r_{0},\phi_{0})=\sum_{m=1}^{n}B_{m}\,\cos\left(m(\phi_{0}-\phi_{m})\right).$
(11)

The magnetic field at position $(r,\phi,z)$ is given by

$\bm{B}_{\text{d}}=\left(\sin\alpha,\cos\alpha,0\right)\,\frac{r_{0}}{r}\,B(r_{0},\phi_{0})\,h_{\text{d}}(z)\,g_{\text{d}}(r),$
(12)

where the corresponding $B(r_{0},\phi_{0})$ can be found by following the
field line along the logarithmic spiral from $(r,\phi)$ to the reference
radius $r_{0}$ at which

$\phi_{0}=\phi-\log(r/r_{0})/\tan\alpha.$ (13)

Due to the expansion of $B(r_{0},\phi_{0})$ in terms of a cosine series, we
will refer to this model of the disk field as “Fourier spiral” in the
following. The solenoidality of this field model is assured because for each
mode in Eq. 11 the same amount of flux enters and exits along the circle at
$r_{0}$ and, due to the $\frac{r_{0}}{r}$ factor in Eq. 12, the magnetic flux
in each arm is constant as a function of radius.

The functions $h_{\text{d}}(z)$ and $g_{\text{d}}(r)$ describe the fade-in and
fade-out of the field in the vertical and radial direction, respectively. We
choose the ansatz

$h_{\text{d}}(z)=1-\operatorname{\mbox{$\sigma$}}\left(\frac{|z|-z_{\text{d}}}{w_{\text{d}}}\right)$
(14)

and

$g_{\text{d}}(r)=\left(1-\operatorname{\mbox{$\sigma$}}\Big{(}\frac{r-r_{\text{2}}}{w_{\text{2}}}\Big{)}\right)\,\operatorname{\mbox{$\sigma$}}\Big{(}\frac{r-r_{\text{1}}}{w_{\text{1}}}\Big{)}\left(1-e^{-r^{2}}\right),$
(15)

where $\operatorname{\mbox{$\sigma$}}(x)$ denotes the logistic sigmoid
function,

$\operatorname{\mbox{$\sigma$}}(x)=\frac{1}{1-e^{-x}},$ (16)

such that the disk field is suppressed to half of its value at $z_{\text{d}}$,
$r_{\text{1}}$ and $r_{\text{2}}$, and the suppression rate is given by the
corresponding transition widths. The additional factor $(1-e^{-r^{2}})$ is
needed to assure that the factor $g_{\text{d}}(r)/r$ in Eq. 12 goes to 0 at
$r=0$.

Note that mathematically $g_{\text{d}}(r)$ violates the solenoidality of the
field and is thus to be understood as an effective modeling of the behavior of
the disk field at small and large radii. We conjecture that the coherent
magnetic flux is dissipated into a random field when it reaches the
circumgalactic medium at a distance $r_{\text{2}}$. In the inner Galaxy, at
around $r_{\text{1}}$, it is plausible that the coherent structured disk field
(whose net radial flux is zero) is randomized by the turbulent motions of the
ISM in the molecular ring and Galactic bar such that locally the coherent
field disappears.777Further possibilities include that the magnetic flux of
the disk field exits vertically in the inner Galaxy (Ferrière & Terral, 2014)
or that the inward- and outward-going field lines of the disk field connect at
$r_{\text{1}}$ and $r_{\text{2}}$ as investigated by Kleimann et al. (2019).
These entail regions of very high coherent field which on energetic grounds
would be disfavored, so we did not pursue them.

In this analysis, we set $r_{\text{1}}$ to 5 kpc, corresponding to the outer
radius of the molecular ring in the Ymw16 model, and fix $r_{\text{2}}$ to 20
kpc. For the transition widths, we adopt $w_{\text{1}}=w_{\text{2}}=0.5$ kpc.
Without loss of generality, we set the reference radius to $r_{0}=5$ kpc, such
that the coefficients $B_{m}$ in Eq. 11 denote the amplitude of the modes in
the inner Galaxy.

An example of a Fourier spiral disk field is shown in the right panel of Fig.
7 and for comparison a fixed-arm spiral field is displayed in the left panel.
Apart from obvious differences in the arm topologies and field strengths,
which are mostly due to different data sets used to fit these models, it can
be seen that the Fourier spiral results in a smooth disk field without the
discontinuities inherent in the previous fixed-arm models.

####  5.1.2 Spiral Spur

As an alternative model to the grand-design magnetic spiral implied by the
Fourier spiral, we investigate if a localized spiral segment can describe the
data equally well. Such a magnetic “spiral spur” can be modeled as a Gaussian
of width $w_{\text{S}}$ at reference radius $r_{0}$

$B(r_{0},\phi_{0})=B_{\text{1}}\,\exp\left(-\frac{1}{2}(\Delta(\phi_{0},\phi_{\text{1}})/w_{\text{S}})^{2}\right),$
(17)

where $\Delta$ denotes the azimuthal difference between $\phi_{0}$ and the
center of the spur $\phi_{\text{1}}$, i.e. the angle between the unit vectors
$\bm{u}_{0}=(\cos\phi_{0},\,\sin\phi_{0})$ and
$\bm{u}_{\text{S}}=(\cos\phi_{\text{1}},\,\sin\phi_{\text{1}})$ with
$\cos\Delta=\bm{u}_{0}\,\bm{u}_{\text{S}}$. For a point at $(r,\phi)$, the
angle $\phi_{0}$ follows again from the logarithmic spiral via Eq. 13 and the
field is similar to Eq. 12 as

$\bm{B}_{\text{S}}=\left(\sin\alpha,\cos\alpha,0\right)\,\frac{r_{0}}{r}\,B(r_{0},\phi_{0})\,h_{\text{d}}(z)\,g_{\text{S}}(\phi),$
(18)

where instead of using the radial attenuation $g_{\text{d}}(r)$ the size of
the spur is determined from its angular center $\phi_{\text{C}}$ and angular
half-length $L_{\text{C}}$. We attenuate the magnetic field at
$\phi=\phi_{\text{C}}\pm L_{\text{C}}$ via

$g_{\text{S}}(\phi)=1-\operatorname{\mbox{$\sigma$}}\left(\frac{\Delta(\phi,\phi_{\text{C}})-L_{\text{C}}}{w_{\text{C}}}\right).$
(19)

In the following, we use a fixed attenuation width of $w_{\text{C}}=5^{\circ}$
and without loss of generality $r_{0}=8.2$ kpc, such that the parameter
$B_{0,\text{S}}$ is close to the magnetic field strength at the center of the
spur, if it is located in the proximity of the solar radius, which is the main
focus of the spiral spur model. The superficial lack of flux conservation is
understood as a transition from coherent to turbulent field and broadly-
distributed flux conservation as discussed in the previous subsection.

####  5.1.3 Further Considerations Regarding the Disk Field

Further variations to the Fourier spiral and spiral spur models were studied,
but not included in the fiducial models presented in this paper. We
investigated a circular “ring field” at small galacto-centric radii. Such a
field was introduced in JF12 at $3~{}\text{kpc}<r<5~{}\text{kpc}$ with an
estimated field strength of $0.1\pm 0.1$ $\upmu$G. Here we confirmed that a
ring field does not significantly improve our fits and therefore we did not
include it in our fiducial models.

We also investigated the sensitivity of the data to the particular choice of
$h_{\text{d}}(z)$ and found no significant changes when replacing the logistic
sigmoid with a Gaussian ($\Delta\chi^{2}=-3$) and a slightly worse description
when using an exponential ($\Delta\chi^{2}=+30$).

Furthermore, we studied a flaring disk field, i.e. an increase of the vertical
extent of the disk with galacto-centric radius as e.g. observed for the atomic
hydrogen of the disk (Kalberla & Kerp, 2009). We also added the Galactic warp
as determined by Levine et al. (2006) to the disk field model. Neither of the
two variations improved our preliminary fits reported in (Unger & Farrar,
2019) and we therefore do not include them in this analysis.

###  5.2 Halo Field

Our knowledge of the global structure of the magnetic halo of Milky Way-like
galaxies relies on a combination of edge-on observations of external galaxies
and high latitude rotation measures and radio emission (e.g. Haverkorn &
Heesen, 2012).

The large-scale antisymmetric features of extragalactic rotation measures, cf.
Eq. 10 and Fig. 2a, follow naturally from a large-scale toroidal888A toroidal
field is a vector field with zero $z$-component,
$\bm{B}=(B_{r},B_{\varphi},0)$ and a poloidal field is a vector field with
zero $\varphi$-component, $\bm{B}=(B_{r},0,B_{z})$ field of opposite sign
above and below the plane, leading to the observed pattern in RM when
superimposed with the local spiral arm of the disk field, schematically

$\underset{\text{RM}_{\text{tot}}}{\mathchoice{\leavevmode\hbox
to25.23pt{\vbox to22.08pt{\pgfpicture\makeatletter\hbox{\hskip
12.6138pt\lower-11.04227pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{12.4138pt}{0.0pt}\pgfsys@curveto{12.4138pt}{5.98808pt}{6.85603pt}{10.84227pt}{0.0pt}{10.84227pt}\pgfsys@curveto{-6.85603pt}{10.84227pt}{-12.4138pt}{5.98808pt}{-12.4138pt}{0.0pt}\pgfsys@curveto{-12.4138pt}{-5.98808pt}{-6.85603pt}{-10.84227pt}{0.0pt}{-10.84227pt}\pgfsys@curveto{6.85603pt}{-10.84227pt}{12.4138pt}{-5.98808pt}{12.4138pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-7.77788pt}{-3.33333pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\\!\\!\frac{-+-+}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to15.8pt{\vbox to16.43pt{\pgfpicture\makeatletter\hbox{\hskip
7.89973pt\lower-8.21384pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{7.69974pt}{0.0pt}\pgfsys@curveto{7.69974pt}{4.42596pt}{4.2525pt}{8.01384pt}{0.0pt}{8.01384pt}\pgfsys@curveto{-4.2525pt}{8.01384pt}{-7.69974pt}{4.42596pt}{-7.69974pt}{0.0pt}\pgfsys@curveto{-7.69974pt}{-4.42596pt}{-4.2525pt}{-8.01384pt}{0.0pt}{-8.01384pt}\pgfsys@curveto{4.2525pt}{-8.01384pt}{7.69974pt}{-4.42596pt}{7.69974pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.44453pt}{-2.33333pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\\!\\!\frac{-+-+}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to9.51pt{\vbox to12.66pt{\pgfpicture\makeatletter\hbox{\hskip
4.75702pt\lower-6.32823pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{4.55702pt}{0.0pt}\pgfsys@curveto{4.55702pt}{3.38457pt}{2.5168pt}{6.12823pt}{0.0pt}{6.12823pt}\pgfsys@curveto{-2.5168pt}{6.12823pt}{-4.55702pt}{3.38457pt}{-4.55702pt}{0.0pt}\pgfsys@curveto{-4.55702pt}{-3.38457pt}{-2.5168pt}{-6.12823pt}{0.0pt}{-6.12823pt}\pgfsys@curveto{2.5168pt}{-6.12823pt}{4.55702pt}{-3.38457pt}{4.55702pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.2223pt}{-1.66666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\\!\\!\frac{-+-+}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to9.51pt{\vbox to12.66pt{\pgfpicture\makeatletter\hbox{\hskip
4.75702pt\lower-6.32823pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{4.55702pt}{0.0pt}\pgfsys@curveto{4.55702pt}{3.38457pt}{2.5168pt}{6.12823pt}{0.0pt}{6.12823pt}\pgfsys@curveto{-2.5168pt}{6.12823pt}{-4.55702pt}{3.38457pt}{-4.55702pt}{0.0pt}\pgfsys@curveto{-4.55702pt}{-3.38457pt}{-2.5168pt}{-6.12823pt}{0.0pt}{-6.12823pt}\pgfsys@curveto{2.5168pt}{-6.12823pt}{4.55702pt}{-3.38457pt}{4.55702pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.2223pt}{-1.66666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\\!\\!\frac{-+-+}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}=\underset{\text{toroidal
halo}}{\mathchoice{\leavevmode\hbox to25.23pt{\vbox
to22.08pt{\pgfpicture\makeatletter\hbox{\hskip 12.6138pt\lower-11.04227pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{12.4138pt}{0.0pt}\pgfsys@curveto{12.4138pt}{5.98808pt}{6.85603pt}{10.84227pt}{0.0pt}{10.84227pt}\pgfsys@curveto{-6.85603pt}{10.84227pt}{-12.4138pt}{5.98808pt}{-12.4138pt}{0.0pt}\pgfsys@curveto{-12.4138pt}{-5.98808pt}{-6.85603pt}{-10.84227pt}{0.0pt}{-10.84227pt}\pgfsys@curveto{6.85603pt}{-10.84227pt}{12.4138pt}{-5.98808pt}{12.4138pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-7.77788pt}{-3.33333pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\\!\\!\frac{++--}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to15.8pt{\vbox to16.43pt{\pgfpicture\makeatletter\hbox{\hskip
7.89973pt\lower-8.21384pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{7.69974pt}{0.0pt}\pgfsys@curveto{7.69974pt}{4.42596pt}{4.2525pt}{8.01384pt}{0.0pt}{8.01384pt}\pgfsys@curveto{-4.2525pt}{8.01384pt}{-7.69974pt}{4.42596pt}{-7.69974pt}{0.0pt}\pgfsys@curveto{-7.69974pt}{-4.42596pt}{-4.2525pt}{-8.01384pt}{0.0pt}{-8.01384pt}\pgfsys@curveto{4.2525pt}{-8.01384pt}{7.69974pt}{-4.42596pt}{7.69974pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.44453pt}{-2.33333pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\\!\\!\frac{++--}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to9.51pt{\vbox to12.66pt{\pgfpicture\makeatletter\hbox{\hskip
4.75702pt\lower-6.32823pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{4.55702pt}{0.0pt}\pgfsys@curveto{4.55702pt}{3.38457pt}{2.5168pt}{6.12823pt}{0.0pt}{6.12823pt}\pgfsys@curveto{-2.5168pt}{6.12823pt}{-4.55702pt}{3.38457pt}{-4.55702pt}{0.0pt}\pgfsys@curveto{-4.55702pt}{-3.38457pt}{-2.5168pt}{-6.12823pt}{0.0pt}{-6.12823pt}\pgfsys@curveto{2.5168pt}{-6.12823pt}{4.55702pt}{-3.38457pt}{4.55702pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.2223pt}{-1.66666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\\!\\!\frac{++--}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to9.51pt{\vbox to12.66pt{\pgfpicture\makeatletter\hbox{\hskip
4.75702pt\lower-6.32823pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{4.55702pt}{0.0pt}\pgfsys@curveto{4.55702pt}{3.38457pt}{2.5168pt}{6.12823pt}{0.0pt}{6.12823pt}\pgfsys@curveto{-2.5168pt}{6.12823pt}{-4.55702pt}{3.38457pt}{-4.55702pt}{0.0pt}\pgfsys@curveto{-4.55702pt}{-3.38457pt}{-2.5168pt}{-6.12823pt}{0.0pt}{-6.12823pt}\pgfsys@curveto{2.5168pt}{-6.12823pt}{4.55702pt}{-3.38457pt}{4.55702pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.2223pt}{-1.66666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\\!\\!\frac{++--}{--++}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}+\underset{\text{local
arm}}{\mathchoice{\leavevmode\hbox to25.23pt{\vbox
to22.08pt{\pgfpicture\makeatletter\hbox{\hskip 12.6138pt\lower-11.04227pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{12.4138pt}{0.0pt}\pgfsys@curveto{12.4138pt}{5.98808pt}{6.85603pt}{10.84227pt}{0.0pt}{10.84227pt}\pgfsys@curveto{-6.85603pt}{10.84227pt}{-12.4138pt}{5.98808pt}{-12.4138pt}{0.0pt}\pgfsys@curveto{-12.4138pt}{-5.98808pt}{-6.85603pt}{-10.84227pt}{0.0pt}{-10.84227pt}\pgfsys@curveto{6.85603pt}{-10.84227pt}{12.4138pt}{-5.98808pt}{12.4138pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-7.77788pt}{-3.33333pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\displaystyle\\!\\!\frac{-\phantom{+-}+}{-\phantom{-+}+}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to15.8pt{\vbox to16.43pt{\pgfpicture\makeatletter\hbox{\hskip
7.89973pt\lower-8.21384pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{7.69974pt}{0.0pt}\pgfsys@curveto{7.69974pt}{4.42596pt}{4.2525pt}{8.01384pt}{0.0pt}{8.01384pt}\pgfsys@curveto{-4.2525pt}{8.01384pt}{-7.69974pt}{4.42596pt}{-7.69974pt}{0.0pt}\pgfsys@curveto{-7.69974pt}{-4.42596pt}{-4.2525pt}{-8.01384pt}{0.0pt}{-8.01384pt}\pgfsys@curveto{4.2525pt}{-8.01384pt}{7.69974pt}{-4.42596pt}{7.69974pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.44453pt}{-2.33333pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\textstyle\\!\\!\frac{-\phantom{+-}+}{-\phantom{-+}+}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to9.51pt{\vbox to12.66pt{\pgfpicture\makeatletter\hbox{\hskip
4.75702pt\lower-6.32823pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{4.55702pt}{0.0pt}\pgfsys@curveto{4.55702pt}{3.38457pt}{2.5168pt}{6.12823pt}{0.0pt}{6.12823pt}\pgfsys@curveto{-2.5168pt}{6.12823pt}{-4.55702pt}{3.38457pt}{-4.55702pt}{0.0pt}\pgfsys@curveto{-4.55702pt}{-3.38457pt}{-2.5168pt}{-6.12823pt}{0.0pt}{-6.12823pt}\pgfsys@curveto{2.5168pt}{-6.12823pt}{4.55702pt}{-3.38457pt}{4.55702pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.2223pt}{-1.66666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptstyle\\!\\!\frac{-\phantom{+-}+}{-\phantom{-+}+}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}{\leavevmode\hbox
to9.51pt{\vbox to12.66pt{\pgfpicture\makeatletter\hbox{\hskip
4.75702pt\lower-6.32823pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{}{}{}{}{}{}{}{}{}{{}\pgfsys@moveto{4.55702pt}{0.0pt}\pgfsys@curveto{4.55702pt}{3.38457pt}{2.5168pt}{6.12823pt}{0.0pt}{6.12823pt}\pgfsys@curveto{-2.5168pt}{6.12823pt}{-4.55702pt}{3.38457pt}{-4.55702pt}{0.0pt}\pgfsys@curveto{-4.55702pt}{-3.38457pt}{-2.5168pt}{-6.12823pt}{0.0pt}{-6.12823pt}\pgfsys@curveto{2.5168pt}{-6.12823pt}{4.55702pt}{-3.38457pt}{4.55702pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.2223pt}{-1.66666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\scriptscriptstyle\\!\\!\frac{-\phantom{+-}+}{-\phantom{-+}+}\\!\\!$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\,.$ (20)

In addition, the large-scale patterns of the $Q$ and $U$ parameters, in
particular the tilted nature of the corresponding polarization vectors, can be
explained by the presence of an addition poloidal halo component, as suggested
by Jansson & Farrar (2012a) who introduced a coherent “X-field” halo component
inspired by the X-shaped radio polarization halos observed in edge-on spiral
galaxies (e.g. Krause et al., 2020).

Here we will model the magnetic halo of the Galaxy either by a) a
superposition of a separate toroidal and poloidal field, or b) by a unified
halo model, as described below.

Figure 8: X-field models: JF12 (left), power-function (middle) and coasting
X-field (right).

####  5.2.1 Toroidal Halo

For the toroidal field, we adopt the logistic-exponential ansatz
(“${\operatorname{\mbox{$\sigma$}}}_{r}\,{\exp}_{z}$”) of JF12. It is a purely
azimuthal field,

$\bm{B}_{\text{t}}=(0,B_{\phi},0)$ (21)

in cylindrical coordinates, with

$B_{\phi}=\\\
\left(1-h_{\text{d}}(z)\right)\,e^{-\frac{|z|}{z_{\text{t}}}}\,\operatorname{\mbox{$\sigma$}}\left(\frac{r-r_{\text{t}}}{w_{\text{t}}}\right)\begin{cases}B_{\text{N}}&z>0\\\
B_{\text{S}}&\text{otherwise}\\\ \end{cases}$ (22)

where $z_{\text{t}}$ is the exponential scale height of the toroidal field and
the radial cutoff of the toroidal halo is modeled by a logistic sigmoid
function with a transition at $r_{\text{t}}$ and a width $w_{\text{t}}$. The
maximum field strength above and below the Galactic plane are given by
$B_{\text{N}}$ and $B_{\text{S}}$, respectively. $h_{\text{d}}$ is defined in
Eq. 14, i.e. the toroidal halo field is phased in by the complement of the
function that phases out the disk field.

Whereas in JF12 the radial extent of the Northern and Southern halo was
allowed to be different, we use only one value $r_{\text{t}}$, since our
preliminary fits showed only a minor deterioration of the fit quality when
enforcing a symmetric toroidal (Unger & Farrar, 2018). We also checked whether
the data is fit better with independent functions $h_{\text{d}}$ (z) for the
disk and the halo separately, but found similar transition heights for the
disk and halo, even if they were allowed to take on different values.

####  5.2.2 Poloidal Halo

Jansson & Farrar (2012a) introduced a coherent “X-field” model that is purely
vertical at $r=0$ and becomes increasingly tilted with increasing radius until
reaching a constant asymptotic angle with respect to the Galactic plane,
$\theta_{X}\sim 50^{\circ}$, at a radius of $r_{X}\sim 5$ kpc. With this model
a good description of the polarized synchrotron data could be achieved.
However, this X-field parametrization has a three types of discontinuities:
one in the inner Galaxy, one at the Galactic plane at $z=0$ and one at
$r=r_{X}$, as can be seen in the left panel of Fig. 8. These discontinuities
are avoided by the improved X-field models described in the following.

##### A) Power-function X-field

Ferrière & Terral (2014) employed a useful method to construct poloidal field
models using Euler potentials. Given the equation of field lines in the form
of $r=f(a,z)$ starting midplane at a radius $a=f^{-1}(r,z=0)$, the divergence-
free cylindrical components of the poloidal field ($B_{\phi}=0$) are given by

$B_{r}=-\frac{a}{r}\,B_{0}(a)\,\frac{\partial a}{\partial z}$ (23)

and

$B_{z}=\frac{a}{r}\,B_{0}(a)\,\frac{\partial a}{\partial r}$ (24)

where any function $B_{0}(a)$ of the radial dependence of the $z$-component of
the field at $z=0$ will preserve the solenoidality of the field.

We extend their parabolic “model C” for the $z$-evolution of the field lines
to arbitrary powers $p$,

$r=a\,(1+|z/z_{\text{p}}|^{p}),$ (25)

corresponding to a midplane radius of

$a(r,z)=\frac{r}{1+|z/z_{\text{p}}|^{p}}.$ (26)

Using Eqs. (23) and (23) the radial and vertical field components are

$B_{r}=p\,\frac{z\,a^{3}}{r^{2}\,z_{\text{p}}^{2}}\,|z/z_{\text{p}}|^{p-2}\,B_{0}(a)$
(27)

and

$B_{z}=\frac{a^{2}}{r^{2}}\,B_{0}(a)$ (28)

for $p\geq 1$. In Sec. 7 we will study several possibilities for the radial
midplane dependence of the field strength,

$B_{0}(a)=B_{\text{p}}\,f_{\text{X}}(a),$ (29)

with normalization constant $B_{\text{p}}$ and one of the following radial
functions:

  * •

exp: $f_{\text{X}}(a)=e^{-\frac{a}{r_{\text{p}}}}$, (30)

  * •

gauss: $f_{\text{X}}(a)={e^{-\frac{1}{2}(\frac{a}{r_{\text{p}}})^{2}}}$, (31)

  * •

sech2: $f_{\text{X}}(a)=\operatorname{sech}(a/r_{\text{p}})^{2}$, (32)

  * •

logistic:
$f_{\text{X}}(a)=\operatorname{\mbox{$\sigma$}}((a-r_{\text{p}})/w_{\text{p}})$.
(33)

where $\operatorname{sech}$ denotes the hyperbolic secant,
$\operatorname{sech}(x)=1/\cosh(x)=\nicefrac{{2}}{{(e^{x}+e^{-x})}}$ and
$\sigma$ is the logistic sigmoid function defined in Eq. 16.

An example of a power-function X-field is shown in the middle panel of Fig. 8
(for $p=2.2$ and an exponential radial dependence with $r_{\text{p}}=2.8$
kpc). As can be seen, the discontinuities present in the original JF12 X-field
are avoided, but it differs qualitatively from the JF12 X-field in that the
field lines become more and more parallel to the Galactic plane as the radius
increases rather than reaching an asymptotic angle.

##### B) Coasting X-field

A “coasting X-field” with parallel field lines beyond a certain reference
radius $a_{\text{c}}$ can be achieved by choosing the field-line equation

$r=\left(a^{n}+c\,\frac{a^{n}}{a^{n}+a_{\text{c}}^{n}}\,|z|^{p}\right)^{\frac{1}{n}}.$
(34)

The positive solution for $a^{n}$ is999Note that $\sqrt{x^{2}+k}-x$ is of low
numerical accuracy if $k\ll x$. Instead we use
$(\sqrt{x^{2}+k}-x)\frac{\sqrt{x^{2}+k}+x}{\sqrt{x^{2}+k}+x}=\frac{k}{\sqrt{x^{2}+k}+x}$
in our numerical implementation.

$\displaystyle a^{n}=\frac{1}{2}\Big{(}$
$\displaystyle\sqrt{\Delta^{2}+4a_{\text{c}}^{n}r^{n}}-\Delta\Big{)},$ (35)

where $\Delta=a_{\text{c}}^{n}+c|z|^{p}-r^{n}$. Inserting into Eqs. (23) and
(23), the radial and vertical field components are obtained as

$B_{r}=B_{0}(a)\,\frac{c\,a^{2}\,p\,z\,|z|^{p-2}}{n\,r\,\sqrt{\Delta^{2}+4a_{\text{c}}^{n}r^{n}}}$
(36)

and

$B_{z}=B_{0}(a)\,\left(\frac{r}{a}\right)^{n-2}\frac{a^{n}+a_{\text{c}}^{n}}{\sqrt{\Delta^{2}+4a_{\text{c}}^{n}r^{n}}}.$
(37)

In this paper we study X-fields for the special case $p=n$ for which Eq. 34
becomes

$r=a\left(1+\frac{1}{1+(a/a_{\text{c}})^{p}}\,\left(\frac{|z|}{z_{\text{p}}}\right)^{p}\right)^{\frac{1}{p}},$
(38)

and Eqs. (36) and (37) simplify accordingly. The new parameter $z_{\text{p}}$
is given by $z_{\text{p}}=\frac{a_{\text{c}}}{c^{1/p}}$. For large values of
the coasting radius $a_{\text{c}}$ this three-parameter function simplifies
further to the two-parameter equation

$r=a\left(1+\left(\frac{|z|}{z_{\text{p}}}\right)^{p}\right)^{\frac{1}{p}}.$
(39)

An example of a coasting X-field with $p=n=2.2$ and $a_{\text{c}}=7$ kpc is
shown in the right panel of Fig. 8.

(a) $t=0$ Myr

(b) $t=25$ Myr

(c) $t=50$ Myr

Figure 9: Illustration of the unified halo model at different times $t=0$, 15
and 25 Myr.

####  5.2.3 Unified Halo Model

It has been speculated for a long time that a toroidal halo field with
different directions in the Northern and Southern hemisphere could be the
result of a differential rotation of a poloidal halo field (e.g. Andreasyan &
Makarov, 1988; Men & Han, 2003; Farrar, 2014). Consider a poloidal field that
is dragged along with the rotation of the Galaxy and evolves via the induction
equation,

$\partial_{t}\bm{B}=\bm{\nabla}\times(\bm{v}\times\bm{B})-\underbrace{\bm{\nabla}\times\eta(\bm{\nabla}\times\bm{B})}_{\rightarrow
0\;{\rm for}\;\sigma\rightarrow\infty},$ (40)

where the “frozen-in condition” applies for perfect conductivity,
$\sigma\rightarrow\infty$, for which the magnetic diffusivity vanishes,
$\eta\varpropto 1/\sigma\rightarrow 0$ (e.g. Parker, 1979). Under these
conditions and for a purely azimuthal rotation velocity $v$ the induction
equation simplifies to

$\partial_{t}\bm{B}=\bm{\nabla}\times(\bm{v}\times\bm{B})=\begin{pmatrix}-\frac{v}{r}\,\partial_{\phi}B_{r}\\\
\partial_{z}(vB_{z})+\partial_{r}(vB_{r})\\\
-\frac{v}{r}\,\partial_{\phi}B_{z}\end{pmatrix}.$ (41)

Thus, for a magnetic field that is poloidal and azimuthally symmetric at $t=0$
($\partial_{\phi}B_{r}=\partial_{\phi}B_{z}=0$), the poloidal components are
constant and only $B_{\phi}$ evolves with time,

$B_{\phi}(t)=\left(B_{z}\,\partial_{z}v+rB_{r}\,\partial_{r}\omega\right)\,t,$
(42)

where we introduced the angular velocity
$\bm{\omega}=\omega\,\bm{e}_{\phi}=\frac{v}{r}\,\bm{e}_{\phi}$ and used the
solenoidality of the poloidal field101010Terral & Ferrière (2017) explored a
twisted halo model, but their Eq. (27) is different from Eq. (42) since their
azimuthal field is an ad-hoc “winding function” rather than the physically-
motivated velocity field investigated here..

Eq. (42) can be applied to evolve any type of poloidal field analytically and
it describes the linear increase of $B_{\phi}$ due to the differential motion
of the plasma the magnetic field is embedded in. The large-scale motion of the
plasma in the Galaxy follows the Galactic rotation curve, which we take to
approximately factorize as

$v(r,z)=v_{0}\,f(r)\,g(z).$ (43)

Using

$f(r)=1-e^{-r/r_{v}}\equiv f_{r},$ (44)

we adjust $v_{0}$ and $r_{v}$ to match the velocities of high-mass star-
forming regions with parallax measurements from Reid et al. (2014), leading to
$v_{0}=-240$ km/s (negative because the Galaxy rotates clock-wise) and
$r_{v}=1.6$ kpc. For the vertical velocity gradient we use

$g(z)=2/(1+e^{\,2\,|z|/z_{v}})\equiv g_{z}$ (45)

with $z_{v}=10$ kpc to reproduce $\partial_{z}v=(-22\pm 6)$ (km/s)/kpc as
measured within 100 pc of the Galactic midplane by Levine et al. (2008).

With this particular choice of $f_{r}$ and $g_{z}$, the evolution of the
azimuthal field via Eq. 42 is

$B_{\phi}(t)=\left(B_{z}\,\Delta_{z}+B_{r}\,\Delta_{r}\right)\,t,$ (46)

with the two shear terms

$\displaystyle\Delta_{z}=\partial_{z}v=-\operatorname{sgn}(z)\,\frac{v_{0}}{z_{v}}\,f_{r}\,g_{z}^{2}\,e^{\,2\,|z|/z_{v}}$
(47)

and

$\Delta_{r}=r\,\partial_{r}\omega=v_{0}\,\left(\frac{1-f_{r}}{r_{v}}-\frac{f_{r}}{r}\right)g_{z}~{}.$
(48)

An example of a unified halo field is shown in Fig. 9 starting at $t=0$ with a
poloidal coasting X-field with $B_{\text{p}}>0$ and a logistic sigmoid cutoff
at $r=8$ kpc. Due to the radial and vertical shear of the rotation curve, Eq.
43, an azimuthal field is created at $t>0$ that has a different sign in the
Northern and Southern hemisphere to which both terms in Eq. (42) contribute
constructively. Obviously, this process cannot be at work on Galactic time
scales or otherwise it would overproduce the desired azimuthal field strength.
Including a suitable dissipative term in the induction equation would lead to
a steady-state azimuthal field, but for the purpose of this paper we interpret
the best-fit twisting time as an effective parameter of this simplest version
of a unified halo model.

In summary, in this section we introduced two variants of the disk field
(grand spiral or spiral spur), two variants of the poloidal field (power-
function or coasting) and two variants for the toroidal halo (explicit or from
twisting). These magnetic field sub-models as well as their parameters are
listed in Table 1 together with a few other model parameters.

Table 1: List of parameter names of different model components. name | explanation | unit | value
---|---|---|---
disk field
common parameters
| $\alpha$ | pitch angle | deg. | free
| $z_{\text{d}}$ | transition height | kpc | free
| $w_{\text{d}}$ | vertical transition width | kpc | free
a) grand-design spiral, Eqs. (12) and (11)
| $B_{m}$ | magnetic field strength of mode $m$ |  $\upmu$G | free
| $\phi_{m}$ | phase of mode $m$ | deg. | free
| $r_{0}$ | reference radius | kpc | 5
| $r_{\text{1}}$ | inner radius | kpc | 5
| $w_{\text{1}}$ | inner radial transition width | kpc | 0.5
| $r_{\text{2}}$ | outer radius | kpc | 20
| $w_{\text{2}}$ | outer radial transition width | kpc | 0.5
b) spiral spur, Eq. (17)
| $B_{\text{1}}$ | magnetic field strength at $r_{0}$ |  $\upmu$G | free
| $\phi_{\text{1}}$ | azimuth at $r_{0}$ | deg. | free
| $w_{\text{S}}$ | Gaussian width | deg. | free
| $\phi_{\text{C}}$ | central azimuth | deg. | free
| $L_{\text{C}}$ | half angular length | deg. | free
| $r_{0}$ | reference radius | kpc | 8.2
| $w_{\text{C}}$ | transition width | deg. | 5
toroidal halo
a) explicit, Eq. (22)
| $B_{\text{N}}$ | Northern magnetic field strength |  $\upmu$G | free
| $B_{\text{S}}$ | Southern magnetic field strength |  $\upmu$G | free
| $z_{\text{t}}$ | vertical scale height | kpc | free
| $r_{\text{t}}$ | transition radius | kpc | free
| $w_{\text{t}}$ | radial transition width | kpc | free
b) twisted, Eq. (42)
| $t$ | twisting time | Myr | free
| $v_{0}$ | Galactic rotation velocity | km/s | -240
| $r_{v}$ | scale radius of rotation curve | kpc | 1.6
| $z_{v}$ | scale height of rotation curve | kpc | 10
poloidal halo
common parameters
| $B_{\text{p}}$ | magnetic field strength |  $\upmu$G | free
| $p$ | field line exponent | – | free
| $z_{\text{p}}$ | scale height | kpc | free
| $r_{\text{p}}$ | radial scale or transition radius | kpc | free
| $w_{\text{p}}$ | transition width | kpc | free
a) power-function, Eqs. (23) and (24)
b) coasting, Eqs. (36) and (37)
| $a_{\text{c}}$ | coasting radius | kpc | free
other model parameters
| $\kappa$ |  $n_{e}$-$B$ correlation coefficient, Eq. 3 | – | free
| $\xi$ | magnetic striation factor Eq. 9 | – | free

##  6 Model Optimization

The parameters $\bm{p}$ of a GMF model are optimized by minimizing the sum of
the variance-weighted squared difference between the predicted (m) and
measured (d) observables,

$\displaystyle\chi^{2}=$ $\displaystyle\sum_{i=\scalebox{0.6}{$\text{RM}$,$Q$,
$U$}}\;\sum_{j=1}^{N_{i}}\frac{(d_{ij}-m_{ij}(\bm{p}))^{2}}{\sigma_{ij}^{2}},$
(49)

where the second sum runs over the $N_{j}$ lines of sight available for the RM
and synchrotron data. Each line of sight data $d_{ij}$ is obtained by
averaging available measurements over a finite solid angle around the line-of-
sight direction. These angular pixels are defined by the HEALPix resolution of
the data. The predicted values $m_{ij}$ of a GMF model,
$\bm{B}(\bm{x};\bm{p})$, are given by the numerical evaluation of the line-of-
sight integrals for RM, $Q$ and $U$. Even for a perfect GMF model, the
differences $d_{ij}-m_{ij}$ are expected to be distributed with a variance of
$\sigma_{ij}^{2}$, since the data is subject to experimental uncertainties and
the model prediction do not include “Galactic variance” originating from
random magnetic fields and from fluctuations of the densities of thermal
electrons and cosmic-ray electrons. In general, Galactic variance will
introduce correlations between adjacent pixels if the size of the
perturbations is larger than the angular size of one pixel. In that case, the
full covariance matrix needs to be included in the calculation of the
$\chi^{2}$. It is, however, non-trivial to calculate the coefficients of the
matrix, because one needs to know the spatial distribution of the turbulence
and coherence length for the magnetic field and the thermal and cosmic-ray
densities. Different approximations have been applied in previous analyses.
Jaffe et al. (2010) and Adam et al. (2016) used a model of the three-
dimensional random field strength with a constant coherence length and
unperturbed electron densities to calculate the diagonal elements of the
covariance matrix. In a more data-driven approach, one can measure the sub-
pixel variance $\hat{\sigma}^{2}$ of the measurements within one pixel and use
it as a weight in the fit. Naively, the average over $N$ observations within
one pixel would then reduce the fluctuations of the mean by a factor of $1/N$,
but due to the aforementioned coherent effects, the effective variance will
typically be larger. Based on a toy model of coherent cells along the line of
sight, Pshirkov et al. (2011) suggested to use $\hat{\sigma}^{2}/\delta^{2}$
as the pixel variance, with $\delta\approx 3$, but the exact value of $\delta$
depends on the integration distance in units of coherence length over which
the GMF contributes to the observations (see also (Terral & Ferrière, 2017)).
Keeping in mind that the overall normalization of the $\chi^{2}$ could be off
by a factor of $\delta$, we follow the procedure of JF12 and simply use
$\hat{\sigma}^{2}$ itself to weight the data points.

The best set of $n$ model parameters $\bm{p}$ are found performing a multi-
dimensional optimization of Eq. 49 with the Minuit program (James & Roos,
1975) using its Migrad method that implements a variable metric gradient
descent (Fletcher & Powell, 1963; Fletcher, 1970). Deterministic gradient
methods like Migrad descend quickly to a local minimum which may not always
coincide with the global minimum. To find this global optimum we use the
heuristic multi-start method, i.e. we perform a number of gradient-descent
minimization runs starting at different positions $\bm{p}_{\text{start}}$
distributed uniformly in the $n$-dimensional hypercube of the parameter space.
The run with the smallest local $\chi^{2}$ is then considered to be the global
minimum. For the model fits performed in this paper we run typically
$\mathcal{O}(100)$ minimizations. The fact that usually many runs from very
different starting positions converge to the same minimum local $\chi^{2}$,
increases our confidence that the $n$-dimensional likelihood contour of the
optimization problem at hand is well-behaved and that we have indeed
identified the global minimum of each model variation investigated. The main
advantage of this simple multi-start method is that such a repeated gradient-
descent can be trivially parallelized on a computing cluster and thus it
provides an efficient method to search for the global minimum of a model,
since each gradient search is very fast111111The optimization of the
20-parameter base model discussed in Sec. 7 needs on average 4000 $\chi^{2}$
evaluations to converge. Each evaluation takes 4 seconds on an
Intel${}^{\text{\textregistered}}$ Xeon${}^{\text{\textregistered}}$ E5
processor at 2.4 GHz. For each evaluation, $\bm{B}$, $n_{e}$, $n_{\text{cre}}$
and $j_{\nu}$ are calculated at $\mathcal{O}(10^{6})$ positions throughout the
Galaxy to calculate the adaptive line-of-sight integrals of RM, $Q$ and $U$. .

For each minimization run, a good first approximation of the covariance matrix
of the best-fit parameters can be obtained by the parabolic estimates derived
from the Hesse-matrix of the second derivatives of the $\chi^{2}$ with respect
to the parameters at the minimum. These estimates come without any additional
computational cost, since the Hesse matrix is anyway calculated by Minuit
during the minimization. For the final fiducial models presented in Sec. 7,
more precise confidence intervals of the parameters at the minimum are derived
with the profile likelihood method, i.e. by finding for each parameter the two
values $p_{\text{up/low}}$ at which
$\smash{\chi^{2}(p_{\text{up/low}})=\chi_{\text{min}}^{2}+m}$, while
marginalizing over the other $n-1$ parameters. These estimates are obtained
using the Minos algorithm of Minuit and each evaluation of one
$p_{\text{up/low}}$ pair is about as computationally expensive as the overall
minimization itself. Unless stated otherwise, we quote the 68% (“one
$\sigma$”) intervals obtained for $m=1$.

It is worthwhile noting that our approach differs from the one of previous GMF
studies, which used Markov-Chain Monte Carlo (MCMC) to explore the model
parameters, (e.g. Jansson et al., 2009; Jaffe et al., 2010; Jansson & Farrar,
2012a; Terral & Ferrière, 2017). MCMC is the preferred method to sample the
posterior distribution of a model given the data, but it is not an efficient
optimizer (e.g. Hogg & Foreman-Mackey, 2018). Moreover, for most applications,
the approximate covariance matrix and the confidence intervals derived with
the profile likelihood method provide information equivalent to the MCMC
samples; see, e.g. (Ade et al., 2014), for a comparison of the two methods in
the context of cosmological parameter estimation. Most importantly, the
ability to efficiently optimize the parameters of many different models is of
paramount importance for the GMF inference, since, as will be shown in Sec. 7,
the systematic differences resulting from different model assumptions are
typically much larger than the precision estimated for the parameter
uncertainties of a particular model.

##  7 Results

disk

toroidal halo

poloidal halo

full model

masked model

data

pull

Figure 10: Data and base model prediction for RM, $Q$ and $U$. The predictions
from the individual model components are shown in the top rows, followed by
full model prediction. The bottom three rows show the masked model, the data
and the pull, i.e. the difference between model and data in units of the
standard deviation of the data.

###  7.1 Base Model

For the results presented in the following, we have optimized the parameters
of more than 200 combinations of magnetic field models and auxiliary models to
obtain an overview of the range of GMF models attainable in the full set of
possible combinations121212The combination of two disk field models, eight
poloidal field models, two toroidal field models, two thermal electron models,
two synchrotron products and six cosmic-ray electron models, results in 1536
possible model variations.. Many of these model variations are performed with
respect to our fiducial “base” model which consists a three-mode grand-design
spiral disk field, an explicit toroidal halo and a coasting X-field with a
logistic sigmoid radial dependence; there is no correlation between $n_{e}$
and $B$ and the striation is a free parameter. The skymaps of RM, $Q$ and $U$
resulting from this model are calculated with the Ymw16 thermal electron model
and a cosmic-ray electron density derived for $h_{D}=6$ kpc. The model
parameters are adjusted to fit the data yielding an optimum with an acceptable
goodness of fit of $\chi^{2}/\text{ndf}=7923/6500=1.22$, where ndf denotes the
number of degrees of freedom, i.e. the number of data points minus the number
of free parameters, which is $n_{\text{par}}=20$ in this case. The
contribution from RM pixels to the $\chi^{2}$ is 4354 ($n_{\text{RM}}=2838$)
and from $Q$ and $U$ pixels it is 3569 ($n_{Q}+n_{U}=3682$), i.e. the model
describes the polarized synchrotron data slightly better than the RMs.

The contributions of each model component (i.e. disk field and toroidal and
poloidal halo) to RM, $Q$ and $U$ is shown separately in the three top rows of
Fig. 10 and the predicted sky maps for the sum of all components is displayed
in the fourth row. Here it is interesting to see how the interplay of
different model components create the large-scale features observed in data.
As pointed out above, the sum of the disk and toroidal field produces the
large-scale structure of RMs, cf. Eq. 20. The polarized intensity on the other
hand, is mostly generated by the toroidal and poloidal halo. And, whereas the
overall RM can be obtained by summing the RM of each component
($\text{RM}\varpropto\bm{B}_{\parallel}=(\sum\bm{B}_{i})_{\parallel}=\sum\bm{B}_{i\parallel}$),
the sky pattern of the Stokes parameters of the full model is not the sum of
its components since e.g.
$Q,U\varpropto\bm{B}_{\perp}^{2}=\left((\sum\bm{B}_{i})_{\perp}\right)^{2}\neq\left(\sum\bm{B}_{i\perp}\right)^{2}$.

The masked model, data and “pull” are shown in the three bottom rows. The pull
is the difference between the data and model in units of standard deviation of
the data and the sum of all squared pulls yields the $\chi^{2}$, Eq. 49. For a
perfect match of data and model, the pull should fluctuate randomly around
zero with a standard-normal distribution, but here several large regions exist
in which the pull has consistently negative or positive values beyond $\pm
1\sigma$. Some of these are close to masked regions, e.g. at the edges of the
masks for the North Polar spur and loop III for $Q$ and $U$ and close to the
mask for the Gum nebula for RM, implying a certain amount of leakage of these
features beyond the mask. Others, like the large region of negative pull for
RM below the Galactic plane at longitudes of $60^{\circ}<\ell<150^{\circ}$
could indicate deficiencies in the modeling of the global structure of the
GMF, e.g. an azimuthal variation of the scale height of the toroidal field. It
is, however, plausible that most regions with large-scale deviations can be
attributed to local perturbations that appear as structures of large angular
scale due to their proximity.

Apart from these potentially local structures, the base model successfully
describes all the large-scale features of the data, in particular the anti-
symmetric structures of the RM sky and the tilted large-scale “lobes” of
negative and positive $Q$ and $U$. The parameters of the base model can be
found in the first column of Table 3.

(a) Half-height, $h_{D}$, of the cosmic-ray diffusion volume.

(b) Correlation coefficient, $\kappa$, between the magnetic field and thermal
electron density.

Figure 11: Goodness of fit, Eq. 49, for two different models of the thermal
electron density, Ne2001 (red) and Ymw16 (blue), and as a function of $h_{D}$
(left) and $\kappa$ (right). Figure 12: Best-fit parameter values of the disk,
toroidal and poloidal GMF components for different values of the half-height
of the cosmic-ray diffusion volume, $h_{D}$, and for two different models of
the thermal electron density, Ne2001 (red) and Ymw16 (blue). See Table 1 for a
short explanation of each parameter. Some parameter values have been
multiplied by a scale factor, as indicated in the axis labels, to fit in one
panel with a single $y$-axis. One-sigma uncertainties are shown as error bars
and arrows indicate the 84% CL lower limits on the parameters Figure 13: Best-
fit parameter values of the disk, toroidal and poloidal GMF components for
different values of the correlation coefficient between the magnetic field and
thermal electron density, $\kappa$, and for two different models of the
thermal electron density, Ne2001 (red) and Ymw16 (blue). See Tab. 1 for a
short explanation of each parameter. Some parameter values have been
multiplied by a scale factor, as indicated in the axis labels, to fit in one
panel with a single $y$-axis. No parameter uncertainties are shown here, see
Fig. 12.

###  7.2 Thermal and Cosmic-Ray Electrons

To study the dependence of the inferred GMF on the auxiliary models of the
thermal and cosmic-ray electron density, we repeat our fits for the model
variations discussed in Sec. 4. The best-fit $\chi^{2}$ values for different
half-heights of the cosmic-ray diffusion volume, $h_{D}$, and for the Ne2001
and Ymw16 model are shown in Fig. 11a. Here the functional forms for the GMF
model are identical to the one used for the base model. As can be seen, the
fit quality deteriorates rapidly for small values of the height of the
diffusion volume and it reaches a near-constant value at $h_{D}\gtrsim 4$ kpc.
For small values of $h_{D}$, the decreasing vertical extent of the cosmic-ray
electron halo cannot be compensated for by a larger magnetic halo131313For
illustration, assume an exponential height dependence of $n_{\text{cre}}$ and
$B$ with scale heights $z_{\text{cre}}$ and $z_{\text{B}}$. Then, using Eq. 4,
the height-dependence of the synchrotron emissivity is $j_{\nu}\varpropto
e^{-|z|/z_{\text{syn}}}$ with
$z_{\text{syn}}=z_{\text{cre}}/(2+z_{\text{cre}}/z_{B})$. It follows that for
a given $z_{\text{syn}}$ of the data, $z_{\text{cre}}$ has to be at least
$\geq 2\,z_{\text{syn}}$ ($z_{B}\rightarrow\infty$).. Interpolating between
the fit results at discrete values of $h_{D}$, we estimate a 5-$\sigma$ lower
limit on the size of the diffusion volume of

$h_{D}\geq 2.9~{}\text{and}~{}3.5~{}\text{kpc},$ (50)

where the two values correspond to the fits with Ne2001 and Ymw16
respectively. This lower limit is, however, only indicative since in this
analysis, there is no feedback between the derived magnetic field and the
propagation of cosmic-ray electrons, see Sec. 4.2.

The interplay of diffusion height and derived magnetic field parameters can be
seen in Fig. 12, where the GMF model parameters of the different fits for
$h_{D}\geq 4$ kpc are displayed. Especially the fitted scale height of the
toroidal field, $z_{\text{t}}$, depends strongly on $h_{D}$. For $h_{D}=4$ kpc
only a lower limit on $z_{\text{t}}$ can be estimated. In the following, we
will therefore only investigate GMF models derived for $n_{\text{cre}}$
densities within $6\leq h_{D}/\text{kpc}\leq 10$ kpc, for which the best-fit
values of $z_{\text{t}}$ are finite and compatible with current estimates of
the size of the diffusive halo estimated from the analysis of unstable
secondary cosmic-ray nuclei.

Concerning the models of the density of thermal electrons in the Galaxy, it
can be seen in Fig. 11a that the fits with the old Ne2001 model result in a
consistently better fit quality than the ones performed with the recent Ymw16
model. Even though the difference is statistically significant
($\Delta\chi^{2}=-167$ at $h_{D}=6$) we consider both models in our analysis,
since the latter $n_{e}$ model gives a better description of Galactic pulsars
DMs than the former. Switching between these two models has a larger
systematic impact on the values of most GMF parameters than a change in
$h_{D}$, cf. Fig. 12.

(a) Grand-design spiral

(b) Spiral spur.

Figure 14: Magnetic field of the disk component for the grand-design spiral
(left) and spiral spur (right) GMF model. The horizontal magnetic field
strength in the $x$-$y$ plane at $z=0$ is displayed with colors ranging from
blue (counter-clockwise field) to red (clock-wise field). Superimposed points
with distance error bars are the locations of tracers of the spiral structure
of the matter density of the Milky Way (high-mass star-forming regions with
parallax distances) and the curved colored lines are the inferred location of
spiral arm segments, both from Reid et al. (2019).

###  7.3 Striation or Correlation?

It is well known that the magnetic field strength derived from the observed
rotation measures is smaller than the one derived from the observed Galactic
polarized synchrotron emission. One way to reconcile the two observables is to
postulate the existence of an anisotropic or striated random field that
fluctuates along the coherent field lines, which leads to an increase in
polarized intensity without changing the rotation measures. For our base
model, the best-fit value of the striation factor $\xi$, cf. Eq. 9, is

$\xi_{\text{{{base}}}}=0.35\pm 0.03,$ (51)

implying that the energy density in the coherent and striated field component
are almost equal (energy density $u_{B}\varpropto B^{2}$,
$B^{\prime}=(1+\xi)B$, $u_{B^{\prime}}/u_{B}=(1+\xi)^{2}=1.8\pm 0.1$).

Another reason for the mismatch between magnetic fields inferred from RM and
PI observations could be an anti-correlation of the thermal electrons and the
magnetic field, leading to smaller Faraday rotation than in the uncorrelated
case. Here we perform, to our knowledge for the first time, an analysis of RM
and PI allowing for a modified RM due to an $n_{e}$-$B$ correlation with
coefficient $\kappa$, as derived by Beck et al. (2003) (cf. Sec. 2). Here we
use the Planck-tune of the JF12 random field $b$ (Jansson & Farrar, 2012b;
Adam et al., 2016) to evaluate Eq. 3 for a fixed value of $b$, but a coherent
field $B$ that is allowed to float freely during the fit.

The fits are performed for different values of $\kappa$ and the resulting
dependence of the fit quality on the correlation coefficient is displayed in
Fig. 11b. As can be seen, a large anti-correlation of $n_{e}$ and $B$ is
disfavored and the five-sigma lower limit is found to be

$\kappa\geq-0.52~{}\text{and}~{}-0.49,$ (52)

where the two values correspond to the fits using the Ne2001 and Ymw16 models
of $n_{e}$, respectively. The optimal fit to the data is at
$\kappa\approx-0.4$, but the $\chi^{2}$ minimum is very broad and all fits
above this value fit the data similarly well. The reason for this degeneracy
can be understood by examining the change of fit parameters with $\kappa$, as
shown in Fig. 13. Given the two degrees of freedom of rescaling the magnetic
field strengths of each model component and the striation factor, it is always
possible to match the RM and PI data as long as $\xi>0$. In the following, we
will not further investigate the fits with $\kappa>0$, since our fiducial fits
with $\kappa=0$ exhibit already a large degree of striation close to
equipartition with the coherent field (see above). Instead, we will
concentrate on the best fit at $\kappa=-0.4$ where the striation parameter is
zero as an alternative to the fiducial model. At this value of $\kappa$, the
magnetic field scale is set by the PI data, whereas at $\kappa=0$ the scale is
set by the RM data. In this extreme scenario, where no striated random fields
contribute to the polarized intensity and the Faraday rotation is diminished
by the anti-correlation of $B$ and $n_{e}$, the fitted magnetic field
strengths ($B_{m}$, $B_{\text{N/S}}$ and $B_{\text{p}}$) are about a factor
1.4 larger than in the fiducial case.

Figure 15: Comparison of data (black points with error bars) and models (red
open circles: base model with a grand-design spiral disk field, blue open
squares: spiral spur model) in five latitude slices around the Galactic plane.
Every point is the RM value detected/predicted in one HEALPix pixel
($N_{\text{side}}=16$). Gray bands are drawn at the longitudes of tangent
points of the spiral arms of the Galaxy. The width indicates the range of
tangent positions in different observables catalogued by Vallée (2022).

###  7.4 Synchrotron Data Products

The choice of data product for the Stokes $Q$ and $U$ parameters affects the
GMF fit due to the differences in the large-scale structure of derived sky
maps, see discussion in Sec. 3.2 and the comparisons shown in Fig. 22 in the
Appendix. Moreover, and perhaps more importantly, since these differences are
relatively small, different products differ by the variance of $Q$ and $U$
within a pixel over which we average these values. Since we use these
variances as a weight of the data when optimizing Eq. 49, they also affect the
outcome of the fit.

For our simple arithmetic average of the $Q$ and $U$ values of Planck and
Wmap, the variance is calculated for all data points before averaging,
therefore we expect it to be the arithmetic average of the Planck and Wmap
variances plus an additional contribution from the systematic difference
between the two products. On average, the variance derived from the Planck and
CosmoGlobe maps is about 0.8 smaller than the one of our default synchrotron
map, and about the same for the Wmap data. Keeping in mind that Wmap could
measure at lower frequencies where the synchrotron intensity is larger by a
factor of $(22.5~{}\text{GHz}/30~{}\text{GHz})^{-3.1}=2.4$, it is not clear to
us to which extent the smaller variance of Planck (and by extension
CosmoGlobe) is actually due to an improved measurement of the synchrotron
radiation or due to the use of more aggressive foreground smoothing priors in
the analysis.

As a consequence of smaller variances, the fits using Planck and CosmoGlobe
have a $\chi^{2}$ that is worse than the one of the base model by a factor of
1.4 and 1.2, respectively. We have checked that the worse fit is indeed mostly
due to the different variances (and not the different values of $Q$ and $U$),
by performing a fit of the Planck data with our default variances, in which
case the fit quality was only worse by a factor of 1.03.

The GMF models obtained with alternative data products are qualitatively
similar, but the parameters obtained in these fits differ from the one of the
base model. The differences are mostly not driven by the systematic
differences between the sky maps, but rather caused by the larger relative
“pull” on the $\chi^{2}$ of $Q$ and $U$ data with respect RM data. The
sensitivity of the parameter values on the data-driven weights are intrinsic
to the method used in this analysis and needs to be considered when estimating
the GMF model uncertainties, see Sec. 7.8 below.

###  7.5 Disk Field

####  7.5.1 Grand-Design Spiral

Many of the previous attempts to model the global structure of the disk field
enforced a certain number of “magnetic arms”, often following the inferred
large-scale structure of matter density in the Galaxy. Due to our more
flexible decomposition of the spiral field into azimuthal Fourier modes, we
can let the data decide how much spiral structure is needed. We have fit the
RM and PI data using a different number of Fourier modes ranging from 1 to 5.
As expected, the fit quality continuously improves with the number of modes,
since each adds two more free parameters (amplitude and phase) to the fit.
Using Wilks’ theorem (Wilks, 1938), we find that the improvements of the fit
quality of the base model are significant up to $n_{\text{mode}}=3$, with a
$\Delta\chi^{2}=153.2$ going from 2 to 3 modes and $\Delta\chi^{2}=12.7$ from
3 to 4 modes. The same conclusion is reached using the Bayesian information
criterion, $\text{BIC}=\chi^{2}+n_{\text{par}}\,\ln(n_{\text{data}})$
(Schwarz, 1978), which is minimal at $n_{\text{mode}}=3$. We have also tested
that the required number of modes remains three when changing the height of
the diffusion volume, $h_{D}$, used to calculate $n_{\text{cre}}$ and when
using a twisted X-field.

Three azimuthal modes correspond to six “magnetic arms” of alternating
polarity. The grand-design configuration of our neCL fit variant (using the
Ne2001 thermal electron model) has already been shown above in the right panel
of Fig. 7. The best-fit disk field of the base model is displayed in the left
panel of Fig. 14, where we also show the location of tracers of the spiral
structure of the matter density of the Milky Way, given here by the
measurements of trigonometric parallaxes of high-mass star forming regions
from Reid et al. (2019). As can be seen, we find a remarkable alignment of the
fitted magnetic spiral structure and these tracers. Large coherent field
strengths are present in the interarm regions, but the coherent field strength
is close to zero at the location of the spiral matter segments derived by Reid
et al. (2019), shown as lines in Fig. 14. This result is similar to what is
observed in external galaxies, where the strongest ordered fields are detected
in the interarm regions (Beck, 2016).

The fitted pitch angle, $\alpha$, of the disk field is found to be nearly
independent of the assumed cosmic-ray halo size, $n_{e}$-$B$ correlation
coefficient ($\kappa\geq-0.4$) and functional form of the magnetic halo.
However, the fits with different models of the thermal electron density result
in pitch angles that are systematically different with respect to each other
by $\pm 1^{\circ}$, cf. Fig. 12. The average of the value obtained using the
Ne2001 and Ymw16 thermal electron models is

$\alpha=(11.0\pm 0.3\;(\text{stat.})\pm 1.0\;(\text{$n_{e}$}))^{\circ}.$ (53)

This value is in good agreement with the pitch angle of the local
(Orion–Cygnus) spiral arm of $(11.4\pm 1.9)^{\circ}$ (Reid et al., 2019) and
the pitch angles of $\alpha=(9.87\dots 11.43)^{\circ}$ of the grand-design
logarithmic spiral model fitted by Hou & Han (2014) using HII regions as
spiral tracers. Thus, the pitch angles of the spiral structure of the magnetic
field and the matter density in the Milky Way are about equal, similar to what
is observed for external spiral galaxies (Van Eck et al., 2015).

####  7.5.2 Local Spur

What drives the fit of the grand-design spiral disk field? Is it the RM values
at longitudes where the line of sight is tangential to a spiral arm? These
magnetic tangents are clearly visible in the predicted RMs in the top left
plot of Fig. 10, but they are not so obvious in the RM of the data, possibly
due to the large variance present in the data at low latitudes. Alternatively,
the fit of the disk field might be mostly determined by the large-scale
“butterfly pattern” of RMs, cf. Eq. 20, and therefore mostly from the disk
field in the vicinity of the Sun.

To test this possibility, we fit the data without a grand-design magnetic
pattern, but include only one local magnetic spur, as introduced in Sec.
5.1.2. The fitted pitch angle of the local spur is $(12.1\pm 0.6)^{\circ}$ and
does not depend on the thermal electron model used. The best-fit local spur is
displayed in the right panel of Fig. 14. It is located at the edge of the
local Orion–Cygnus spur and has a magnetic field strength of
$4.30~{}\upmu\textrm{G}$ at the reference radius of 8.2 kpc. The quality of
this fit is found to be slightly worse than the one of the base model (7991
instead of 7923), but close enough to conclude that both models are
approximately equivalent.

For a closer look at the differences between the spur and base model, we show
the RM values in bands of latitude for data and the two models in Fig. 15. As
can be seen, both models give the same overall good description of the
longitude-dependence of RM of the data. The base model (shown as red open
circles) exhibits distinct RM features at low latitudes, shown in the middle
panel for $|b|<5^{\circ}$, but none of these are visible in the data.

The near-equivalence of these two radically different disk field models is a
consequence of the fact that the extragalactic RMs used in the fit constrain
only the integrated Faraday rotation to the edge of the Galaxy. The model
degeneracy could, in principle, be broken by including RMs from Galactic
pulsars into the fit to provide cumulative RMs at different distances.
However, pulsar RMs are currently of limited use, since the distances to most
pulsars are not known with sufficient precision. Some previous GMF studies
used the “DM-distance” (i.e. the distance $d$ at which
$\text{DM}^{\prime}=\int_{0}^{d}n_{e}(\bm{x}(r))\,{\text{d}}r$ equals the
observed DM) to calculate the model RMs, but this introduces an additional
dependence on thermal electron density models which is difficult to account
for in the GMF model optimization. Notwithstanding these caveats, it is
interesting to note that Han et al. (2018) derived a six-arm grand-design
spiral model from the RMs of pulsars and extragalactic sources, which is
qualitatively very similar to the disk field of our base model. More studies
are however needed to unequivocally prove the existence of a grand-design
magnetic spiral in our Galaxy. Until then, a conservative approach is to
consider the two models introduced in this section as extreme possibilities
for the coherent magnetic field in the disk of the Milky Way.

Figure 16: Best-fit radial dependence of the $z$-component of the poloidal
magnetic field at $z=0$, $B_{0}(a)$, cf. Eqs. (23) and (24). Different radial
functions are shown in different colors, cf. Eqs. (• ‣ 5.2.2) - (• ‣ 5.2.2).
Variations in the thermal electron density and toroidal field type used during
the fit are shown with different line styles, as indicated in the figure
legend.

###  7.6 Poloidal Field

We begin by reassessing the need for a poloidal component to describe the
data, performing a fit without including such a component. The resulting fit
quality deteriorates tremendously, by $\Delta\chi^{2}>+2100$ with respect to
the base model, irrespective of the thermal electron model used in the fit. We
also tried to fit the data with a N-S reflection-symmetric X-field
(characteristic of a S1 dynamo), instead of a dipolar one (A0 dynamo), without
success ($\Delta\chi^{2}>+1700$).

Inspecting the pixel-by-pixel contributions to the $\chi^{2}$ we find that the
poloidal fit is mainly driven by the $Q$ and $U$ data at longitudes
$|\ell|<60^{\circ}$. Unless a very peculiar141414We performed a few initial
tests with a magnetic field compressed at the edge of a local foreground
bubble (e.g. Ferriere et al., 1991) using a simple spherical model with
magnetic flux conservation, but could not find a configuration that described
the data. Further, more systematic studies are needed to definitely exclude
this possibility. foreground is responsible for these large-scale features in
$Q$ and $U$, we confirm with high significance the conclusions of Jansson &
Farrar (2012a) that the polarized synchrotron intensity cannot be described
without the presence of a dipolar X-field in the Galaxy.

We find that the data is described equally well by a power-function and a
coasting X-field ($\chi^{2}=7926$ and 7923, respectively), since the major
differences between the field lines of the two models are at large radii,
where the field is small, see Fig. 8. There is, however, a clear preference
regarding the choice of the radial dependence of the mid-plane vertical field,
Eqs. (• ‣ 5.2.2)-(• ‣ 5.2.2). The best fit to the data is obtained for the
logistic sigmoid function, i.e. for a vertical field strength that is constant
with galactocentric radius and then vanishes to zero at a certain radius
$r_{\text{p}}$ with a transition width $w_{\text{p}}$. The tested alternatives
of a Gaussian, hyperbolic secant or exponential radial dependence result in
fits with are worse by $\Delta\chi^{2}=+370$, $+452$ and $+508$. This clear
preference for a logistic radial cutoff remains for all model variations
tested (power-function or coasting X-field, twisted X-field and using
different thermal electron models).

The improvement of $\chi^{2}$ when using a logistic instead of an exponential
function originates mostly from the RM ($\Delta\chi^{2}_{\text{RM}}=315$) and
$Q$ data ($\Delta\chi^{2}_{Q}=155$ and only to a lesser extent from the $U$
data ($\Delta\chi^{2}_{U}=36$). Much of the improvement in
$\chi^{2}_{\text{RM}}$ is concentrated in one swathe of the RM sky located at
$\ell\approx(-30\pm 15)^{\circ}$ and $90^{\circ}<b<150^{\circ}$ and most of
the decrease in $\chi^{2}_{Q}$ is due to a better description of the data at
$15^{\circ}<\ell<30^{\circ}$ and $b>0^{\circ}$. Since we cannot be certain
that these features at intermediate angular scales are global, we
conservatively keep the exponential radial function as an extreme variation in
our model ensemble.

The best-fit mid-plane poloidal magnetic field strengths, derived for
different radial functions, are shown in Fig. 16. The largest differences of
$B_{0}(a)$ are at small radii, $a\lesssim 4$ kpc, where the fit of the
poloidal field strength is unconstrained due to the large amount of $Q$ and
$U$ pixels masked at $|\ell|<30^{\circ}$. The logistic model used in our base
model has a mid-plane vertical field strength of $B_{\text{p}}=(0.98\pm 0.03)$
$\upmu$G and cuts off inside the solar circle at $r_{\text{p}}=(7.29\pm
0.06)~{}\text{kpc}$ with a transition width of $w_{\text{p}}=(0.112\pm
0.029)~{}\text{kpc}$. In the case of the twisted X-field model, the transition
is fitted to be at a somewhat large radius and the transition is broader
($r_{\text{p}}\sim 7.9$ kpc, $w_{\text{p}}\sim 0.3$ kpc). The expX model leads
to a large vertical magnetic field at the Galactic Center,
$B_{\text{p}}=(5.8\pm 0.4)$ $\upmu$G that falls off exponentially with scale
length of $r_{\text{p}}=(2.5\pm 0.1)$ kpc. For comparison, the poloidal field
of the JF12 model is of similar strength with $B_{\text{p}}=(4.6\pm 0.3)$
$\upmu$G and $r_{\text{p}}=(2.9\pm 0.1)$ kpc.

Thus, the expX and base models bracket the possible range of the radial
dependence of the strength of the poloidal field. The corresponding
uncertainty of the GMF, quantified here for the first time, must be taken into
account when interpreting the arrival directions of ultrahigh-energy cosmic
rays and when discussing the radial dependence of cosmic-ray energy spectra
derived from gamma-ray observations (e.g. Gabici et al., 2019). A better
understanding of the polarized foregrounds is needed to be able to further
constrain the GMF fits by including PI data towards the Galactic center.

###  7.7 Toroidal Field

Fitting for an “explicit” toroidal halo field, we obtain a radial extent of
$r_{\text{t}}=(10.2\pm 0.2)$ kpc and a vertical scale height of
$z_{\text{t}}=(4.0\pm 0.7)$ kpc within our base model. The estimate of the
radial extent of the halo is very stable with respect to any of the model
variations we have studied. The vertical scale is less certain, ranging from
2.9 to 6.1 kpc, where the lowest value is obtained assuming a large cosmic-ray
diffusion volume with $h_{D}=10$ kpc.

The magnitude of azimuthal field strength in the Northern and Southern
hemispheres was found to vary between about 2 and 5 $\upmu$G, depending on the
model variation under study. For the base model, the best-fit values are
$B_{\text{N}}=3.3\pm 0.3$ $\upmu$G and $B_{\text{S}}=-3.1\pm 0.3$ $\upmu$G.
Remarkably, in all the studied variations, the magnitude of the Northern and
Southern magnetic field strength is found to be compatible within the
estimated uncertainties,

$B_{\text{N}}=-B_{\text{S}},$ (54)

which strongly suggests a common origin of the toroidal field in the two
hemispheres.

Our twisted X-field model naturally explains this equality, since the magnetic
field in both hemispheres originates from the shearing caused by the same
velocity field. However, the fit quality of twistX is considerably worse than
the one of the explicit fit, $\Delta\chi^{2}=+462$. On a technical level, this
can be explained by the freedom of two different radial cutoffs values,
$r_{\text{p}}$ and $r_{\text{t}}$, for the toroidal and poloidal components
that are available in the fit with an explicit toroidal halo. A future
implementation of a less simplistic unified halo model would naturally
introduce additional radial and vertical dependencies, e.g. from the spatial
dependence of the dissipation rate, which affects the toroidal and poloidal
components differently.

A practical virtue of the unified halo model is that its built-in connection
between the radial extent of the poloidal and toroidal components should make
it less prone to over-fitting local structures in the data. In that case, a
fitted common radial scale with a large transition width (see the dash-dotted
blue line in Fig. 16) might prove to be a more accurate description of the
global structure of the GMF.

In summary, the unified halo model describes the data very well using only 6
free parameters for the halo, instead of the 10 parameters needed for the base
model. This is the first quantitative demonstration of earlier conjectures,
that the toroidal halo of the Milky Way likely arose dynamically from
differential rotation.