A Uncertainty estimation

# Constraints from Faraday rotation on the magnetic field structure in the Galactic halo

###### Key Words.:
ISM: magnetic fields – Galaxy: halo – galaxies: halos – galaxies: magnetic fields – galaxies: spiral

## Abstract

Context:

Aims:We examine the constraints imposed by Faraday rotation measures of extragalactic point sources on the structure of the magnetic field in the halo of our Galaxy. Guided by radio polarization observations of external spiral galaxies, we look in particular into the possibility that field lines in the Galactic halo have an X shape.

Methods:We employ the analytical models of spiraling, possibly X-shape magnetic fields derived in a previous paper to generate synthetic all-sky maps of the Galactic Faraday depth, which we fit to an observational reference map with the help of Markov Chain Monte Carlo simulations.

Results:We find that the magnetic field in the Galactic halo is slightly more likely to be bisymmetric (azimuthal wavenumber, ) than axisymmetric (). If it is indeed bisymmetric, it must appear as X-shaped in radio polarization maps of our Galaxy seen edge-on from outside, but if it is actually axisymmetric, it must instead appear as nearly parallel to the Galactic plane.

Conclusions:

## 1 Introduction

Interstellar magnetic fields are an important component of the interstellar medium (ISM) of galaxies. They play a crucial role in a variety of physical processes, including cosmic-ray acceleration and propagation, gas distribution and dynamics, star formation However, their properties remain poorly understood. The main difficulty is the lack of direct measurements, apart from Zeeman-splitting measurements in dense, neutral clouds. In addition, all existing observational methods provide only partial information, be it the field strength, the field direction/orientation, or the field component parallel or perpendicular to the line of sight.

In the case of our own Galaxy, magnetic field observations are particularly difficult to interpret, because the observed emission along any line of sight is generally produced by a number of structures, whose contributions are often hard to disentangle and locate along the line of sight. In contrast, observations of external galaxies can give us a bird’s eye view of the large-scale structure of their magnetic fields. High-resolution radio polarization observations of spiral galaxies have shown that face-on galaxies have spiral field lines, whereas edge-on galaxies have field lines that are parallel to the galactic plane in the disk (e.g., Wielebinski & Krause, 1993; Dumke et al., 1995) and inclined to the galactic plane in the halo, with an inclination angle increasing outward in the four quandrants; these halo fields have been referred to as X-shape magnetic fields (Tüllmann et al., 2000; Soida, 2005; Krause et al., 2006; Krause, 2009; Heesen et al., 2009; Braun et al., 2010; Soida et al., 2011; Haverkorn & Heesen, 2012).

In a recent paper, Ferrière & Terral (2014) (Paper 1) presented a general method to construct analytical models of divergence-free magnetic fields that possess field lines of a prescribed shape, and they used their method to obtain four models of spiraling, possibly X-shape magnetic fields in galactic halos as well as two models of spiraling, mainly horizontal (i.e., parallel to the galactic plane) magnetic fields in galactic disks. Their X-shape models were meant to be quite versatile and to have broad applicability; in particular, they were designed to span the whole range of field orientation, from purely horizontal to purely vertical.

Our next purpose is to resort to the galactic magnetic field models derived in Paper 1 to explore the structure of the magnetic field in the halo of our own Galaxy. The idea is to adjust the free parameters of the different field models such as to achieve the best possible fits to the existing observational data – which include mainly Faraday rotation measures (RMs) and synchrotron (total and polarized) intensities – and to determine how good the different fits are. Evidently, a good field model must be able to fit both Faraday-rotation and synchrotron data simultaneously. However, we feel it is important to first consider both types of data separately and examine independently the specific constraints imposed by each kind of observations. This is arguably the best way to understand either why a given model is ultimately acceptable or on what grounds it must be rejected.

In the present paper, we focus on fitting our galactic magnetic field models to the existing Faraday-rotation data. In practice, since we need to probe through the entire Galactic halo, we retain exclusively the RMs of extragalactic point sources (as opposed to Galactic pulsars). For each considered field model, we simulate an all-sky map of the Galactic Faraday depth (FD), which we confront to an observational reference map based on the reconstructed Galactic FD map of Oppermann et al. (2015). The fitting procedure relies on standard minimization, performed through a Markov Chain Monte Carlo (MCMC) analysis.

Our paper is organized as follows: In Sect. 2, we present the all-sky map of the Galactic FD that will serve as our observational reference. In Sect. 3, we review the magnetic field models derived in Paper 1 and employed in this study. In Sect. 4, we describe the fitting procedure. In Sect. 5, we present our results and compare them with previous halo-field models. In Sect. 6, we summarize our work and provide a few concluding remarks.

