The Galactic Faraday depth sky revisited

Context. The Galactic Faraday depth sky is a tracer for both the Galactic magnetic ﬁeld and the thermal electron distribution. It was previously reconstructed from polarimetric measurements of extra-Galactic point sources. Aims. Here we improve on these works by using an updated inference algorithm and by taking into account the electron emission measure as traced by free–free emission measured by the Planck survey. In the future the data situation will improve drastically thanks to the next generation Faraday rotation measurements from the SKA and its pathﬁnders. Anticipating this, a further aim of this paper is to update the map reconstruction method with some of the latest developments in Bayesian imaging. Methods. To this end we made use of information ﬁeld theory, an inference scheme that is particularly powerful in cases of noisy and incomplete data. Results. We demonstrate the validity of the new algorithm by applying it to an existing data compilation. Even though we used exactly the same data set, a number of novel ﬁndings are made; for example, a non-parametric reconstruction of an overall amplitude ﬁeld resembles the free–free emission measure map of the Galaxy. Folding this emission measure map into the analysis provides more detailed predictions. The joint inference enables us to identify regions with deviations from the assumed correlations between the emission measure and Faraday data, thereby pointing us to Galactic structures with distinguishably different physics. We ﬁnd evidence for an alignment of the magnetic ﬁeld within the lines of sight along both directions of the Orion arm.


Introduction
The Faraday rotation effect is one of the primary sources of information on astrophysical and cosmological magnetic fields.This includes the fields of planets, stars, other galaxies and galaxy clusters, as well as more curious objects such as radio jets, lobes, and relics.In particular, the study of the Faraday rotation induced by the Milky Way magnetic field is of twofold importance.First of all it is an interesting research topic on its own, due its connection to the formation and structure of our home galaxy, but furthermore it also constitutes a non-negligible foreground component.All polarized light stemming from cosmological sources passes through the Galaxy and interacts with the Galactic magnetic field, which affects its polarization direction.Accurate and highly resolved Galactic foreground templates are therefore a necessary condition for any precision measurement of extra-galactic magnetic fields via Faraday rotation.Multiple efforts have been made to to map the Galactic Faraday sky, quite a few of them already to remarkable accuracy, e.g. in Frick et al. (2001); Johnston-Hollitt et al. (2004); Dineen & Coles (2005); Xu et al. (2006); Oppermann et al. (2012) and also in Oppermann et al. (2015).Especially the approach in Oppermann et al. (2012) is important to us, as we will rely on the same theoretical framework for our inference algorithms.Our goal in this work is to further sharpen our knowledge on the Faraday sky by exploiting correlations to bremsstrahlung measured from electron proton interaction in the interstellar medium, commonly known as the Galactic free-free emission.This is well motivated by observation as well as physical considerations, as we will outline later on.In a more general context, we hope to demonstrate an interesting test case of multi-messenger astronomy, and show that the synthesis of long existing data sets still can yield undiscovered information, under the condition that our a priori physical knowledge is employed and the data is treated in a consistent way.Our main result, the revised Faraday sky including the free-free data is shown in Fig. 1a and the revised map without free-free data in Fig. 1b.The results of our inference are publicly available 1 .The paper is structured as follows.In Sec. 2 we develop our updated reconstruction method.This involves an amplitude field aiming at representing the overall Galactic structure and a sign field representing magnetic reversals and smaller scaled structures.It turns out that the reconstructed amplitude field resembles the free-free emission measure (EM) of the Galaxy, which is physically plausible.In order to take advantage of this, we extend the method to fold in the measurements of the free-free sky in Sec. 3. We conclude in Sec. 4. Mathematical details of the method can be found in the Appendix.

The physics
The rotation angle ∆ λ of the polarization plane of linearly polarized light due to Faraday rotation can be written in most astrophysical settings as (see e.g.Oppermann et al. 2012Oppermann et al. , 2015) where λ is the wavelength at which the light is observed.
The Faraday depth φ depends on the thermal electron den-  (Oppermann et al. 2015).We show the differences of the three maps in Fig. 2.
sity n th and the line of sight (LOS) component of the magnetic field B LOS .More specifically for a polarized extragalactic source shining through the interstellar medium, the Galactic part of the Faraday depth φ gal can be written as (Oppermann et al. 2012)  where B gal is the galactic magnetic field, e is the electron charge, m e is the electron mass and c the speed of light.Unfortunately, contributions to the observed rotation angle not only come from the Galactic magnetic field but also contain extra-galactic, potentially source internal contribu-tions.Similar to the approach taken by Oppermann et al. (2012), it will be helpful to separate the Faraday depth signal into φ = φ gal + φ etc . (3) The terms φ gal and φ etc describe the aforementioned Galactic and other e.g.extra-galactic, source intrinsic, or ionospheric contributions, respectively.Our goal in this work will be the estimation of φ gal as a correlated field on the sky.This implies the necessity for an at least stochastic estimation of the other contributions.The connection of φ to the observed data d φ is the topic of the next section.

