11institutetext: Université Libre de Bruxelles, CP225 Boulevard du Triomphe,
1050 Brussels, Belgium 11email<EMAIL_ADDRESS>22institutetext:
Université de Paris Cite, CNRS, Astroparticule et Cosmologie, F-75013 Paris,
France

# The coherent magnetic halo of Milky Way

A.Korochkin1 D.Semikoz2 and P.Tinyakov1

###### Abstract

Context. Recent catalog of Faraday rotation measures (RM) of extragalactic
sources together with the synchrotron polarization data from WMAP and Planck
provide us with the wealth of information on magnetic fields of the Galaxy.
However, the integral character of these observables together with our
position inside the Galaxy makes the inference of the coherent Galactic
magnetic field (GMF) complicated and ambiguous.

Aims. We combine several phenomenological components of the GMF — the spiral
arms, the toroidal halo, the X-shaped field and the field of the Local Bubble
— to construct a new model of the regular GMF outside of the thin disk.

Methods. To have control over the relative contributions of the RM and
polarization data to the fit we pay special attention to the estimation of
errors in data bins. To this end we develop a systematic method which is
uniformly applicable to different data sets. This method takes into account
individual measurement errors, the variance in the bin as well as fluctuations
in the data at angular scales larger than the bin size. This leads to decrease
of the errors and, as a result, to better sensitivity of the data to the model
content. We cross checked the stability of our method with the new LOFAR data.

Results. We found that the four components listed above are sufficient to fit
both the RM and polarization data over the whole sky with only a small
fraction masked out. Moreover, we have achieved several important improvements
compared to previous approaches. Due to account of our location inside of the
Local Bubble our model does not require introduction of striated fields. For
the first time we showed that the Fan Region can be modeled as a Galactic-
scale feature. The pitch angle of the magnetic field in our fit converged to
the value around 20 degrees. Interestingly, with value is very close to the
direction of arms inferred recently from Gaia data on upper main sequence
stars.

###### Key Words.:

Galaxies: magnetic fields, ISM: magnetic fields, cosmic rays

##  1 Introduction

It is well known that our Galaxy is permeated with magnetic field. Such a
field is intricately linked to the evolution of the galaxy, affecting gas
dynamics and, in particular, influencing the star formation, propagation of
cosmic rays, and shaping the Galactic emission at radio frequencies. Accurate
account of the Galactic magnetic field (GMF) is one of the keys for solving
the long standing puzzle of sources of ultra-high-energy cosmic rays (UHECR).
However, the structure, origin and evolution of the GMF is still not well
understood.

Apart from the short-scale turbulent field, radio observations of external
galaxies indicate the presence of a large-scale coherent field ordered on
scales comparable to the size of the galaxy (Beck, 2015). Interestingly, two
essentially different coherent components can be identified. For the first
component, the magnetic field lies in the plane of the galactic disk and forms
so called magnetic arms, which are not necessarily coincident with the gaseous
or optical spiral arms. These spiral magnetic structures are found in almost
every spiral galaxy observed so far (e.g., magnetic arms of M51 (Fletcher et
al., 2011) or the disk field of M101 (Berkhuijsen et al., 2016)). The second
component, on the contrary, extends out of the disk into the galactic halo and
forms an X-shaped pattern which can be traced up to the distances of several
kpc away from the galactic plane (e.g., the halo field of NGC 4631 (Mora &
Krause, 2013) or NGC 5775 (Soida et al., 2011)).

In the standard paradigm, the galactic coherent field is explained as a result
of action of large-scale dynamo (for a review see Beck et al. (2019) and
Brandenburg & Ntormousi (2023)) which can also link together the coherent
field of the disk and the X-shaped field of the halo (Moss & Sokoloff, 2008).
Alternatively, it can be the result of the gravitational compression of a
primordial magnetic field of nanogauss strength occurred during structure
formation (Howard & Kulsrud, 1997).

While other galaxies are observed from outside and their global magnetic field
structure can be relatively easily recognized, our location inside the
Galactic disk makes such measurements for the Milky Way a non-trivial problem.
Typically used tracers of magnetic field such as Faraday rotation measures
(RM) of polarized sources or synchrotron emission of cosmic ray electrons
(CRE) give information integrated over the line of sight. Thus, one has to
solve the inverse problem to reconstruct the global 3-dimensional GMF
structure from sky projections of several GMF-sensitive observables. This is
not possible without assumptions about the global GMF structure.

The situation is further complicated by the fact that these integral tracers
get contributions from magnetized and ionized bubbles located in the Galactic
disk (Wolleben et al., 2010) typically within $100-200$ pc from the Galactic
plane. In fact, the solar system itself is located inside one of such bubbles
usually referred to as the Local Bubble (Pelgrims et al., 2020), which is
generally believed to be remnant of supernova explosions. Several other nearby
bubbles seen from the exterior can also be identified. They produce large-
scale anomalies in the RM and synchrotron polarization data at high Galactic
latitudes that have to be subtracted or masked. The problem gets worse as more
bubbles contribute for the lines of sight closer to the Galactic plane, and
appears untractable without invoking distance-sensitive (tomographical)
methods. As a result, the information on the global structure of the magnetic
field in the disk, particularly in the direction of central Galaxy, is very
uncertain.

Several phenomenological models of the coherent GMF have been proposed
previously which differ by complexity, underlying data used, and fitting
techniques. First generation of models used synchrotron data and earlier
compilations of rotation measures of pulsars and extragalactic sources (Han et
al., 1997; Tinyakov & Tkachev, 2002; Beck, 2001). After the arrival of the
first all-sky rotation measure catalog of extragalactic objects (Taylor et
al., 2009) the next generation of models included it as the main data source
for reconstruction of the magnetic field in the halo (Pshirkov et al., 2011).
The combination of the rotation measures with the synchrotron polarisation
data gives an additional powerful tool for model building (Jansson & Farrar,
2012a). Recent models using various combinations of the existing data include
Han et al. (2018); Xu & Han (2019); Unger & Farrar (2023); Xu & Han (2024).

The purpose of this paper is to present a new model of the GMF with a focus on
the halo field, by which we understand in what follows the field outside of
the thin disk. There are several motivations to update the existing models. On
one hand, a significant amount of new RM data has recently appeared. For
instance, the S-PASS/ATCA survey Schnitzeler et al. (2019a) filled the blind
spot of the RM catalog Taylor et al. (2009) used in previous studies. The
Wilkinson Microwave Anisotropy Probe (WMAP) synchrotron polarization data have
been updated in Bennett et al. (2013) and more recently in Watts et al.
(2023). On the other hand, even the most recent existing models leave a lot of
room for improvement. The model of Pshirkov et al. (2011) did not include the
synchrotron polarization data. The model of Jansson & Farrar (2012a) assumed
the magnetic field to be to some degree ”striated” (flipping direction in
random domains) which is essentially a way to boost the synchrotron
polarization signal with respect to the rotation measures by an arbitrary
factor treated as a phenomenological parameter. The last but not the least,
there is a large freedom in choosing and parameterizing the main components of
the GMF. Comparing the results of independently constructed fits gives some
understanding of the systematic uncertainties involved. While this paper was
in preparation there appeared a new model of the GMF by Unger & Farrar (2023);
their results and those presented here should be considered complementary in
this sense.

In this paper we use the WMAP synchrotron polarization measurements at 23 GHz
(Bennett et al., 2013) and the catalog of about 59000 RMs of extragalactic
sources (Van Eck et al., 2023) to fit the parameters of our model. Even though
mostly the same data were used in previous models including the recent model
by Unger & Farrar (2023), we have made several essential improvements.

First, we model and take into account the contribution of the Local Bubble.
Its importance was first pointed out by Alves et al. (2018) when considering
polarized dust emission observed by Planck at 353 GHz. In previous models that
used the RM and synchrotron data together it was found that the magnetic field
preferred by the RM data alone is insufficient by a factor of 1.5-3 (Unger &
Farrar, 2023) to explain the observed polarized intensity at high Galactic
latitudes. To reconcile the two sets of data the ”striation” parameter was
introduced. Here we show that a sizeable contribution of the Local Bubble to
the polarization signal at 23 GHz solves this ”synchrotron deficit” problem
while being fully compatible with the RM data. Using the simplest model of the
Local Bubble, we achieved a significant improvement in the quality of the fit,
especially at high Galactic latitudes $|b|>60^{\circ}$, thus completely
eliminating the need for the ”striated” magnetic fields.

Second, we have developed a uniform approach to the error estimation which,
although technical, is a crucial step in fitting together data of a different
nature, the RM and synchrotron polarization in our case. We therefore have
control of a relative importance of the RM and synchrotron polarization data
in the global fit. The key observation here is that one of the main origins of
the ”noise” in the data is common in both cases — the magnetized bubbles in
the disk. In general, our error bars are smaller than obtained by previously
used methods.

Third, we found from our fit that the pitch angle of the local magnetic arms
is $20^{\circ}$. This pitch angle is significantly larger than found in
previous studies where the pitch angle did not exceed $\sim 10^{\circ}$. This
large value is in perfect agreement with the pitch of the Galactic arms
inferred from Gaia data by Poggio et al. (2021), as is illustrated below in
Fig. 3.

Finally, combining magnetic field in the Local arm with that of the Perseus
arm we, for the first time, obtained a good fit of the Fan region — a large
bright region on the synchrotron polarization maps roughly in the direction of
the Galactic anti-center. According to recent studies (Hill et al., 2017) the
polarization signal in this region is produced at $1-2$ kpc from the solar
system and thus makes part of the global structure (for the model of the Fan
Region as a local feature see West et al. (2021)). In previous GMF models this
region was masked as a local anomaly.

The rest to of the paper is organized as follows. In the Sect. 2 we describe
our method for calculation of errorbars of rotation measure and synchrotron
data and discuss data we used. In Sect. 3 we describe our model for Galactic
Magnetic Field in the halo. In Sect. 4 we present the results of the global
fit of our model to the data. In Sect. 5 we discuss the performance of our
model. In particular, in Sect. 5.1 we go in details of new feature introduced
in the present model for the first time, namely the Local Bubble. Finally we
sum up our results in the Sect. 6. The Appendices A and B are devoted to the
method of error assignment. In the Appendix C we compare our results to the
previous models.

Throughout the paper and in the numerical code, we fix the coordinate system
such that the center of the Galaxy is at the origin, the $z$-axis is directed
towards the Galactic North pole, and the solar system coordinates are {-8.2
kpc, 0, 0}. For the observer at Earth the direction towards positive $y$ is
$(l,b)=(90^{\circ},0^{\circ})$, towards positive $z$ –
$(l,b)=(0^{\circ},90^{\circ})$.

##  2 Data preparation

Galactic magnetic field can be probed with different tracers. Among them are
polarized thermal emission of dust, Faraday rotation measures of Galactic and
extragalactic sources and synchrotron emission of cosmic ray electrons at
different frequencies. Each of these tracers is sensitive to different
properties of magnetic field but all of them provide us with valuable
information about coherent GMF structure. Below we will briefly review each
tracer and justify the choice of datasets for our analysis.

#### Dust.

Observations of dust polarized emission reveal the average orientation of the
component $B_{\perp}$ of GMF perpendicular to the line of sight, without
providing any information on its strength. The dust is concentrated in the
thin disk with the height $\sim 100$ pc, thus making the method mainly
sensitive to our local neighborhood. Significant progress has been achieved
thanks to the Planck measurements at a frequency of 353 GHz, where dust
emission dominates. Firstly, it was shown that none of the existing models can
accurately predict dust emission (Planck Collaboration et al., 2016).
Secondly, it was demonstrated that the emission at high Galactic latitudes is
dominated by the Local Bubble wall (Skalidis & Pelgrims, 2019; Pelgrims et
al., 2020). Moreover, it was found that a model where the field in the wall is
obtained by compressing an initially uniform field successfully explains the
structure of polarization after fitting the direction of the initial field. We
adopted this approach when building our model of the Local Bubble, see Sect.
3.

#### Rotation measures: Galactic pulsars.

Faraday rotation is the phenomenon of rotation of the polarization plane of a
linearly polarized radiowave when propagating through the magnetized plasma.
The rotation angle $\Delta\theta$ is proportional to the square of the
wavelength $\lambda$; the coefficient is called the rotation measure (RM):