Two different coordinate systems are used in the paper. The all-sky maps are plotted in Galactic coordinates , with longitude increasing eastward (to the left) and latitude increasing northward (upward).1 In contrast, the magnetic field models are described in Galactocentric cylindrical coordinates, , with azimuthal angle increasing in the direction of Galactic rotation, i.e., clockwise about the -axis, from in the azimuthal plane through the Sun. As a result, the coordinate system is left-handed. For the Galactocentric cylindrical coordinates of the Sun, we adopt .

## 2 Observational map of the Galactic Faraday depth

### 2.1 Our observational reference map

Faraday rotation is the rotation of the polarization direction of a linearly-polarized radio wave that passes through a magneto-ionic medium. The polarization angle rotates by the angle , where is the observing wavelength and RM is the rotation measure given by

 RM=(0.81 rad m−2) ∫L0ne B∥ ds , (1)

with the free-electron density in cm, the line-of-sight component of the magnetic field in (positive [negative] for a magnetic field pointing toward [away from] the observer) and the path length from the observer to the source measured in pc.2 Clearly, RM is a purely observational quantity, which can be defined only for a linearly-polarized radio source located behind the Faraday-rotating medium.

More generally, one may use the concept of Faraday depth (Burn, 1966; Brentjens & de Bruyn, 2005),

 FD(d)=(0.81 rad m−2) ∫d0ne B∥ ds , (2)

a truly physical quantity, which has the same formal expression as RM but can be defined at any point of the ISM, independent of any background source. FD simply corresponds to the line-of-sight depth of the considered point, , measured in terms of Faraday rotation. Here, we are only interested in the FD of the Galaxy, in which case represents the distance from the observer to the outer surface of the Galaxy along the considered line of sight.

As an observational reference for our modeling work, we adopt the all-sky map of the Galactic FD reconstructed by Oppermann et al. (2015), denoting the observational Galactic FD by (see Fig. 1). To create their map, Oppermann et al. compiled the existing catalogs of extragalactic RMs, for a total of 41 632 data points (37 543 of which are from the NVSS catalog of Taylor et al. (2009), which covers fairly homogeneously the sky at declination ), and they employed a sophisticated signal reconstruction algorithm that takes spatial correlations into account.

The all-sky map of in Fig. 1, like previous all-sky RM maps, shows some coherent structure on large scales. This large-scale structure in the RM sky has often been assumed to reflect, at least to some extent, the large-scale organization of the Galactic magnetic field (e.g., Simard-Normandin & Kronberg, 1980; Han et al., 1997; Taylor et al., 2009). In reality, however, nearby small-scale perturbations in the magneto-ionic ISM can also leave a large-scale imprint in the RM sky (see Frick et al., 2001; Mitra et al., 2003; Wolleben et al., 2010; Mao et al., 2010; Stil et al., 2011; Sun et al., 2015, and references therein). The most prominent such perturbation identified to date is the nearby magnetized bubble uncovered by Wolleben et al. (2010) through RM synthesis (a.k.a. Faraday tomography) on polarization data from the Global Magneto-Ionic Medium Survey (GMIMS). This bubble, estimated to lie at distances around , is centered at and extends over in longitude and in latitude, thereby covering of the sky.

Another potential source of strong contamination in the RM sky is the North Polar Spur (NPS), which extends from the Galactic plane at nearly all the way up to the north Galactic pole. However, Sun et al. (2015) showed, through Faraday tomography, that the actual Faraday thickness (FD) of the NPS is most likely close to zero. They also showed that the FD of the Galactic ISM behind the NPS cannot account for the entire Galactic FD toward the section of the NPS around . From this they concluded that if the NPS is local (such that the FD of the ISM in front of the NPS is very small), the Galactic FD must be dominated by the FD of Wolleben et al.’s (2010) magnetized bubble, which must then be larger than estimated by Wolleben et al. Clearly, the argument can be taken the other way around: if the FD of the magnetized bubble is as estimated by Wolleben et al. (an assumption we will make here), the Galactic FD must have a large contribution from the ISM in front of the NPS, which implies that the section of the NPS around is not local. This last conclusion is consistent with the results of several recent studies, which systematically place the lower part of the NPS beyond a few 100 pc – more specifically: beyond the polarization horizon at 2.3 GHz, , for (from Faraday depolarization; Sun et al., 2014), behind the Aquila Rift, at a distance , for (from X-ray absorption; Sofue, 2015), and beyond the cloud complex distributed between and and probably much farther away (again from X-ray absorption; Lallement et al., 2016).

Since the focus of our paper is on the large-scale magnetic field, we will subtract from the map of Fig. 1 the estimated contribution from Wolleben et al.’s (2010) magnetized bubble.3 We will then proceed on the assumption that the resulting map can be relied on to constrain the large-scale magnetic field structure. As it turns out, the removal of Wolleben et al.’s bubble will systematically lead to an improvement in the quality of the fits.

