DERIVING THE GLOBAL STRUCTURE OF THE GALACTIC MAGNETIC FIELD FROM FARADAY ROTATION MEASURES OF EXTRAGALACTIC SOURCES

A realistic model of the Galactic magnetic field (GMF) is needed for various applications, such as direct dark matter searches, studies of cosmic rays, and others. In particular, the GMF plays a crucial role in the propagation of ultrahigh energy cosmic rays (UHECRs)—those with energies in excess of 10¹⁹ eV. Their deflections (assuming these are charged particles) in the GMF are large enough to prevent the identification of some, perhaps most, sources directly from the UHECR data. However, in the case of protons these deflections are sufficiently small to be corrected for, provided the GMF is known with enough accuracy. Thus, the possibility of astronomy with charged particles at ultrahigh energies may depend crucially on our knowledge of the GMF.

The magnetic field (MF) of our Galaxy is thought to contain both regular and turbulent components (see, e.g., Beck 2001, 2008). While the regular component is subdominant in strength by a factor of a few, it is expected to give a dominant contribution in integral quantities such as total deflections of ultrahigh energy particles. Unlike the turbulent component, the regular component adds up coherently and comes to dominate (Tinyakov & Tkachev 2005) when integrated over many coherence lengths. Estimating the strength and global structure of the regular MF of our Galaxy is the main purpose of this paper. Our primary interest is in the MF outside the thin Galactic disk, since that is the field that determines the deflections for most arrival directions of UHECRs.

The regular component of the GMF has been studied previously (e.g., Simard-Normandin & Kronberg 1980; Sofue & Fujimoto 1983; Han & Qiao 1994; Frick et al. 2001; Vallée 2005; Han et al. 2006; Brown et al. 2007; Sun et al. 2008; Jansson et al. 2009), and in particular in the context of UHECR propagation (Stanev 1997; Tinyakov & Tkachev 2002; Prouza & Šmída 2003). The advance of this paper is that it is based on improved data on rotation measures (RMs) and the ionized gas distribution. We also develop improved numerical models commensurate with these improved data.

Observations of other spiral galaxies reveal MFs of a few μG strength, coherent over distances of several kiloparsecs (see, e.g., Beck 2008; Krause 2009; Braun et al. 2010). Observations from the face-on and edge-on perspectives indicate that these fields consist of at least two components: the "disk" component which is concentrated toward the galactic plane and has largely a spiral pattern and the "halo" component that resides at some distance from the disk plane and has radial size comparable with the size of the galaxy. As the Milky Way seems to be a typical member of the entire class of spiral galaxies, MFs of similar strength and configuration should be expected in our Galaxy.

There are several methods to estimate the MF in the Galaxy, such as Zeeman splitting measurements (Crutcher 1999), infrared polarization studies (Nishiyama et al. 2010), synchrotron polarization and intensity surveys (Wielebinski 2005; Sun et al. 2008; Jansson et al. 2009; Jaffe et al. 2010; Ruiz-Granados et al. 2010), starlight polarization studies (Heiles 1996), and, finally, observation of the Faraday RMs of different Galactic pulsars and extragalactic radio sources (Simard-Normandin et al. 1981; Han et al. 2006; Men et al. 2008; Mao et al. 2010; Kronberg & Newton-McGee 2011). Except for the Zeeman measurements all these methods are somewhat indirect as they involve additional assumptions about the density of the warm ionized medium or concentrations of cosmic ray electrons. In this paper we use recent data on the Faraday RMs.

When propagating through a magnetized plasma, the polarization plane of a linearly polarized electromagnetic wave of wavelength λ rotates by the angle Δψ proportional to the square of the wavelength,

$$\begin{equation} \Delta \psi = {\rm RM}\cdot \lambda ^2. \end{equation} \tag{ 1 }$$

Multi-frequency observations thus allow one to measure the value of RM. The RM can be expressed in terms of the properties of the intergalactic/interstellar medium (IGM/ISM) and permeating MFs:

$$\begin{equation} {\rm RM} = 0.81\int _{0}^{D} n_{\rm e}B_{||} dl, \end{equation} \tag{ 2 }$$

where n_e is the density of free electrons measured in cm⁻³, B_|| is the component of the MF parallel to the line of sight measured in μG (positive when directed toward the observer), and D is the distance from the observer to the source in pc. Combined with data on the Galactic free electron distribution, Equations (1) and (2) can be used to estimate the GMF from the Faraday RMs of a large number of sources.

Two types of sources are normally used to measure Galactic Faraday rotation: Galactic pulsars and extragalactic radio sources. Pulsars have an advantage that their intrinsic RM values are negligibly small. However, the use of pulsars requires the knowledge of their distances (usually inferred from their dispersion measures), which introduces additional assumptions about the ISM density distribution. In addition, pulsars are concentrated around the Galactic plane and their total number (∼500) is currently much smaller than the number of extragalactic radio sources with measured RMs. For these reasons, in this paper we concentrate on the RMs of the extragalactic radio sources.

Our general strategy is the following. We consider several types of analytical GMF models, each containing a finite number of parameters. For each model we calculate the expected RMs and fit them to the observations by the binned χ² method. The quality of the fit is then used to determine the best values of the GMF parameters and their allowed ranges.