$\Delta\theta={e^{3}\lambda^{2}\over 2\pi
m_{e}^{2}}\int_{0}^{l}n_{e}B_{||}ds\equiv\mathrm{RM}\cdot\lambda^{2},$ (1)

where $e$ and $m_{e}$ are the electron charge and mass, and the integral is
taken along the line of sight between the observer and the source. The RM is
proportional to the density of free electrons $n_{e}$ and the component of the
magnetic field parallel to the line of sight $B_{||}$,

$\mathrm{RM}\approx
0.812\,\,\int_{0}^{l}\left[\frac{n_{e}(s)}{\mathrm{cm^{-3}}}\right]\left[\frac{B_{||}(s)}{10^{-6}\,\mathrm{G}}\right]\,\left[\frac{\mathrm{d}s}{\mathrm{pc}}\right]\,\,\,\,\mathrm{rad/m^{2}}.$
(2)

Positive RM corresponds to the field pointing towards the observer.
Measurement of the polarization on at least two different wavelengths allows
one to calculate RM.

Figure 1: Compilation of the data used in this paper. Upper panel:
extragalactic RM skymap, prepared by Hutschenreuter et al. (2022). Middle and
bottom panels: WMAP synchrotron Stokes Q and U skymaps at 23 GHz (Bennett et
al., 2013).

Pulsars are bright sources of polarized radio emission. The advantage of
pulsars is their negligibly small intrinsic RM values. Also their location
within the galaxy allows for the investigation of the inner structure of the
GMF. Study of pulsar RMs gives an indication of the magnetic field reversal in
the direction to the Galactic center (Han & Qiao, 1994). However, pulsars are
mostly concentrated in the disk, which makes them not so useful for studying
the GMF halo. Moreover, the distances to pulsars are usually inferred from
their dispersion measures, which introduces additional uncertainties. For
these reasons we do not use pulsar data in this paper and leave the validation
of our model against pulsars for the future.

#### Extragalactic RM.

At the moment more than 59000 extragalactic RMs have been measured which are
summarized in publicly available master catalog by Van Eck et al. (2023). In
our study we use the latest compilation of this catalog, namely CIRADA
consolidated catalog version v1.2.0 111The catalog can be found online at the
link https://github.com/CIRADA-Tools/RMTable. The core of the catalog is the
NRAO VLA Sky Survey (NVSS) survey (Taylor et al., 2009). Additionally, it
comprises RMs from Schnitzeler et al. (2019b); Van Eck et al. (2021); Betti et
al. (2019); Farnes et al. (2014); Mao et al. (2010); Tabara & Inoue (1980);
Simard-Normandin et al. (1981); Broten et al. (1988); Riseley et al. (2020);
Brown et al. (2003); Taylor et al. (2024); Mao et al. (2012b, a); Feain et al.
(2009); Van Eck et al. (2011); Ma et al. (2020); O’Sullivan et al. (2017);
Kaczmarek et al. (2017); Anderson et al. (2015); Brown et al. (2007); Klein et
al. (2003); Heald et al. (2009); Shanahan et al. (2019); Clarke et al. (2001);
Minter & Spangler (1996); Ranchod et al. (2024); Van Eck et al. (2018); Law et
al. (2011); Riseley et al. (2018); Livingston et al. (2022); Mao et al.
(2008); Roy et al. (2005); Livingston et al. (2021); Oren & Wolfe (1995);
Clegg et al. (1992); Kim et al. (2016); Battye et al. (2011); Ma et al.
(2019); Rossetti et al. (2008); Costa & Spangler (2018); Vernstrom et al.
(2018); Gaensler et al. (2001); Costa et al. (2016), including the most recent
extensive Low-Frequency Arrar (LOFAR) survey LoTSS (O’Sullivan et al., 2023)
and the Australian Square Kilometre Array Pathfinder survey POSSUM
(Vanderwoude et al., 2024).

The totality of extragalactic RMs covers the whole sky almost uniformly with
roughly 1 source per 1 deg2 on average, which is very convenient for studying
the global structure of the GMF. The smoothed and cleaned RM skymap
constructed by Hutschenreuter et al. (2022)222Note that the RM skymap of
Hutschenreuter et al. (2022) while cleaned of outliers is not suitable for our
purposes since the correlations between individual pixels of the map do not
allow for the estimation of the errorbars of larger bins, see Sec 2.1. from
the previous version of the catalog (v0.1.8) is shown in Fig. 1, upper panel.
An even earlier version of the catalog was used for the first time to study
magnetic field in the Galactic halo by Pshirkov et al. (2011).

Unlike pulsars, extragalactic sources generally do have intrinsic RMs which
sometimes may be very large, so the RM data have to be cleaned from these
outliers. The exact procedure is described in the Appendix A. The bulk of the
extragalactic sources have intrinsic RMs with the zero mean and the r.m.s. not
exceeding $\sim 7$ rad/m2 (Schnitzeler, 2010; Pshirkov et al., 2013).

#### Synchrotron.

Relativistic electrons of energy $E$ spiraling in the magnetic field emit
radio waves near the critical frequency

$\nu_{c}\approx
1.6\,\left[\frac{B_{\perp}}{10^{-6}\,\mathrm{G}}\right]\left[\frac{E}{10\,\mathrm{GeV}}\right]^{2}\,\,\,\mathrm{GHz}.$
(3)

The resulting synchrotron radiation is polarized in the plane perpendicular to
the direction of the magnetic field. Assuming the distribution of cosmic ray
electrons (CRE) in energy is known, the information about $B_{\perp}$ can be
extracted from synchrotron measurements. This information is complementary to
that contained in RMs.

The CMB experiments WMAP and Planck performed the full-sky synchrotron
polarization measurements with high precision at 23 GHz and 30 GHz,
respectively. The WMAP synchrotron skymap was first combined with RM data to
study the GMF in Jansson et al. (2009); Sun & Reich (2010); Jansson & Farrar
(2012a, b).

At higher frequencies the polarized emission is dominated by the emission from
dust. On the other hand, at lower frequencies the Faraday rotation becomes
important as can be seen from equations 1 and 3. Significant Faraday rotation
leads to depolarization of initially polarized beam. The combination of
Faraday rotation with the synchrotron measurements at large number of nearby
frequencies allows for the so-called Faraday tomography (Burn, 1966), when the
synchrotron emissivity is determined as a function of RM accumulated along the
line of sight. Recently such measurements were performed by GMIMS
collaboration in the Northern sky in the frequency band 1280 - 1750 MHz
(Wolleben et al., 2021), and in the Southern sky in the range 300 - 480 MHz
(Wolleben et al., 2019).

In the current study we use the final nine-year WMAP 23 GHz dataset to fit the
model. We choose it because firstly, it covers the entire sky, and secondly,
it is measured at sufficiently high frequency so that depolarization is
negligible. The use of frequencies at which depolarization is significant
implies an increasingly significant role of local structures which we do not
model in this study. On the other hand WMAP 23 GHz and Planck 30 GHz data are
very similar after rescaling for the difference in frequency. We could use
either one; we choose WMAP because its seven-year version has already been
used in Jansson & Farrar (2012a). The residual systematic difference after
rescaling between WMAP and Planck was studied in detail in the Cosmoglobe
project (Watts et al., 2023). We note that this difference is small compared
to the current precision of the GMF models, and thus is not important for our
purposes. Figure 1 shows the sky maps of Stokes Q and U parameters (middle and
bottom panels, respectively) corresponding to the WMAP 23 GHz dataset.

###  2.1 Binning and error estimation

In order to extract maximum benefit from the combination of the RM and
synchrotron polarization data, we prepare both datasets in a similar way. This
includes four steps: removal of outliers, masking, binning and error
estimation. The procedures for outliers and masking are different for RMs and
synchrotron skymaps and are discussed separately in the subsections 2.2 and
2.3. On the other hand, for both datasets we adopt the same binning scheme and
errors estimation algorithm which are described below.

After getting rid of outliers and applying the masks we are left with a
collection of small pixels (or individual RMs) that can be directly used for
modelling. However, as we are interested in the coherent part of the GMF it is
useful to average the masked RM and synchrotron data over the angular bins of
a relatively large size. Indeed, the signal of coherent GMF is expected to
vary on scales of tens of degrees, and so bins of a smaller, but not much
smaller angular size are sufficient to trace it. Moreover, averaging over
large bins effectively cancels small-scale fluctuations in the data and
reduces errorbars.

The same approach was already used in previous studies. Pshirkov et al. (2011)
and Terral & Ferrière (2017) used the bins with the size $\sim 10^{\circ}$,
Jansson & Farrar (2012a); Unger & Farrar (2023), used smaller bins. Since
there is no strict criterion for choosing the optimal bin size we use the same
bins as Pshirkov et al. (2011) for both RM and synchrotron data, see Fig. 2.
The advantage of this binning scheme is that bins form isolatitude bands. The
binning scheme was chosen independently of the masks. Those bins where more
than 50% of the original content (by area for WMAP and by the number of
sources for RM) has been masked were discarded.

Figure 2: The binning scheme used in this analysis. The bins are rectangular
in Galactic coordinates and are arranged in iso-latitude bands of $10^{\circ}$
width. The two equatorial bands are divided into 36 square $10^{\circ}\times
10^{\circ}$ bins each. The number of bins in other bands is chosen in such a
way as to get the closest approximation to equal-area bins.

After setting the binning scheme the most important step is the estimation of
the mean value in each bin and assigning corresponding error bars. Jansson &
Farrar (2012a) and Unger & Farrar (2023) used the standard deviation within
the bins as an error estimate. We do not see a justification for this choice.
If the variations of individual sources (or pixels in the case of polarization
data) from the mean were all independent, the error in the bin should be the
standard deviation divided by the square root of the number of sources,
typically much smaller than the standard deviation itself. However, bubbles
and other structures of large angular size (e.g. Radcliffe Wave (Panopoulou et
al., 2024)), which are also part of the noise, give coherent contributions to
individual bins, so the errors should be larger, but likely not as large as
the standard deviation. Pshirkov et al. (2011) used a more complex error
estimation procedure which for many bins resulted in smaller errors equal to
the standard deviation divided by three. Such a procedure, however, lacks a
statistical justification and thus is difficult to generalize to the case of
two data sets of a different nature.

Figure 3: Schematic picture of the Galactic arms as viewed from the North
Galactic pole. Black spiral curves shows the axes of the best-fit spiral
magnetic arms of our model. The intensity of the color indicates the magnetic
field strength in the Galactic plane $z=0$. The blue (orange) colored regions
correspond to the magnetic field directed inward/clockwise (outward/counter-
clockwise) as marked by the black arrows at the bottom of the plot. Black dot
marks the Galactic center while an ellipse represents the Galactic bar. Red
star marks the position of the Sun. Grey regions show the spiral arms as
deduced from Gaia observations. The pitch angle in the solar vicinity is
$20^{\circ}$.

In this work we propose a new method of error estimation. The observed
quantities — RMs or polarized intensity maps — are sums of contributions from
the coherent GMF and from fluctuations of different origin, including
inhomogeneities of the interstellar medium (ISM), the turbulent component of
the GMF, magnetized or ionized clouds, intrinsic RMs of sources, etc., and of
different scales including those exceeding the bin size. We do not model such
fluctuations and consider them as noise obscuring the coherent GMF signal. Our
method to estimate the contribution of this noise to the errors assigned to
bins is based on the assumption that its statistical properties only depend on
the latitude $b$. We assess the parameters of the noise at given $b$ by the
Fourier analysis of the data in the isolatitude strips, which we then use to
estimate the errors.

If the data in a strip only contained noise, then the variance associated with
a bin of size $L$ would be given by the sum over modes of the power spectrum
with the window function $\operatorname{sinc}^{2}(kL/2)$, see Appendix B for
the derivation. Note that this variance includes contributions from scales
larger than the bin size. The actual data also contain the coherent part which
is concentrated in the lowest modes and should be subtracted. In fact, the
model that we will fit to the data only has sizeable power in the first three
harmonics $k=0,1,2$, so we subtract these harmonics before calculating the
variance, see Fig. 10. Specifically, for all bins of the same strip we
calculate the variance $\sigma_{\mathrm{L}}$ according to the formula:

$\sigma^{2}_{\mathrm{L}}=2\sum_{k_{0}}^{\infty}\,\operatorname{sinc}^{2}\left(\frac{kL}{2}\right)S_{k}$
(4)

where $S_{k}$ is the power spectrum of the strip of data to which the bin
belongs, and the bin size $L$ is of order 10∘. Eq. 4 results in errors that
are smaller than the simple data variance over the bin by a factor 2-3. Note
also that this approach allows us to treat the RM and synchrotron polarization
data on equal footing.

Figure 4: Section of the best-fit GMF model by the plane $y=0$, which is
perpendicular to Galactic plane and come through Sun (red star) and Galactic
center. Blue and orange regions are the magnetic arms, northern and southern
toroidal field and the Local Bubble. The color map is the same as in Fig. 3.
Black arrows show the direction of the X-shaped field. Red star marks the
position of the Sun.

Eq. 4 assumes that the mean value in a given bin can be calculated exactly,
while in reality this is not the case: the data only sample the underlying
function in a number of points, so the mean calculated from these values has
an associated statistical error $\sigma_{\mu}$. We determine this error by
fitting the distribution of values in the bin by a Gaussian, which gives both
the mean and its statistical error: $\mu\pm\sigma_{\mu}$. The total error in
the bin is thus

$\sigma_{\text{tot}}=\sqrt{\sigma^{2}_{\mathrm{L}}+\sigma_{\mu}^{2}}.$ (5)

For the polarization data, and for most (but not all) of the RM bins the
contribution of $\sigma_{\mu}$ is negligible.

We note that the errors $\sigma_{\mathrm{L}}$ and $\sigma_{\mu}$ have
different behavior when increasing the number of observations $N$. The error
on the mean value of the bin $\sigma_{\mu}$ tends to zero approximately as
$\sim$ 1/$\sqrt{N}$, while $\sigma_{\mathrm{L}}$ remains constant.

The details of the calculation of $\sigma_{\mu}$ and $\sigma_{\mathrm{L}}$ are
described in the Appendices A and B, respectively. The procedure is the same
for both RM and synchrotron datasets.

###  2.2 RM mask

The value of the RM in some directions may be dominated by local features of
the ISM. Thus, when fitting the coherent GMF they should be masked out. An
illustrative example is the well-known Sh2-27 cloud in the direction
$(l,b)\approx(10^{\circ},25^{\circ})$, ionized by the massive star Zeta
Ophiuchi passing through it. This cloud is located at the distance of $\sim$
100 pc and is clearly visible as a $\sim 10^{\circ}$ bright blue spot in the
RM map, see Fig. 1. The average RM in the direction of Sh2-27 goes below
$-100$ while RMs in its vicinity are close to zero (Thomson et al., 2019).

Since RM is proportional to the density of thermal electrons, potentially
biased regions are expected to correlate with highly ionized nearby bubbles.
In order to mask them out it is convenient to use maps of H$\alpha$ emission.
Similar approach was recently used by Unger & Farrar (2023). Using the full-
sky H$\alpha$ map compiled by Finkbeiner (2003) we masked pixels where the
intensity of H$\alpha$ emission is at least twice as strong as the average
intensity at the same latitude.

Additionally, based on visual inspection of the RM map (see Fig. 1) we masked
out the arc-like segment with positive RMs around
$(l,b)\approx(100^{\circ},0^{\circ}-50^{\circ})$ which is probably related to
the Loop III. Also we removed from the analysis bright bubble-like structure
at $(l,b)\approx(160^{\circ},-15^{\circ})$ and a region with mostly zero RM
surrounded by strongly negative RMs near
$(l,b)\approx(110^{\circ},-10^{\circ})$.

Finally, we excluded the Galactic plane ($|b|<10^{\circ}$) because its
modelling depends on the enhancement of thermal electrons density in the
spiral arms. Unfortunately, existing models of thermal electrons rely on pre-
Gaia measurements of the spiral arms whose pitch was believed to be
significantly smaller than that inferred from the Gaia data (Poggio et al.,
2021). Our preliminary fits, as well of the final result, favor the same pitch
as from Gaia, thus making these models incompatible with our setup. The
description of the adopted model of thermal electrons is given in Sect. 3.2.

The resulting mask is shown in Fig. 6. In total the RM dataset mask covers 26%
of the sky.

###  2.3 Synchrotron mask

The synchrotron mask is built independently of the RM mask. First, we
discarded all pixels of the WMAP Stokes Q and U maps marked as corrupted. Then
we masked out the brightest, most probably nearby features. Namely, we
excluded Loop I and Loop II by cutting out corresponding ring segments. In
contrast with previous studies we do not mask the Galactic plane and the Fan
region. We apply this mask for both Stokes Q and U skymaps, see Fig. 6. The
mask covers 11% of the sky.

##  3 Model of the Galactic Magnetic Field

Sky maps of both Faraday rotation and synchrotron polarization data reveal
large regions of constant sign that are aligned with the Galactic plane and
the direction to the Galactic center, leaving no doubt of the existence of the
coherent magnetic field in the Galaxy. These regions form a particular pattern
which has to be reproduced by the magnetic field model. Even though we only
measure the projections onto the sky plane of 3-dimensional structures, this
requirement turns out to be rather restrictive if one assumes that the model
consists of a few simple components. We found that 4 major components are
sufficient to describe the global structure of the Faraday rotation and
polarization maps: the thick disk, the toroidal halo, the X-shape halo field
and the field of the local bubble. Different parts of the data are more
sensitive to different components, but all four are needed for a good global
fit. In this section we overview these components and their role in
reproducing the key features of the data, and present the results of the
global fit.

###  3.1 Magnetic field components

#### Thick disk.

While most of Galactic stars are concentrated within a couple hundred parsecs
from the Galactic plane where they populate the Galactic arms, thermal
electrons and coherent magnetic fields are usually assumed to form a larger
structure often referred to as the thick disk which generally traces the
spiral arms and extends up to 1-2 kpc from the Galactic plane. We assume in
this work that the magnetic field associated with this structure is directed
along the arms either towards or away from the Galactic center. The solar
system is situated at the inner edge of the Local arm; the ”local” magnetic
field is directed clockwise as viewed from the Galactic North pole, with a
pitch towards the center. The detailed structure of the spiral arms in the
solar vicinity has been mapped out by Gaia (Poggio et al., 2021) and will be
discussed further in Sect. 4. There are indications from the pulsar data of
the field reversal at $\sim 0.5$ kpc in the direction to the Galactic center
(Han & Qiao, 1994). We therefore assume that the field changes direction to
counter-clockwise in the inner Sagittarius-Carina arm.

The typical thickness of the disk is $\sim 1$ kpc, which means that only the
region within $\sim 6$ kpc from the Sun contributes to the observables at
Galactic latitudes $|b|>10^{\circ}$. Correspondingly, only the arms that pass
in the relative vicinity of the Sun are constrained from our fits to the data.
These are, in the first place, the Local and Perseus arms in the outer Galaxy
and the Sagittarius-Carina and Scutum arms in the inner Galaxy. Our fits are
not sensitive to the field in the Norma and Outer arms. The general shape of
the thick disk field is shown in Fig. 3 where the colored region corresponds
to distances within $6$ kpc from the Sun.

In our model each arm is an independent GMF component. The position of an arm
is given by its axis, represented by the logarithmic spiral:

$x=x_{0}+be^{k\phi}\cos{\phi},\quad y=y_{0}+be^{k\phi}\sin{\phi},$ (6)

where $k=\tan\alpha$ is the tangent of the pitch angle $\alpha$, $b$ is an
initial phase, and $x_{0}$ and $y_{0}$ are the origin of the spiral. The cross
section of an arm in the plane perpendicular to its axis is a squircle: the
point of the cross section plane belongs to an arm if its coordinates in that
plane $\\{\tilde{x},z\\}$ satisfy

$\left|\frac{\tilde{x}}{r_{\mathrm{disk}}}\right|^{n}+\left|\frac{z}{r_{z}}\right|^{n}\leq
1.$ (7)

Here $r_{\mathrm{disk}}$ and $r_{z}$ are two semi-axes of the squircle. The
continuous parameter $n$ controls the shape of the cross section. If $n=2$ the
squircle reduces to an ellipse, while larger $n$ makes it more similar to a
rectangle.

The magnetic field of an arm is directed along the arm and varies in inverse
proportion to the arm cross section as to guarantee its divergence-free
structure. It is parameterized by the value $B$ at the intersection of the arm
axis and the line passing through the Sun and the Galactic center. To minimize
the number of fitting parameters we take the field to be constant over the arm
cross section.

In total each arm is described by 8 parameters: $x_{0}$, $y_{0}$, $b$,
$\alpha$, $B$, $r_{\mathrm{disk}}$, $r_{z}$, $n$. The pitch angles $\alpha$
are taken the same for all arms and are controlled by the common pitch
$\alpha_{0}$ which is allowed to vary freely while fitting. The
$y$-coordinates of the arm origins are zero $y_{0}=0$ for all arms. To the
contrary, the $x$-coordinates $x_{0}$ are individual for each arm and are set
by hand to make our magnetic arms better match with the spiral design of the
Galaxy. The adjustments of $x_{0}$ do not affect the result as long as $x_{0}$
is small compared to the distance from the Sun to the Galactic center. Indeed,
in this case the change in $x_{0}$ is degenerate with the initial phase of the
spiral arm $b$ which is a fitting parameter. We make use of this fact to
roughly connect major Perseus and Scutum arms with the ends of the Galactic
bar. The widths of the arms in the Galactic plane are also set to predefined
values close to 1 kpc, see Table 5.

After these adjustments we are left with 4 fitting parameters per arm. For the
nearby Local arm and the Sagittarius-Carina arm we fit all 4 parameters, while
for the more distant Scutum arm only the strength $B$.

The Perseus arm is a special case. It contributes to the observables in the
direction of the outer Galaxy at low Galactic latitudes where the Fan Region
is located, which is one of the most prominent features on the synchrotron
polarization maps. In order to fit the Fan Region we had to assume that the
inner part of the Perseus arm close to the disk has stronger magnetic field as
shown in Fig. 4. We assumed for simplicity that the field takes two different
constant values across the inner and outer sections of the arm (this can be
viewed, of course, as a crude approximation for the field decreasing away from
the disk). In total, there are three fitting parameters for the Perseus arm:
the field strength $B$, the additional field in the inner part $B^{\prime}$
(so that the total field in this part is $B+B^{\prime}$) and the width of the
arm $r_{\mathrm{z}}$. The width of the inner part $r_{\mathrm{z}}^{\prime}$ is
fixed to 0.4 kpc. The phase $b$ and shift $x_{0}$ are chosen in such a way as
to make the Perseus and Scutum arms symmetric with respect to the Galactic
$z$-axis: $x_{0}(\mathrm{Perseus})=-x_{0}(\mathrm{Scutum})$,
$b(\mathrm{Perseus})=b(\mathrm{Scutum})+180^{\circ}$, see Fig. 4 and Table 5.

#### Toroidal halo.

The field in the thick disk is not sufficient to account for the observed
Faraday rotations — a component extending further in the Galactic halo is
needed. Previous models have shown the necessity for a toroidal component
above and below the Galactic disk with the azimuthal field directions opposite
in the North and South. Such a field may arise as a result of the Galactic
dynamo (Wielebinski & Krause, 1993). This component is schematically shown in
Fig. 4.

It has been first noted in Han et al. (1997, 1999) and later confirmed by
Pshirkov et al. (2011) that the combination of the thick disk and the toroidal
halo reproduces well the characteristic all-sky $-+-+/-+$ North/South pattern
of the Faraday rotations, provided the directions of the disk field around the
Sun location and the toroidal component are opposite in the North and
coincident in the South. As a result, the magnetic field in the Galactic
Northern sky is somewhat weaker than in the South. If only Faraday rotation
data are included, these two components give a reasonable fit to the data
(Pshirkov et al., 2011).

We model the toroidal components as cylindrical sections whose axes coincide
with the Galactic axis, filled with purely azimuthal magnetic field. Each
cylinder is described by 4 parameters: the outer radius $r_{\mathrm{out}}$,
lower bound $z_{\mathrm{min}}$, upper bound $z_{\mathrm{max}}$ and the
strength of the field $B_{\mathrm{tor}}$.