For convenience, we bin the data and average them within the different bins, both before (Fig. 2) and after (Fig. 3) removal of Wolleben et al.’s (2010) bubble. Following a binning procedure similar to that proposed by Pshirkov et al. (2011), we divide the sky area into 18 longitudinal bands with latitudinal width , and we divide every longitudinal band into a number of bins chosen such that each of the two lowest-latitude bands contains 36 bins with longitudinal width and all the bins have roughly the same area . This leads to a total of 428 bins. For each of these bins, we compute the average before removal of Wolleben et al.’s bubble, and we plot it in Fig. 2, with a red circle or a blue square according to whether it is positive or negative. In both cases, the color intensity increases with the absolute value of the average .

We then repeat the binning and averaging steps after removal of Wolleben et al.’s bubble, and for visual purposes, we wash out the obviously anomalous bin at . The underlying RM anomaly can be identified as the H ii region Sh 2-27 around Oph (Harvey-Smith et al., 2011). Our exact treatment of the anomalous bin is of little importance, because its weight in the fitting procedure will be drastically reduced through an artificial tenfold increase of its estimated uncertainty, , in the expression of (Eq. 38). In practice, we just blend the anomalous bin into the background by replacing its average with that of the super-bin enclosing its 8 direct neighbors. The bin-averaged after removal of Wolleben et al.’s bubble and blending of the anomalous bin is denoted by .

Finally, since we are primarily interested in the magnetic field in the Galactic halo, we exclude the 72 bins with , which leaves us with a total of 356 bins. The of each of these bins is plotted in the top panel of Fig. 3, again with a red circle if and a blue square if . The map thus obtained provides the observational reference against which we will test our large-scale magnetic field models in Sect. 4.

### 2.2 General trends

A few trends emerge from the observational map of the average in Fig. 2:

• a rough antisymmetry with respect to the prime () meridian,

• a rough antisymmetry with respect to the midplane in the inner Galactic quadrants away from the plane ( and ),

• a rough symmetry with respect to the midplane in the inner quadrants close to the plane ( and ) and in the outer quadrants ().

These three trends appear to persist after removal of Wolleben et al.’s (2010) bubble (see top panel of Fig. 3, for ). We will naturally try to reproduce them with our magnetic field models, but before getting to the models, a few comments are in order.

First, the rough symmetry with respect to the midplane at and all suggests that the disk magnetic field is symmetric (or quadrupolar),4 as already pointed out by many authors (e.g., Rand & Lyne, 1994; Frick et al., 2001).

Second, the rough antisymmetry [symmetry] with respect to the midplane at and [] suggests one of the two following possibilities: either the halo magnetic field is antisymmetric in the inner Galaxy and symmetric in the outer Galaxy, or the halo magnetic field is everywhere antisymmetric, but only toward the inner Galaxy does its contribution to the Galactic FD at exceed the contribution from the disk magnetic field. The first possibility is probably not very realistic (although it cannot be completely ruled out): while the disk and halo fields could possibly have different vertical parities (see, e.g., Moss & Sokoloff, 2008; Moss et al., 2010), it seems likely that each field has by now evolved toward a single parity. The second possibility may sound a little counter-intuitive at first, but one has to remember that the halo field contribution to the Galactic FD is weighted by a lower free-electron density than the disk field contribution (especially toward the outer Galaxy); moreover, the halo field can very well be confined inside a smaller radius than the disk field. Such is the case in the double-torus picture originally sketched by Han (2002) and later modeled by Prouza & Šmída (2003); Sun et al. (2008); Jansson & Farrar (2012a): in these three models, the toroidal field of the halo falls off exponentially with at large radii, i.e., much faster than the disk field, which falls off approximately as .

Third, the rough antisymmetry with respect to the prime meridian (east-west antisymmetry) can a priori be explained by an axisymmetric, predominantly azimuthal magnetic field. However, the situation is a little more subtle.

For the disk, RM studies (mainly of Galactic pulsars, and also of extragalactic point sources) converge to show that the magnetic field is indeed predominantly azimuthal, though not with the same sign everywhere throughout the disk, i.e., the field must reverse direction with Galactic radius (e.g., Rand & Lyne, 1994; Han et al., 1999, 2006; Brown et al., 2007). The exact number and the radial locations of these reversals are still a matter of debate, which we do not whish to enter. Here, we are content to note that field reversals naturally arise in a bisymmetric (azimuthal modulation , with ) or higher-order () configuration, but can also be produced in the axisymmetric case (Ferrière & Schmitt, 2000).

For the halo, most existing models have settled on a purely azimuthal, and hence axisymmetric, magnetic field. Here, we inquire into the other possibilities.

Let us first ask whether a purely poloidal (possibly, but not necessarily, X-shape) magnetic field would be a viable alternative. Clearly, a purely poloidal field that is axisymmetric automatically leads to an FD pattern with east-west symmetry, inconsistent with the observed FD distribution. In contrast, a purely poloidal field that is bisymmetric can lead to a variety of FD patterns, depending on the azimuthal angle of maximum amplitude: if the maximum amplitude occurs in the azimuthal plane parallel to the plane of the sky () [in the azimuthal plane through the Sun ()], the FD pattern has east-west antisymmetry [symmetry], in agreement [disagreement] with the observed FD distribution. Hence, a purely poloidal magnetic field can possibly match the FD data, provided it is bisymmetric (or higher-order) and favorably oriented.