The paper is organized as follows. Section 2 describes the data. Sections 3 and 4 describe the models used and our method. The results are presented in Section 5. Section 6 contains the results and conclusions.

In the present work we use two sets of RM data. The first is the compilation of RMs recently obtained by Taylor et al. (2009) who reanalyzed NRAO VLA Sky Survey (NVSS) data. NVSS is the largest—by number survey—of polarized radio sources at declinations > − 40° (Condon et al. 1998). It was made in two nearby bands, 1364.9 and 1435.1 MHz, each having a width of 42 MHz. Simultaneous observations in these two different bands give estimates of the RMs of the sources. The total number of sources observed was 37,543 with a mean error, estimated by the authors, of ∼11 rad m⁻². The sky map of the NVSS⁸ set is shown in Figure 1, which displays the location and RM sign of each source. The data have a large blind spot corresponding to declinations below ∼ − 40°.

Zoom In Zoom Out Reset image size

An advantage of the NVSS RM catalog is the large number of sources. However, observations being made in two close bands lead to relatively large individual RM errors. Also, because in general the RM is not always constant in λ², it can have modest variations over a wider range of λ. Thus an NVSS RM, measured in just one small wavelength range at 21 cm (λ = 20.9–22.0 cm), can differ from an RM measured over a wider wavelength range which tends to average out fluctuations in RM(λ²). Such systematic differences will add to the mean error of ∼11 rad m⁻² quoted above by up to ∼13 rad m⁻² and will contribute to the widths of the RM comparison plots with the results from the Kronberg & Newton-McGee (2011) compilation (KNM11) shown in Figure 3 below. On the other hand, because λ²₂ − λ₁² is small, ambiguities of nπ in the NVSS RMs are unlikely to occur, especially at |b| > 10°.

It can be easily seen that the NVSS data set should be treated with care when considering the sources with large observed RMs, which mainly reside near the Galactic plane and in the region of the inner Galaxy. In order to avoid this problem, we use the NVSS data only at Galactic latitudes |b| > 10°.

The second all-sky data set we used is the compilation by Kronberg & Newton-McGee (2011) (2257 sources). It consists of ∼1500 revised and statistically more accurate extragalactic radio source RMs of P. P. Kronberg et al. (2011, in preparation). These are derived from polarization measurements made over a large "baseline"of λ². The remaining RMs are taken from the Canadian Galactic Plane Survey (Brown et al. 2003), the Southern Galactic Plane Survey (Brown et al. 2007), and two smaller published lists from Klein et al. (2003) (108 RMs) and Mao et al. (2008) (68 RMs). The extragalactic radio sources in the KNM11 set were carefully observed at multiple frequencies and give RMs largely free from ambiguities. Errors in the RM measurements in the KNM11 set are typically of order ∼3 rad m⁻². We use this set for control.

The sky map of the KNM11 RMs is shown in Figure 2, upper panel. The lower panel shows the same data smoothed over the circle of 21° centered at the location of each source as described by Kronberg & Newton-McGee (2011). One may observe a generally good agreement in the RM sky distributions of the NVSS and KNM11. This is particularly visible when comparing Figure 1 with the lower panel of Figure 2.

Zoom In Zoom Out Reset image size

The RM sky structure is rather peculiar. The regions of positive and negative RMs are roughly symmetric with respect to the Galactic plane in the outer Galaxy 90° < l < 270° and antisymmetric in the inner Galaxy. Note, however, that the boundary between positive and negative RMs in the inner Galaxy is displaced away from the Galactic plane (toward negative b for 0 < l < 90°), that is, RMs do not change sign across the Galactic plane itself. This fact was observed, e.g., in Kronberg & Newton-McGee (2011). This overall structure of RMs suggests, qualitatively, a Galactic field structure that consists of two components, with the disk field being symmetric in b with respect to the Galactic plane and the halo field antisymmetric.

Although the second data set is significantly smaller, it has much smaller individual RM errors and can be used in two ways. First, it allows for an independent evaluation of the quality of the NVSS data. In order to assess the accuracy of the NVSS RMs, we identified, by positional cross-correlation, the sources that are common to both catalogs. We did this separately for the two regions |b| > 10° and |b| < 10°. We found 1338 pairs and 306 pairs of sources, respectively, separated by less than 005. For these sources, we built distributions of differences between RMs given in the NVSS catalog and in the set of KNM11. The difference distributions are shown in Figure 3.

Zoom In Zoom Out Reset image size

In both regions, the distributions of ΔRM have peaks centered at zero and are well fitted by Gaussians with dispersions σ ≃ 13 rad m⁻² and σ ≃ 36 rad m⁻², respectively. In the region |b| > 10° the tails contain ∼25% of points and are confined within ±100 rad m⁻². The effective error introduced by these tails is of the order of the Gaussian part. Thus, in this region the NVSS data can be used for our purposes.

On the contrary, in the region |b| < 10° the tails (containing a similar fraction of points) are spread over a very wide range, ±900 rad m⁻². Such tails would introduce errors in the low-b bins which are too large to constrain the fit. For this reason we exclude the bins with |b| < 10° from the NVSS fit, as will be described in Section 4.

The KNM11 RM set provides statistically more accurate RMs in the Galactic plane, including areas where the NVSS data do not give sky coverage. This allows us to check whether the models which fit the NVSS data also fit the RM data in the Galactic plane.