#### X-shape field.

When the synchrotron polarization data are included, we could not fit both RM
and synchrotron maps with the thick disk and toroidal halo. With only these
two components one cannot reproduce the change of the sign of the Q-parameter
around the meridian line $l=0$, which is particularly visible in the data in
the Northern hemisphere, see Fig. 1. As was first noted in Jansson & Farrar
(2012a), this can be achieved by adding a second halo component — the X-shape
field — in the inner part of the Galaxy. Topologically this field corresponds
to the Model C of Ferrière & Terral (2014). Similar fields are observed in
some other galaxies viewed edge-on (Beck et al., 2019). To reduce the number
of parameters, we take the X-field lines to be symmetric with respect to the
Galactic plane, axially symmetric and having everywhere the same inclination
angle $\theta$ with respect to the Galactic axis. The footprint of the X-field
in the disk extends up to $r_{X}\sim 6$ kpc from the center. The X-field is
schematically shown in Fig.4.

The fitting parameters of the X-field are the strength $B_{X}$ at $z=0$, the
radius $r_{X}$ of the X-field footprint at $z=0$ and the inclination angle
$\theta$. The radial profile of the X-field was chosen to be constant. Thus,
the footprint of the X-field on a plane above or below the Galactic plane is a
ring of inner radius $|z|/\cos\theta$ and outer radius $|z|/\cos\theta+r_{X}$.
The area of this ring grows with $|z|$; the strength of the field decreases in
inverse proportion so as to keep the field divergence-free.

#### Local Bubble.

The three components described so far reproduce well the overall pattern of
both RM and synchrotron data. However, when normalized to the rotation
measures, they underestimate the overall strength of the synchrotron
polarization signal at high latitudes by about factor 2. In the previous
models of Jansson & Farrar (2012a); Unger & Farrar (2023) where the
synchrotron data were included in the fit this problem was solved by
introducing the ’striation’ which allows one to boost the synchrotron
polarization signal relative to the rotation measures. We found that inclusion
of the contribution from the Local Bubble may solve the problem without the
striation factor.

Figure 5: Schematic drawing of the Local Bubble and surrounding magnetic field
lines. The radius of the bubble in our model is fixed to 200 pc, while the
best-fit thickness of its wall was found to be 30 pc.

The Local Bubble is a cavity in the interstellar medium created, presumably,
by several local supernova. It is surrounded by clouds of cold dust compressed
by hot gas resulted from supernova explosion(s). It has an irregular shape,
but in the zeroth approximation can be fitted by an ellipsoid with half-axes
of the size of $x=100$ pc by $y=200$ pc by $z=300$ pc (Alves et al., 2018;
Pelgrims et al., 2020). The Solar system is located in the inner part of the
bubble relatively close to its center.

It turns out that even in a crude approximation of a spherical bubble, the
compressed magnetic field on the bubble walls can explain the missing part of
the polarization intensity while having approximately right sky distribution
and not spoiling the RM signal. The latter can even be improved by adjusting
the parameters of the bubble as discussed in more detail below in Sect.5.1.
This is possible because the polarization parameters $Q$ and $U$ are quadratic
in the field strength, while the RMs are linear.

In our model the bubble is a spherical shell with the inner radius fixed to
$r_{\mathrm{LB}}=200$ pc. While fitting we tuned the thickness of the bubble
wall $\delta r_{\mathrm{LB}}$, the coordinates of the center of the bubble
$x_{\mathrm{LB}}$, $y_{\mathrm{LB}}$, $z_{\mathrm{LB}}$, and the direction
$l_{\mathrm{LB}}$, $b_{\mathrm{LB}}$ of the magnetic field before compression.
The field strength in the bubble wall is determined as a result of the
compression of magnetic lines of the initially uniform field by the expanding
spherically symmetric bubble. We assume that the field outside of the bubble
is unperturbed, the field inside is zero, and the field in the wall is
tangential to its surface and uniformly distributed across the wall in the
radial direction. We also assume that it is symmetric with respect to the axis
passing through the center of the bubble and parallel to the original
uncompressed field. The field configuration is schematically shown in Fig. 5.
Strictly speaking, such field is not divergence-free. However, by adjusting
the dependence of the field strength on the distance from the axis one can
make it divergence-free at scales larger than the thickness of the wall. The
conservation of the magnetic flux at these scales implies the following
relation:

$B_{\mathrm{LB}}(\theta)=B_{\mathrm{LB,0}}\left[1+\frac{r_{\mathrm{LB}}^{2}}{2r_{\mathrm{LB}}\delta
r_{\mathrm{LB}}+\delta r_{\mathrm{LB}}^{2}}\right]\sin\theta,$ (8)

where $\theta$ is the angle between the direction of the initial field and the
radius vector from the center of the bubble to the point in which the field is
calculated. The strength of the initial field $B_{\mathrm{LB,0}}$ is taken to
be the same as in the Local arm.

###  3.2 Thermal electrons

The calculation of the RMs requires the knowledge of the distribution of free
electrons in the Galaxy. In our study we adopt the simple plane-parallel model
of the density of free electrons $n_{e}$:

$n_{e}(z)=n_{0}\,e^{-|z|/z_{0}}$ (9)

where ${n_{0}=0.015}$ cm-3 is the mid-plane density and ${z_{0}=1.57}$ kpc is
the scale height. This corresponds to the best-fit model of Ocker et al.
(2020) where it was shown to correctly predict the dispersion measures (DM) of
pulsars located above the Galactic plane with ${|b|>20^{\circ}}$.

Similar values of $n_{0}$ and $z_{0}$ are used in more sophisticated models
such as NE2001 (Cordes & Lazio, 2002) and YMW16 (Yao et al., 2017).
Additionally, these models take into account the increase of the density of
electrons in the spiral arms whose effect is most pronounced near the Galactic
plane $|b|<10^{\circ}$. Since in our analysis we treat the positions of
magnetic arms as free parameters, the fixed positions of spiral arms of free
electrons may significantly affect the results. For this reason we exclude the
Galactic plane $|b|<10^{\circ}$ from the analysis of RMs, as has already been
stated in Sect. 2.2. At higher Galactic latitudes the predictions of NE2001,
YMW16 and plane-parallel models are approximately the same (Ocker et al.,
2020) and thus the use of the latter model does not lead to a loss of
accuracy.

In addition, we take into account the modifications of the thermal electron
density caused by the Local Bubble. On the wall of the Local Bubble the
density is higher and depends on its thickness $\delta r_{\mathrm{LB}}$ as

$n_{\mathrm{LB}}=n_{0}\,\frac{(r_{\mathrm{LB}}+\delta
r_{\mathrm{LB}})^{3}}{(r_{\mathrm{LB}}+\delta
r_{\mathrm{LB}})^{3}-r_{\mathrm{LB}}^{3}}.$ (10)

This expression follows from the assumption that the explosions which formed
the Local Bubble swept the ambient electrons of original constant density
$n_{0}$ to its walls. The value of the electron density inside the bubble is
not important because there we assume the magnetic field to be zero.

###  3.3 Cosmic ray electrons

Cosmic ray electrons (CRE)333Throughout the paper by CRE we mean both
relativistic electrons and positrons. are responsible for the Galactic
synchrotron emission and are believed to be mostly accelerated by Galactic
supernovae. Their propagation in the Galaxy is of a diffusive nature and can
be described by combined diffusion in the regular and turbulent magnetic
fields.

In a simplified picture, CRE propagation can be describe as an isotropic
diffusion determined by the single diffusion coefficient $D_{0}$. We modeled
the distribution of the CRE in the Galaxy using the numerical code DRAGON
(Evoli et al., 2017), taking into account CRE energy losses by the inverse
Compton scattering and the synchrotron emission. As a template for synchrotron
losses we used the modified version of Jansson & Farrar (2012b) turbulent
magnetic field with the normalization of all arms set to 5 $\mu$G.

An isotropic diffusion is assumed to be within the simulation volume of radius
${R_{\mathrm{CRE}}=16}$ kpc and the halo height of ${H=4}$ kpc. We neglected
anisotropic diffusion and assumed also that the diffusion coefficient is
spatially independent. The dependence on the CRE rigidity $R$ was taken to be
the following:

$D(R)=D_{0}\,\left(\frac{R}{R_{0}}\right)^{\delta},$ (11)

where $D_{0}=3.6\times 10^{28}$ cm2/s is the diffusion coefficient at rigidity
$R_{0}=4$ GV and $\delta=0.47$ is the spectral index. The propagation
parameters $H$, $D_{0}$ and $\delta$ were chosen based on the analysis of
secondary-to-primary ratios by De la Torre Luque et al. (2024). Similar values
were found in other studies, see Korsmeier & Cuoco (2021). The sources of
cosmic rays were assumed to be distributed in the Galaxy according to Lorimer
et al. (2006). The dynamic range of the simulation goes from 1 GeV up to 1
TeV. Our resulting CRE distribution is the sum of the distribution of primary
electrons and secondary electrons and positrons produced in interactions of
cosmic ray protons. Note that the contribution of secondary leptons to the
total CRE density does not exceed $10\%$ at energies of interest. Finally, we
normalize the distribution to the local measurements by the Alpha Magnetic
Spectrometer (AMS-02) experiment (Aguilar et al., 2019), as shown in Fig. 11.

Modelling the Fan Region required an additional CRE component in the Perseus
arm compared to the DRAGON output. Denoting as $L$ the length of the arm axis
starting from its origin we increase the CRE density in the Perseus arm
according to the bump-like profile:

$n(L)=\frac{f\,n_{0}}{(1+e^{-k(L-L_{1})})(1+e^{-k(L_{2}-L)})}$ (12)

where $n_{0}$ is the CRE density on the Solar vicinity, $f=0.7$ is the
relative amplitude of the additional bump, the parameter ${k=2}$ kpc controls
the sharpness of the bump boundaries, and $L_{1}$ and $L_{2}$ determine their
coordinates. As a result, the CRE density is enhanced roughly in the direction
from $l\simeq 110^{\circ}$ to $l\simeq 180^{\circ}$. The parameters $L_{1,2}$
were fixed from the geometry of the Fan Region. The amplitude of the
enhancement was also fixed; its value in the fit is largely degenerate with
the amplitude of the magnetic field.

To the contrary, we do not consider the modification of the CRE density on the
wall of the Local Bubble. Studying propagation of cosmic rays with full
account of the details of the local magnetic field constitutes a separate
problem and is beyond the scope of this paper.

Ideally, the CRE propagation should be fitted together with the GMF and
include the effects of anisotropic diffusion. In practice, diffusion of CRE in
space-dependent realistic GMF is a difficult problem which has not been solved
so far Giacinti et al. (2014); Kachelriess & Semikoz (2019). In this sense the
additional CRE component required for the Fan Region in our model may indicate
the necessity for deviation from the simplified isotropic picture.

Given the space and energy distributions of CRE in the Galaxy we calculate
synchrotron Stokes $Q$ and $U$ parameters starting from the synchrotron
emission of a single electron (Landau & Lifschits, 1975). We thus do not adopt
the commonly used power-law approximation of the electron spectrum for the
volume emissivity and instead integrate numerically over the CRE spectrum.
This allow us to take into account the curvature of the CRE spectrum as well
as it’s variation over the line of sight.

##  4 Results of the global fit

Figure 6: RM (top row) and polarized synchrotron skymaps for Stokes Q (middle
row) and Stokes U (bottom row) produced with the best-fit model (left column)
in comparison with the data (right column). Masks on the data maps are
discussed in the Sect. 2.2.

The fitting of the model is performed by minimizing the $\chi^{2}$. We first
adjusted individual GMF components, as well as their various combinations
until we reproduced qualitatively the key features of the data, and only then
performed the combined fit. The model is fitted to the binned data as
described in Sect. 2.1 with the binning scheme shown in Fig. 2. We used only
those bins, which are left after removing outliers and masking. All the data
bins and corresponding errorbars used for fitting are shown with black points
in Figs. 12,13, and 14. For each data bin we evaluated the model up to the
distance of 30 kpc from the position of the Sun in the direction of the line
of sight (LOS) passing through the center of the bin. We checked that if we
take 10 evenly distributed LOSes per bin, the resulting $\chi^{2}$ changes
very slightly, by 1-2%. Therefore, for the sake of saving computational time
during fitting, we used only the central LOS.