The data
A measurement of the Galactic Faraday depth can be modeled via where we ignore for the moment the potential non-galactic contributions mentioned in Eq. ( 3).The operator R φ is called the response and describes the (noiseless) measurement process translating the signal into data space.In this case, the evaluation of the response at the location of source i is sufficiently described via with e being the radial unit vector on the 2-sphere.The quantity n φ in Eq. ( 4) is the measurement noise, which in the easiest case can be assumed to be independent and Gaussian distributed, according to with zero mean and known noise covariance N φ .The noise will play an important role for our estimation of nongalactic contributions to φ, which will be discussed in more detail in Sec.2.3.2.For the moment we assume N φ to be simply given by the observational measurement uncertainties, therefore being a diagonal matrix in data space.
As we are conducting a Bayesian inference, our main task is the evaluation of the posterior probability distribution.
An integral part of this distribution is the likelihood, which in our case is best described by the following Gaussian: The likelihood will be further determined by the the specific modeling of φ, which we will present in Sec. 4.There we will also specify the different prior terms.
The data set used in this work is exactly the same compilation of Faraday rotation measures as used in Oppermann et al. (2012) and is mostly publicly accessible 2 .This compilation contains rotation measures, corresponding standard deviations and source positions for 41632 sources3 .Similar to Oppermann et al. (2015) we multiply the uncertainties of the Taylor et al. (sub-)data set (Taylor et al. 2009) with 1.22 according to Stil et al. (2011).No additional polarimetric information such as e.g.Faraday spectra was used in this work.The data and the corresponding standard deviations are shown in Figs.3a and 3b.Next, we specify the modeling of φ in Eq. ( 4).We will furthermore explain the treatment of the non-galactic components of d φ .

The sky model
The results of the previous reconstruction (Oppermann et al. 2015) indicates that the Faraday sky has two characteristics that necessarily need to be modeled for an optimal inference.Firstly, the sign of the sky field changes rapidly, which reflects the rather abrupt directional changes of the LOS magnetic field component.Secondly, the absolute value of the signal varies over two orders of magnitude, most notably the Galactic disk will have a much stronger signal than the Galactic poles.For these reasons, we chose to parameterize φ : S 2 → R as a point-wise product between two fields living on the sky, or explicitly where χ : S 2 → R will be called the sign field, as it is able to capture the information on the sign of the Faraday depth, and ρ : S 2 → R will be called the amplitude field, as e ρ is strictly positive after the exponentiation and is able to model the large amplitude variations of the Faraday sky over orders of magnitude.Both fields are unitless.This is a more generic approach as in Oppermann et al. (2012), where instead of the amplitude field a profile function was used to capture the latitudinal dependence of the overall Galactic Faraday dispersion profile.Note that the sign field is in no way constrained to only contain information on the sign of the signal, but also can capture morphological features, leading to a degeneracy between the two fields.This could be broken by either imposing constraints on the correlation structure for (one of) the fields or by introducing new data that informs the algorithm on one of the two fields.However, it turns out that the symmetry is sufficiently broken by the Gaussian process priors we impose on both fields separately in combination with the specific functional form of φ as given by Eq. ( 8).Sign variations can only be captured by the sign field, ensuring that it is structured.The The resolution of the sky pixelisation is lower in this image than in the reconstruction for illustrational purposes.We note that the region corresponding to the terrestrial south pole is only weakly constrained by data, except for the Galactic disk.
overall amplitudes of these fluctuations, which change as a function of position on the sky, are, however, more easily represented by the amplitude field, which therefore preferentially absorbs those.Both fields are assumed to have independent and Gaussian isotropic statistics.Their joint prior probability function is However, in contrast to the distribution in Eq. ( 6), the covariances S χ and S ρ are unknown.If this were not the case, the above model in a Bayesian setting would directly lead to a specific application of the well known (non-linear) Wiener filter (Wiener 1966;Enßlin et al. 2009).This is not the case unfortunately, as in our case we have little a priori knowledge on the correlation structure.Similar inference problems have been solved in the past, e.g. in the previous reconstructions of the Faraday sky by (Oppermann et al. 2012(Oppermann et al. , 2015)), but also e.g. in (Selig et al. 2015;Pumpe et al. 2018;Leike & Enßlin 2019).
All of these inferences relied on a framework laid out by information field theory (IFT) (Enßlin et al. 2009;Enßlin 2018).IFT connects Bayesian statistics with methods from Fig. 4: Hierarchical Bayesian model for the initially revised Faraday sky.We decompose the Faraday depth φ into the fields χ and ρ, which are supposed to capture the sign and the overall amplitude of the signal, respectively (see Eq. ( 8)).The Faraday depth field φ together with the measurement noise n φ determine the observed data d φ via Eq.( 4).The noise n φ of every measurement i is assumed to be drawn from a Gaussian with variance η i σ 2 i , where σ i is the reported uncertainty and η i a unknown uncertainty correction fudge factor.statistical and quantum field theory, joining them into a inference scheme that connects noisy, incomplete data with the underlying continuous field(s).Considering the problem mentioned above, the requirement of simultaneously inferring the map and the correlation structure of a field leads to the critical filter formalism, which was first formulated in IFT by Enßlin & Weig (2010) and Enßlin & Frommert (2011).We will follow the most modern formulation, as outlined in Knollmüller et al. (2017).We will shortly summarize the modeling of the respective covariance matrices in the Appendix.