We adopt a general model of the GMF consisting of two different components: a disk field and a halo field (Prouza & Šmída 2003; Sun et al. 2008). Each of these two components is parameterized independently. According to the qualitative features of the RM distribution discussed above, for most of this work we take the disk field to be symmetric with respect to the Galactic plane and the halo field to be antisymmetric (other possibilities will be discussed in Section 5). The general layout of the model is shown in Figure 4.

Zoom In Zoom Out Reset image size

3.1. The Disk Field

We consider several disk GMF models. The first one is the widely used logarithmic spiral model (e.g., Simard-Normandin & Kronberg 1980; Han & Qiao 1994; Stanev 1997; Tinyakov & Tkachev 2002). There are two versions of this model depending on whether the direction of the field in two different arms is the same⁹ (axisymmetric or ASS model) or opposite (bisymmetric or BSS model). They are defined as follows:

$$B_{\theta } = B\cos {p},\quad B_{r} = B\sin {p},$$

$$\begin{equation} B(r,\theta, z)=B(r)\left|\cos {\left(\theta - b\ln {\frac{r}{R_{\odot }}} +\phi \right)}\right|\exp {(-|z|/z_0)}, \end{equation} \tag{ 3 }$$

or

$$\begin{equation} B(r,\theta,z)=B(r)\cos {\left(\theta - b \ln {\frac{r}{R_{\odot }}} +\phi \right)}\exp {(-|z|/z_0)}, \end{equation} \tag{ 4 }$$

for ASS and BSS models, respectively. Here

$$b\equiv \frac{1}{{\rm tan}\, {p}}$$

$$\phi =b \ln \left(1+\frac{d}{R_{\odot }}\right)-\frac{\pi }{2},$$

where p is the pitch angle and d is the distance to the first field reversal. Negative d means that the nearest reversal occurs in the direction to the Galactic center; positive d corresponds to the opposite direction. R_☉ is the distance to the Galactic center (we adopt R_☉ = 8.5 kpc in this paper). The amplitude of the GMF is a function of the radial coordinate r,

$$\begin{equation} B(r)= \left\lbrace \begin{array}{@{}ll}\displaystyle B_0\frac{R_{\odot }}{r\cos {\phi }}, \quad &r>R_c,\\\ms \displaystyle B_0\frac{R_{\odot }}{R_c\cos {\phi }}, & r<R_c, \end{array} \right. \end{equation} \tag{ 5 }$$

where R_c is the radius of the central region where the disk field is assumed to have constant magnitude. In the course of the simulations, the parameters were varied within the following ranges:

Parameter	Minimum	Maximum
Rc	3 kpc	6 kpc
z0	0.5 kpc	1.5 kpc
d	−1.4 kpc	1.4 kpc
p	−15°	15°

Download table as:  ASCIITypeset image

The second model of the disk field was the axisymmetric field with reversals in concentric rings. It was taken from Sun et al. (2008):

$$B(r,\theta,z)=D_1(r,z)D_2(r)$$

$$\begin{equation} D_1(r,z)=\left\lbrace \begin{array}{@{}ll}\displaystyle B_0\exp \left(-\frac{r-R_{\odot }}{R_0} -\frac{|z|}{z_0}\right),\; & r>R_c\\ \displaystyle B_0\exp \left(-\frac{|z|}{z_0}\right), &r\le R_c \end{array} \right. \end{equation} \tag{ 6 }$$

$$\begin{equation} D_2(r)= \left\lbrace \begin{array}{@{}ll}+1,\quad &r>7.5\, \mathrm{kpc},\\ -1,&6\, \mathrm{kpc}<r\le 7.5\, \mathrm{kpc},\\ +1,&5\, \mathrm{kpc}<r\le 6\, \mathrm{kpc},\\ -1,&r\le 5\, \mathrm{kpc}, \end{array} \right. \end{equation} \tag{ 7 }$$

where R₀ is radial scale of the MF. We studied the following parameter space:

Parameter	Minimum	Maximum
Rc	4.5 kpc	5.5 kpc
R0	8 kpc	15 kpc
z0	0.5 kpc	1.5 kpc

Download table as:  ASCIITypeset image

We also varied slightly the radial positions of the inner MF reversals (the boundaries of regions of different sign in Equation (7)) in the ranges 4.8–5.2 kpc, 5.8–6.2 kpc, and 7.3–8.0 kpc.

The magnitude of the disk MF in the vicinity of the solar system, B₀, was taken to be 2 μG. Note that, in the total GMF, the disk and halo components will be combined with arbitrary coefficients (to be determined by fitting), so the choice B₀ = 2 μG at this point is merely a matter of definition.

3.2. The Halo Model

The basic halo model was taken from Prouza & Šmída (2003) and Sun et al. (2008):

$$\begin{equation} B_{\theta }^H(r,z)= B_0^H\left[1+\left(\frac{|z|-z_0^H}{z_1^H}\right)^2\right]^{-1} \frac{r}{R_0^H} \exp {\left(1- {r\over R_0^H} \right)}, \end{equation} \tag{ 8 }$$