The results of the global fit are presented in Fig. 6 which shows the
comparison between the model (left column) and the masked data (right column)
for three quantities: the RMs and the $Q$ and $U$ polarization parameters. One
may observe an overall qualitative agreement between the model and data maps.
The best-fit total $\chi^{2}/\text{ndf}$ is 1.37. The breakup of $\chi^{2}$
for the corresponding three contributions is given in Table 4, first three
columns. The last two columns show what these quantities would be if the
errors were calculate as mere variances in the bins, after removal of
outliers. The over-estimation of errors in that case is obvious. The bin-by-
bin comparison between the model and the data is given in the Appendix C.

Table 1: The result of the global fit.

| $\chi^{2}$ |  $\chi^{2}$/ndf | ndf | $\chi_{\mathrm{var}}^{2}$ |  $\chi_{\mathrm{var}}^{2}$/ndf
---|---|---|---|---|---
RM | 537 | 1.90 | 283 | 145 | 0.51
Q | 397 | 1.14 | 348 | 238 | 0.68
U | 492 | 1.41 | 348 | 251 | 0.72
total | 1426 | 1.37 | 1037 | 634 | 0.61
444First three columns show the performance of the best-fit model: the overall
$\chi^{2}$ and $\chi^{2}$ per degree of freedom, as well as their breakup
between RM and polarization Q and U maps. The last two columns show the same
quantities when the standard deviation within the bins (calculated after
removal of outliers) is taken as an error estimate.

The best-fit parameters of the model are summarized in Table 5. The strength
of the magnetic field in the Local Arm was found to be 3.5 $\mu$G which is
higher than previously thought. The increased strength of the field is a
consequence of making the shape of the cross section of the Local Arm a free
parameter (the squircle parameter $n$). In previous models a constant height
of the disk across the arms was assumed. To the contrary, in our model the
best-fit value is $n=1.45$ causing the Local magnetic arm to be thicker in the
center than above the Sun, see Fig. 4. This shape of the cross section
significantly improves the fit allowing for stronger field in the direction of
the outer Galaxy. As a result, the model better describes the wave-like
pattern in the RM maps of the outer Galaxy between $l\sim 90^{\circ}$ and
$l\sim 270^{\circ}$, Fig. 12.

In the neighboring Sagittarius-Carina arm in the direction of the Galactic
center the magnetic field was found to be about 1.3 $\mu$G, which is weaker
than the field in the Local Arm. In addition, its direction is opposite to
that of the field in the Local Arm. Reversing the field direction in the
Sagittarius-Carina arm (so that it aligns with the direction of the field in
the Local Arm) while keeping the same strength worsens the RM $\chi^{2}$ by
$\Delta\chi^{2}\approx 110$. Thus, we confirm the field reversal in accordance
with Han & Qiao (1994) and Pshirkov et al. (2011).

Similarly, we confirm the necessity of the X-field. Setting it to zero
increases the total chi-square by $\Delta\chi^{2}\approx 260$. The best-fit
strength of the X-field near the Galactic plane was found to be 2 $\mu$G which
is twice as large as in the base model of Unger & Farrar (2023). Our model
prefers a sharp drop of the X-field strength at the Galactic radius of about 6
kpc. We discuss below in Sect. 5.2 how much it costs in terms of $\chi^{2}$ to
extend the X-type field to the solar vicinity as suggested by the arguments
based on the escape of Galactic cosmic rays.

The toroidal halo field is one of the key components of the model. If one sets
the field in the toroidal halo to zero, the total chi-square nearly doubles:
$\Delta\chi\approx 1010$. In the best-fit configuration the radial extensions
of the norther and southern tori were found to be different, with the southern
torus being larger than the northern one. However, the significance of this
difference is not strong. When both tori are forced to have the same radius,
the best-fit value of this common radius is $r_{\mathrm{max}}=13.5$ kpc, which
worsens the fit by only $\Delta\chi^{2}\approx 20$.

Inclusion of the Local Bubble is a distinctive and novel feature of our model.
As explained in Sect. 5.1, we approximated it as a spherical shell of finite
thickness with a compressed magnetic field on its wall. Notably, the fitted
parameters of the Local Bubble were found to be close to those known from
other observations. The best-fit thickness of the Local Bubble wall in our
model is $\delta r_{\mathrm{LB}}=30$ pc which is in perfect agreement with the
typical thickness of $35$ pc determined in the recent analysis of O’Neill et
al. (2024) based on the distribution of the local dust. Given that the radius
of the bubble was fixed to $200$ pc the field strength in the most compressed
regions of the wall reaches $\sim 10$–$14$ $\mu$G in accordance with Andersson
& Potter (2006), Medan & Andersson (2019) and Kalberla & Haud (2023). To be
more precise, using the Chandrasekhar-Fermi method Andersson & Potter (2006)
found the plane-of-the-sky strength of the magnetic field
$B_{\perp}=8^{+5}_{-3}$ $\mu$G in the wall of the Local Bubble toward
$(l,b)\approx(300^{\circ},0^{\circ})$. This is in good agreement with the
$B_{\perp}=11$ $\mu$G predicted by our model in the same direction.

The best-fit direction of the bubble field before compression in our model is
$(l_{\mathrm{LB}}$, $b_{\mathrm{LB}})=(230^{\circ}$, $-2^{\circ})$, which is
nearly opposite to the field of the Local Arm which has
$(l,b)=(70^{\circ},0^{\circ})$. Remarkably, Pelgrims et al. (2020) found that
the direction of the pre-compressed bubble field lies in the range
$(l,b)\approx(72^{\circ}\pm 1^{\circ},15^{\circ}\pm 2^{\circ})$ by analysing
the polarized emission of the dust in the Galactic polar caps with
$|b|>60^{\circ}$. In view of the sign ambiguity inherent in the synchrotron
polarization data, this is in perfect agreement with our fit. In fact, in our
fit the reversed (as compared to the Local Arm) direction is preferred only by
the RM data. It is worth noting, however, that the Local Bubble contribution
to the RM strongly depends on the shape of the bubble wall, so the problem of
the opposite field orientation should be re-examined together with a more
accurate modelling of the bubble.

The best-fit position of the center of the Local Bubble was found to be
shifted by $\sim 100$ pc from the position of the Solar System mostly in the
direction of the positive $y$-axis. This is also generally consistent with the
observations of the dust (Pelgrims et al., 2020; O’Neill et al., 2024).
Overall, the contribution from the Local Bubble is strongly favored by the
data. Its removal worsens the fit by $\Delta\chi^{2}\approx 660$.

The shape and overall brightness of the Fan Region in our model generally
matches the observations, as one can see from Fig. 6. The main contributions
come from the Local Arm and the Perseus Arm. The parameters of the former are
strongly constrained by the RM data. Adjusting mainly the parameters of the
Perseus Arm allowed us to fit the Fan Region. The correct vertical extension
of the Fan Region was achieved by fitting the (outer) vertical height of the
Perseus Arm which converged to $r_{z}=1.2$ kpc. The field strength in the
inner part around the Galactic plane was found to be $5.5$ $\mu$G in order to
reproduce the brightest part of the Fan Region at $|b|<10^{\circ}$. In the
outer part of the Perseus Arm the best-fit field strength is about
$3.5\,\mu$G.

The common pitch angle of the spiral magnetic field arms was an independent
parameter of the fit and converged to the value around
$\alpha_{0}=20^{\circ}$. As one can see from Fig. 3 this is close to the pitch
angle of spiral arms segments as inferred from the Gaia data (Poggio et al.,
2021).

Table 2: List of the parameters of the model.

Component | Parameter | Value | Units | Fit
---|---|---|---|---
Local Arm | $B$ | -3.5 |  $\mu$G | free
| $b$ | -2.2 | deg | free
| $n$ | 1.45 |  | free
| $r_{\mathrm{z}}$ | 0.73 | kpc | free
| $r_{\mathrm{disk}}$ | 1.0 | kpc |
| $x_{0}$ | -0.15 | kpc |
| $y_{0}$ | 0.0 | kpc |
Sagittarius- | $B$ | 1.3 |  $\mu$G | free
Carina Arm | $b$ | -80.0 | deg | free
| $n$ | 2.3 |  | free
| $r_{\mathrm{z}}$ | 1.0 | kpc | free
| $r_{\mathrm{disk}}$ | 0.8 | kpc |
| $x_{0}$ | 1.37 | kpc |
| $y_{0}$ | 0.0 | kpc |
Perseus Arm | $B$ | -3.5 |  $\mu$G | free
| $B^{\prime}$ | -2.0 |  $\mu$G | free
| $b$ | 46.0 | deg |
| $n$ | 2.0 |  |
| $r_{\mathrm{z}}$ | 1.2 | kpc | free
| $r_{\mathrm{z}}^{\prime}$ | 0.4 | kpc |
| $r_{\mathrm{disk}}$ | 1.1 | kpc |
| $x_{0}$ | -1.0 | kpc |
| $y_{0}$ | 0.0 | kpc |
Scutum Arm | $B$ | 4.9 |  $\mu$G | free
| $b$ | -134.0 | deg |
| $n$ | 2.0 |  |
| $r_{\mathrm{z}}$ | 0.8 | kpc |
| $r_{\mathrm{disk}}$ | 1.0 | kpc |
| $x_{0}$ | 1.0 | kpc |
| $y_{0}$ | 0.0 | kpc |
common pitch | $\alpha_{0}$ | 20.0 | deg | free
Northern | $B$ | 3.2 |  $\mu$G | free
Torus | $z_{\mathrm{min}}$ | 1.18 | kpc | free
| $z_{\mathrm{max}}$ | 1.9 | kpc | free
| $r_{\mathrm{max}}$ | 9.0 | kpc | free
Southern | $B$ | -3.5 |  $\mu$G | free
Torus | $z_{\mathrm{min}}$ | -2.5 | kpc | free
| $z_{\mathrm{max}}$ | -1.22 | kpc | free
| $r_{\mathrm{max}}$ | 14.0 | kpc | free
X-field | $B$ | 1.9 |  $\mu$G | free
| $r_{\mathrm{max}}$ | 6.2 | kpc | free
| $\theta$ | 28.0 | deg | free
Local | $x_{\mathrm{LB}}$ | -8.2 | kpc |
Bubble | $y_{\mathrm{LB}}$ | 0.095 | kpc | free
| $z_{\mathrm{LB}}$ | -0.05 | kpc | free
| $r_{\mathrm{LB}}$ | 0.2 | kpc |
| $\delta r_{\mathrm{LB}}$ | 0.03 | kpc | free
| $l_{\mathrm{LB}}$ | 230.0 | deg | free
| $b_{\mathrm{LB}}$ | -2.0 | deg | free
555See Sect. 3 for the description of each parameter.

##  5 Discussion

###  5.1 Local Bubble and nearby loops

The Local Bubble introduced in our model significantly improves the fit for
the polarization data and, to a lesser extent, for the RMs. Below we discuss
in more detail its role in these fits.

On the RM skymap the effect of the Local Bubble is clearly visible around the
northern polar cap. In this region between $0^{\circ}\lesssim l\lesssim
180^{\circ}$ the Local Bubble gives a positive contribution to the RM which
appears as a nearly uniform orange blur in Fig. 6. In the binned
representation one can see this area as three data bands with
$50^{\circ}<b<80^{\circ}$ with mostly positive RMs, see Fig. 12. This is the
region of the RM sky that forces the field of the bubble to take direction
opposite of the Local Arm. We should stress, however, that the sign and the
amplitude of the bubble contribution to the RM depends significantly on the
shape of its wall and the relative position of the solar system within the
bubble. For instance, in our model where the bubble is spherical, this
contribution can be set to zero by placing the solar system in the center as
the field on the wall, which is purely tangential, is then perpendicular to
the line of sight for all directions. It is clear, therefore, that a more
detailed modelling of the bubble shape is required to firmly determine whether
the direction of the field is parallel or anti-parallel to the field in the
Local Arm.