Let us now turn to the more general and more realistic situation where the magnetic field has both azimuthal and poloidal components, as required by dynamo theory. If this general field is axisymmetric, the observed east-west antisymmetry in FD implies that it must be mainly azimuthal. In contrast, if the general field is bisymmetric, it could in principle cover a broad range of orientation, from mostly azimuthal to mostly poloidal. What happens here is that, in the presence of an azimuthal field component, the azimuthal modulation, which is presumably carried along field lines, crosses azimuthal planes. The resulting FD pattern becomes more difficult to predict: it may become more complex and fluctuating than in the axisymmetric case, and large-scale longitudinal trends may be partly washed out, but a priori the predominance of neither the azimuthal nor the poloidal field component can be ruled out.

To sum up, the halo magnetic field could either be axisymmetric and mainly azimuthal or bisymmetric (or higher-order) with no clear constraint on its azimuthal-versus-poloidal status.

## 3 Magnetic field models

### 3.1 Magnetic fields in galactic halos

In Paper 1, we derived four different models of spiraling, possibly X-shape magnetic fields in galactic halos, which can be applied to the halo of our own Galaxy. Each of these models was initially defined by the shape of its field lines together with the distribution of the magnetic flux density on a given reference surface. Using the Euler formalism (e.g., Northrop, 1963; Stern, 1966), we were able to work out the corresponding analytical expression of the magnetic field vector, , as a function of Galactocentric cylindrical coordinates, . The main characteristics of the four halo-field models are summarized below, and a few representative field lines are plotted in Fig. 4.

In models A and B, we introduced a fixed reference radius, , and we labeled field lines by the height, , and the azimuthal angle, , at which they cross the centered vertical cylinder of radius . In models C and D, we introduced a fixed reference height, (more exactly, one reference height, , in model C where all field lines cross the galactic midplane, and two reference heights, , in model D where field lines do not cross the midplane), and we labeled field lines by the radius, , and the azimuthal angle, , at which they cross the horizontal plane (or one of the two horizontal planes) of height . Thus, in all models, the point of a given field line can be regarded as its footpoint on the reference surface (vertical cylinder of radius in models A and B and horizontal plane(s) of height in models C and D).

#### Poloidal field

The four models are distinguished by the shape of field lines associated with the poloidal field (hereafter referred to as the poloidal field lines), and hence by the expressions of the radial and vertical field components.

In model A, the shape of poloidal field lines is described by the quadratic function

 z=z1 1+ar21+ar21 , (3)

where is a strictly positive free parameter governing the opening of field lines away from the -axis, is the prescribed reference radius, and is the vertical label of the considered field line. Conversely, the vertical label of the field line passing through is given by

 z1=z 1+ar211+ar2 ⋅ (4)

It then follows (see Paper 1 for the detailed derivation) that the poloidal field components can be written as

 Br = r1r z1z Br(r1,φ1,z1) (5) Bz = 2ar1z11+ar2 Br(r1,φ1,z1) , (6)

with, for instance, obeying Eq. (11).

In model B, the corresponding equations read

 z=1n+1 z1 [(rr1)−n+nrr1] , (7)

with a power-law index satisfying the constraint ,

 z1=(n+1) z [(rr1)−n+nrr1]−1 , (8)

and

 Br = r1r z1z Br(r1,φ1,z1) (9) Bz = −nn+1 r1z21r2z [(rr1)−n−rr1] Br(r1,φ1,z1) ⋅ (10)

In both models A and B, the radial field component on the vertical cylinder of radius is chosen to have a linear-exponential variation with and a sinusoidal variation with :

 Br(r1,φ1,z1) = B1 fsym [|z1|H exp(−|z1|−HH)] (11) × cos(m(φ1−gφ(r1,z1)−φ⋆)) ,

where is the normalization field strength, is a factor setting the vertical parity of the magnetic field ( for symmetric fields and for antisymmetric fields; see footnote 2.2), is the exponential scale height, is the azimuthal wavenumber, is the shifted winding function defined by Eq. (23), and is the orientation angle of the azimuthal pattern. The term , which will be discussed in more detail in Sect. 3.1.2, ensures that the phase of the sinusoidal modulation remains constant on the winding surfaces defined by Eq. (26) (see comment following Eq. (26)).

Models C and D are the direct counterparts of models A and B, respectively, with the roles of the coordinates and inverted in the equations of field lines.

In model C, the reference height is set to , and the shape of poloidal field lines is described by the quadratic function

 r=r1 (1+az2) , (12)