where the direction of the field is reversed in the southern hemisphere (Brandenburg et al. 1992; Han et al. 1997). We repeat here the following definitions from Sun et al. (2008): B^H₀ is the halo MF strength, R^H₀ is its radial scale, the parameter z^H₀ defines the vertical position of the halo, and z^H₁ is its vertical scale; this scale could differ in directions toward and away from the Galactic plane: |z| < z^H₀ and |z| > z^H₀, respectively. We will refer to these parameters as z^H₁₍₁₎ and z^H₁₍₂₎. Here and below the index "H" denotes parameters belonging to the halo. The original values of the parameters in the model of Sun et al. (2008) were B^H₀ = 10 μG, R^H₀ = 4 kpc, z^H₀ = 1.5 kpc, and z^H₁ = 0.2 kpc for |z| < z^H₀; otherwise z^H₁ = 0.4 kpc(z^H₁₍₁₎ = 0.2 kpc and z^H₁₍₂₎ = 0.4 kpc). We also tested halo models with slightly different radial profiles:

$$\begin{equation} B_{\theta }^H(r,z)= B_0^H \left[ 1+\left(\frac{|z| -z_0^H}{z_1^H}\right)^2\right]^{-1} \exp {\left(-\left(\frac{r -R_0^H}{R_0^H}\right)^2\right)}, \end{equation} \tag{ 9 }$$

and

$$\begin{equation} B_{\theta }^H(r,z)= B_0^H \left[ 1+\left(\frac{|z|- z_0^H}{z_1^H}\right)^2\right]^{-1} \exp {\left(-\left|\frac{r -R_0^H}{R_T}\right|\right)}. \end{equation} \tag{ 10 }$$

We investigated the following parameter space:

Parameter	Minimum	Maximum
RH0	3.5 kpc	15 kpc
RT	1 kpc	5 kpc
zH0	0.8 kpc	3.5 kpc
zH1(1)	0.2 kpc	0.4 kpc
zH1(2)	0.3 kpc	0.5 kpc

Download table as:  ASCIITypeset image

Given that the halo field has not been firmly established as entirely azimuthal, we further studied the basic halo model, varying the pitch angle within the range −10° < p^H < 10°.

3.3. The Electron Density

4.1. Binning

RMs of individual sources are affected by uncertainties due to the intrinsic source contributions to the RM and by small-scale fluctuations of the GMF. The effect of the coherent component of the GMF can be revealed by averaging over a large number of nearby sources. The angular bin size should be small enough so that variations of the coherent component of GMF within the bin can be neglected, and it should include a sufficient number of sources to suppress source-to-source fluctuations. Applying these requirements, we chose the following binning: the celestial sphere is divided into 18 bands of 10° width, from −90° to +90° in Galactic latitude b. Each band is further divided into bins so that individual bins would have roughly equal area, about 100 deg², as shown in Figure 5. So, there are 36 bins in the 9th band (0° < b < 10°), but only four in the first one (80° < b < 90°). The total number of bins is equal to 422. After averaging, all the NVSS data are reduced to 422 values of average RM in each bin. No outliers were removed at this stage. Figures 6 and 7 show the averaged RMs in the eight bands with 10° < b < 50° and −50° < b < −10°, respectively. The coherent structure is clearly visible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The GMF models were treated in a similar way. For each model we simulated 20,000 values¹⁰ of RMs at randomly chosen locations and then applied the same binning and averaging procedure to these sets. The resulting averaged RMs were compared to observations.

A somewhat different procedure was used for the set of sources from KNM11. To check the disk field, we used the sources residing within 1° of the Galactic plane. As their number is more limited (about 240) and they are located irregularly within their respective bins (10 deg² in this case), we simulated model RM values in their exact locations to avoid an additional error that would be produced by different distributions of simulated and real sources within the bin.

4.2. Errors

4.3. Fitting Procedure

The best values of the GMF parameters were identified by the maximum likelihood method. This amounts to minimizing the value of the χ² defined as follows:

$$\begin{equation} \chi ^2 = \sum _{i=1}^{{N_{\rm tot}}} \frac{({\rm RM}_{\rm obs}-{\rm RM}_{\rm sim})_{i}^2}{\delta _{i}^2}, \end{equation} \tag{ 11 }$$

where N_tot is the total number of bins that we have included in the fit, (RM_obs − RM_sim)_i is the difference between the observed and simulated average RM values for the ith bin, and δ_i is the adopted error of the observed average.

Since the RMs (Equation (2)) are linear in the MF, the total RM of a given source is a simple sum of the disk and halo contributions. Making use of this property, we separately simulated more than 30,000 halo RM maps and 10,000 disk RM maps, corresponding to different parameters varied within the ranges described in Section 3. At this stage the strengths of the disk and halo components were held fixed and set to values given in Section 3. Then each disk configuration was combined with each halo configuration with two arbitrary coefficients α and β. For a global GMF model thus obtained, the value of χ² in Equation (11) is a quadratic polynomial in α and β, whose minimum with respect to α and β we found analytically. The resulting values of α and β determined the strength of the disk and halo fields.

The parameters of the GMF field were varied, within the ranges given in Section 3, by two different methods: by uniformly covering a given parameter range with a regular grid (grid-type simulation) and by randomly picking parameter values from a uniform distribution over the desired range (Monte Carlo type simulations). The disk field was simulated using the grid method, whereas for the halo we used both methods in roughly equal proportion.