The sign flip of the initial bubble field has no effect on the polarization
maps. For both initial directions the Local Bubble produces the correct shape
of the signal in $Q$ and $U$ maps near the Galactic poles, provided the bubble
field is aligned to that in the Local Arm. The bubble contribution then almost
completely compensates the deficit of the polarized intensity resulting from
the other components when they are normalized to the RM data. It thus works
similarly or even better than the ’striation’ of the magnetic fields assumed
in previous works. We think, therefore, that the alignment of the fields in
the bubble and the Local Arm is a more robust result of our fit than the
flipped direction of the bubble field which may be an artifact of our over-
simplified modeling of the bubble geometry.

Figure 7: Polarized intensity of the halo and halo+Local Bubble above given
latitude $b$ normalized to the polarized intensity of the full model in the
same region of the sky. In this plot ”Halo” means everything above $|z|>1$
kpc. See text for the discussion.  Figure 8: Profile of the Stokes U signal,
calculated with the full model and the contribution from the Local Bubble.

The role of the Local Bubble in fitting the polarization data is illustrated
in Fig. 7 which shows the polarized intensity $\text{PI}=\sqrt{Q^{2}+U^{2}}$
integrated over the polar caps above given Galactic latitude $b$, as a
function of $b$. The thick dashed line shows the fractional contribution of
the ”halo” part of the field defined here as everything at $|z|>1$ kpc,
normalized to the full model. The thick solid line shows the contribution of
the ”halo” and the Local Bubble, without other ”disk” fields at $|z|<1$ kpc.
One can see that these two contributions together practically saturate the PI
at high latitudes in roughly equal shares. For comparison, the same quantity
is plotted in the model of Unger & Farrar (2023) with (thin solid curve) and
without (thin dashed curve) the striation. The effects of the Local Bubble and
striation are quite similar. The contribution of the Local Bubble to the
Stokes U is explicitly shown in the Fig. 8.

The significance of the Local Bubble can also be understood with the simple
analytical estimate. The brightness of the PI is proportional to the square of
the magnetic field $B^{2}$, the CR electron density $n_{\mathrm{CRE}}(z)$ and
the distance $d$. At high Galactic latitudes in the model of Unger & Farrar
(2023) the main contribution to the PI comes from the the toroidal halo, which
starts at $z\sim 1$ kpc and has a width of $d\sim 1$ kpc. Thus
${\mathrm{PI}_{\mathrm{halo}}\propto(3\mu\mathrm{G})^{2}\,n_{\mathrm{CRE}}(1\,\mathrm{kpc})\cdot\mathrm{kpc}}$.
At the same time the contribution from the Local Bubble reads:
${\mathrm{PI}_{\mathrm{LB}}\sim(10\mu\mathrm{G})^{2}\,n_{\mathrm{CRE}}(0)\cdot
0.03\,\mathrm{kpc}}$. Noting that ${n_{\mathrm{CRE}}(1\,\mathrm{kpc})\approx
n_{\mathrm{CRE}}(0)/2}$ one comes to the conclusion that both contributions
are approximately equal:

$\mathrm{PI}_{\mathrm{halo}}\approx\mathrm{PI}_{\mathrm{LB}}.$ (13)

The Local Bubble also affects propagation of cosmic rays from nearby sources.
In particular, it helps to explain the cosmic ray flux at the knee by a
contribution from the Vela supernova remnant (Bouyahiaoui et al., 2019). Also,
interactions of cosmic rays with the gas on the walls of the Local bubble can
give a significant contribution to the astrophysical neutrino flux (Andersen
et al., 2018; Bouyahiaoui et al., 2020).

There are other large-scale anomalies in the data which are also, presumably,
due to supernova bubbles, but viewed from outside, such as Loop I, Loop II,
Loop III, etc. The fact that these anomalies are prominent on the polarization
sky suggests that the contribution of the Local Bubble is also sizeable.
Inversely, it is possible that these anomalies can also be modelled in a
similar way and included in the fit, rather than masked (Thomson et al.,
2021). We leave this interesting question for future work.

###  5.2 Local X-shape field and escape of cosmic rays

The ratio of secondary to primary fluxes of cosmic rays indicates that locally
most of the cosmic rays escape from the Galaxy before interaction with the
interstellar gas (for the review see Kachelriess & Semikoz (2019) and
references therein). In the simplified models of cosmic ray propagation the
escape is described by diffusion. However, for pure turbulent field the
required diffusion coefficient turns out to be too high, corresponding to
unrealistically small magnetic field of the order of $10^{-5}$ $\mu$G. In
Giacinti et al. (2018) it was shown that local $z$-component of the coherent
GMF with the strength of a fraction of $\mu$G can resolve this cosmic ray
escape problem. In the recent collection of GMF models of Unger & Farrar
(2023) the local $z$-component is present in most of the models (see Fig. 16
in Unger & Farrar (2023)).

In our best-fit model the X-field vanishes beyond the galactocentric radius of
6.2 kpc. Adding the local $z$-component beyond 6.2 kpc with the strength of
$0.2$ $\mu$G slightly worsens the fit. If the direction of this additional
component is the same as the direction of the X-field (i.e., towards positive
$z$) the increase in the chi-square is $\Delta\chi^{2}\approx 90$. For the
opposite direction the worsening is smaller, $\Delta\chi^{2}\approx 30$. In
both cases all the increase comes from the RM fit since this field is too weak
to contribute to the polarized intensity which is quadratic in the magnetic
field strength. On the other hand, the contribution to the RM is of the same
order as the effect of the Local Bubble discussed in Sect. 5.1. Thus, the
search for such a weak additional component can only be performed together or
after the accurate modeling of the Local Bubble.

###  5.3 General remarks

Our Galaxy is believed to be a typical representative of the class of spiral
galaxies. Therefore, it appears likely that the structure of its magnetic
field should be similar to that of other spirals. Below we discuss how well
our model fits into the bulk of knowledge about other galaxies.

Our fit requires the presence of the disk and halo fields in our Galaxy, which
were also found in external spiral galaxies observed so far. We found that the
pitch angle of the disk field is close to that of spiral gaseous arms. The
rough similarity of the magnetic and stellar pitch angles were also observed
in other spirals (Beck et al., 2019).

From the polarization observations of other galaxies one finds that the
strongest field of the disk often resides in interarm regions, while the
regular field in the gaseous arms is somewhat suppressed. As it is clear from
Fig. 4, if our Galaxy was observed from outside, its spacial distribution of
the polarized intensity would be the opposite: the brightest parts would
rather coincide with the gaseous arms. In order to see whether this is an
unavoidable feature of our model we replace the constant field in the arms by
hollow magnetic profiles such that the magnetic field in the inner 70% of the
arm section is zero. Notably, this change of the profile does not spoil the RM
fit while shifting the strongest polarized emission regions to the interarm
space. The polarization skymaps also remain practically unchanged at latitudes
above $\sim 30^{\circ}$ where the polarized emission is mostly produced by the
Local Bubble and the magnetic tori. However, one should adjust the strength of
the field in the arms to have the correct emission near the Galactic plane. We
do not see any fundamental contradictions that would prevent this modification
of the model. We leave the detailed analysis of this possibility for future
work.

Another feature that is observed in some other galaxies is field reversals. As
one can see from Fig. 3, our model has one field reversal in the disk. The
galaxies with one reversal in the disk have also been observed, see Beck et
al. (2019).

##  6 Conclusions

In this paper we presented the new model of the regular magnetic field in the
Galactic halo, by which we understand the field everywhere except the thin
disk. The model makes use of the most recent catalog of the Faraday rotation
measures and the synchrotron polarization data by WMAP at 23 GHz. It builds on
the previously existing models but, apart from using the most recent data, has
several novel features as compared to previous works.

For the first time we took into account the contribution of the magnetic field
in the wall of the Local Bubble. By constructing an explicit, albeit rough,
model of the Local Bubble as a sphere with the compressed magnetic field
concentrated in the thick wall, we have shown that its contribution to the
polarized intensity can explain the observations at high Galactic latitudes
without contradiction with the RM data. This is in contrast with the previous
works where the deficit of model polarized intensity in these regions had to
be compensated by assuming the ”striation” of the field. The best-fit
parameters of the bubble were found to be in good qualitative agreement with
the existing independent measurements.

For the first time we have fitted the Fan Region — the large bright region in
the outer galaxy which was masked out in the previous works that were using
the polarization data. There recently appeared arguments that a large part of
the polarized emission in this region comes from distances of $1-2$ kpc and
thus cannot be explained by a local anomaly. In order to fit the overall
brightness and morphology of the polarized emission we had to assume that the
magnetic field in the Perseus arm varies along the arm cross section, being
stronger in the inner part close to the disk.

With these new features included, the best-fit value of the pitch angle of the
magnetic arms comes out close to $\sim 20^{\circ}$ which is about 2 times
larger than found in previous works. This value, however, is in good agreement
with the recent measurements by Gaia of the pitch angle of stellar arms
neighbouring the solar system.

Finally, we have developed a systematic method of error estimation which gives
control over the relative contribution of the data of different nature (the RM
and polarization data in our case) into the global fit. This method takes into
account the errors of individual measurements, the statistics of measurements
in the bins as well as the contribution of fluctuations at length scales
exceeding the bin size. We found that the errors estimated by this method are
typically smaller than were used in previous works, thus increasing the
constraining power of the data.

###### Acknowledgements.

We would like to thank Ioana Codrina Maris for illuminating discussions. Our
special thanks to Michael Unger, Maxim Pshirkov and Vincent Pelgrims for
thorough reading of the manuscript and numerous constructive comments. The
work of AK and PT is supported by the IISN project No. 4.4501.18.