with a strictly positive free parameter governing the opening of field lines away from the -axis and the radial label of the considered field line. Conversely, the radial label of the field line passing through is given by

 r1=r1+az2 ⋅ (13)

The poloidal field components can then be written as

 Br = 2ar31zr2 Bz(r1,φ1,z1) (14) Bz = r21r2 Bz(r1,φ1,z1) , (15)

with, for instance, obeying Eq. (20).

In model D, where every field line remains confined to one side of the galactic midplane, a reference height, , is prescribed on each side of the midplane, such that the ratio is always positive. We then have

 r=1n+1 r1 [(zz1)−n+nzz1] , (16)

with ,

 r1=(n+1) r [(zz1)−n+nzz1]−1 , (17)

and

 Br = −nn+1 r31r2z [(zz1)−n−zz1] Bz(r1,φ1,z1) (18) Bz = r21r2 Bz(r1,φ1,z1) ⋅ (19)

In both models C and D, the vertical field component on the horizontal plane(s) of height is chosen to have an exponential variation with and a sinusoidal variation with :

 Bz(r1,φ1,z1)=B1 ¯fsym exp(−r1L) cos(m(φ1−gφ(r1,z1)−φ⋆)) , (20)

with the normalization field strength, a factor setting the vertical parity of the magnetic field ( in model C, which is always antisymmetric, and in the antisymmetric version of model D, and in the symmetric version of model D; see footnote 2.2), the exponential scale length, the azimuthal wavenumber, the shifted winding function (Eq. 23), and the orientation angle of the azimuthal pattern.

#### Azimuthal field

In all four models, field lines are assumed to spiral up or down according to the equation

 φ=φ1+fφ(r,z) , (21)

where is a winding function starting from the field line’s footpoint on the reference surface and, therefore, satisfying . If we consider that field lines are somehow anchored in the external intergalactic medium,5 it proves more convenient to shift the starting point of the winding function to infinity. This can be done by letting

 fφ(r,z)=gφ(r,z)−gφ(r1,z1) , (22)

where is a shifted winding function starting from the field line’s anchor point at infinity6 and, therefore, satisfying for and for . A reasonably simple choice for the shifted winding function is

 gφ(r,z)=cotp0 ln[1−exp(−rLp)]1+(|z|Hp)2 , (23)

with the pitch angle (i.e., the angle between the horizontal projection of a field line and the local azimuthal direction) at the origin (), the scale height, and the exponential scale length. With this choice, defines, in horizontal planes (), spirals that are logarithmic (constant pitch angle) at small and become increasingly loose (increasing pitch angle) at large . The spirals also loosen up with increasing .

Combining Eqs. (21) and (22), we can rewrite the equation describing the spiraling of field lines in the form

 φ=gφ(r,z)+[φ1−gφ(r1,z1)] , (24)

with constant along field lines. Conversely, the field line passing through can be traced back to the azimuthal label

 φ1=φ−gφ(r,z)+gφ(r1,z1) , (25)

where and either have prescribed values or are given functions of (see Eqs. (4), (8), (13), and (17) in Sect. 3.1.1). It then follows from Eqs. (11) and (20) that the sinusoidal modulation of the magnetic field goes as .

If we now consider the anchor point at infinity of the field line passing through , denote its azimuthal angle by , and recall that for , we find that the spiraling of field line can also be described by

 φ=gφ(r,z)+φ∞ ⋅ (26)

Eq. (26) defines a so-called winding surface, formed by all the field lines with the same azimuthal angle at infinity, , and thus the same phase in the sinusoidal modulation, . The crest of the modulation occurs at phase , i.e., at , which means that the free parameter can be interpreted as the azimuthal angle at infinity of the crest surface.

Since the magnetic field is by definition tangent to field lines, Eq. (24) implies that its azimuthal component is related to its radial and vertical components through

 Bφ=(r ∂gφ∂r) Br+(r ∂gφ∂z) Bz ⋅ (27)

Roughly speaking, the two terms on the right-hand side of Eq. (27) represent the azimuthal field generated through shearing of radial and vertical fields by radial and vertical gradients in the galactic rotation rate, respectively. With the choice of Eq. (23), Eq. (27) becomes

 Bφ = cotp0 11+(|z|Hp)2 rLp exp(−rLp)1−exp(−rLp) Br (28) − 2 cotp0 zH2p⎡⎣1+(|z|Hp)2⎤⎦2 r ln[1−exp(−rLp)] Bz ⋅

It is easily seen that the factor of has the same sign as (usually negative) everywhere, while the factor of has the same sign as above the midplane () and the opposite sign below the midplane ().

#### Total field