Because of the uncertainty in the absolute normalization of χ² resulting from uncertainties in the error estimation (see Section 4.2), we did not combine the northern and southern Galactic hemispheres in a single fit. The two hemispheres are not on an equal footing because of the blind spot in the southern sky. Instead, our general fitting strategy was as follows. First we performed the fit using the NVSS data for the northern hemisphere only, thus obtaining the preliminary set of models. Next, we fitted the NVSS data in the southern hemisphere. We chose only the models that performed well in both tests and had the same value of the disk field. Note that this allowed for different halo fields in the northern and southern hemispheres. In these tests we imposed the additional constraint that $|\textbf {B}|$ be greater than 1.8 μG in both hemispheres to conform with the results of other observations. Models were selected as acceptable if they satisfied the following requirement:

$$\begin{equation} \chi ^2\big/\chi ^2_{\rm min}=\chi ^2_{\rm red}\big/\chi ^2_{\rm red, min}<1+\sqrt{ 2/N_{\rm {\rm dof}}}, \end{equation} \tag{ 12 }$$

where χ²_min (χ²_{red, min}) is a minimal value of χ² (χ²_red). These are the models having parameters within ∼1σ from the best-fit ones. Note that the criterion (12) is insensitive to the absolute normalization of the errors. We choose such a criterion to reduce the effect of the ambiguities in the error estimation discussed in Section 4.2.

Independently, we fitted the sources from the compilation of KNM11 in the narrow strip of the Galactic plane |b| < 1°. This imposed additional constraints on the structure of the field in the thin disk to which our fit of NVSS data is insensitive. Requiring that fits to both the NVSS and KNM11 data be good, we performed the final selection of the models.

Overall, we required that a successful model passes three separate fits: to the NVSS data in the northern and southern hemispheres, and to the KNM11 data in the disk. Note that this is, in general, more restrictive than performing a single combined fit.

Figure 8 shows a comparison between one of the best-fit models and the binned NVSS data. Positive average RMs are shown by red circles and negative average RMs by blue squares. The intensity of the color indicates the absolute value of RM. It is seen that the general structure of the field is reproduced quite well.

Zoom In Zoom Out Reset image size

The results of the fits are as follows.

North. The absolute minimum value of χ² was equal to 232 for 104 degrees of freedom (dof); χ²_red = 2.23. A bin-by-bin comparison between the data and the best-fit model is shown in Figure 9. Parameter ranges that were obtained from the fits are presented in Table 4 (see the Appendix). The fit was most sensitive to the pitch angle p. The vertical scale of the disk MF z₀ and the distance to the nearest reversal d affect the fit less. There was no significant influence of the R_c parameter, so we excluded it from further studies, fixing its value at R_c = 5 kpc. The same concerns the parameters z^H₁₍₁₎ and z^H₁₍₂₎ which were fixed at z^H₁₍₁₎ = 0.25 kpc and z^H₁₍₂₎ = 0.4 kpc.

Zoom In Zoom Out Reset image size

We could not make a successful fit using the disk component only. The addition of the halo component decreases χ²_min by a factor of 1.5. Due to the weaker effect of the halo MF, its parameters are estimated with larger uncertainties than the disk ones. In addition, there is a degeneracy between the height of the halo z^H₀, the amplitude of the halo field and, to some extent, the radial size R^H₀. Even unrealistically large values of the halo field of order tens of μG could give an acceptable fit if the halo is placed at a large distance from the Galactic plane where the concentration of the electrons greatly decreases (cf. Figure 10). There is also a degeneracy between the vertical scales of the GMF model (both halo and disk) and the vertical scale of the electron density (recall that the latter was held fixed during the fits).

Zoom In Zoom Out Reset image size

The strength of the disk field corresponding to the allowed set of models resides in the 1.8–2.8 μG range, largely concentrated around 2.0 μG. The amplitude of the allowed halo field is in the range 2–12 μG. Note that the halo field strength is constrained by the synchrotron radio emission from the population of Galactic relativistic electrons (Haslam et al. 1982). These constraints favor halo fields with the magnitude in the range 2–5 μG.

Ring models of the disk field performed rather poorly in the northern hemisphere no matter what the halo field was, giving a best value of χ²_min = 288 for 101 dof, which corresponds to χ²_red = 2.85 as compared to χ²_red = 2.23 for the best-fit model with the spiral disk. Thus, these models are disfavored by the NVSS data, confirming the earlier finding of Simard-Normandin & Kronberg (1980).

Finally, our attempts to improve the quality of the fit by introducing a non-zero halo pitch angle were unsuccessful: there was no significant change of χ² when a non-zero pitch angle of the halo field was varied within the range [ − 10°, 10°]. We also could not improve the fit by choosing the alternative model of the halo described by Equations (9) and (10).

South. The data from the southern hemisphere are more sparse because of the blind spot, and show more scatter, so the fits are generally worse. The best-fit ASS-type model gives χ² = 243 for 65 dof, which corresponds to χ²_red = 3.74. The results for the best-fit BSS model are χ² = 235 for 65 dof (χ²_red = 3.62). The bin-by-bin comparison between the data and the best-fit model in the southern hemisphere is shown in Figure 11.

Zoom In Zoom Out Reset image size