## References

  * Aguilar et al. (2019) Aguilar, M., Ali Cavasonza, L., Alpat, B., et al. 2019, Phys. Rev. Lett., 122, 101101
  * Alves et al. (2018) Alves, M. I. R., Boulanger, F., Ferrière, K., & Montier, L. 2018, A&A, 611, L5
  * Andersen et al. (2018) Andersen, K. J., Kachelrieß, M., & Semikoz, D. V. 2018, Astrophys. J. Lett., 861, L19
  * Anderson et al. (2015) Anderson, C. S., Gaensler, B. M., Feain, I. J., & Franzen, T. M. O. 2015, ApJ, 815, 49
  * Andersson & Potter (2006) Andersson, B. G. & Potter, S. B. 2006, ApJ, 640, L51
  * Battye et al. (2011) Battye, R. A., Browne, I. W. A., Peel, M. W., Jackson, N. J., & Dickinson, C. 2011, MNRAS, 413, 132
  * Beck (2001) Beck, R. 2001, Space Sci. Rev., 99, 243
  * Beck (2015) Beck, R. 2015, A&A Rev., 24, 4
  * Beck et al. (2019) Beck, R., Chamandy, L., Elson, E., & Blackman, E. G. 2019, Galaxies, 8, 4
  * Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20
  * Berkhuijsen et al. (2016) Berkhuijsen, E. M., Urbanik, M., Beck, R., & Han, J. L. 2016, A&A, 588, A114
  * Betti et al. (2019) Betti, S. K., Hill, A. S., Mao, S. A., et al. 2019, ApJ, 871, 215
  * Bouyahiaoui et al. (2019) Bouyahiaoui, M., Kachelriess, M., & Semikoz, D. V. 2019, JCAP, 01, 046
  * Bouyahiaoui et al. (2020) Bouyahiaoui, M., Kachelrieß, M., & Semikoz, D. V. 2020, Phys. Rev. D, 101, 123023
  * Brandenburg & Ntormousi (2023) Brandenburg, A. & Ntormousi, E. 2023, ARA&A, 61, 561
  * Broten et al. (1988) Broten, N. W., MacLeod, J. M., & Vallee, J. P. 1988, Ap&SS, 141, 303
  * Brown et al. (2007) Brown, J. C., Haverkorn, M., Gaensler, B. M., et al. 2007, ApJ, 663, 258
  * Brown et al. (2003) Brown, J. C., Taylor, A. R., & Jackel, B. J. 2003, ApJS, 145, 213
  * Burn (1966) Burn, B. J. 1966, MNRAS, 133, 67
  * Clarke et al. (2001) Clarke, T. E., Kronberg, P. P., & Böhringer, H. 2001, ApJ, 547, L111
  * Clegg et al. (1992) Clegg, A. W., Cordes, J. M., Simonetti, J. M., & Kulkarni, S. R. 1992, ApJ, 386, 143
  * Cordes & Lazio (2002) Cordes, J. M. & Lazio, T. J. W. 2002, arXiv e-prints, astro
  * Costa & Spangler (2018) Costa, A. H. & Spangler, S. R. 2018, ApJ, 865, 65
  * Costa et al. (2016) Costa, A. H., Spangler, S. R., Sink, J. R., Brown, S., & Mao, S. A. 2016, ApJ, 821, 92
  * De la Torre Luque et al. (2024) De la Torre Luque, P., Winkler, M. W., & Linden, T. 2024, JCAP, 05, 104
  * Evoli et al. (2017) Evoli, C., Gaggero, D., Vittino, A., et al. 2017, JCAP, 02, 015
  * Farnes et al. (2014) Farnes, J. S., Gaensler, B. M., & Carretti, E. 2014, ApJS, 212, 15
  * Feain et al. (2009) Feain, I. J., Ekers, R. D., Murphy, T., et al. 2009, ApJ, 707, 114
  * Ferrière & Terral (2014) Ferrière, K. & Terral, P. 2014, A&A, 561, A100
  * Finkbeiner (2003) Finkbeiner, D. P. 2003, Astrophys. J. Suppl., 146, 407
  * Fletcher et al. (2011) Fletcher, A., Beck, R., Shukurov, A., Berkhuijsen, E. M., & Horellou, C. 2011, MNRAS, 412, 2396
  * Gaensler et al. (2001) Gaensler, B. M., Dickey, J. M., McClure-Griffiths, N. M., et al. 2001, ApJ, 549, 959
  * Giacinti et al. (2014) Giacinti, G., Kachelrieß, M., & Semikoz, D. V. 2014, Phys. Rev. D, 90, 041302
  * Giacinti et al. (2018) Giacinti, G., Kachelriess, M., & Semikoz, D. V. 2018, JCAP, 07, 051
  * Han et al. (1997) Han, J. L., Manchester, R. N., Berkhuijsen, E. M., & Beck, R. 1997, A&A, 322, 98
  * Han et al. (1999) Han, J. L., Manchester, R. N., & Qiao, G. J. 1999, Mon. Not. Roy. Astron. Soc., 306, 371
  * Han et al. (2018) Han, J. L., Manchester, R. N., van Straten, W., & Demorest, P. 2018, ApJS, 234, 11
  * Han & Qiao (1994) Han, J. L. & Qiao, G. J. 1994, A&A, 288, 759
  * Heald et al. (2009) Heald, G., Braun, R., & Edmonds, R. 2009, A&A, 503, 409
  * Hill et al. (2017) Hill, A. S., Landecker, T. L., Carretti, E., et al. 2017, MNRAS, 467, 4631
  * Howard & Kulsrud (1997) Howard, A. M. & Kulsrud, R. M. 1997, Astrophys. J., 483, 648
  * Hutschenreuter et al. (2022) Hutschenreuter, S., Anderson, C. S., Betti, S., et al. 2022, A&A, 657, A43
  * Jansson & Farrar (2012a) Jansson, R. & Farrar, G. R. 2012a, Astrophys. J., 757, 14
  * Jansson & Farrar (2012b) Jansson, R. & Farrar, G. R. 2012b, Astrophys. J. Lett., 761, L11
  * Jansson et al. (2009) Jansson, R., Farrar, G. R., Waelkens, A. H., & Enßlin, T. A. 2009, J. Cosmology Astropart. Phys., 2009, 021
  * Kachelriess & Semikoz (2019) Kachelriess, M. & Semikoz, D. V. 2019, Prog. Part. Nucl. Phys., 109, 103710
  * Kaczmarek et al. (2017) Kaczmarek, J. F., Purcell, C. R., Gaensler, B. M., McClure-Griffiths, N. M., & Stevens, J. 2017, MNRAS, 467, 1776
  * Kalberla & Haud (2023) Kalberla, P. M. W. & Haud, U. 2023, A&A, 673, A101
  * Kim et al. (2016) Kim, K. S., Lilly, S. J., Miniati, F., et al. 2016, ApJ, 829, 133
  * Klein et al. (2003) Klein, U., Mack, K. H., Gregorini, L., & Vigotti, M. 2003, A&A, 406, 579
  * Korsmeier & Cuoco (2021) Korsmeier, M. & Cuoco, A. 2021, Phys. Rev. D, 103, 103016
  * Landau & Lifschits (1975) Landau, L. D. & Lifschits, E. M. 1975, Course of Theoretical Physics, Vol. Volume 2, The Classical Theory of Fields (Oxford: Pergamon Press)
  * Law et al. (2011) Law, C. J., Gaensler, B. M., Bower, G. C., et al. 2011, ApJ, 728, 57
  * Livingston et al. (2021) Livingston, J. D., McClure-Griffiths, N. M., Gaensler, B. M., Seta, A., & Alger, M. J. 2021, MNRAS, 502, 3814
  * Livingston et al. (2022) Livingston, J. D., McClure-Griffiths, N. M., Mao, S. A., et al. 2022, MNRAS, 510, 260
  * Lorimer et al. (2006) Lorimer, D. R. et al. 2006, Mon. Not. Roy. Astron. Soc., 372, 777
  * Ma et al. (2020) Ma, Y. K., Mao, S. A., Ordog, A., & Brown, J. C. 2020, MNRAS, 497, 3097
  * Ma et al. (2019) Ma, Y. K., Mao, S. A., Stil, J., et al. 2019, MNRAS, 487, 3432
  * Mao et al. (2010) Mao, S. A., Gaensler, B. M., Haverkorn, M., et al. 2010, ApJ, 714, 1170
  * Mao et al. (2008) Mao, S. A., Gaensler, B. M., Stanimirović, S., et al. 2008, ApJ, 688, 1029
  * Mao et al. (2012a) Mao, S. A., McClure-Griffiths, N. M., Gaensler, B. M., et al. 2012a, ApJ, 755, 21
  * Mao et al. (2012b) Mao, S. A., McClure-Griffiths, N. M., Gaensler, B. M., et al. 2012b, ApJ, 759, 25
  * Medan & Andersson (2019) Medan, I. & Andersson, B. G. 2019, ApJ, 873, 87
  * Minter & Spangler (1996) Minter, A. H. & Spangler, S. R. 1996, ApJ, 458, 194
  * Mora & Krause (2013) Mora, S. C. & Krause, M. 2013, A&A, 560, A42
  * Moss & Sokoloff (2008) Moss, D. & Sokoloff, D. 2008, A&A, 487, 197
  * Ocker et al. (2020) Ocker, S. K., Cordes, J. M., & Chatterjee, S. 2020, Astrophys. J., 897, 124
  * O’Neill et al. (2024) O’Neill, T. J., Zucker, C., Goodman, A. A., & Edenhofer, G. 2024, arXiv e-prints, arXiv:2403.04961
  * Oren & Wolfe (1995) Oren, A. L. & Wolfe, A. M. 1995, ApJ, 445, 624
  * O’Sullivan et al. (2017) O’Sullivan, S. P., Purcell, C. R., Anderson, C. S., et al. 2017, MNRAS, 469, 4034
  * O’Sullivan et al. (2023) O’Sullivan, S. P., Shimwell, T. W., Hardcastle, M. J., et al. 2023, MNRAS, 519, 5723
  * Panopoulou et al. (2024) Panopoulou, G. V., Zucker, C., Clemens, D., et al. 2024, arXiv e-prints, arXiv:2406.03765
  * Pelgrims et al. (2020) Pelgrims, V., Ferrière, K., Boulanger, F., Lallement, R., & Montier, L. 2020, A&A, 636, A17
  * Planck Collaboration et al. (2016) Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016, A&A, 596, A103
  * Poggio et al. (2021) Poggio, E., Drimmel, R., Cantat-Gaudin, T., et al. 2021, A&A, 651, A104
  * Pshirkov et al. (2011) Pshirkov, M. S., Tinyakov, P. G., Kronberg, P. P., & Newton-McGee, K. J. 2011, Astrophys. J., 738, 192
  * Pshirkov et al. (2013) Pshirkov, M. S., Tinyakov, P. G., & Urban, F. R. 2013, Mon. Not. Roy. Astron. Soc., 436, 2326
  * Ranchod et al. (2024) Ranchod, S., Mao, S. A., Deane, R., et al. 2024, arXiv e-prints, arXiv:2403.13500
  * Riseley et al. (2020) Riseley, C. J., Galvin, T. J., Sobey, C., et al. 2020, PASA, 37, e029
  * Riseley et al. (2018) Riseley, C. J., Lenc, E., Van Eck, C. L., et al. 2018, PASA, 35, 43
  * Rossetti et al. (2008) Rossetti, A., Dallacasa, D., Fanti, C., Fanti, R., & Mack, K. H. 2008, A&A, 487, 865
  * Roy et al. (2005) Roy, S., Rao, A. P., & Subrahmanyan, R. 2005, MNRAS, 360, 1305
  * Schnitzeler (2010) Schnitzeler, D. H. F. M. 2010, Mon. Not. Roy. Astron. Soc., 409, L97
  * Schnitzeler et al. (2019a) Schnitzeler, D. H. F. M., Carretti, E., Wieringa, M. H., et al. 2019a, MNRAS, 485, 1293
  * Schnitzeler et al. (2019b) Schnitzeler, D. H. F. M., Carretti, E., Wieringa, M. H., et al. 2019b, MNRAS, 485, 1293
  * Shanahan et al. (2019) Shanahan, R., Lemmer, S. J., Stil, J. M., et al. 2019, ApJ, 887, L7
  * Simard-Normandin et al. (1981) Simard-Normandin, M., Kronberg, P. P., & Button, S. 1981, ApJS, 45, 97
  * Skalidis & Pelgrims (2019) Skalidis, R. & Pelgrims, V. 2019, A&A, 631, L11
  * Soida et al. (2011) Soida, M., Krause, M., Dettmar, R. J., & Urbanik, M. 2011, A&A, 531, A127
  * Sun & Reich (2010) Sun, X.-H. & Reich, W. 2010, Research in Astronomy and Astrophysics, 10, 1287
  * Tabara & Inoue (1980) Tabara, H. & Inoue, M. 1980, A&AS, 39, 379
  * Taylor et al. (2024) Taylor, A. R., Sekhar, S., Heino, L., et al. 2024, MNRAS, 528, 2511
  * Taylor et al. (2009) Taylor, A. R., Stil, J. M., & Sunstrum, C. 2009, The Astrophysical Journal, 702, 1230
  * Terral & Ferrière (2017) Terral, P. & Ferrière, K. 2017, A&A, 600, A29
  * Thomson et al. (2019) Thomson, A. J. M., Landecker, T. L., Dickey, J. M., et al. 2019, MNRAS, 487, 4751
  * Thomson et al. (2021) Thomson, A. J. M., Landecker, T. L., McClure-Griffiths, N. M., et al. 2021, MNRAS, 507, 3495
  * Tinyakov & Tkachev (2002) Tinyakov, P. G. & Tkachev, I. I. 2002, Astropart. Phys., 18, 165
  * Unger & Farrar (2023) Unger, M. & Farrar, G. 2023 [arXiv:2311.12120]
  * Van Eck et al. (2021) Van Eck, C. L., Brown, J. C., Ordog, A., et al. 2021, ApJS, 253, 48
  * Van Eck et al. (2011) Van Eck, C. L., Brown, J. C., Stil, J. M., et al. 2011, ApJ, 728, 97
  * Van Eck et al. (2023) Van Eck, C. L., Gaensler, B. M., Hutschenreuter, S., et al. 2023, ApJS, 267, 28
  * Van Eck et al. (2018) Van Eck, C. L., Haverkorn, M., Alves, M. I. R., et al. 2018, A&A, 613, A58
  * Vanderwoude et al. (2024) Vanderwoude, S., West, J. L., Gaensler, B. M., et al. 2024, arXiv e-prints, arXiv:2403.15668
  * Vernstrom et al. (2018) Vernstrom, T., Gaensler, B. M., Vacca, V., et al. 2018, MNRAS, 475, 1736
  * Watts et al. (2023) Watts, D. J. et al. 2023, Astron. Astrophys., 679, A143
  * West et al. (2021) West, J. L., Landecker, T. L., Gaensler, B. M., Jaffe, T., & Hill, A. S. 2021, ApJ, 923, 58
  * Wielebinski & Krause (1993) Wielebinski, R. & Krause, F. 1993, The Astronomy and Astrophysics Review, 4, 449 – 485
  * Wolleben et al. (2010) Wolleben, M., Fletcher, A., Landecker, T. L., et al. 2010, ApJ, 724, L48
  * Wolleben et al. (2019) Wolleben, M., Landecker, T. L., Carretti, E., et al. 2019, AJ, 158, 44
  * Wolleben et al. (2021) Wolleben, M., Landecker, T. L., Douglas, K. A., et al. 2021, AJ, 162, 35
  * Xu & Han (2019) Xu, J. & Han, J. L. 2019, MNRAS, 486, 4275
  * Xu & Han (2024) Xu, J. & Han, J. L. 2024, arXiv e-prints, arXiv:2404.02038
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29