To sum up, the poloidal field is described by Eqs. (5) – (6) in model A and Eqs. (9) – (10) in model B, with given by Eq. (11), and by Eqs. (14) – (15) in model C and Eqs. (18) – (19) in model D, with given by Eq. (20). The azimuthal field is described by Eq. (28) in all models. Together, the above equations link the magnetic field at an arbitrary point to the normal field component at the footpoint on the reference surface of the field line passing through . The footpoint, in turn, is determined by the reference coordinate ( in models A and B, and in models C and D, with in model C and in model D), the poloidal label of the field line ( given by Eq. (4) in model A and Eq. (8) in model B, and given by Eq. (13) in model C and Eq. (17) in model D), and its azimuthal label ( given by Eq. (25) in all models).

#### Regularization of model A

An important remark should be made regarding model A. Eq. (5) together with Eq. (4) imply that for , which reflects the fact that all field lines, from all azimuthal planes, converge to the -axis. The problem does not arise in model B, as shown by Eq. (9) together with Eq. (8), because all field lines are deflected vertically before reaching the -axis (see Fig. 4b). The singularity in model A is inescapable in the axisymmetric case ( in Eq. (11)), where has the same sign all around the -axis, so that, at every height, a non-zero magnetic flux reaches the -axis. However, the singularity can be removed in non-axisymmetric () configurations, where changes sign sinusoidally around the -axis, so that, at every height, no net magnetic flux actually reaches the -axis. It is then conceivable to detach field lines from the -axis, spread them apart inside a centered vertical cylinder, whose radius can be chosen to be (remember that is still a free parameter at this stage), and connect up two by two field lines with the same , opposite (with respect to a node of the sinusoidal modulation in Eq. (11)), and hence opposite , such as to ensure conservation of the magnetic flux as well as continutity in its sign.

For instance, in the bisymmetric () case, field lines with the same and opposite (with respect to ) can be connected up inside the vertical cylinder of radius through straight-line segments parallel to the direction . In this case, magnetic flux conservation (or continuity of ) across the surface of the cylinder gives for the magnetic field at an arbitrary point inside the cylinder:

 \@vec\boldmath{B} = Br(r1,gφ(r1,z1)+φ⋆,z1) ^e⋆(z1) (29) = B1 fsym [|z1|H exp(−|z1|−HH)] ^e⋆(z1) ,r≤r1 ,

with , the unit vector in the direction , and the other symbols having the same meaning as in Eq. (11).

In the following, the regularized bisymmetric version of model A is referred to as model A1, where the number 1 indicates the value of the azimuthal wavenumber.

### 3.2 Magnetic fields in galactic disks

We now present three different models of spiraling, mainly horizontal magnetic fields in galactic disks, which can be applied to the disk of our own Galaxy. These models are directly inspired from the three models for magnetic fields in galactic halos that have nearly horizontal field lines at low , namely, models A1, B, and D in Sect. 3.1, and they are named models Ad1, Bd, and Dd, respectively. Models Bd and Dd are directly taken up from Paper 1, with a minor amendment in model Dd (see below), while model Ad1 is a new addition.

Model Ad1 has nearly the same descriptive equations as model A1. Outside the reference radius, , the shape of field lines is described by Eq. (3), with the parameter assuming a much smaller value than for halo fields, and by Eq. (24); the vertical and azimuthal labels of the field line passing through are given by Eqs. (4) and (25); and the three magnetic field components obey Eqs. (5), (6), and (28), with

 Br(r1,φ1,z1)=B1 fsym exp(−|z1|H) cos(m(φ1−gφ(r1,z1)−φ⋆)) ⋅ (30)

All the symbols entering Eq. (30) have the same meaning as in Eq. (11); the difference between both equations resides in the -dependence of , which leads to a peak at (appropriate for halo fields) in Eq. (11) and a peak at (appropriate for disk fields) in Eq. (30). Inside , field lines are straight, horizontal (), and parallel to the unit vector in the direction , and the magnetic field vector is given by Eq. (29) adjusted to Eq. (30):

 \@vec\boldmath{B} = Br(r1,gφ(r1,z1)+φ⋆,z1) ^e⋆(z1) (31) = B1 fsym exp(−|z1|H) ^e⋆(z1) ,r≤r1 ⋅

Model Bd is a slight variant of model B. Given a reference radius, , the shape of field lines is described by

 z=12n+1 z1 [(rr1)−n+2n√rr1] , (32)

with , and by Eq. (24); the vertical and azimuthal labels of the field line passing through are given by

 z1=(2n+1) z [(rr1)−n+2n√rr1]−1 (33)

and by Eq. (25), respectively; and the three magnetic field components obey

 Br = r1r z1z Br(r1,φ1,z1) (34) Bz = −n2n+1 r1z21r2z [(rr1)−n−√rr1] Br(r1,φ1,z1) , (35)

and Eq. (28), together with Eq. (30).

Model Dd is nearly the same as model D. Given a reference height on each side of the midplane, , the shape of field lines is described by Eq. (16), with , and by Eq. (24); the radial and azimuthal labels of the field line passing through are given by Eqs. (17) and (25); and the three magnetic field components obey Eqs. (18), (19), and (28), with

 Bz(r1,φ1,z1)=B1 ¯fsym F(r1) cos(m(φ1−gφ(r1,z1)−φ⋆)) (36)