Disk. Finally, we fitted the RMs from KNM11 at very low latitudes |b| < 1°. The best-fit ASS-type model gives χ² = 30 for 29 dof (χ²_red = 1.03), while the best-fit BSS model gives χ² = 37 for 29 dof (χ²_red = 1.28). The comparison between these data and the best-fit model is shown in Figure 12. The results are fully compatible with the other fits. They also agree with the results of Kronberg & Newton-McGee (2011). Unfortunately, there is not enough information about the region around the direction to the Galactic center where the discrepancy between the predictions of different models is the largest. This reduces the discriminative power of the fit. Nevertheless, the fit was very sensitive to the field strength in the disk, so its results impose most stringent constraints on this parameter. Also, it was shown by Kronberg & Newton-McGee (2011) that the KNM11 data favor a pitch angle near the disk of −55. Taking this constraint into account allows us to exclude positive pitch angles, which is not possible from the fits to the NVSS data alone.

Zoom In Zoom Out Reset image size

In both north and south fits of the NVSS data, the reduced χ² values significantly exceed unity (especially in the southern hemisphere). This could result either from our underestimation of errors or from features in the data that cannot be accounted for by our models. Visual examination of the fits (Figures 9 and 10) shows that there exists a systematic discrepancy between the data and the best-fit models in the region l ≈ 130°–220° which corresponds to directions away from the Galactic center. The discrepancy is present in both northern and southern hemispheres.

In fact, the observed variation of RM in the outer Galaxy is too strong to be accounted for in any model that has a thick disk or halo with a nearly tangential field in the outer Galaxy (i.e., small pitch angle) if the field magnitude decays away from the Galactic center. So, either these models are completely wrong in this region or the large RMs are caused by some nearby anomaly in the MF or the electron density, which is sufficiently close to subtend a large angle of about ∼90°. The latter option is supported by the fact that the variation of RM which is responsible for this feature does not visibly decrease away from the Galactic plane, as it would if it were due to an additional Galactic arm or another global (and more distant) feature. Interestingly, the anomaly in the nearby region (the so-called region A) was pointed out previously (Simard-Normandin & Kronberg 1980).

For an additional check of this assumption we made a separate fit of RMs of the sources located in the direction of the inner Galaxy, i.e., at −90° < l < 90°. Excluding this apparent anomaly should improve the fit and show the extent to which the best-fit parameters were affected by its presence. The results of this additional test are as follows.

North-inner. The best fit gives χ²_min = 43 for 46 dof, considerably better than for the full sky. The allowed values of parameters obtained from the restricted fit largely coincide with those obtained from the full fit. This indicates that the full fit is not dominated by the anomalous region.

South-inner. χ²_min = 35 for 22 dof. The quality of the fit has considerably increased as compared to the fit in the full range of longitudes, but is still worse than in the northern hemisphere.

All our attempts to fit the outer Galaxy separately (90° < l < 270°), even those including an extra disk component, were unsuccessful in the sense that we could not get χ²_red smaller than five. This confirms our conclusion that it is not possible to fit the anomalous region with GMF models of the type considered here.

The results of individual fits are presented in Table 4 (see the Appendix) in the form of acceptable ranges of model parameters. It should be stressed that there are dependencies and degeneracies among them; Table 4 only shows the maximum variation of each parameter within the whole acceptable set of models. In other words, not every combination of acceptable parameters gives an acceptable model. Parameters of the models which fit the data best and are in agreement with all fits are summarized in Tables 1 and 2.

Finally, let us comment briefly on the other tests that we have made. First, as has been already mentioned in Section 3.3, we have run preliminary simulations with the original, lower vertical scale of electron density. For the realistic strength of the disk field B₀ > 1.8 μG and the original electron density, we found only χ² values larger than 390 as compared to 232 for the best-fit case. If we vary the magnitude of the disk field, the value of χ² becomes close to the best-fit value for field strengths ∼0.7 μG, which is clearly unacceptable. In conclusion, our simulations disfavor the original vertical scale of the electron density.

Second, we have tested different combinations of symmetries (with respect to the Galactic plane) of the halo and disk fields in order to check whether the NVSS data alone favor the symmetric disk and antisymmetric halo as we have chosen in our model. With our fitting procedure, which treats the northern and southern hemispheres independently, no additional simulations are required for this test. For instance, for an antisymmetric disk field we simply have to require that the coefficients with which the disk contribution enters the northern and southern hemispheres are opposite in sign. We could not find a satisfactory fit with either the ASS or with the BSS antisymmetric disk field, obtaining values of χ² more than four times larger than the best-fit value. In the case of the symmetric disk and symmetric halo, the best fit gives χ² 2.5 times larger than the best-fit model. Thus, the combination of the symmetric (with respect to b) disk and antisymmetric (with respect to b) halo is favored by the NVSS data.

The latest Faraday rotation data show an unmistakable coherent component of the RMs over the sky. The RMs are correlated in sign and value over the largest angular scales. Moreover, the distribution of RMs is correlated with the Galactic plane and the direction to the Galactic center. This strongly indicates the Galactic origin of this coherence and provides firm evidence for the existence of a regular GMF.

Apart from the Galactic plane region |b| ≲ 10°, the global pattern of the RMs has different symmetry with respect to the Galactic plane in the inner and outer Galaxy: symmetric outside and antisymmetric inside (see Figure 8). This change of symmetry cannot be achieved with a single (symmetric or antisymmetric) component of the GMF. At least two components—disk and halo—are required (Figure 4). Since the study of the MF in the disk shows that it does not change sign across the Galactic plane itself, the combination of the symmetric disk field and the antisymmetric halo field is favored. We have demonstrated that such a combination, with properly chosen amplitudes of components, reproduces the observed pattern of RMs very well.