##  Appendix A Removal of outliers

Each data bin that we consider (a rectangular region on the sphere limited by
lines of constant $l$ and $b$), we fit with a Gaussian. Namely, we assume that
the measurements contained within the bin are samples from a Gaussian
distribution with unknown parameters. By fitting data with a Gaussian we
determine the mean $\mu$ and width $\sigma$ of this distribution, taking into
account errors of individual measurements.

To do this we construct a special likelihood function $\mathcal{L}$. Denoting
as $\mu_{\mathrm{i}}$ and $\sigma_{\mathrm{i}}$ the mean and error of each
measurement in the bin it reads:

$\mathcal{L}=\prod_{i=1}^{N}\int\limits_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi\sigma^{2}}\sqrt{2\pi\sigma_{\mathrm{i}}^{2}}}\exp\left({-\frac{(\mu-x^{\prime})^{2}}{2\sigma^{2}}}{-\frac{(\mu_{\mathrm{i}}-x^{\prime})^{2}}{2\sigma_{\mathrm{i}}^{2}}}\right)\,\mathrm{d}x^{\prime}.$
(14)

where $N$ is the total number of measurements within the bin. Physically each
integral under the product means the probability of observing the value
$\mu_{\mathrm{i}}$ given the parameters of the distribution $\mu$ and $\sigma$
and an observational error $\sigma_{\mathrm{i}}$. This likelihood explicitly
relies on the assumption of Gaussianity of the errors and in the limit
$\sigma_{\mathrm{i}}\to 0$ reduces to a standard likelihood function.

Figure 9: Example of RM bin fitted with the gaussian. Points at the bottom of
the figure indicate individual measurements within the bin together with their
observational errors. Black points were discarded as outlier during the fit,
blue points mark the LOFAR LoTSS data and orange points mark all other data.
Blue boxes show the histogram of all sources except for outliers. Orange curve
is the gaussian obtained as a result of the fit to all data (including LOFAR)
and the pale green band is 68% confidence interval of its $\mu$. The
parameters of this gaussian are indicated on the plot. The blue curve is the
best-fit gaussian if only LOFAR LoTSS data is used.

Evaluating the integral we find

$\mathcal{L}=\prod_{i=1}^{N}\frac{1}{\sqrt{2\pi(\sigma^{2}+\sigma_{\mathrm{i}}^{2})}}\exp\left({-\frac{(\mu-\mu_{\mathrm{i}})^{2}}{2(\sigma^{2}+\sigma_{\mathrm{i}}^{2})}}\right).$
(15)

The best fit values of $\sigma$ and $\mu$ are obtained by minimizing the log
likelihood $\chi^{2}=-2\log\mathcal{L}$. When applied to the real data we use
an iterative procedure which also allow us to self-consistently remove
outliers. After fitting the bin we remove those measurements whose values
deviate from the mean $\mu$ by more than
$3\sqrt{\sigma^{2}+\sigma_{\mathrm{i}}^{2}}$ and then repeat the fitting. The
algorithm terminates when no outliers are detected. As a final step we
determine $1\sigma$ confidence intervals ${\mu\pm\sigma_{\mu}}$ and
${\sigma\pm\sigma_{\sigma}}$ by exploring the contour of
${\chi^{2}=\chi^{2}_{\mathrm{min}}+1}$.

An example of the bin is shown in Fig 9. The best fit values as well as its
uncertainties are shown in the legend. The resulting Gaussian curve is
narrower than the histogram since it shows how the distribution would look
like assuming zero measurement errors. In total the bin initially contained
245 sources, 8 of which were detected as outliers (only 6 are visible in the
plot). Since the algorithm takes into account errors of individual
measurements, it can be noticed that a source with an RM of about -30 rad/m2
and a large error is included in the fit, while a source with a smaller RM of
$\sim$ 25 rad/m2 and a much smaller error is marked as an outlier.

Figure 10: Fourier decomposition of the RM strip $40^{\circ}<b<50^{\circ}$.
The upper panel shows the RM signal binned into 360 bins by $l$. The middle
panel shows the first three Fourier modes (namely, the constant and
$360^{\circ}$ and $180^{\circ}$ wavelength modes) of the RM signal, while the
lower panel shows the contribution from all other modes except for the first
three.

##  Appendix B Derivation of the bin variance

Here we derive the variances of the bins associated with the coherent
fluctuations in the ISM. First we derive the formula and then discuss how it
is applied to the data.

Assume that the observable $f$ depends on the discrete positions $x^{i}$ and
denote $f^{i}=f(x^{i})$. The positions are evenly spaced so that
$x_{i}-x_{i-1}=\Delta x$ for any $i$. Assume also that $f$ can be represented
as a sum of the signal $f_{0}$ and a random noise $\epsilon$:

$f^{i}=f^{i}_{0}+\epsilon^{i}.$ (16)

Consider then the quantity $f_{\mathrm{L}}^{i}$ which is the value of $f$
averaged over $N$ nearby positions, or in other words averaged over the bin of
the size $L=N\,\Delta x$

$f_{\mathrm{L}}^{i}=\frac{1}{N}\sum\limits_{j=i}^{i+N-1}f^{j}.$ (17)

We are interested in its variance $\sigma_{\mathrm{L}}^{2}$:

$\sigma_{L}^{2}=\langle\left(f_{\mathrm{L}}^{i}-\langle
f_{\mathrm{L}}^{i}\rangle\right)^{2}\rangle.$ (18)

The angle brackets here denote averaging over different realizations of the
noise $\epsilon$ while keeping the same signal $f_{0}$. We further assume that
the statistical properties of the noise do not depend on the position $x^{i}$
and that the noise has zero mean $\langle\epsilon\rangle=0$. Substituting 16
and 17 into 18 one gets

$\sigma_{L}^{2}=\langle\left(\frac{1}{N}\sum\limits_{j=i}^{i+N-1}\epsilon^{j}\right)^{2}\rangle=\frac{1}{N^{2}}\sum\limits_{j=0}^{N-1}\sum\limits_{q=0}^{N-1}\langle\epsilon^{j}\epsilon^{q}\rangle.$
(19)

As we assumed that the noise does not depend on $x$, the correlator
$\langle\epsilon^{j}\epsilon^{q}\rangle$ is the function only of the
difference $|x_{j}-x_{q}|$ but not the position itself. This allows one to
evaluate the double sum in the Eq. 19 by collecting and counting equal terms.
Denoting the correlator
$\langle\epsilon^{j}\epsilon^{q}\rangle=r(|x_{j}-x_{q}|)$ one obtains the
following expression:

$\sigma_{L}^{2}=\frac{1}{N}\sum\limits_{j=-N+1}^{N-1}\left(1-\frac{|j|}{N}\right)\,r(|j\,\Delta
x|),$ (20)

where $r(x)$ is the correlation function of the noise. If the noise can be
extracted from the measurements, this formula allows one to directly calculate
the variance for the extended bin of the size $L$.

In order to clarify the physical meaning of 20 it is useful to consider the
continuous limit $\Delta x\to 0$ in which the sum in 20 can be replaced with
an integral:

$\sigma_{L}^{2}=\frac{2}{L}\int\limits_{0}^{L}\mathrm{d}x\left(1-\frac{x}{L}\right)\,r(x).$
(21)

Using Wiener-Khinchin theorem the correlation function $r(x)$ can be expressed
in terms of the noise power spectrum $S(k)$

$r(x)=\int\mathrm{d}k\,\,S(k)e^{ikx}=\int\mathrm{d}k\,\,S(k)\cos(kx).$ (22)

Substituting 22 into 21 and evaluating the integral over $x$ one obtains

$\sigma^{2}_{L}=\int\mathrm{d}k\,\,\operatorname{sinc}^{2}\left(\frac{kL}{2}\right)S(k).$
(23)

When applied to the real data this equation transforms into a sum over the
finite number of modes, where $S_{k}$ can be calculated for example with the
help of the Fast Fourier Transform (FFT). Thus, we arrive at the equation 24:

$\sigma^{2}_{\mathrm{L}}=2\sum_{k=k_{\mathrm{min}}}^{\infty}\,\operatorname{sinc}^{2}\left(\frac{kL}{2}\right)S_{k}.$
(24)

In practice we calculate the power spectrum $S_{k}$ by making a FFT of the
data strips of constant latitude binned into 360 bins by longitude $l$. The
FFT is made after removal of outliers, as described in Appendix A. An example
of such transformation is shown in Fig. 10 for the RM data strip. In order to
exclude the contribution of the regular field we take into account only the
modes with $k\geq k_{\mathrm{min}}$, where $k_{\mathrm{min}}=3$ corresponds to
the wavelength of $90^{\circ}$.

##  Appendix C Details of the final fit and comparison with existing models

In this section we present a comparison of our model with the base model of
Unger & Farrar (2023). In order to reproduce the results of Unger & Farrar
(2023) we switched to the thermal electron model of Yao et al. (2017),
recomputed cosmic ray electron distribution following the prescriptions given
in Unger & Farrar (2023) and included the effect of the striated field.

The bin by bin performances of the models are shown in Figs. 12 13 14. Our
model and the model of Unger & Farrar (2023) were fitted assuming different
masks on the data. Apart from small scale features, the main difference
between the masks comes from the Galactic plane and the Fan region. The
Galactic plane in the RM data is excluded in our analysis, while present in
the analysis of Unger & Farrar (2023). To the contrary, in the polarization
data the Galactic plane is masked by Unger & Farrar (2023) while included in
our fits. Additionally, Unger & Farrar (2023) masked the Fan region on the
Stokes $Q$ and $U$ maps. Another two bright radio features Loop I and Loop II
were excluded in both analysis in Stokes $Q$ and $U$.

Summarizing the above, our mask is larger for the RM dataset which can be seen
by the absence of data points in the lower panel of Fig. 12. On the other
hand, the mask on the polarization data is larger in Unger & Farrar (2023).
The additional regions masked out by Unger & Farrar (2023) are roughly shown
with grey shading in the Figs. 13 and 14. Note that the model by Unger &
Farrar (2023) does not fit the data in these regions.

As one can see from Figs. 12 13 14, in general the predictions of the models
are quite similar. However, the direct comparison of the models in terms of
$\chi^{2}$ is not straightforward since they were fitter assuming different
masks on the data. For cross-check purposes, the conservative approach would
be to compare the models only in those regions where both were fitted.
Reproducing roughly the mask of Unger & Farrar (2023) and taking into account
only those bins which were included in both analyses we find that the total
chi-squared/ndf for our model is $\chi^{2}=1070/717$, while for the base model
of Unger & Farrar (2023) one finds $\chi^{2}=1436/717$. Most of the difference
comes from the RM data.

Figure 11: Spectrum of relativistic leptons. Black dots show the measurements
by the AMS-02 experiment. The black curve corresponds to the CRE spectrum
calculated with the DRAGON code as described in Sec. 3.3 and used in the
analysis. The CRE within the energies indicated with the yellow band produce
more than 95% of the synchrotron emission at 23 GHz. The deviation of the
AMS-02 measurements below 10 GeV is due to the Solar modulation which we do
not take into account. Figure 12: Our best fit model (black line) and UF23
base model (red line) in comparison with the binned RM data (black points with
errorbars).  Figure 13: Same as Fig.12 but for the Stokes Q parameter. The
gray shading roughly indicates the regions which were excluded in the analysis
of Unger & Farrar (2023). Figure 14: Same as Fig.12 but for the Stokes U
parameter.