The noise model
As mentioned before, we also need to find an estimate for the non-galactic contributions to the Faraday rotation.In this part of the inference, we follow the approach in Oppermann et al. (2012) very closely.We will start by inserting Eq. (3) into the measurement equation (4) and define an effective noise term ñφ : We have no reason to drop the assumption of Gaussianity for the new noise term ñφ , but have to adapt the covariance to include the increase in uncertainty.As we have no a priori guesses on the specific systematics of the different sources, we will infer this increase in uncertainty as well.This is implemented by introducing the following model for the new noise covariance N where η φ are the parameters that need to be inferred and σ 2 φ is the reported measurement uncertainty.We will assume the η φ to follow an inverse gamma distribution: The hyper-parameters α φ and β φ need to be specified for the inference.In contrast to Oppermann et al. (2012) we do not set α φ = 1 and β φ such that the prior expectation value of ln(η φ ) is zero, but demand a steeper slope with α φ = 2 and the prior expectation value of η φ to be unity.The reason for this is that we are allowing for more degrees of freedom in our new sky model as compared to Oppermann et al. (2012), as we here have a two dimenensional amplitude field, where as Oppermann et al. (2012) only had a one dimensional profile function that was empirically determined without any prior.The field prior we have to impose on this amplitude field works against the data, increasing the tendency to classify data structures as noise.The larger value of α φ just counter balances this tendency to a degree that our resulting initially revised Faraday map (Fig. 1b) looks similar to that of Oppermann et al. (2015) in Fig. 1c.

The inference
The hierarchical Bayesian model described in the previous sections can be symbolically depicted as a hierarchical tree, of which the lower branches are shown in Fig. 4. The higher branches of the tree contain the hyper-parameters and are discussed in the appendix together with the more explicit likelihood and prior functions.The evaluation of the posterior is a non-trivial task in our case, as the non-linearities in the models lead to highly non-Gaussian distributions for which no analytical description is known to exist.Furthermore, the fact that we aim for an simultaneous inference of the fields determining the Faraday map as well as their correlation structures leads to a strong interdependence between all degrees of freedoms.Together with the high dimensionality of the problem with of the order of 10 6 unknown parameters for the desired resolution of ≈ 10 arcmin this leads to considerable numerical complexity.IFT solves such problems by providing efficient and well tested algorithms for posterior evaluation.In contrast to Oppermann et al. (2012), we employ here the Metric Gaussian Variational Bayes method (MGVI), which is better suited for simultaneous map and correlation structure inference compared to the maximum a posteriori algorithm used in Oppermann et al. (2012), see Knollmüller et al. (2017) for details.We furthermore implement the aforementioned hierarchical tree as a forward model, which aims among other things for a more efficient numerical coupling between different parts of the model and allows the employment of sampling methods (Knollmüller et al. 2017).On a technical side, the inference was conducted using the newest version (version 5) of the NIFTy python package (Selig et al. 2013;Steininger et al. 2017;The NIFTy5 team et al. 2019) for numerical information field theory.NIFTy is specifically designed for IFT applications and provides libraries for implementing the above mentioned models as well as for inferring their parameters and the uncertainties of those.

Comparison to older results
The initially revised Faraday depth map relying on rotation measure (RM) data only and inferred according to the model depicted in Fig. 4 is shown in Fig. 1b, the previous one of Oppermann et al. (2015) in Fig. 1c, with corresponding standard deviations in Figs.5b and 5a, respectively.The difference of both maps are shown in Fig. 2b.A visual com-parison of the two maps reveals satisfying agreement.The difference map, however, shows that there are indeed significant small scaled deviations of the maps in the area of the Galactic disk.These may be a result of the different noise estimation parameters in the updated inference and/or the different treatment of the amplitude field (see the discussion in Sec.2.3.2).One should note that the revised reconstruction was conducted with a doubled resolution, which may also be responsible for some of the very small scaled differences.The uncertainty maps nicely demonstrate the impact of the sky modeling on the inference.Fig. 5a shows that the uncertainty of the previous reconstruction, which exhibits imprints of the longitudinal galactic profile used in the amplitude modeling.Our updated uncertainty map in Fig. 5b shows a decrease in uncertainty in certain regions, especially in the galactic disk.This simply reflects the fact that our updated inference scheme can use the correlation information encoded in the data much more efficiently.
A further quantity which needs to be inspected are the power spectra depicted in Fig. 6, as they have been inferred with different algorithms in the previous and initially revised approaches.We note that we can confirm the slope of the spectrum.The revised map has significantly less power on small scales.Again, the parameter change in the noise estimation may be the reason for this discrepancy, as we know that it can have a considerable impact on the smaller scales, as we change the ability of the algorithm to increase uncertainties to accommodate outliers in the data.With respect to the error estimation, we show in Fig. 7a the comparison of the reported observational noise standard deviation σ and our estimate for the noise σ according to Eq. ( 10).The uncertainties of some sources have been multiplied by factors of up to 1000, indicating strong internal or extra-galactic contributions to these particular RM's.Nonetheless, we confirm or even decrease the uncertainties for most measurements, as shown in Fig. 7b, where we show a normalized density plot in the σ − σ plane.With this we conclude the comparison between the previous and the initially revised results and turn our face to the analysis of the underlying components of our sky model defined in Eq. (8).

The sign and amplitude fields
The component fields χ and ρ are shown in Figs.8b and 9b, respectively.The sign field captured most of the small scale structure, but has more or less completely lost all information on the Galactic disc profile, which was absorbed by the amplitude field, as intended by our modeling.χ shows some similarity to the corresponding signal field inferred by Oppermann et al. (2015), which is shown in Fig. 8a.The exponentiated amplitude field appears to be relatively smooth with few distinguished features apart from the large scaled disk profile.These features, however, show remarkable similarity to galactic free-free emission maps, which was precisely measured by missions investigating the Cosmic Microwave Background (CMB) such as Planck (Planck Collaboration 2016a), as it is a important foreground component in the microwave sky.We show the Planck free-free EM map in Fig. 9a.The color scales were chosen to highlight the structures in the maps of which we think that they originate from the same galactic structures.This apparent resemblances in both maps have motivated us in including  the free-free map in our inference.A investigation of the physical plausibility of this correlation is given in the next section.2015) (dashed black) to the spectrum of the revised Faraday map including the free-free EM data (red) and the results of the initially revised inference in Sec. 2 (blue).Both of the latter inference algorithms also produce uncertainties in the power spectrum, which are depicted as shaded areas in the respective colors.