Several further conclusions can be made.

• 1.  
  The cross-check between the NVSS and the smaller KNM11 set shows that the quality of the NVSS set is sufficient for our purposes over all the sky except in the |b| < 10° band around the Galactic plane.
• 2.  
  There is a feature in the observed RMs near the direction to the Galactic anticenter that cannot be accounted for by models of the type considered in this paper (nearly tangential field decaying away from the Galactic center). This feature may be due to some unspecified local structure.
• 3.  
  A ring model of the GMF is disfavored as compared to the model with logarithmic spiral arms. In the former case the best value of χ²_red is a factor of ⩾1.3 larger than in the latter.
• 4.  
  Our fits favor the increased vertical scale of the electron density distribution as proposed by Gaensler et al. (2008). We were unable to find any fit with the original vertical scale of the NE2001 model and strong disk field. This is in agreement with the results of Sun & Reich (2010).
• 5.  
  Our fits favor somewhat smaller values of the pitch angle p than have previously been discussed in the literature. These smaller values agree with the recent results of Kronberg & Newton-McGee (2011) where the reference value of p = −55 was obtained by a totally different method.

Our fitting procedure could not discriminate between the ASS and BSS models because a major contribution to the extragalactic source RMs (especially at |b| > 15°) comes from the near vicinity of the solar system; in this region both models behave quite similarly. A more detailed analysis of the Galactic plane region may be more sensitive to the differences between the ASS and BSS models.

Two "benchmark" models (one with the ASS and one with the BSS disk field) are presented in Table 3. These models fit the RM data well and have values of parameters favored by other studies (pitch angle p not too small, halo field not too strong, and the magnitude in the vicinity of the Earth close to 2 μG). They can be used, e.g., for propagation of cosmic rays through the Galaxy. At ultrahigh energies, an accurate three-dimensional description of the GMF may be crucial for identification of extragalactic cosmic ray sources.

In this paper we used only the Faraday RMs of extragalactic radio sources. More precise understanding of the GMF can be achieved by taking into account other data, also sensitive to the MF, such as  the RMs of pulsars, measurements of galactic synchrotron radiation, etc. This would permit better constraints on the unknown quantities such as the three-dimensional electron density profiles and other properties of the ISM. It will also help to resolve some degeneracies inherent in our models.

The authors are grateful to D. D. Sokoloff and JinLin Han for fruitful discussions. Also, we thank an anonymous referee for careful reading of the manuscript and valuable remarks. This work is supported by IISN project no. 4.4509.10 and the Natural Sciences and Engineering Research Council of Canada (NSERC), the Australian Research Council, the U.S. Department of Energy, and the use of NASA's Astrophysics Data System.

In our fitting procedure, each bin included in the analysis was characterized by an average RM and its error δ. Here we describe the algorithm by which these errors were assigned on a bin-by-bin basis.

First, we estimated the intrinsic dispersion of RMs of extragalactic sources. To this end we investigated the regions on the sky where the influence of the GMF is expected to be suppressed. Since RMs are proportional to the MF component that is parallel to the line of sight, these are the regions around the Galactic poles. Considering the sources with |b| > 80° we found that the dispersion of RMs of individual sources in these regions was about σ_qso = 15 rad m⁻². Comparing to Figure 3, we see that the dominant contribution to this number comes from the measurement errors. In the error assignment algorithm, we adopted the value of σ_qso as an indicator of whether a given bin is perturbed by effects other than the random scatter of sources.

If a bin showed the dispersion of RMs σ_i smaller than σ_qso we considered it unperturbed and assigned its error on a statistical basis. Since the average of a sample of N elements varies around the true average with the width $\sigma /\sqrt{N}$, where σ² is the variance of the distribution, we assigned to such bins an error of $\delta _i = \sigma _{\rm qso}/\sqrt{N_i}$, N_i being the number of sources in the bin. In practice, this situation occurs only for high-latitude bins.

In the opposite case, i.e., when the dispersion σ_i was larger than σ_qso, we considered the bin as perturbed by coherence effects. These may be due to, e.g., random MF with the coherence length comparable to the bin size, fluctuations of electron density, strong variation of the regular component across the bin, etc. In this case the statistical approach is not correct and a precise determination of the error is difficult. We have assumed that the error in such bins is some fraction of the bin variance σ_i; we have adopted 1/3 as the fraction, so that the error in perturbed bins was taken as δ_i = σ_i/3. We cannot rigorously derive this coefficient (although it can be motivated by a toy model presented below). Since this uncertainty propagates directly into the value of χ², it results in some uncertainty of the absolute normalization of the latter. Note, however, that the ratios of χ² are less affected by this uncertainty.