and

 Unknown environment '%

Eq. (36) is the counterpart of Eq. (20), with the radial exponential cut off for . The effect of this cutoff is to reduce the crowding of field lines along the midplane.

## 4 Method

### 4.1 Simulated Galactic Faraday depth

As explained earlier, our purpose is to study the structure of the magnetic field in the Galactic halo. We are not directly interested in the exact attributes of the disk magnetic field, but we may not ignore the disk field altogether: since all sightlines originating from the Sun pass through the disk (if only for a short distance before entering the halo), we need to have a good description of the disk field in the interstellar vicinity of the Sun.

Therefore, all our Galactic magnetic field models must include a disk field and a halo field. Moreover, in view of the general trends discussed in Sect. 2.2, we expect the disk field to be roughly symmetric and the halo field roughly antisymmetric with respect to the midplane. Alternatively, we may consider that the Galactic magnetic field is composed of a symmetric field, which dominates in the disk, and an antisymmetric field, which dominates in the halo. In other words, we may write the complete Galactic magnetic field as

 \@vec\boldmath{B}=\@vec\boldmath{B}S+\@vec\boldmath{B}A , (37)

where is the symmetric field, described by one of the disk-field models (model Ad1, Bd, or Dd; see Sect. 3.2), and is the antisymmetric field, described by one of the halo-field models (model A1, B, C, or D; see Sect. 3.1).

Our knowledge of the large-scale magnetic field at the Sun, , imposes two constraints on our magnetic field models. Since the Sun is assumed to lie in the midplane (), where vanishes, these constraints apply only to . From the RMs of nearby pulsars, Rand & Kulkarni (1989) derived a local field strength , while Rand & Lyne (1994); Han & Qiao (1994) obtained . All three studies found to be running clockwise about the rotation axis (). Han & Qiao (1994) also derived a local pitch angle (implying ), close to the value inferred from the polarization of nearby stars (Heiles, 1996). In the following, we restrict the range of to and the range of to . The implications of these restrictions are discussed at the end of Appendix C.1.

To calculate the Galactic FD (Eq. 2), we also need a model for the free-electron density, . Here, we adopt the NE2001 model of Cordes & Lazio (2002), with revised values for the parameters of the thick disk component. Gaensler et al. (2008) found that the scale height of the thick disk in the NE2001 model was underestimated by almost a factor of 2; they derived a new scale height of 1.83 kpc (up from 0.95 kpc) and a new midplane density of 0.014 cm (down from 0.035 cm). This revision prompted several authors to adopt the NE2001 model and simply replace the parameters of the thick disk with the new values derived by Gaensler et al., while keeping all the other components unchanged (Sun & Reich, 2010; Pshirkov et al., 2011; Jansson & Farrar, 2012a). However, Schnitzeler (2012) showed that this simple modification entails internal inconsistencies, as the thick disk of Gaensler et al. overlaps with the local ISM and local spiral arm components of the NE2001 model. By correctly accounting for all the components of the NE2001 model, Schnitzeler (2012) updated the thick disk with a scale height of 1.31 kpc and a midplane density of 0.016 cm. These are the values that we adopt for our FD calculation. For future reference, the scale length of the thick disk is .

For any considered magnetic field model, we can now use Eq. (2) to calculate the modeled Galactic FD, , as a function of position in the sky. To produce a modeled FD map that is directly comparable to the observational map of displayed in the top panel of Fig. 3, we divide the sky area into the same 428 bins as in the map, and we only retain the same bins with . For each of these bins, we compute along 20 well-separated sightlines (which is sufficient, in view of the smoothness of our large-scale field models, to reach the required accuracy), and we average the computed to obtain . The resulting modeled map of can then be compared to the observational map of .

To disentangle the impact of the different model parameters and to guide (and ultimately speed up) the global fitting process, it helps to first consider and separately, denoting the associated by and , respectively. If the slight vertical asymmetry in the free-electron density distribution (which essentially arises from the local ISM component in the NE2001 model) is ignored, and represent the symmetric and antisymmetric parts of . Accordingly, we will first adjust to and to , where and are the symmetric and antisymmetric parts of (displayed in the middle and bottom panels of Fig. 3). We will then combine the separate best-fit and into a starting-guess for the actual fit of to .

### 4.2 Fitting procedure

For each magnetic field model, our fitting procedure relies on the least-square method, which consists of minimizing the parameter defined as

 χ2=Nbin∑i=1(¯¯¯¯¯¯¯¯FDobs,i−¯¯¯¯¯¯¯¯FDmod,i)2σ2i , (38)

where is the total number of retained bins, subscript is the bin identifier, and are the average values of the observational and modeled Galactic FDs, respectively, within bin , and is the uncertainty in (estimated in Appendix A). We also define the reduced parameter,

 χ2red=χ2Nbin−Npar , (39)