The physics
In addition to the visual similarities between the free-free map and our inferred amplitude field, we want to further motivate this apparent congruence (and deviations from it) by investigating the physics of the free-free emission and the Faraday rotation a bit more closely.The EM traced by the bremsstrahlung emitted by thermal electrons on protons in the warm interstellar medium can be quantified via (see e.g.Smoot (1998)) where we defined em ff ≡ n2 th as a source term for EM ff .If we turn to the Faraday depth and want to study Eq.
(2) with the intention of an comparison to EM ff , we must quantify the dependence of B LOS on the thermal electron density.In general, we can assume the magnetic flux to be frozen into the ISM (de Gouveia Dal Pino 2006).The exact nature of the B − n th relationship, however, is complicated as it strongly depends on the morphology and dynamics of the plasma under consideration.In the simplest case of an isotropically collapsing structure, flux freezing leads to an |B| ∝ n shows the estimated noise standard deviations compared to the observed ones.This plot is based the same information as Fig. 7a, but shows the underlying probability density function estimated via Gaussian kernels.In both plots, equality between the two quantities is marked by the dashed line.a short time span, the B − n th correlation may be subject to large variability in all regimes of measured electron density (Vázquez-Semadeni 2015).We therefore assume the dependence of the magnetic field strength on n th to follow a power law with an unknown and spatially dependent coefficient |B| ∝ n p B th .We can then rewrite the proportionality of the absolute value of integrand of the Faraday depth in Eq. (2) on n th Here we introduced the angle θ indicating the orientation of the magnetic field with respect to the LOS.Its impact will be discussed later on.If p B ≈ 1 for all locations along a line of sight, the absolute value of the Faraday integrand does strongly depend on the free-free emissivity and the additional electron density dependence has completely canceled.
In the unrealistic case that p B = −1, the absolute value of the Faraday integrand would be completely uncorrelated to the free-free source term.Furthermore as discussed before, in some cases the electron density will not be a good tracer for the magnetic field strength and the LOS projection might partly mask the relation in the observables, which again limits the above relationship.Nonetheless, we think that the above considerations motivate the inclusion of the free-free map as a proxy for the Faraday amplitude field in our inference, as for most realistic cases of the power law index p B the integral over both sides of Eq. ( 14) will lead to a correlation between the free-free and the Faraday sky.On the other hand this analysis also indicates that the effects of ISM dynamics on the magnetic field strength have to be considered for a reliable inference.
A second effect that needs to be taken into account are sign reversals of the magnetic field along the line of sight, which could strongly affect the value of φ in contrast to EM ff .To consider these, we approximate the line of sight integral via where L is the length of the line of sight and . . .L ≡ 1 L LOS dl . . .indicates an averaging process along the line of sight.The last approximation implies statistical independence between the magnetic field strength and θ.This shows that to good approximation, the geometry of the magnetic field can lead to a multiplicative term, which is not constrained by the free-free data alone and has to be considered independently.In our setup, it can only be further determined by discrepancies between the Faraday and free-free data sets.Similar to the first part of this paper, we will first discuss properties of the EM data and then show the modeling of the respective sky maps in Sec.3.3.

The data
In of the true free-free EM with known uncertainties.Unfortunately, upon inspecting the Planck data, we found some pixels to have zero variances.As this is rather surprising given the large degeneracy between the components, we expect this to be a numerical artifact of the sampling procedure or of the marginalization process that results in the uncertainties for specific components.Our algorithm relies on accurate uncertainty measures for all pixels, as in case of a severe underestimation of the errors in some specific area, this region is over-weighted by the algorithm in the inference of the correlation structure, thereby immediately affecting other parts of the sky and subsequently all other maps of the inference.This has to be accounted for if one is not completely sure that the noise values are free of systematics.We will therefore perform a similar noise estimation procedure as for the Faraday data.With this in mind, we write the measurement equation for the free-free data as In case of Planck, the response R ff is just a unit matrix, as the data is already provided as a full sky map on the reso-lution we desire.Here, we decompose the noise n ff into two components namely a observed one obtained by COMMANDER and an unobserved one, which contains systematics.As we do not consider the correlation structure of the noise in for any data set, the procedure of noise estimation is the same from here on as in 2.3.2,apart from different hyper priors.
In this case, we choose α ff = 1 and β ff such that the expectation value over ln(η ff ) is zero.This is Jeffreys prior for the inverse gamma distribution (Oppermann et al. 2012).The lower branches of the free free sky model are shown in Fig. 12.The full extended model for the hierarchical Bayesian model is shown in Fig. 13.Using the considerations above, the likelihood for the free-free sky is again Gaussian in analogy to Eq. ( 7) with the noise covariance N ff defined in analogy to Eq. ( 10).Again, we still need to explain the detailed model for EM ff and its priors.This is presented in the next section.
3.3.The free free model