The following toy model shows how the proportionality between the error and the dispersion in the bin may come about. Let the random part of the MF consist of cells of a typical size d ∼ 100 pc, such that their contributions to RM have zero mean and the dispersion σ₀. If the cells were statistically independent, individual sources would have an RM dispersion given by $\sigma = \sigma _0 \sqrt{N_{\rm tot}}$, where N_tot is the total number of cells along the line of sight (we have neglected the internal scatter). However, for sources in one bin a number N_coh of nearby cells contributes coherently. If N_coh ≪ N_tot this will not change the dispersion σ over the bin. The error, on the contrary, will be dominated by the coherent contribution, $\delta = \sigma _0\sqrt{N_{\rm coh}} = \sigma \sqrt{N_{\rm coh}/N_{\rm tot}}$. Thus, the error and the dispersion are proportional to one another with the coefficient $\sqrt{N_{\rm coh}/N_{\rm tot}}$. The latter coefficient is purely geometrical. Taking N_coh ∼ 1/θ, where θ is the bin size, and N_tot ∼ L/d, where L is the length of the trajectory inside the GMF, we find δ ∼ σ/3. Interestingly, similar conclusions (proportionality and a similar coefficient) are obtained by considering a perturbation caused by a close group of a few sources having similar RMs substantially different from the rest of the sources in the bin.

As a last step, all bins were subjected to visual scrutiny, with the result that some of the bins were completely removed from the analysis. These are the bins showing obvious pathologies like small number of sources or non-uniform distribution across the bin (e.g., some bins at the boundary of the NVSS blind spot), compact aggregates of sources with RMs drastically different from the surrounding value, etc. Finally, some bins with such pathologies were assigned large errors (2–3 times larger than prescribed by the above algorithm) by hand. This was done in cases where the corresponding bin might serve as a good discriminator between models, even with a large error. A constant error of 40 rad m⁻² was used in this case. This value was chosen in such a way that these bins do not contribute significantly to the value of χ², at least near the best-fit point. We have checked that the results of the fits do not depend on the precise choice of these errors. The number of such bins is about 10%; they are easily identifiable in Figures 6 and 7.

Here is the complete list of excluded bins and bins with manually assigned errors.

Northern hemisphere. An error of 40 rad m⁻² was assigned to the bins in the region 30° < l < 60°, 10° < b < 40° (local structures described by Wolleben et al. 2010); to the bins in the area −115° < l < −90°, 20° < b < 30° (a prominent structure with large positive RMs, possibly related with the Gum nebula); to the bin located at (15°, 25°) (strong non-uniformity: only 12 sources are located at 10° < b < 15°, while there are 38 at 15° < b < 20°). Excluded from the analysis were the bins in the area −10° < l < 30°, 20° < b < 30° (unresolved structure in RM and a large positive spot at (− 10°, 25°)); bins in the area 70° < l < 90°, 10° < b < 20° (large negative spot centered at (80°, 15°); bin affected by Gum nebula at (l = −100°, b < 20°); the bin centered at (− 30°, 15°) (strongly non-uniform distribution).

Southern hemisphere. An error of 40 rad m⁻² was assigned to the bin (105°, −15°) (local structure resulting in positive excess); to the bin (165°, −15°) (positive spot); to the bin centered at (180°, −15°) (negative spot near the anticenter direction); and to the bin (170°, −25°) (large positive spot).

SUPPLEMENTARY FIGURES AND TABLES

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Ranges of Disk and Halo Parameters of the Acceptable Models

Disk	χ2	χ2red	p	z0	d	α	zH0	RH0	β
Models			(deg)	(kpc)	(kpc)		(kpc)	(kpc)
North
A	240	2.30	[−10, −3]	[0.5, 1.3]	[−1.1, −0.3]	[0.9, 2.0]	[1.0, 1.4]	[4, 15]	[0.4, 0.7]
B	232	2.23	[−6, −10]	[0.5, 1.3]	[−1.1, −0.5]	[0.9, 1.5]	[1.0, 1.4]	[6, 15]	[0.3, 1.0]
South
A	243	3.74	[−10, 10]	[0.5, 1.5]	[−1.1, 1.1]	[0.9, 2.0]	[1.0, 3.5]	[3.5, 15]	[0.2, 2.0]
B	235	3.62	[−10, 10]	[0.5, 1.1]	[−1.0, 1.0]	[0.2, 2.0]	[1.0, 3.5]	[3.5, 15]	[0.2, 2.0]
Plane
A	30	1.03	[−7, 10]	[0.5, 1.4]	[−1.4, 1.4]	[0.9, 1.1]
B	37	1.27	[−8, −3] ∪ [5, 8]	[0.5, 1.5]	[−0.8, 1.2]	[1.0, 1.2]
North-inner
A	43	0.93	[−5, 6]	[0.5, 1.2]	[−1.0, 0.9]	[0.9, 3.0]	[1.2, 2.4]	[6.0, 8.0]	[0.4, 1.2]
B	45	0.98	[−7, 9]	[0.5, 1.4]	[−0.6, 1.1]	[0.9, 1.1]	[1.0, 2.0]	[8.0, 15]	[0.4–1.2]
South-inner
A	40	1.81	[−7, 11]	[0.5, 1.3]	[−0.2, −0.9]	[0.9, 3.5]	[1.0, 2.5]	[8.5, 15]	[0.1, 0.6]
B	35	1.59	[−5, 9]	[0.5, 1.2]	[−0.3, −0.8]	[0.9, 1.2]	[1.2, 2.2]	[3.5, 5.5]	[0.2, 0.5]

Notes. A and B denote ASS and BSS models, respectively. The parameters α and β determine the strengths of the disk and halo components as explained in Section 4.3.

Download table as:  ASCIITypeset image