where is the number of free parameters entering the considered magnetic field model. The advantage of using instead of is that gives a more direct idea of how good a fit is. In principle, a fit with [] is good [bad]. Here, however, because of the large uncertainties in (see Appendix A), we do not ascribe too much reality to the absolute , but we rely on the relative to rank the different fits by relative likelihood. In the following, the minimum values of and are denoted by and , respectively.

Because of the relatively large number of free parameters, a systematic grid search over the -dimensional parameter space would be prohibitive in terms of computing time. Therefore, we opt for the much more efficient MCMC sampling method. We employ two different versions of this method, corresponding to two different implementations of the Metropolis-Hastings algorithm (Metropolis et al., 1953; Hastings, 1970): in the first version, the proposal distribution is centered on an educated guess for the best-fit solution, whereas in the second version (random walk Metropolis-Hastings algorithm), the proposal distribution is centered on the current point of the Markov chain (see, e.g., Robert, 2015, and references therein). After verifying that both versions lead to the same results (within the Monte Carlo uncertainties) and remarking that the second version is systematically faster, we use only the second version in most of our numerous trial runs. However, in our last series of runs (leading to the final results listed in Table 2), we use again both versions as a reliability check.

For every minimization, we run several, gradually more focused simulations, whose sequence is optimized by trial and error to converge as fast as possible to the correct target distribution. The key is to find a good compromise between sufficiently broad coverage of the parameter space (to make sure the absolute minimum is within reach) and sufficiently narrow bracketing of the minimum- solution (to recover the target distribution in a manageable time). We begin with a simulation having broad proposal () and target () distributions. For the former, we generally adopt a -dimensional Gaussian with adjustable widths, which we take initially large enough to encompass all plausible parameter values. For the latter, we adopt , with tuned in parallel with the Gaussian widths to achieve the optimal acceptance rate for the next proposed point in the Markov chain (see, e.g., Roberts et al., 1997; Jaffe et al., 2010). Once a likelihood peak clearly emerges, we start a new simulation with narrower (smaller Gaussian widths) and (smaller ) and with the starting point of the Markov chain close to the emerging peak. We repeat the procedure a few times until , as ultimately required. Note that we prefer to keep a handle on how the proposal and target distributions are gradually narrowed down, rather than implement one of the existing adaptive algorithms that can be found in the literature.

We monitor convergence to the target distribution by visually inspecting the evolution of the histogram, the running mean, and the minimum- (i.e., best-fit) value of each parameter as well as the histogram of itself. We stop sampling when a stationary state appears to have been reached, and after checking that the best-fit value of each parameter has achieved an accuracy better than 10% of the half-width of its confidence interval.7 The latter is taken to be the projection onto the parameter’s axis of the -dimensional region containing the points of the Markov chain (or second half thereof) with the lowest . For our final results (listed in Table 2), we simulate multiple Markov chains in parallel, with overdispersed starting points, we discard the first halves of the chains to skip the burn-in phase, and relying on the Gelman-Rubin convergence diagnostic (Gelman & Rubin, 1992), we consider that our chains have converged when, for each parameter, the variance of the chain means has dropped below 10% of the mean of the chain variances, or nearly equivalently, when the potential scale reduction factor, , has dropped below 1.05.

As mentioned at the end of Sect. 4.1, we apply our fitting procedure first to the symmetric and antisymmetric components of the Galactic FD separately, then to the total (symmetric + antisymmetric) Galactic FD. The uncertainties in and , and , are related to the uncertainties in and , and , through , where subscript is used to denote the bin that is the symmetric of bin with respect to the midplane.

## 5 Results

### 5.1 Preliminary remarks

Before considering any particular magnetic field model, let us make a few general remarks that will guide us in our exploration of the different parameter spaces and help us optimize the Metropolis-Hastings MCMC algorithm used for our fitting procedure (see Sect. 4.2) – in particular, by making an informed initial choice for the proposal distribution, . These remarks will also enable us to check a posteriori that the best-fit solutions do indeed conform to our physical expectation.

First, the all-sky maps of the symmetric and antisymmetric components of the bin-averaged observational Galactic Faraday depth, and , displayed in the middle and bottom panels of Fig. 3, show that is broadly distributed, both in latitude and in longitude, while the distribution of is thinner in latitude10 and more heavily weighted toward the outer Galaxy. This suggests that the symmetric field, , has a flatter distribution, which extends out to larger radii, while the antisymmetric field, , has a more spherical distribution, which extends up to larger heights. Clearly, these different distributions support our choosing disk-field models (Ad1, Bd, Dd) to describe and halo-field models (A1, B, C, D) to describe . Furthermore, out of our four halo-field models, we may already anticipate that model C, where the field strength has the most spherical distribution (see Fig. 4c), is the best candidate to represent the antisymmetric field. This expectation is in line with the notion that models A, B, and D, where low- field lines are nearly horizontal (see Figs.