The sky model
We now turn our head to the modeling of the free-free sky and its connection to the amplitude field of the Faraday sky under the constraint of the physical limitations considered in Sec.3.1.We model EM ff with a log-normal approach, again to enforce positivity and to capture the expected large variability of EM ff on the sky: If the true free-free sky would be a perfect estimator for the amplitude field of the Faraday sky, we could simply replace ρ in Eq. ( 8) with to simultaneously infer both quantities from both data sets.However, as discussed in Section 3.1 a necessary condition for a combination of the information of free-free EM's and Faraday rotation measures is to account for magnetic field reversals, different scalings with the thermal electron density and potential systematics in the data sets.All issues can be modeled by introducing further terms to the model.
We start with our approach to model the differences in the electron density scaling and the geometry factor introduced by the sign reversals.As we can only work with the LOS integrated quantities, the most natural way of inferring the three dimensional effective scaling factor 1+p B 2 in Eq. ( 14) is not possible, as this would necessarily imply setting up the full 3D problem, for which we are lacking the necessary radial information.We therefore choose to introduce two scaling fields γ : S 2 → R and ψ : S 2 → R into the modeling according to Eq. ( 19) which essentially sets up a linear fit for the logarithm of the Faraday amplitude field in terms of the reconstructed logarithmic free-free sky: φ = e γ +ψ χ. ( This model somewhat mimics the approximation derived in Eq. ( 15), with γ approximating the role of the power law coefficients and vice versa e ψ the geometry factor.Unfortunately, we loose the interpretation of the γ field in terms of 1+p B 2 in Eq.( 14) due to the aforementioned lack of consistent 3D modeling.Furthermore the large scale modes of the ψ and γ fields are to a degree degenerated, making a quantitative interpretation of the ψ field very hard, although the smaller scaled morphology may be informative.Altogether, these terms should be able to capture sign reversal effects as well as strong deviations from the proportionality of the magnetic field strength on the thermal electron density presented in Eq. ( 14).We further have to consider the case where the Faraday sky contains amplitude structure not present in the EM ff sky.This is done by expanding the model according to φ = e γ +ψ + e δ χ. (20) Here, the second amplitude field δ : S 2 → R is not connected to the free-free map, giving it the possibility to capture structures that are in the Faraday data but not in the free-free data.All newly introduced fields are assumed to follow Gaussian statistics, so their prior is similar to Eq. ( 9), again with unknown correlation structure.One should note at this point that the obvious a priori degeneracies between the newly introduced helper fields are usually broken during the inference, if sufficiently supported by data.This is possible due to the joint correlation structure inference, which constrains the respective fields just to the degree it is demanded by the inference problem.
Although there is motivation for each term, the above modeling might be viewed as ad hoc, as one may certainly argue for the addition or omission of certain terms or a completely different modeling altogether.There are certain limitations for potential models, as they all have to be able to capture the physical characteristics (say e.g. the sign reversal of the Faraday map), and should be kept as simple as possible.Unfortunately these requirements already allow for a range of different parametrizations.We have varied the model during our analysis by e.g.adding or omitting further terms in Eq. ( 20).We found little variation in the overall morphology of the Faraday sky, apart form the disk, where e.g omission of the γ and ψ field lead to an extremely pronounced features, with RM values reaching regimes a magnitude higher than previous reconstructions.As we think that the above presented physical considerations are too important to be neglected, we deem the resulting maps in these cases as unrealistic.Given that all newly introduced helper fields have captured distinct structures as we will show in Sec.3.4, we view the model in Eq. ( 20) as one of the simplest cases which can capture all relevant effects, but nonetheless we deem the systematic errors higher than the statistical ones in this inference.
To gain some insights on the robustness of the obtained free-free map, we will further perform an inference just for the free free sky with the model in Eq. ( 18), unconstrained by Faraday data.The comparison of the result to the outcome of the joint Faraday and free-free inference will be presented in the results section.We can implement this extended model numerically with similar algorithms as described in Sec.2.4.

The revised Faraday map
In the following we show the results of the model depicted in Fig. 13.We show the mean, the differences to the previous and initially revised reconstructions and the power spectrum in Figs.1a, 2a, 2c and 6, respectively.The inclusion of the free-free sky has lead to a notable difference in our newly revised estimate of the Faraday sky.Driven by the strong disc in the free-free data, the Faraday sky now has a similarly strongly pronounced disc feature as well.Other features mostly farther away from the Galactic plane remained relatively stable, such as the northern arc of the Gum nebula, which is also strongly visible in the freefree data.In Fig. 5c, we show the updated uncertainties.These have again narrowed considerably compared to our initial free-free data free reconstruction as shown in Fig. 5b.The innermost part of the disc, however, is still quite uncertain.Furthermore one should be aware that any model uncertainties are not considered in this plot.The new power spectrum of the Faraday map is shown as the red line in Fig. 6.It is very similar to that of the previous reconstruction, with an notable offset towards smaller scales, implying a slightly steeper power law.

The components
We now show some of the components of the Faraday model in Eq. ( 20), apart from the field in Fig. 10a, which will be more closely discussed in the following section dedicated to the inference results of the free-free sky.
The revised Faraday amplitude is shown in Fig. 9d.The comparison to the initially revised amplitude field in Fig. 9b and the free-free data in Fig. 9a demonstrated the influence of both the Faraday and the free-free data on the new field.The field is enhanced by roughly a factor of 4 compared to the old amplitude field.This may be a result of the symmetry breaking process described in Sec.2.3, as now the amplitude field is partially constrained by the freefree sky.The comparison of the sign fields will demonstrate that this factor was mostly absorbed there.
The second amplitude field δ shown in Fig. 10b, its influence on the total logarithmic Faraday amplitude is demonstrated by Figs.10c and 10d.The field is dominated by two diffuse structures of enhanced δ, centered in the Galactic plane and separated by roughly 180 • longitude.Remembering the discussion in Sec.3.1, the δ field is supposed predominately capture amplitude variations between the free-free and the Faraday sky.Larger values of δ can be caused by systematic errors in either the free-free or the Faraday data and/or a high alignment of the magnetic field with the LOS in this direction.The approximated longitudinal difference of the two δ enhancements of about 180 • as well as their positions seem to correspond to the angular positions of the Orion arm on the sky, which are depicted in the plot as red crosses.In the first quadrant of the Milky Way, at about 60 • longitude, this minor arm is known to extend over several kiloparsec up to the Sagittarius arm, with which it might merge at last (Xu et al. 2009).In the other direction towards the third quadrant of the Galaxy, the structure of the arm is much more complex.Even a bifurcation or a crossing with the Perseus arm seems possible (Vázquez et al. 2008).We have indicated possible continuations of the arm in the plot at 240 • and 260 • longitude, respectively.The δ field seems to reflect not only the position of the arm, but also its morphology, as the structure in the third quadrant is more diffuse, as opposed to the relatively pronounced structure at 60 • longitude.All in all, this would argue for a high magnetic field alignment as a viable possibility for the interpretation of the second amplitude field.
As the Sun is positioned directly in the Orion arm, we expect more coherent magnetic fields along the LOS directed pointed directly to the dominant direction of the arm.Re-  We model EM ff via the exponentiated field (see Eq. ( 18)).
The field EM ff is together with the noise n ff connected to the observed data via Eq.( 16) .
membering the discussion in Sec.3.1, this then may lead to a disproportionate increase in the observed Faraday measures as the free-free map would indicate at these positions.Furthermore the Faraday sky in Fig. 1a indeed shows opposite magnetic signs for those regions, as one would expect for a field traversing the solar location.We show the field in Fig. 11a, which corresponds to our estimate of the enhancement of the Orion Faraday depth, should the above considerations be true and if the χ field really only captured the sign.The (in terms of amplitude contributions) complementary Faraday component is shown in Fig. 11b.The considerable differences to Fig. 1a in the sky regions corresponding to the Orion arm reveal again the influence of the δ field there.Of course we can not rule out any of the other explanations for the observed δ enhancements, the plot in Fig. 11a should therefore be taken with care.
The comparison of the amplitude field in Fig. 10c without the δ contributions with the field alone (see Fig. 10a) reveal the strong rescaling and reweighting of the free-free contribution due to the γ and ψ fields.At last we discuss the sign field χ in Fig. 8c.Little has changed in this component morphologically compared to Fig. 8b.The absolute value of the field has decreased by approximately a factor of 4. This factor was mostly absorbed into the new amplitude field, as discussed before.Apart from that, the field has again fully served its intended purpose, namely capturing the sign of the Faraday sky.
We now discuss the results of the inference of the free-free EM sky.

The revised free-free map
The revised free-free EM sky constrained by both the Planck and the Faraday data is shown in Fig. 9c, the logarithm of it is shown in Fig. 10a.The difference to the map inferred using only Planck data is shown in Fig. 14.There is little difference between the maps, indicating that the systematic differences of the data sets were mostly absorbed by the second amplitude field δ and the rescaling fields γ and .
Notable differences appear only towards the Galactic poles, where Planck has measured little to no free-free component, while considerable Faraday rotation exists.The visual comparison of Fig. 9c to the data shown in Fig. 9a demonstrates good accordance.The main difference is the lack of small structures far away from the disc, which where mostly absorbed into the noise.Closer inspection reveals some small scaled deviations in the disc, which we want to investigate by showing an excerpt of the disc ranging from 0 • to 60 • longitude in Fig. 15.This shows approximately the same region as the first plot of Fig. 7 in Planck Collaboration (2016b), where the Planck team compares their free-free temperature map with independent measurements via the radio recombination line survey (RRL).While our reconstruction reproduced the overall morphology, some peaks in the picture experience deviations around 20 %.In general our results seem to underestimate the amplitudes compared to Planck.It is noteworthy that the comparison between Planck and RRL resulted in a similar statement, as the Planck features were generally much more pronounced.We will refrain from drawing physical conclusions from these discrepancies as these can only be resolved by a detailed discussion of the Planck data analysis, which is beyond the scope of this work.Note again that the Faraday data has little influence on this result.The changes are mainly driven by the noise estimation and the full consideration of the correlation structure of the free-free EM sky.We think that this deems the inclusion of the correlation structure in the inference an important step for any component separation algorithm, as also shown in Knollmüller & Enßlin (2017).

Summary and conclusion
We infer the Galactic Faraday sky map using the same data as Oppermann et al. ( 2012) using an updated algorithm.
Fig. 13: The full hierarchical model excluding the power spectrum hyper priors for the fields.The lowest layer contains the data sets d ff and d φ .These are connected by the equations ( 4) and ( 16) to the sky maps EM ff and φ.The sky maps in turn are connected to the Gaussian fields in the uppermost layer via the respective models defined in Eqs. ( 18) and (20).From there on, the higher branches for the respective correlation structure inference follow.These are not depicted here, but shown and explained in the Appendix.
Fig. 14: Difference between the field defined in Eq. ( 18), inferred only by free-free data and the same field resulting from the joint inference with Faraday data.This involves the separation of the Faraday sky into the point wise product of a sign and amplitude field, which both are supposed to mainly capture the sign and the Galactic profile of the Faraday sky, respectively.In the first part of the paper, we can largely confirm the old results in map and power spectrum.Smaller differences can be explained by the algorithmic advancements.The biggest improvement here is the updated uncertainty map, where we were able to significantly reduce the error bars over the whole sky, especially in some locations in the galactic disc.
In the further proceedings we notice that the amplitude field has notable similarities to the Galactic free-free EM map.In the second part of the paper we therefore incorporate free free data as a potential proxy for the Faraday amplitude field into our inference.This leads to a new Faraday map which now is notably different to previous results.The Galactic disc is much more pronounced and uncertainties are largely reduced.Since a pronounced Galactic is very sensible, as well as required by the data, which now is more satisfied as in any of the other approaches, we regard this map as the most reliable and therefore to be the main result of this work.Moreover, the component fields introduced in our model to resolve discrepancies between the two data sets reveal indicators for more LOS aligned magnetic field structures in the direction of the Orion arm.We also produce a denoised free-free map, which agrees with the Planck data to large extent.
Appendix B: Correlation structure model Again we will not follow the details in Oppermann et al. (2012), but use an updated algorithm for the correlation structure inference.A short discussion on the differences of the two approaches and potential implications on the spectrum follows at the end of this section.We have verified our results via a comparison to the previous map only using the data used in Oppermann et al. (2012) in Sec.2.5.
The assumption of statistical homogeneity and isotropy together with Gaussianity for a field f implies via the Wiener-Khintchin theorem that the respective covariance is diagonal in harmonic space.We will make use of this theorem, as it allows us to describe the correlation structure via a field with at most N pix dimensions, as opposed to the N 2 pix dimensions in the non-harmonic space.The correlation structure of a Gaussian field f is then fully described via its power spectrum defined as where P is a projection operator onto spectral bands of the power spectrum.This operator simply gives us the possibility to reduce the degrees of freedom of the spectrum, it may very well be set to a unit matrix in case one e.g. expects rather sharp features such as characteristic frequencies.
For the inference, we follow the reparametrization trick laid out in Knollmüller et al. (2017).We decompose the sky fields into Y is the spherical harmonic transform.The amplitude fields A f contain the information on the correlation structure on the respective field, while the ξ f are a priori white Gaussian fields with zero mean and unit covariance which contain the information of the specific field configuration of the field f given the correlation structure in A f .Both fields live in harmonic space.One can show that the amplitude fields are related to the covariance via This formulation has immediate consequences for the Hamiltonian, as the prior terms are now independent of the correlation structure as (B.4) using the above equations.For a more general discussion of the advantages of this approach we refer the reader to Knollmüller et al. (2017).Now we are left with the inference of the amplitude fields.Those are modeled via a log-normal field τ f to ensure positivity: It is decomposed into a amplitude field A f containing the correlation structure and a a priori white Gaussian field ξ f connecting this correlation information to the actual spatial representation.A f on its own again depends on hyper priors l f , σ f and another white field ζ f .This tree wraps up Section B and is the core building block of our correlation structure inference.
The inference then follows the approach of Arras et al. (2019), in which τ f is again transformed according to the reparametrisation trick above in to a white field ξ f and hyper parameters σ f and l f which steer the expected variability in the spectrum.In Oppermann et al. (2012) the power spectrum was modeled via a log-normal field together with the product of an inverse gamma and a regularizing smoothness likelihood.The above demonstrated flattening of the hierarchical structure was not performed.We show a comparison in the section 2.5 to demonstrate that the revised algorithm reproduces the inference of the correlation structure in Oppermann et al. (2012), and also comment on the minor differences occurring between the two spectra.The full Hamiltonian for the simultaneous inference of Faraday sky, the free-free sky, the correlation structure of their respective component fields and noise corrections terms is then: where the last term contains the hyper priors for the τ f inference and is described in detail in Arras et al. (2019).

Appendix C: Evaluating the posterior
The Hamiltonian in Eq. (B.6) fully describes our information on the problem.Unfortunately, the non-linearities in this expression imply that this posterior distribution is non-Gaussian, making the evaluation a non-trivial task as a general analytic solution is not known to exist.Therefore, we need to choose an appropriate numerical evaluation scheme.For this end, variational Bayes approaches have proven to combine satisfying accurateness in their uncertainty estimates with speed and stability during the minimization (Knollmüller et al. 2017).In a mathematical language, this implies that we will minimize the Kullback-Leibler Divergence (KL) where P is a Gaussian probability distribution and H is the respective Hamiltonian.... P indicates an averaging procedure of the quantity in the brackets with respect to P.
Minimizing the above expression can be interpreted as a fit of P to P. This approximated posterior is then easily evaluated via sampling due to its Gaussian nature.
Both the Hamiltonian and the KL functional can be conveniently implemented into NIFTy (Selig et al. 2013;Steininger et al. 2017;The NIFTy5 team et al. 2019), which in its newest version then takes care of calculating the correct gradients via auto-differentiation and provides a library of suitable minimization schemes for optimizing Eq. (C.1).
Convergence is assessed via predefined convergence criteria and the stability is tested via starting the inference from different (both deliberately and randomly chosen) initial conditions.For further details on the here used MGVI technique see Knollmüller & Enßlin (2019).

Fig. 1 :
Fig. 1: Reconstructions of the Faraday sky. Figure (a) shows our main result.This map uses the free-free EM map as a proxy for the Faraday amplitude.The model furthermore contains fields that translate the free-free map into a Faraday amplitude, thereby balancing between effects that support and such that disturb a direct relation of these two quantities.Figure (b) contains the posterior mean of the initially revised reconstruction.Here we did not use the additional freefree data.Figure (c) shows the posterior mean of the previous reconstruction(Oppermann et al. 2015).We show the differences of the three maps in Fig.2.

Fig. 2 :
Fig. 2: Differences of the reconstructions of the Faraday sky shown in Fig. 1. Figure (a) shows the difference between the previous map and our revised reconstruction including the free free data.Figure (b) shows the difference between the previous map and our initially revised reconstruction.Figure (c) finally shows the difference between the Faraday sky resulting from our revised reconstruction including the free free data and the initially revised sky.
Fig. 2: Differences of the reconstructions of the Faraday sky shown in Fig. 1. Figure (a) shows the difference between the previous map and our revised reconstruction including the free free data.Figure (b) shows the difference between the previous map and our initially revised reconstruction.Figure (c) finally shows the difference between the Faraday sky resulting from our revised reconstruction including the free free data and the initially revised sky.

Fig. 3 :
Fig. 3: Sky map of the RM data set (figure (a)) used in this work and the corresponding uncertainties in figure (b).The resolution of the sky pixelisation is lower in this image than in the reconstruction for illustrational purposes.We note that the region corresponding to the terrestrial south pole is only weakly constrained by data, except for the Galactic disk.

Fig. 5 :
Fig. 5: Uncertainties of the different reconstructions.Figure (a) is the uncertainty corresponding to the previous reconstruction of Oppermann et al. (2015) shown in Fig. 1c. Figure (b) shows the uncertainty corresponding to the initially revised reconstruction shown in Fig. 1b.Figure (c) finally shows the uncertainty corresponding to the revised reconstruction including the free free data shown in Fig. 1a.
Fig. 5: Uncertainties of the different reconstructions.Figure (a) is the uncertainty corresponding to the previous reconstruction of Oppermann et al. (2015) shown in Fig. 1c. Figure (b) shows the uncertainty corresponding to the initially revised reconstruction shown in Fig. 1b.Figure (c) finally shows the uncertainty corresponding to the revised reconstruction including the free free data shown in Fig. 1a.

Fig. 6 :
Fig. 6: Comparison of the power spectrum inferred by Oppermann et al. (2015) (dashed black) to the spectrum of the revised Faraday map including the free-free EM data (red) and the results of the initially revised inference in Sec. 2 (blue).Both of the latter inference algorithms also produce uncertainties in the power spectrum, which are depicted as shaded areas in the respective colors.

Fig. 7 :
Fig. 7: Plots describing the impact of the noise estimation procedure.Figure (a) shows the estimated noise standard deviations compared to the observed ones.Note that the visual impression is dominated by the outliers.Figure (b)shows the estimated noise standard deviations compared to the observed ones.This plot is based the same information as Fig.7a, but shows the underlying probability density function estimated via Gaussian kernels.In both plots, equality between the two quantities is marked by the dashed line.

Fig. 8 :
Fig. 8: Sign fields of the different reconstructions.Figure (a) shows the field s of the previous reconstruction by Oppermann et al. (2015).Figure (b) shows the field χ of the initially revised reconstruction.Figure (c) finally shows the field χ of the revised reconstruction including the free free data.Note that the difference in the scales of these fields is compensated by corresponding differences in the respective Faraday amplitude fields.

Fig. 9 :
Fig. 9: Free-free data, reconstructions and Faraday amplitude fields.Figure (a) shows the Galactic free-free EM map as obtained by the Planck Collaboration (Planck Collaboration 2016a).Figure (b) shows the exponentiated amplitude field ρ of the initially revised reconstruction.This field is not constrained by free-free data.Figure (c) shows the exponentiated field.This field is part of the revised Faraday reconstruction as well as the denoised free-free sky.A logarithmic version of this plot is shown in Fig. 10a. Figure (d) shows the Faraday amplitude field of the revised reconstruction including free-free data.This field is part of the result of the model in Eq. (20).A logarithmic version of this plot is shown in Fig. 10d.

Fig. 10 :
Fig. 10: Logarithmic amplitude fields.Figure (a) shows the field defined in Eq. (20), constraint by free-free and Faraday data.This is also the reconstructed log free free map. Figure (b) shows the second amplitude field δ defined in Eq. (20).The red crosses indicate the approximate angular positions of the Orion arm, as given by Vázquez et al. (2008) and Xu et al. (2009).Figure (c) shows the logarithmic amplitude of the Faraday sky as defined in Eq. (19) without the additional δ contribution, while this contribution is included in figure (d).
Fig. 10: Logarithmic amplitude fields.Figure (a) shows the field defined in Eq. (20), constraint by free-free and Faraday data.This is also the reconstructed log free free map. Figure (b) shows the second amplitude field δ defined in Eq. (20).The red crosses indicate the approximate angular positions of the Orion arm, as given by Vázquez et al. (2008) and Xu et al. (2009).Figure (c) shows the logarithmic amplitude of the Faraday sky as defined in Eq. (19) without the additional δ contribution, while this contribution is included in figure (d).

Fig. 12 :
Fig.12: Hierarchical tree model for the free free EM sky.We model EM ff via the exponentiated field (see Eq. (18)).The field EM ff is together with the noise n ff connected to the observed data via Eq.(16) .

Fig. 15 :
Fig. 15: Excerpt of the free-free reconstruction in the first slice, followed by the same excerpt in the Planck data.The difference of the two slices is shown in the last row.
Fig. B.1: Hierarchical tree model for a generic Gaussian field f living on the sky.It is decomposed into a amplitude field A f containing the correlation structure and a a priori white Gaussian field ξ f connecting this correlation information to the actual spatial representation.A f on its own again depends on hyper priors l f , σ f and another white field ζ f .This tree wraps up Section B and is the core building block of our correlation structure inference.