© ESO, 2016

1. Introduction

Regions of the interstellar medium that are (partly) ionized play an important role in a number of effects such as pulse dispersion and scattering, and Faraday rotation. Additionally, ionized parts of the interstellar medium emit radiation through free-free emission and H_α emission. The magnitude of these effects depends on the distribution of free electrons, the free electron density. It is therefore of great interest to model or reconstruct the free electron density as accurately as possible.

Reconstruction and modeling of the Milky Way has been an ongoing topic of research for many years. The free electron density has been modeled by Taylor & Cordes (1993), Cordes & Lazio (2002), and Gaensler et al. (2008), among others. For a comparison and discussion of various existing models see Schnitzeler (2012) and for a review of the mapping of HI regions see Kalberla & Kerp (2009). The interstellar magnetic field has been modeled by Sun et al. (2008), Sun & Reich (2010) and Jansson & Farrar (2012b,a). The dust distribution has been modeled for example by Berry et al. (2012), and even nonparametric tomography has been performed by Lallement et al. (2014) and Sale & Magorrian (2014).

We plan to use the dispersion measures (DM) of pulsar signals together with accurate pulsar distances to map the distribution of ionized gas in the Milky Way. The dispersion measure is defined as the line-of-sight integral over the free electron density between the observer and the pulsar, $\begin{array}{ccc}\mathit{DM}\mathrm{=}\underset{\mathrm{pulsar}}{\overset{\mathrm{observer}}{\mathrm{\int }}}\mathrm{d}\mathit{r}\hspace{0.17em}{\mathit{n}}_{\mathrm{e}}\mathit{,}& & \end{array}$(1)where n_e is the three-dimensional free electron density. DM can be estimated by measuring the arrival time of a pulse at different frequencies, since the time delay is proportional to DM/ν². While there is a vast number of known dispersion measures, very few of them are complemented by an independent distance estimate. The NE2001 model by Cordes & Lazio (2002) is currently the most popular model for the free electron density of the Milky Way. It uses 1143 DM measurements, 112 of which were complemented by distance estimates of varying quality. Additionally, it uses 269 pulsar scattering measurements, which only provide very indirect distance constraints.

We here perform nonparametric tomography of a simulation of the Galactic free electron density from pulsar dispersion measures complemented by independent distance estimates. By nonparametric tomography we mean a reconstruction with a virtually infinite¹ number of degrees of freedom using a close to minimal set of prior assumptions that only resolves structures that are supported by the data. Our assumptions are that the electron density is positive and spatially correlated and that the large-scale electron distribution only shows a variation with distance from the Galactic center and height above the Galactic plane. Both the correlation structure and the scaling behavior have to be inferred from the data. As a consequence, our reconstruction is focused on the large (kpc) scales of the Galactic free electron density. Small-scale structures such as HII regions and supernova remnants as well as spiral arms are only recovered if they are sufficiently deeply probed and constrained by the data.

Our tomography algorithm is derived from first principles in a Bayesian setting. This has the advantage that all assumptions are clearly stated as priors. Additionally, it allows us to provide uncertainty maps of our reconstructions, which are important for any subsequent scientific analysis.

To obtain a meaningful map with minimal assumptions, a data set of high quality is needed, of course. Currently, there are around 100 pulsars known with reliable (independent) distance estimates. This only allows for a nonparametric reconstruction of the largest features in the Milky Way. New measurements with the Very Long Baseline Array will soon double the number of pulsars with accurate distances (see Deller et al. 2011). However, with the planned Square Kilometer Array radio interferometer (SKA) the number of pulsars with parallax distance estimates might increase to around 10 000 (see Smits et al. 2011). In this paper we therefore investigate the feasibility of nonparametric tomography of the free electron density and demonstrate the performance of our algorithm by applying it to mock data sets similar to what the SKA might deliver. To this end, we create four Galaxy models from the NE2001 code by Cordes & Lazio (2002) with varying degrees of fluctuations and contrast as well as observational mock data sets for up to 10 000 pulsars with distance estimates of varying quality and apply our algorithm to these data sets.

The remainder of this paper is structured as follows: first, we derive our tomography algorithm in Sect. 2, explaining our notation, our underlying assumptions, and all probability density functions (PDFs) involved. Second, we explain our Galaxy models and mock observations in detail in Sect. 3. In Sect. 4 we compare the electron density distributions reconstructed from mock observations with those from the Galaxy models used to produce them. We summarize our discussion in Sect. 5.

2. Reconstruction algorithm

The reconstruction algorithm applied in this work was derived within the framework of information field theory introduced by Enßlin et al. (2009). We also follow for most parts the notation used by them. To reconstruct the Galactic free electron density from pulsar dispersion measurements, we used a filter formalism very similar to the one presented by Junklewitz et al. (2016), which in turn is based on the critical filter formalism developed by Enßlin & Weig (2010), Enßlin & Frommert (2011), and refined by Oppermann et al. (2013).

2.1. Signal model

In the inference formalism we aim to reconstruct the free electron density field ρ, a three-dimensional scalar field. We assume it is related to the observed dispersion measure data DM by a linear measurement equation subject to additive and signal-independent measurement noise, $\begin{array}{ccc}\mathit{DM}\mathrm{=}\mathit{R\rho }\mathrm{+}\mathit{n,}& & \end{array}$(2)where n is the measurement noise and Rρ is the application of the linear response operator R on the field ρ, $\begin{array}{ccc}{\left(\mathit{R\rho }\right)}_{\mathit{i}}\mathrm{\equiv }\mathrm{\int }{\mathrm{d}}^{\mathrm{3}}\mathit{x}\mathit{R}\mathrm{\left(}\mathit{i,}{x}\mathrm{\right)}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}{x}\mathrm{\right)}\mathit{.}& & \end{array}$(3)The response operator R describes line-of-sight integrals through the density. It can be defined as $\begin{array}{ccc}\mathit{R}\mathrm{\left(}\mathit{i,}{x}\mathrm{\right)}\mathrm{=}\underset{\mathrm{0}}{\overset{\mathrm{|}{{d}}_{\mathit{i}}\mathrm{|}}{\mathrm{\int }}}\mathrm{d}\mathit{r}\mathit{\delta }\left({x}\mathrm{-}\mathit{r}{{d̂}}_{\mathit{i}}\right)\mathit{,}& & \end{array}$(4)where d_i is the position of pulsar i in a coordinate system centered on the Sun, and δ(·) is the three-dimensional Dirac delta-distribution and .

Formally, the free electron density is a continuous field. In practice, we reconstructed a discretized version of this field, for instance, a three-dimensional map with some pixel size. The discretized density field can be thought of as a vector of dimension N_pix with each component containing the field value in a specific pixel. The dispersion data DM and the noise n can be regarded as vectors of dimension N_data, where each component of DM contains a specific measurement result and the corresponding component of n the noise contribution to it. Thus, the response operator becomes a matrix with N_pix columns and N_data rows.

We parametrized the density as $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\mathrm{\Delta }\mathrm{\left(}{x}\mathrm{\right)}\begin{array}{c}\mathrm{˜}\\ \mathit{\rho }\end{array}\mathrm{\left(}{x}\mathrm{\right)}\mathit{,}& & \end{array}$(5)where Δ is the Galactic profile field that describes the disk shape of the Milky Way. All deviations from the Galactic profile are described by $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\rho ,}\end{array}$ for which we assumed no distinguished direction or position a priori. To ensure positivity of the density, these fields were in turn parametrized as $\begin{array}{ccc}\mathrm{\Delta }\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}& \mathrm{=}& \exp \left(\mathit{\alpha }\left(\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}\right)\mathrm{+}\mathit{\beta }\left(\mathrm{\left|}\mathit{z}\mathrm{\right|}\right)\right)\mathit{,}\\ \begin{array}{c}\mathrm{˜}\\ \mathit{\rho }\end{array}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}& \mathrm{=}& \end{array}$(6)Thus, Δ can only represent the vertical and radial scaling behavior of the density and has the degrees of freedom of two one-dimensional functions. On the other hand, $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\rho }\end{array}$ retains all degrees of freedom of a three-dimensional field and can represent arbitrary structures. Neither Δ nor $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\rho }\end{array}$ are known a priori and will be inferred from the data.

We summarize our modeling in Fig. 1. The logarithmic density ρ is parametrized by three additive components, one three-dimensional field and two one-dimensional fields. As we outline in Sects. 2.2.2 and 2.3, all three fields are assumed to follow Gaussian statistics a priori. For the 1D fields a specific correlation structure was assumed, while the correlation structure of the three-dimensional field was unknown, but assumed to be homogeneous and isotropic. Therefore, our modeling prefers smooth structures, fluctuations that scale with the density, and exponential scaling in radial and vertical directions. Of course, this is a strong simplification of the Galaxy, where the behavior of the fluctuations can depend on the phase of the interstellar medium or the position within the Galaxy, for instance. However, all of these properties can be recovered if the data demand it, since all degrees of freedom are retained. They are just not part of the prior knowledge entering our inference.

2.2. Necessary probability density functions

Our goal is to derive an algorithm that yields an estimate of the logarithm of the Galactic free electron density. Hence, we constructed the posterior PDF , which is the PDF for the signal given the data set { DM,d_obs(erved) }, using Bayes’ theorem, $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\mathit{s}\mathrm{\right|}\mathrm{data}\mathrm{)}\mathrm{=}\frac{\mathrm{𝒫}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}{{d}}_{\mathrm{obs}}\mathrm{)}}{\mathrm{𝒫}\mathrm{\left(}\mathit{DM}\mathrm{\right|}{{d}}_{\mathrm{obs}}\mathrm{)}}\mathrm{=}\frac{\mathrm{𝒫}\mathrm{\left(}\mathit{s}\mathrm{\right|}{{d}}_{\mathrm{obs}}\mathrm{\left)}\mathrm{𝒫}\mathrm{\right(}\mathit{DM}\mathrm{\left|}\mathit{s,}{{d}}_{\mathrm{obs}}\mathrm{\right)}}{\mathrm{𝒫}\mathrm{\left(}\mathit{DM}\mathrm{\right|}{{d}}_{\mathrm{obs}}\mathrm{)}}\mathrm{·}& & \end{array}$(7)On the right-hand side, we have three PDFs: the prior , the likelihood , and the evidence . The evidence is independent of the signal and therefore automatically determined by the normalization of the posterior. The prior and the likelihood are addressed in the following sections. For notational convenience we drop the dependence on the observed pulsar positions d_obs throughout the rest of this paper.

Throughout this section we assume the Galactic profile field to be given. We adress its inference in Sect. 2.3.

2.2.1. Likelihood

The likelihood is the PDF that an observation yields dispersion measures DM assuming a specific realization of the underlying signal field s. If both the noise n and the pulsar distances d_i ≡ | d_i | were known, the relation between the dispersion measure data and signal would be deterministic, $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\mathit{DM}\mathrm{\right|}\mathit{s,n,d}\mathrm{)}\mathrm{=}\mathit{\delta }\mathrm{(}\mathit{DM}\mathrm{-}\mathit{R\rho }\mathrm{-}\mathit{n}\mathrm{)}\mathit{,}& & \end{array}$(8)with ρ(x) = Δ(x)e^s(x). We do not know the realization of the noise, nor do we aim to reconstruct it. It was assumed to follow Gaussian statistics with zero mean and known covariance structure², $\begin{array}{ccc}{\mathrm{⟨}{\mathit{n}}_{\mathit{i}}{\mathit{n}}_{\mathit{j}}\mathrm{⟩}}_{\mathrm{𝒫}\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{=}{\mathit{N}}_{\mathit{ij}}\mathrm{=}{\mathit{\delta }}_{\mathit{ij}}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}\mathit{,}& & \end{array}$(9)where σ_i is the root mean square error of the observation i, and we assumed independent measurements. Distance information is usually given in the form of parallaxes from which distance estimates can be derived. As all observables, these are subject to uncertainties, which is why the information about the distances of the pulsars is described by a PDF³, , which can be non-Gaussian. Since we are doing inference on s, we need the noise and distance⁴ marginalized likelihood $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\mathit{DM}\mathrm{\right|}\mathit{s}\mathrm{)}\mathrm{=}\mathrm{\int }\mathrm{𝒟}\mathit{n}\mathrm{𝒟d}\mathrm{𝒫}\mathrm{(}\mathit{DM}\mathrm{\left|}\mathit{s,n,d}\mathrm{\right)}\mathrm{𝒫}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{𝒫}\mathrm{\left(}\mathit{d}\mathrm{\right)}\mathit{,}& & \end{array}$(10)where we assumed n and d to be independent of s and each other. The symbols and denote integration over the full configuration space of n and d, that is, the space of all possible configurations ().

Integration over n in Eq. (10)is trivial and yields $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\mathit{DM}\mathrm{\right|}\mathit{s}\mathrm{)}\mathrm{=}\mathrm{\int }\mathrm{𝒟d}\mathrm{𝒢}\mathrm{(}\mathit{DM}\mathrm{-}\mathit{R\rho ,N}\mathrm{\left)}\mathrm{𝒫}\mathrm{\right(}\mathit{d}\mathrm{)}\mathit{,}& & \end{array}$(11)where indicates a Gaussian PDF, $\mathrm{𝒢}\mathrm{\left(}\mathit{x,X}\mathrm{\right)}\mathrm{:}\mathrm{=}\mathrm{|}\mathrm{2}\mathit{\pi X}{\mathrm{|}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{†}{\mathit{X}}^{-1}\mathit{x}}$. Integration over d, however, cannot be done analytically, but the marginalized likelihood can be approximated by a Gaussian characterized by its first two moments in DM. The first moment is $\begin{array}{ccc}{⟨\mathit{DM}⟩}_{\mathrm{𝒫}\mathrm{\left(}\mathit{DM}\mathrm{\right|}\mathit{s}\mathrm{)}}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{\rho ,}& & \end{array}$(12)with $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}_{\mathit{i}}\end{array}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}{⟨{\mathit{R}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}⟩}_{\mathrm{𝒫}\mathrm{\left(}\mathit{d}\mathrm{\right)}}\mathrm{=}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\mathrm{\int }}}\mathrm{d}\mathit{r}\mathit{\delta }\left({x}\mathrm{-}\mathit{r}{{d̂}}_{\mathit{i}}\right)\hspace{0.17em}\mathit{P}\mathrm{[}{\mathit{d}}_{\mathit{i}}\mathit{>}\mathit{r}\mathrm{]}\mathit{,}& & \end{array}$(13)where P [ d_i>r ] is the probability that the pulsar distance d_i is larger than r. The second moment is $\begin{array}{ccc}{⟨\mathit{DM}\hspace{0.17em}\mathit{D}{\mathit{M}}^{†}⟩}_{\mathrm{𝒫}\mathrm{\left(}\mathit{DM}\mathrm{\right|}\mathit{s}\mathrm{)}}\mathrm{=}\mathit{N}\mathrm{+}{⟨\left(\mathit{R\rho }\right){\left(\mathit{R\rho }\right)}^{†}⟩}_{\mathrm{𝒫}\mathrm{\left(}\mathit{d}\mathrm{\right)}}\mathit{.}& & \end{array}$(14)For non-diagonal elements the second term on the right-hand side decouples, $\begin{array}{ccc}{\mathrm{⟨}{\left(\mathit{R\rho }\right)}_{\mathit{i}}{\left(\mathit{R\rho }\right)}_{\mathit{j}}\mathrm{⟩}}_{\mathrm{𝒫}\mathrm{\left(}\mathit{d}\mathrm{\right)}}& \mathrm{=}& {\mathrm{⟨}{\left(\mathit{R\rho }\right)}_{\mathit{i}}\mathrm{⟩}}_{\mathrm{𝒫}\mathrm{\left(}\mathit{d}\mathrm{\right)}}{\mathrm{⟨}{\left(\mathit{R\rho }\right)}_{\mathit{j}}\mathrm{⟩}}_{\mathrm{𝒫}\mathrm{\left(}\mathit{d}\mathrm{\right)}}\\ & \mathrm{=}& {\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}{\mathit{\rho }}^{\mathrm{)}}}_{\mathit{i}}{\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}{\mathit{\rho }}^{\mathrm{)}}}_{\mathit{j}} \mathrm{for} \mathit{i}\ne \mathit{j}\mathit{.}\end{array}$(15)Diagonal elements yield $\begin{array}{ccc}{⟨{\left(\mathit{R\rho }\right)}_{\mathit{i}}{\left(\mathit{R\rho }\right)}_{\mathit{i}}⟩}_{\mathrm{𝒫}\mathrm{\left(}\mathit{d}\mathrm{\right)}}& \mathrm{=}& \underset{}{{\mathrm{\int }}^{R^{\mathrm{3}}}}{\mathrm{d}}^{\mathrm{3}}\mathit{x}\underset{}{{\mathrm{\int }}^{R^{\mathrm{3}}}}{\mathrm{d}}^{\mathrm{3}}\mathit{y}\mathit{\rho }\mathrm{\left(}{x}\mathrm{\right)}\mathit{\rho }\mathrm{\left(}{y}\mathrm{\right)}\\ & & \mathrm{\times }{⟨{\mathit{R}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}{\mathit{R}}_{\mathit{i}}\mathrm{\left(}{y}\mathrm{\right)}⟩}_{\mathrm{𝒫}\mathrm{\left(}\mathit{d}\mathrm{\right)}}\mathit{,}\end{array}$(16)with $\begin{array}{ccc}{⟨{\mathit{R}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}{\mathit{R}}_{\mathit{i}}\mathrm{\left(}{y}\mathrm{\right)}⟩}_{\mathrm{𝒫}\mathrm{\left(}\mathit{d}\mathrm{\right)}}& \mathrm{=}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\mathrm{\int }}}\mathrm{d}\mathit{r}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\mathrm{\int }}}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathit{\delta }\mathrm{(}{x}\mathrm{-}\mathit{r}{{d̂}}_{\mathit{i}}\mathrm{)}\mathit{\delta }\mathrm{(}{y}\mathrm{-}{\mathit{r}}^{\mathrm{\prime }}{{d̂}}_{\mathit{i}}\mathrm{)}\\ & & \mathrm{\times }\mathit{P}\mathrm{[}{\mathit{d}}_{\mathit{i}}\mathit{>}\begin{array}{c}\max \end{array}\mathrm{(}\mathit{r,}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right]}\mathit{.}\end{array}$(17)Using these first two moments, we can approximate⁵ the likelihood by a Gaussian $\mathrm{𝒢}\mathrm{(}\mathit{DM}\mathrm{-}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathit{\rho ,}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{N}\end{array}\mathrm{)}$ with $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ {\mathit{N}}_{\mathit{ii}}\end{array}\mathrm{=}{\mathit{N}}_{\mathit{ii}}\mathrm{+}{\mathit{\rho }}^{†}{\mathit{F}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathit{\rho ,}& & \end{array}$(18)where⁶$\begin{array}{ccc}{\mathit{F}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{\left(}{x}\mathit{,}{y}\mathrm{\right)}& \mathrm{:}\mathrm{=}& {⟨{\mathit{R}}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}{\mathit{R}}_{\mathit{i}}\mathrm{\left(}{y}\mathrm{\right)}⟩}_{\mathrm{𝒫}\mathrm{\left(}{\mathit{d}}_{\mathit{i}}\mathrm{\right)}}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}_{\mathit{i}}\end{array}\mathrm{\left(}{x}\mathrm{\right)}\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}_{\mathit{i}}\end{array}\mathrm{\left(}{y}\mathrm{\right)}\\ & \mathrm{=}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\mathrm{\int }}}\mathrm{d}\mathit{r}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\mathrm{\int }}}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}\mathit{\delta }\mathrm{(}{x}\mathrm{-}\mathit{r}{{d̂}}_{\mathit{i}}\mathrm{)}\mathit{\delta }\mathrm{(}{y}\mathrm{-}{\mathit{r}}^{\mathrm{\prime }}{{d̂}}_{\mathit{i}}\mathrm{)}\\ & & \mathrm{\times }\mathit{P}\mathrm{[}{\mathit{d}}_{\mathit{i}}\mathit{>}\begin{array}{c}\max \end{array}\mathrm{(}\mathit{r,}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right]}\mathit{P}\mathrm{[}{\mathit{d}}_{\mathit{i}}\mathit{<}\begin{array}{c}\min \end{array}\mathrm{(}\mathit{r,}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right]}\mathit{.}\end{array}$(19)The noise covariance matrix of this effective likelihood is signal dependent, which increases the complexity of the reconstruction problem. Therefore, we approximated the density in Eq. (18)by its posterior mean, $\begin{array}{ccc}{\mathit{\rho }}^{†}{\mathit{F}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathit{\rho }\mathrm{=}\mathrm{tr}\mathrm{(}\mathit{\rho }{\mathit{\rho }}^{†}{\mathit{F}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{)}\mathrm{\approx }\mathrm{tr}\mathrm{(}{⟨\mathit{\rho }⟩}_{\mathrm{𝒫}\mathrm{\left(}\mathit{\rho }\mathrm{\right|}\mathit{DM}\mathrm{)}}⟨\mathit{\rho }⟩\begin{array}{c}†\\ \mathrm{𝒫}\mathrm{\left(}\mathit{\rho }\mathrm{\right|}\mathit{DM}\mathrm{)}\end{array}{\mathit{F}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{)}\mathit{.}& & \end{array}$(20)Since depends on $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{N,}\end{array}$ this yields a set of equations that need to be solved self-consistently (see Sect. 2.4).

2.2.2. Priors

The signal field s is unknown a priori, but we assumed that it has some correlation structure. We describe this correlation structure by moments up to second order in s. The principle of maximum entropy therefore requires that our prior probability distribution has a Gaussian form, $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\mathit{s}\mathrm{\right|}\mathit{S}\mathrm{)}\mathrm{=}\mathrm{𝒢}\mathrm{(}\mathit{s,S}\mathrm{)}\mathrm{:}\mathrm{=}{\left|\mathrm{2}\mathit{\pi S}\right|}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{s}}^{†}{\mathit{S}}^{-1}\mathit{s}\right)\mathit{,}& & \end{array}$(21)with some unknown correlation structure, $\begin{array}{ccc}\mathit{S}\mathrm{\left(}{x}\mathit{,}{y}\mathrm{\right)}\mathrm{=}{⟨\mathit{s}\mathrm{\left(}{x}\mathrm{\right)}\mathit{s}\mathrm{\left(}{y}\mathrm{\right)}⟩}_{\mathrm{𝒫}\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathit{.}& & \end{array}$(22)The first moment of s is set to zero because it can be absorbed into Δ(x). This means that the a priori mean of s is contained in Δ(x).

A priori, our algorithm has no preferred direction or position for s. This reduces the number of degrees of freedom of the correlation structure S. It is fully described by a power spectrum p(k), $\begin{array}{ccc}\mathit{S}\mathrm{\left(}{x}\mathit{,}{y}\mathrm{\right)}\mathrm{=}\sum _{\mathit{k}}\hspace{0.17em}{\mathit{S}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}{x}\mathit{,}{y}\mathrm{\right)}\hspace{0.17em}\mathit{p}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}& & \end{array}$(23)where S^(k) is the projection operator onto the spectral band k, with its Fourier transform defined as $\begin{array}{ccc}{\mathit{S}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}{q}\mathit{,}{{q}}^{{\prime }}\mathrm{\right)}\mathrm{=}\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}\mathit{\delta }\mathrm{\right(}{q}\mathrm{-}{{q}}^{{\prime }}\mathrm{)}{1}_{\mathit{k}}\left(\mathrm{\left|}{q}\mathrm{\right|}\right)\mathit{,}& & \end{array}$(24)with $\begin{array}{ccc}{1}_{\mathit{k}}\left(\mathrm{\left|}{q}\mathrm{\right|}\right)\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(25)The power spectrum p(k), however, is still unknown. The prior for the power spectrum is constructed out of two parts: First, an inverse gamma distribution ℐ(p(k);α_k,q_k) for each k-bin (see Appendix A), which is a conjugate prior for a Gaussian PDF, second a Gaussian cost-function that punishes deviations from power-law spectra (see Oppermann et al. 2013), $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{\propto }\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\underset{\mathit{k}}{\mathrm{􏽙}}\mathrm{ℐ}\mathrm{\left(}\mathit{p}\mathrm{\right(}\mathit{k}\mathrm{\left)}\mathrm{;}{\mathit{\alpha }}_{\mathit{k}}\mathit{,}{\mathit{q}}_{\mathit{k}}\mathrm{\right)}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\left(}\log \mathit{p}{\mathrm{)}}^{†}\mathit{T}\mathrm{\right(}\log \mathit{p}\mathrm{)}\right)\mathit{.}& & \end{array}$(26)T is an operator that fulfills $\begin{array}{ccc}\mathrm{\left(}\log \mathit{p}{\mathrm{)}}^{†}\mathit{T}\mathrm{\right(}\log \mathit{p}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{p}}^{\mathrm{2}}}\mathrm{\int }\mathrm{d}\mathrm{(}\log \mathit{k}\mathrm{)}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\log \mathit{p}\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathit{\partial }\mathrm{(}\log \mathit{k}{\mathrm{)}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\mathit{,}& & \end{array}$(27)and σ_p is a parameter that dictates how smooth the power spectrum is expected to be. In our paper log  refers to the natural logarithm. We explain our choice of the parameters α_k, q_k, and σ_p in Appendix A.

2.2.3. Power spectrum posterior

With the signal and power spectrum priors, we can derive a posterior for the power spectrum, $\begin{array}{ccc}& & \mathrm{𝒫}\mathrm{\left(}\mathit{p}\mathrm{\right|}\mathit{DM}\mathrm{)}\mathrm{\propto }\mathrm{\int }\mathrm{𝒟}\mathit{s}\mathrm{𝒫}\mathrm{(}\mathit{DM}\mathrm{\left|}\mathit{s,p}\mathrm{\right)}\hspace{0.17em}\mathrm{𝒫}\mathrm{\left(}\mathit{s}\mathrm{\right|}\mathit{p}\mathrm{\left)}\hspace{0.17em}\mathrm{𝒫}\mathrm{\right(}\mathit{p}\mathrm{)}\\ & & \mathrm{=}\mathrm{\int }\mathrm{𝒟}\mathit{s}\mathrm{𝒫}\mathrm{\left(}\mathit{DM}\mathrm{\right|}\mathit{s}\mathrm{\left)}\hspace{0.17em}\mathrm{𝒢}\mathrm{\right(}\mathit{s,S}\mathrm{\left)}\hspace{0.17em}\mathrm{𝒫}\mathrm{\right(}\mathit{p}\mathrm{)}\mathit{.}\end{array}$(28)We calculated the integral using a saddle point approximation up to second order around the maximum for the s-dependent part, $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\mathit{DM}\mathrm{\right|}\mathit{s}\mathrm{\left)}\hspace{0.17em}\mathrm{𝒢}\mathrm{\right(}\mathit{s,S}\mathrm{)}& & \mathrm{\approx }\mathrm{𝒫}\mathrm{\left(}\mathit{DM}\mathrm{\right|}\mathit{m}\mathrm{\left)}\hspace{0.17em}\mathrm{𝒢}\mathrm{\right(}\mathit{m,S}\mathrm{)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{s}\mathrm{-}\mathit{m}{\mathrm{)}}^{†}{\mathit{D}}^{-1}\mathrm{(}\mathit{s}\mathrm{-}\mathit{m}\mathrm{)}}\\ & & \mathrm{\propto }\mathrm{𝒢}\mathrm{\left(}\mathit{m,S}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{s}\mathrm{-}\mathit{m}{\mathrm{)}}^{†}{\mathit{D}}^{-1}\mathrm{(}\mathit{s}\mathrm{-}\mathit{m}\mathrm{)}}\mathit{,}\end{array}$(29)where m and D are defined as m^(s) and D^(s) in Sect. 2.2.4 and only s − and p-dependent factors are kept after the proportionality sign. With this approximation we arrive at $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\mathit{p}\mathrm{\right|}\mathit{DM}\mathrm{)}\mathrm{\propto }{\left|\mathrm{2}\mathit{\pi D}\right|}^{\frac{\mathrm{1}}{\mathrm{2}}}\hspace{0.17em}{\left|\mathrm{2}\mathit{\pi S}\right|}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{m}}^{†}{\mathit{S}}^{-1}\mathit{m}}\hspace{0.17em}\mathrm{𝒫}\mathrm{(}\mathit{p}\mathrm{)}\mathit{.}& & \end{array}$(30)Maximizing this PDF with respect to log (p) (see Oppermann et al. 2013) leads to $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{{\mathit{q}}_{\mathit{k}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{tr}\mathrm{(}{\mathit{S}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{(}\mathit{m}{\mathit{m}}^{†}\mathrm{+}\mathit{D}{\mathrm{)}}^{\mathrm{)}}}{{\mathit{\alpha }}_{\mathit{k}}\mathrm{-}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\varrho }}_{\mathit{k}}\mathrm{+}\mathrm{(}\mathit{T}\log \mathit{p}{\mathrm{)}}_{\mathit{k}}}\mathit{,}& & \end{array}$(31)where ${\mathit{\varrho }}_{\mathit{k}}\mathrm{=}\mathrm{tr}{}^{\mathrm{(}}\mathit{S}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}{}^{\mathrm{)}}$ is the number of degrees of freedom in the spectral band k. This formula for the power spectrum p(k) should be solved self-consistently because m and D depend on p(k) as well. Thus we arrive at an iterative scheme, where we look for a fixed point of Eq. (31).

2.2.4. Signal posterior

The signal posterior can be expressed as $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\mathit{s}\mathrm{\right|}\mathit{DM}\mathrm{)}\mathrm{=}\mathrm{\int }\mathrm{𝒟}\mathrm{(}\log \mathit{p}\mathrm{\left)}\hspace{0.17em}\mathrm{𝒫}\mathrm{\right(}\log \mathit{p}\mathrm{\left|}\mathit{DM}\mathrm{\right)}\hspace{0.17em}\mathrm{𝒫}\mathrm{\left(}\mathit{s}\mathrm{\right|}\mathit{p,DM}\mathrm{)}\mathit{,}& & \end{array}$(32)where is the signal posterior with a given power spectrum. Instead of calculating the marginalization over log p, we used Eq. (31)for the power spectrum; that is, we approximated by a Dirac peak at its maximum. This procedure is known as the empirical Bayes method. The signal posterior with a given power spectrum is proportional to the product of the signal prior and the likelihood (see Eq. (7)), $\begin{array}{ccc}& & \mathrm{𝒫}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}\mathrm{)}\mathrm{\propto }\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{s}}^{†}{\mathit{S}}^{-1}\mathit{s}\right)\\ & & \mathrm{\times }\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{DM}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{\rho }{\mathrm{)}}^{†}\stackrel{}{{{\mathit{N}}^{-1}}_{\mathrm{˜}}}\mathrm{(}\mathit{DM}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{\rho }\mathrm{)}\right)\mathit{.}\end{array}$(33)However, as has been demonstrated in Sects. 2.2.3 and 2.2.1, S and $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{N}\end{array}$ depend on the mean and the covariance of $\mathrm{𝒫}\mathrm{\left(}\mathit{s}\mathrm{\right|}\mathit{DM,p,}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{N}\end{array}\mathrm{)}$ leading to a circular dependence that needs to be solved self-consistently.

We approximate the mean of the posterior by minimizing the joint Hamiltonian $\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{N}\end{array}\mathrm{)}\mathrm{:}\mathrm{=}\mathrm{-}\log \mathrm{𝒫}\mathrm{(}\mathit{s,DM}\mathrm{\left|}\mathit{p,}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{N}\end{array}\mathrm{\right)}$ with respect to s, $\begin{array}{ccc}{\mathit{m}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{\approx }\underset{\mathit{s}}{\arg \hspace{0.17em}\min }\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}\mathrm{)}\mathit{,}& & \end{array}$(34)and its covariance by the inverse Hessian at that minimum, $\begin{array}{ccc}{\mathit{D}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{\approx }{\left({\frac{{\mathit{\delta }}^{\mathrm{2}}}{\mathit{\delta s\delta }{\mathit{s}}^{†}}\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}\mathrm{)}|}_{\mathit{s}\mathrm{=}\mathit{m}}\right)}^{-1}\mathit{.}& & \end{array}$(35)These estimates are the maximum a posteriori (MAP) estimates of s. Consequently, m^(ρ) and D^(ρ) are estimated as $\begin{array}{ccc}{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{\approx }\mathrm{\Delta }\mathrm{\left(}{x}\mathrm{\right)}{\mathrm{e}}^{{\mathit{m}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{\left(}{x}\mathrm{\right)}}& & \end{array}$(36)and $\begin{array}{ccc}{\mathit{D}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}\mathrm{\left(}{x}\mathit{,}{y}\mathrm{\right)}& & \mathrm{\approx }\mathrm{\Delta }\mathrm{\left(}{x}\mathrm{\right)}{\mathrm{e}}^{{\mathit{m}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{\left(}{x}\mathrm{\right)}}\left({\mathrm{e}}^{{\mathit{D}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{\left(}{x}\mathit{,}{y}\mathrm{\right)}}\mathrm{-}\mathrm{1}\right){\mathrm{e}}^{{\mathit{m}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{\left(}{y}\mathrm{\right)}}\mathrm{\Delta }\mathrm{\left(}{y}\mathrm{\right)}\mathit{.}\\ & & \end{array}$(37)S is then constructed as $\begin{array}{ccc}\mathit{S}\mathrm{\left(}{x}\mathit{,}{y}\mathrm{\right)}\mathrm{=}\sum _{\mathit{k}}\hspace{0.17em}{\mathit{S}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}{x}\mathit{,}{y}\mathrm{\right)}\mathit{p}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}& & \end{array}$(38)with p(k) given by Eq. (31). $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{N}\end{array}$ is constructed using Eqs. (18)and (20)as $\begin{array}{ccc}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}{\mathrm{)}}_{\mathit{ij}}\mathrm{=}\mathrm{(}\mathit{N}{\mathrm{)}}_{\mathit{ij}}\mathrm{+}{\mathit{\delta }}_{\mathit{ij}}\hspace{0.17em}\mathrm{tr}\mathrm{(}{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}†}{\mathit{F}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{)}\mathit{,}& & \end{array}$(39)where δ_ij is the Kronecker delta.

2.3. Galactic profile inference

To infer the Galactic profile field Δ, we introduce $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{s}\end{array}\mathrm{\equiv }\mathit{s}\mathrm{+}\log \mathrm{\left(}\mathrm{\Delta }\mathrm{\right)}\mathrm{\equiv }\log \mathrm{\left(}\mathit{\rho }\mathrm{\right)}$. The Galactic profile is to capture the most prominent symmetries of a disk galaxy, namely its rotational symmetry and the scaling behavior with radial distance from the Galactic center and vertical distance from the Galactic plane. Using Eq. (6), μ ≡ log (Δ) becomes $\begin{array}{ccc}\mathit{\mu }\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathrm{=}\mathit{\alpha }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}\mathit{\beta }\mathrm{\left(}\mathrm{\right|}\mathit{z}\mathrm{\left|}\mathrm{\right)}\mathit{,} \mathrm{with} \mathit{r}\mathrm{\equiv }\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(40)where α and β are one-dimensional functions describing the average behavior with respect to the radial distance and the vertical distance from the Galactic center. Including the shift by μ from s to $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{s}\end{array}$ yields the signal prior $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ \mathit{s}\end{array}\mathrm{\right|}\mathit{\alpha ,\beta }\mathrm{)}\mathrm{=}\mathrm{𝒢}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{s}\end{array}\mathrm{-}\mathit{\mu ,S}\mathrm{)}\mathit{.}& & \end{array}$(41)We do not wish to assume specific functions α and β, but to infer them. To this end we chose a Gaussian prior, $\begin{array}{ccc}\mathrm{𝒫}\mathrm{\left(}\mathit{\alpha ,\beta }\mathrm{\right)}\mathrm{\propto }\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\alpha }}{\mathit{\partial }{\mathit{r}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{\beta }}^{\mathrm{2}}}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\beta }}{\mathit{\partial }\mathrm{|}\mathit{z}{\mathrm{|}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\right)\mathit{,}& & \end{array}$(42)with the second derivative of α (or β, respectively) as the argument. This prior prefers linear functions for α and β and thus Galactic profile fields with an exponential fall-off (or rise). To simplify the notation, we defined ξ(r, | z | ) = (α(r),β( | z | ))^T and introduced the linear operators Ξ and X, where $\begin{array}{ccc}\mathit{X\xi }\mathrm{=}\mathit{\alpha }\mathrm{+}\mathit{\beta }\mathrm{\equiv }\mathit{\mu }& & \end{array}$(43)and $\begin{array}{ccc}{\mathit{\xi }}^{†}\mathrm{\Xi }\mathit{\xi }\mathrm{=}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\alpha }}{\mathit{\partial }{\mathit{r}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{\beta }}^{\mathrm{2}}}{\left(\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\beta }}{\mathit{\partial }\mathrm{|}\mathit{z}{\mathrm{|}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\mathrm{·}& & \end{array}$(44)Now we can write the Hamiltonian of ξ given a specific electron density as $\begin{array}{ccc}\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{\xi }\mathrm{\right|}\begin{array}{c}\mathrm{˜}\\ \mathit{s}\end{array}\mathrm{)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}{\left(\begin{array}{c}\mathrm{˜}\\ \mathit{s}\end{array}\mathrm{-}\mathit{X\xi }\right)}^{†}{\mathit{S}}^{-1}\left(\begin{array}{c}\mathrm{˜}\\ \mathit{s}\end{array}\mathrm{-}\mathit{X\xi }\right)\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\xi }}^{†}\mathrm{\Xi }\mathit{\xi }\mathrm{+}\mathrm{const}\mathit{.}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}{\mathit{\xi }}^{†}\mathrm{(}{\mathit{X}}^{†}{\mathit{S}}^{-1}\mathit{X}\mathrm{+}{\mathrm{\Xi }}^{\mathrm{)}}\mathit{\xi }\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {\mathit{s}}^{†}\end{array}{\mathit{S}}^{-1}\mathit{X\xi }\mathrm{+}\mathrm{const}\mathit{.}\\ & \mathrm{\equiv }& \end{array}$(45)with ${\mathit{D}}_{\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}\mathrm{=}{{}^{\mathrm{(}}\mathit{X}^{†}{\mathit{S}}^{-1}\mathit{X}\mathrm{+}{\mathrm{\Xi }}^{\mathrm{)}}}^{-1}$ and ${\mathit{j}}_{\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}\mathrm{=}{\mathit{X}}^{†}{\mathit{S}}^{-1}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{s}\end{array}$. Since this Hamiltonian is a quadratic form in ξ, the mean of the corresponding Gaussian PDF is $\begin{array}{ccc}{⟨\mathit{\xi }⟩}_{\mathrm{\left(}\mathit{\xi }\mathrm{\right|}\begin{array}{c}\mathrm{˜}\\ \mathit{s}\end{array}\mathrm{)}}\mathrm{=}{\mathit{D}}_{\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}{\mathit{j}}_{\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}\mathit{.}& & \end{array}$(46)

2.4. Filter equations

Using the posterior estimates presented in the previous section, we arrive at the following iterative scheme for reconstructing the density ρ:

• 1.

  Make an initial guess for the power spectrum (e.g., some power law) and the additive term in the noise covariance (e.g., simple relative error propagation).

• 2.

  With the current estimates for p and $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{N,}\end{array}$ the Hamiltonian, $\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{N}\end{array}\mathrm{)}\mathrm{+}\mathrm{const}\mathit{.}\mathrm{\equiv }\log \mathrm{𝒫}\mathrm{(}\mathit{s,DM}\mathrm{\left|}\mathit{p,}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{N}\end{array}\mathrm{\right)}$, is $\begin{array}{ccc}\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}\mathrm{)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}{\mathit{s}}^{†}\left(\sum _{\mathit{k}}{\mathit{S}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathit{p}}_{\mathit{k}}^{-1}\right)\mathit{s}\\ & & \mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\left({\mathrm{e}}^{\mathit{s}}\mathrm{\ast }\mathrm{\Delta }\right)}^{†}\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}^{†}\end{array}\stackrel{}{{{\mathit{N}}^{-1}}_{\mathrm{˜}}}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\left({\mathrm{e}}^{\mathit{s}}\mathrm{\ast }\mathrm{\Delta }\right)\\ & & \mathrm{-}\mathit{D}{\mathit{M}}^{†}\stackrel{}{{{\mathit{N}}^{-1}}_{\mathrm{˜}}}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\left({\mathrm{e}}^{\mathit{s}}\mathrm{\ast }\mathrm{\Delta }\right)\mathit{,}\end{array}$(47)where theasterisk denotes point-wise multiplication in position space.

• 3.

  The MAP estimate of this Hamiltonian is calculated as $\begin{array}{ccc}{\mathit{m}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{=}\underset{\mathit{s}}{\arg \hspace{0.17em}\min }\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}\mathrm{)}\mathit{,}& & \end{array}$(48)with the covariance estimate (see Appendix B) $\begin{array}{ccc}{\mathit{D}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{=}{\left({\frac{{\mathit{\delta }}^{\mathrm{2}}}{\mathit{\delta s\delta }{\mathit{s}}^{†}}\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}\mathrm{)}|}_{\mathit{s}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{m}}\right)}^{-1}\mathit{.}& & \end{array}$(49)

• 4.

  The updated power spectrum is the solution (with respect to p(k)) of the equation $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{{\mathit{q}}_{\mathit{k}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{tr}\mathrm{(}{\mathit{S}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{(}\mathit{m}{\mathit{m}}^{†}\mathrm{+}\mathit{D}{\mathrm{)}}^{\mathrm{)}}}{{\mathit{\alpha }}_{\mathit{k}}\mathrm{-}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\varrho }}_{\mathit{k}}\mathrm{+}\mathrm{(}\mathit{T}\log \mathit{p}{\mathrm{)}}_{\mathit{k}}}\mathrm{·}& & \end{array}$(50)

• 5.

  The updated effective noise covariance is calculated as $\begin{array}{ccc}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}{\mathrm{)}}_{\mathit{ii}}\mathrm{=}\mathrm{(}\mathit{N}{\mathrm{)}}_{\mathit{ii}}\mathrm{+}\mathrm{tr}\mathrm{(}{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}†}{\mathit{F}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{)}\mathit{,}& & \end{array}$(51)with $\begin{array}{ccc}{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\mathrm{\Delta }\mathrm{\left(}{x}\mathrm{\right)}\exp \mathrm{(}{\mathit{m}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{(}{x}{\mathrm{)}}^{\mathrm{)}}\mathit{.}& & \end{array}$(52)

• 6.

  The updated Galactic profile field is $\begin{array}{ccc}& & \mathrm{\Delta }\mathrm{=}\exp \mathrm{(}\mathit{X}{\mathit{m}}^{\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}\mathrm{)} \mathrm{with}\\ & & {\mathit{m}}^{\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}\mathrm{=}{\mathrm{(}{\mathit{X}}^{†}{\mathit{S}}^{-1}\mathit{X}\mathrm{+}{\mathrm{\Xi }}^{\mathrm{)}}}^{-1}{\mathit{X}}^{†}{\mathit{S}}^{-1}\log \mathrm{\left(}{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}\mathrm{\right)}\mathit{.}\end{array}$(53)

• 7.

  Repeat from step 2 until convergence is reached.

When the solution of this set of equations is converged, the estimate of the density ρ is $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}{x}\mathrm{\right)}\mathrm{\approx }{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{\pm }{\mathit{\sigma }}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}\mathrm{\left(}{x}\mathrm{\right)}\mathit{,}& & \end{array}$(54)where the confidence interval σ^(ρ) is defined as $\begin{array}{ccc}{\mathit{\sigma }}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{:}\mathrm{=}\sqrt{{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{(}{\mathrm{e}}^{{\mathit{D}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{\left(}{x}\mathit{,}{x}\mathrm{\right)}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}\mathrm{\left(}{x}\mathrm{\right)}}\mathit{.}& & \end{array}$(55)

3. Application to simulated data

To test the reconstruction of the Galactic free electron density distribution with the SKA, we generated mock data sets of pulsars with various distance uncertainties. We simulated pulsar populations using the PSRPOPpy package by Bates et al. (2014), which is based on the pulsar population model by Lorimer et al. (2006). The generated populations take the observational thresholds of the SKA into account (mid-frequency). These data sets sample modified versions of the NE2001 model by Cordes & Lazio (2002) through dispersion measures.

3.1. Galaxy model

We deactivated⁷ all local ISM components as well as all clumps and voids in the NE2001 model. We kept the clump in the Galactic center, since it is the only one at a distinguished position. We evaluated⁸ the resulting free electron density model in a 512 × 512 × 64 pixel grid centered on the Galactic center with a pixel edge length of 75 pc. This means that our model extends out to 2400 pc from the Galactic plane. We assumed a density of zero outside of this regime when calculating the dispersion measures. The resulting density field is very smooth. We generated three Gaussian random fields that follow a power-law distribution with a spectral index⁹ of −4.66, but have different fluctuation amplitudes. We verified that the Sun is located in an underdensity in these random fields. Then we added these three random field maps to our smooth map of log (n_e) to create three different modified versions of NE2001. In Fig. 2 we depict the power spectra of the smooth NE2001 field (without local features, clumps, and voids) and the power spectra of its three modified versions.

3.1.1. Contrast-enhanced model

The three Galaxy models we generated from NE2001 have relatively little contrast in the sense that under- and overdense regions differ by relatively moderate factors. For example, the density in the region between the Perseus and Carina-Sagittarius arm where the Sun is located is only a factor of three lower than in the Perseus arm itself. Since the Perseus arm is a less than 1 kpc in width, any excess dispersion measure that is due to the arm can also be explained by an underestimated pulsar distance for many lines of sight. In consequence, we expect the reconstruction quality to improve if the input model has higher contrast. Therefore, we prepared one additional model with enhanced contrast. To this end, we took the input model with medium-strength fluctuations as described above. We divided out the scaling behavior in radial and vertical directions using the scale heights from NE2001. We squared the density and divided it by a constant to ensure that the mean density in the Galactic plane remains unchanged¹⁰. Finally we multiplied the resulting density with the scaling functions to restore the original scaling in radial and vertical directions.

This procedure yielded a Galaxy model sharing the same morphology and scaling behavior as the input model. Averaged over the lines of sight, the value of the dispersion measures is roughly unchanged. But the contrast is twice as strong, that is, the previously mentioned factor between the density in the Perseus arm and the inter-arm region is now squared from 3 to 9. We show a picture of the density in the Galactic plane of this model in Sect. 4.1, where we compare it with its reconstruction.

3.2. Simulated population and survey

We used the SKA template in the PSRPOPpy package, but reduced the maximum declination in equatorial coordinates to 50° (this is due to the SKAs position on the Southern Hemisphere, see e.g. Smits et al. 2009). This yielded a detected population of roughly 14 000 pulsars. Out of these, we took the first 1000, 5000, or 10 000 pulsars as our test populations. The population was not ordered in any sense, which means that the first 1000 pulsars, for example, represent a random sample from the whole detected population. In reality, the Malmquist bias preferentially selects pulsars that lie close to the Sun. We chose a random selection, however, to see the effect of the population size on the quality of the reconstruction more clearly. In Fig. 3 we depict the population of 10 000 pulsars projected onto the sky. The pulsars are concentrated toward the center of the Galaxy. The gap in the equatorial Northern Hemisphere is clearly evident in the left part of the plot.

3.3. Simulated dispersion measures and distances

We calculated the line integrals through the Galaxy models from the positions generated by the PSRPOPpy package to the Sun to generate simulated dispersion measures. We added Gaussian random variables to the pulsar distances to simulate measurement uncertainties of the distances. For each pulsar we generated one random number and scale this to 5%, 15%, or 25% of the distance of the pulsar. In reality the distance PDF would be non-Gaussian. The exact form depends on the combination of observables that are used to infer the distance. We used Gaussian PDFs for simplicity. As long as the real distance PDFs are unimodal, we do not expect this choice to have a significant effect on our study. We did not simulate additional measurement noise for the dispersion measures because it is expected to be small compared to the distance uncertainty. This leaves us with a number of data sets that are described in Table 1. In this table we omit the combinations of 1000 pulsars at 25% distance error (because we do not hope for a good reconstruction in that case) and 10 000 pulsars at 5% distance error (because we deem it to be too unrealistic).

The aforementioned measurement scenarios were chosen to see the effect of the population size and the distance error on the reconstruction in isolation. A more realistic setting is of course a mix of distance uncertainties where more distant pulsars have larger distance errors on average. We therefore created one additional measurement scenario for 10 000 pulsars, where we assigned the uncertainty magnitude of each pulsar randomly¹¹. The distance uncertainties are distributed as shown in Fig. 4. In this measurement set, 2969 pulsars have a 5% distance error, 3400 pulsars have a 15% distance error, and 3631 pulsars have a 25% distance error. Throughout the rest of this paper we refer to this data set as the mixed data set.

A very rough estimate of the scales that we can hope to resolve is given by the mean distance between neighboring pulsars and the average misplacement that is due to distance errors. The mean distance between neighboring pulsars is 490 pc for 1000, 290 pc for 5000, and 230 pc for 10 000 pulsars. The average misplacement is 380 pc for 5%, 1100 pc for 15%, and 1900 pc for 25% distance errors and 1300 pc for the mixed data set. Interpreting these distances as independent uncertainties, we can combine them by adding the squares and taking the square root. This provides a rough estimate of sampling distances. In Table 2 we list these distances for each data set.

3.4. Algorithm setup

The algorithm was set up in a 128 × 128 × 48 pixel grid centered on the Galactic center with pixel dimensions¹² of 281.25 pc × 281.25 pc × 250 pc. While the dispersion measures in our data sets were free from instrumental noise, the noise was assumed to be 2% in the algorithm. This provides a lower limit for the effective noise covariance (Eq. (18)) and thus ensures stability of the inference without losing a significant amount of precision. The initial guess for the power spectrum is a broken power law with an exponent of −3.66¹³. For the propagated distance uncertainty it is $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{i}}\mathrm{=}\frac{\sqrt{\mathrm{Var}\mathrm{\left[}{\mathit{d}}_{\mathit{i}}\mathrm{\right]}}}{{\mathit{d}}_{\mathit{i}}}\mathit{D}{\mathit{M}}_{\mathit{i}}\mathit{,}& & \end{array}$(56)where d_i is the distance of the pulsar (see Sect. 2.2.1). The initial guesses of the Galactic profile functions¹⁴ are $\begin{array}{ccc}\mathit{\alpha }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{-}\mathit{r}}{\mathrm{28}\hspace{0.17em}\mathrm{000}\hspace{0.17em}\mathrm{pc}} \mathrm{and} \mathit{\beta }\mathrm{\left(}\mathrm{\right|}\mathit{z}\mathrm{\left|}\mathrm{\right)}\mathrm{=}\frac{\mathrm{-}\mathrm{\left|}\mathit{z}\mathrm{\right|}}{\mathrm{1600}\hspace{0.17em}\mathrm{pc}}\mathrm{·}& & \end{array}$(57)We discuss the convergence and final values of the power spectrum, effective errors, and profile functions in Appendix C.

4. Simulation evaluation

Our algorithm accounts for most of the variance in the data while regularizing the result to avoid overfitting. Most of the reconstructions shown in this section have corresponding reduced χ² values close to 1, indicating that they show all structures that are sufficiently well constained by the data. We discuss the reduced χ² values in detail in Appendix D.

4.1. Density in the midplane

The simulations show that with the amount of pulsars with reliable distance estimates that the SKA should deliver reconstruction of the free electron density in the vicinity of the Sun becomes feasible (see Fig. 5). However, small-scale features are difficult to identify in the reconstruction. Identifying spiral arms remains challenging as well, especially beyond the Galactic center. To resolve the spiral arms in the vicinity of the Sun, between 5000 and 10 000 pulsars with distance accuracies between 5% and 15% are needed. As is evident from the figure, small distance uncertainties increase the quality of the reconstruction significantly. The reconstruction from 5000 pulsars with 5% distance uncertainty is better in quality than the one from 10 000 pulsars with 15% distance uncertainty¹⁵. All reconstructions smooth out small-scale structure in the electron density, for example, at the Galactic center. If an overdensity appears at the wrong location, this indicates that the data do not constrain the overdensity well. For completeness we also show the recovered Galactic profile in the Galactic plane for 5000 pulsars with 5% distance uncertainty in Fig. 6. For other data sets the plot would look very similar. In Appendix E we show a reconstruction where the Galactic profile and the correlation structure a known a priori and in Appendix F we show and discuss the uncertainty estimate of the algorithm.

In Fig. 7 we compare the performance of the reconstruction algorithm for the three input model fluctuation strengths using 5000 pulsars with 5% and 15% distance uncertainty. The strength of the fluctuations clearly does not influence the quality of the reconstructions by much. The reconstructions of the models with stronger fluctuations exhibit stronger fluctuations as well, while all reconstruction omit or smear features to a similar degree. However, it clearly becomes more difficult to reconstruct the Perseus arm toward the Galactic anticenter if the fluctuations in the electron density are strong. This is to be expected because the spiral arm is also harder to recognize in the original model as the fluctuations become stronger.

In Fig. 8 we show the contrast-enhanced Galaxy model and its reconstruction using 5000 pulsars with distance uncertainties of 5%. As is clear from the figure, the algorithm is able to resolve much more detailed structure than the reconstruction of the unenhanced Galaxy model (bottom middle panel in Fig. 5). We stress that the pulsar population and their distance uncertainties are exactly the same for both cases. The increase in quality comes merely from the increased contrast and the resulting stronger imprint of under- and overdensities in the dispersion data. Therefore, we conclude that if the contrast of the real Galaxy is much stronger than in NE2001, our algorithm could resolve the Galaxy much better that the study on NE2001 indicates.

4.2. Vertical fall-off

A quantity of interest in any model of the Galactic free electron density is the drop-off of the average density with respect to distance from the Galactic plane. This behavior is shown in Fig. 9, which displays a vertical cut through the Galactic profile reconstructed using 5000 pulsars with 5% distance uncertainty.

In our parametrization the function β in Eq. (40)describes the average log-density at a certain distance from the Galactic plane. In Fig. 10 we show the estimates for β corresponding to the reconstructions shown in Fig. 5 along with their uncertainties (see Appendix G for their calculation). The uncertainty regions reflect that a vertical fall-off can be explained by a global profile as well as by density fluctuations close to the Sun. This uncertainty is nearly independent of the quality of the data set, but depends on the strength of fluctuations on kpc scales. These are always present unless the data probe a simplistic disk. Therefore, there is a lower bound of precision to which our algorithm can determine the vertical fall-off behavior. We compared the reconstructed vertical scaling to a global and a local estimate generated from the original input model. The global estimate describes vertical fall-off throughout the whole model, whereas the local estimate describes the vertical fall-off close to the Sun¹⁶. For completeness we also provide best-fitting scale heights for exponential fall-offs in the figure, that is, we fit the vertical scaling to n_e ∝ e^{− | z | /H_z}. The uncertainties of these estimates were calculated by performing the fit on multiple posterior samples of β. The local and global estimate both have significantly lower scale heights than the 950 pc from NE2001 (thick disk). This is probably due to the combination of the thick disk with the thin disk of NE2001 (which has a scale height of 140 pc). As is evident from the figure, the reconstructed z-profile is dominated by the local behavior of the density and agrees with it within the error bars throughout all data sets¹⁷ for the regime | z | < 2400 pc. However, the width of the uncertainty region prohibits a clear decision whether the vertical fall-off follows a single exponential function or a thick disk and a thin disk, as is the case for NE2001. In our input model we set n_e to zero for | z | > 2400 pc. In this regime our reconstruction is unreliable.

5. Summary and conclusions

We presented an algorithm that performs nonparametric tomography of the Galactic free electron density using pulsar dispersion measures and distances. The algorithm produces a three-dimensional map and a corresponding uncertainty map. It automatically estimates the correlation structure and the scales of the disk shape, requiring only approximate initial guesses for them. The uncertainties of pulsar distance estimates are consistently propagated.

Using our algorithm, we investigated the feasibility of nonparametric tomography with the upcoming Square Kilometer Array. To this end, we created three Galaxy models with various fluctuation strengths and one with enhanced contrast and simulated mock observations of these models using between 1000 and 10 000 pulsars. Our results indicate that with the amount of pulsars that the SKA is expected to deliver, nonparametric tomography becomes feasible. However, detecting spiral arms in the free electron density from pulsar dispersion measures alone remains challenging if the input model has unenhanced contrast. We find that to distinguish the spiral arms in the vicinity of the Sun, between 5000 and 10 000 pulsars with distance accuracies of between 5% and 15% are needed. The vertical fall-off behavior of the free electron density was recovered for all mock data sets we investigated. However, a clear decision whether the vertical fall-off of free electron density is best described by a single exponential function or a thick disk and a thin disk could not be made by our algorithm.

One way to increase the sensitivity of the algorithm for Galactic features would be to include them into the prior description. By including higher order statistics (non-Gaussian priors), the inference could be made more sensitive for spiral arm structures in the electron density. In cosmology, for example, higher order statistics allowed for a better recovery of cosmological filaments (e.g., Jasche & Wandelt 2013). Modeling HII regions and supernova remnants is also beyond the scope of Gaussian statistics. Another approach would be to include parametrized structures that are known from stellar observations, such as spatial templates for spiral arm locations. This would connect data-driven tomography (with infinite degrees of freedom) with classical model fitting.

The algorithm can also be used for other tomography problems with line-of-sight measurements, such as stellar absorption coefficients. It could also be extended to infer vector fields, enabling inference of the Galactic magnetic fields from pulsar rotation measures. Furthermore, a joint reconstruction of the Galactic free electron density and the magnetic field using pulsar dispersion, measures, pulsar rotation measures, and extragalactic Faraday sources needs to be investigated.

---

Acknowledgments

We wish to thank Niels Oppermann and Marco Selig for fruitful collaboration and advice. We also wish to thank Henrik Junklewitz, Jimi Green, Jens Jasche and Sebastian Dorn for useful discussions. The calculations were realized using the NIFTY¹⁸ package by Selig et al. (2013). Some of the minimizations described in Sect. 2.4 were performed using the L-BFGS-B algorithm (Byrd et al. 1995). We acknowledge the support by the DFG Cluster of Excellence “Origin and Structure of the Universe”. The simulations have been carried out on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP). We are grateful for the support by Frederik Beaujean through C2PAP. This research has been partly supported by the DFG Research Unit 1254 and has made use of NASA’s Astrophysics Data System.

References

Appendix A: Parameters of the power spectrum prior

The inverse gamma distribution is defined as $\begin{array}{ccc}\mathrm{ℐ}\mathrm{\left(}{\mathit{p}}_{\mathit{k}}\mathrm{;}{\mathit{\alpha }}_{\mathit{k}}\mathit{,}{\mathit{q}}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{q}}_{\mathit{k}}\mathrm{\Gamma }\mathrm{(}{\mathit{\alpha }}_{\mathit{k}}\mathrm{-}\mathrm{1}\mathrm{)}}{\left(\frac{{\mathit{p}}_{\mathit{k}}}{{\mathit{q}}_{\mathit{k}}}\right)}^{\mathrm{-}{\mathit{\alpha }}_{\mathit{k}}}\hspace{0.17em}\exp \left(\mathrm{-}\frac{{\mathit{q}}_{\mathit{k}}}{{\mathit{p}}_{\mathit{k}}}\right)\mathrm{·}& & \end{array}$(A.1)The mean and variance of this distribution are $\begin{array}{ccc}& & {⟨{\mathit{p}}_{\mathit{k}}⟩}_{\mathrm{\left(}{\mathit{p}}_{\mathit{k}}\mathrm{\right)}}\mathrm{=}{\mathit{q}}_{\mathit{k}}\mathit{/}\mathrm{(}{\mathit{\alpha }}_{\mathit{k}}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{for}\mathit{\alpha }\mathit{>}\mathrm{2}\\ & & {\mathrm{⟨}{{\mathit{p}}_{\mathit{k}}^{\mathrm{2}}}^{\mathrm{⟩}}}_{\mathrm{\left(}{\mathit{p}}_{\mathit{k}}\mathrm{\right)}}\mathrm{-}⟨{\mathit{p}}_{\mathit{k}}⟩\begin{array}{c}\mathrm{2}\\ \mathrm{\left(}{\mathit{p}}_{\mathit{k}}\mathrm{\right)}\end{array}\mathrm{=}\frac{{\mathit{q}}_{\mathit{k}}^{\mathrm{2}}}{\mathrm{(}{\mathit{\alpha }}_{\mathit{k}}\mathrm{-}\mathrm{3}\mathrm{)}\mathrm{(}{\mathit{\alpha }}_{\mathit{k}}\mathrm{-}\mathrm{2}{\mathrm{)}}^{\mathrm{2}}}\mathrm{for}\mathit{\alpha }\mathit{>}\mathrm{3.}\end{array}$(A.2)There are three properties of the prior that we wish to fulfill by choosing α_k and q_k. The prior of the monopole p₀, which corresponds to the variance of a global prefactor of the density, should be close to Jeffreys prior, that is, the limit α₀ → 1, q₀ → 0. The reason for this is that we need the algorithm to stay consistent if the units are changed, for instance, from pc to kpc. These changes introduce a global prefactor in front of the density. Jeffreys prior has no preferred scale because it is flat for log p₀ and therefore all prefactors are equally likely a priori. For other k-bins the parameters should favor, but not enforce, falling power spectra. Furthermore, since p(k) is the average power of many independent Fourier components, its a priori variance should be inversely proportional to ϱ_k (the amount of degrees of freedom in the respective k-bin), while thea priorimean should be independent of ϱ_k. We therefore set the parameters as $\begin{array}{ccc}{\mathit{q}}_{\mathit{k}}\mathrm{=}\mathit{f}\hspace{0.17em}{\mathit{\varrho }}_{\mathit{k}} \mathrm{and} {\mathit{\alpha }}_{\mathit{k}}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathit{k}}{\mathrm{100}{\mathit{k}}_{\min }}{\mathit{\varrho }}_{\mathit{k}}\mathit{,}& & \end{array}$(A.3)where k_min is the first non-zero k-value and f is a prefactor, which defines a lower cut-off of the power spectrum calculated in Eq. (31). The choice of f does not influence the result as long as it is suitably low, but higher f accelerate the convergence of the algorithm. The denominator of 100k_min before ϱ_k is chosen to introduce a preference for falling power spectra starting two orders of magnitude from the fundamental mode k_min (we note that 1 is subtracted in Eq. (27)). As long as the denominator is not too small, it has very little influence on the result of the algorithm, but smaller denominators increase the convergence speed. We found 100k_min to be a good compromise.

The parameter σ_p in Eq. (27)describes how much the power spectrum is expected to deviate from a power law. We chose σ_p = 1. If the power spectrum is locally described by a power law, σ_p = 1 means that the typical change of the exponent within a factor of e in k should be on the order of 1.

Appendix B: Functional derivatives of the Hamiltonian

To minimize the Hamiltonian in Sect. 2.4, the first derivative with respect to s is needed. It is $\begin{array}{ccc}\frac{\mathit{\delta }}{\mathit{\delta }{\mathit{s}}^{†}}\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}\mathrm{)}\mathrm{=}{\mathit{S}}^{-1}\mathit{s}\mathrm{+}\begin{array}{c}\mathrm{􏽤}\\ \mathrm{(}{\mathrm{e}}^{\mathit{s}}\mathrm{)}\end{array}\mathit{M}\left({\mathrm{e}}^{\mathit{s}}\right)\mathrm{-}\begin{array}{c}\mathrm{􏽤}\\ \mathrm{(}{\mathrm{e}}^{\mathit{s}}\mathrm{)}\end{array}\mathit{j,}& & \end{array}$(B.1)where the hat converts a field to a diagonal operator in position space, for example, $\begin{array}{c}{}^{\mathrm{􏽢}}\\ \mathit{\xi }\end{array}\mathrm{\left(}{x}\mathit{,}{y}\mathrm{\right)}\mathrm{=}\mathit{\xi }\mathrm{\left(}{x}\mathrm{\right)}\mathit{\delta }\mathrm{(}{x}\mathrm{-}{y}\mathrm{)}$, and we used the shorthand notations $\begin{array}{ccc}& & {\mathit{S}}^{-1}\mathrm{\equiv }\sum _{\mathit{k}}{\mathit{S}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathit{p}}_{\mathit{k}}^{-1}\\ & & \mathit{M}\mathrm{\equiv }\begin{array}{c}\mathrm{􏽢}\\ \mathrm{\Delta }\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}^{†}\end{array}\stackrel{}{{{\mathit{N}}^{-1}}_{\mathrm{˜}}}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\begin{array}{c}\mathrm{􏽢}\\ \mathrm{\Delta }\end{array}\\ & & \mathit{j}\mathrm{\equiv }\begin{array}{c}\mathrm{􏽢}\\ \mathrm{\Delta }\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}^{†}\end{array}\stackrel{}{{{\mathit{N}}^{-1}}_{\mathrm{˜}}}\mathit{DM}\mathit{.}\end{array}$(B.2)The second derivative in Eq. (49)is $\begin{array}{ccc}\frac{{\mathit{\delta }}^{\mathrm{2}}}{\mathit{\delta s\delta }{\mathit{s}}^{†}}\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}\mathrm{)}& \mathrm{=}& {\mathit{S}}^{-1}\mathrm{+}\begin{array}{c}\mathrm{􏽤}\\ \mathrm{(}{\mathrm{e}}^{\mathit{s}}\mathrm{)}\end{array}\mathit{M}\begin{array}{c}\mathrm{􏽤}\\ \mathrm{(}{\mathrm{e}}^{\mathit{s}}\mathrm{)}\end{array}\\ & & \mathrm{+}\begin{array}{c}\mathrm{􏽤}\\ \mathrm{(}{\mathrm{e}}^{\mathit{s}}\mathrm{)}\end{array}\begin{array}{c}[\\ \mathit{M}\mathrm{(}{\mathrm{e}}^{\mathit{s}}\mathrm{)}\end{array}\mathrm{-}\begin{array}{c}\mathrm{􏽤}\\ \mathrm{(}{\mathrm{e}}^{\mathit{s}}\mathrm{)}\end{array}\hspace{0.17em}\begin{array}{c}\mathrm{􏽢}\\ \mathit{j}\end{array}\mathit{.}\end{array}$(B.3)The last term in the second derivative can be problematic because it can break the positive definiteness¹⁹ of the second derivative, which is crucial for efficiently applying inversion techniques such as the conjugate gradient method. However, a closer inspection of the last two terms (omitting the hats for readability) $\begin{array}{ccc}\mathit{M}\left({\mathrm{e}}^{\mathit{s}}\right)\mathrm{-}\mathit{j}\mathrm{=}\begin{array}{c}\mathrm{􏽢}\\ \mathrm{\Delta }\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{R}}^{†}\end{array}\stackrel{}{{{\mathit{N}}^{-1}}_{\mathrm{˜}}}\left(\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}\begin{array}{c}\mathrm{􏽢}\\ \mathrm{\Delta }\end{array}{\mathrm{e}}^{\mathit{s}}\mathrm{-}\mathit{DM}\right)\mathrm{\propto }\begin{array}{c}\mathrm{􏽧}\\ \mathit{DM}\end{array}\mathrm{-}\mathit{DM}& & \end{array}$(B.4)shows that their contribution is proportional to the difference between the real dispersion data DM and the idealized data generated by the map, $\begin{array}{c}{}_{\mathrm{􏽧}}\\ \mathit{DM}\end{array}\mathrm{=}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\begin{array}{c}{}^{\mathrm{􏽢}}\\ \mathrm{\Delta }\end{array}{\mathrm{e}}^{\mathit{s}}$. These two terms counteract each other at the minimum, and we therefore omit them to gain numerical stability. Hence the second derivative is approximated as $\begin{array}{ccc}\frac{{\mathit{\delta }}^{\mathrm{2}}}{\mathit{\delta s\delta }{\mathit{s}}^{†}}\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{s,DM}\mathrm{\right|}\mathit{p,}\begin{array}{c}\mathrm{˜}\\ \mathit{N}\end{array}\mathrm{)}\mathrm{\approx }{\mathit{S}}^{-1}\mathrm{+}\begin{array}{c}\mathrm{􏽤}\\ \mathrm{(}{\mathrm{e}}^{\mathit{s}}\mathrm{)}\end{array}\mathit{M}\begin{array}{c}\mathrm{􏽤}\\ \mathrm{(}{\mathrm{e}}^{\mathit{s}}\mathrm{)}\end{array}\mathit{.}& & \end{array}$(B.5)

Appendix C: Convergence

In this section, we display the convergence behavior of the power spectrum, effective errors, and profile functions. For the sake of brevity, we limit the discussion to the reconstruction of the Galaxy model with fluctuations of medium strength using the data set of 5000 pulsars with a distance error of 5%. In Appendix E we show the reconstruction of this model and data set using the real power spectrum and profile functions as a benchmark on how well our iterative estimation of them does.

In Fig. C.1 we show the convergence of the power spectrum. As is evident, the power moves away from the initial guess to a fixed point. Compared to the spectrum of the logarithmic input model, the converged spectrum misses power in the large- and small-scale regime. The loss of power in the large-scale regime arises because the profile field absorbs large features, in the small-scale regime it is due to the general loss of small-scale power.

The loss of small-scale power comes from two effects. First, the dispersion measure data sample the density sparse and irregularly. Without the regularization imposed by the prior, this would lead to severe aliasing, as is commonly known from Fourier analysis because the prior typically suppresses aliasing from large to small scales and the algorithm consequently interprets the power as noise. Aliasing from small to large scales is negligible, since the input model is spatially correlated and therefore has a falling power spectrum. Second, some loss of power is due to the distance uncertainties. These make the likelihood less informative about small-scale structures, which are in consequence suppressed by the prior. This effect yields no aliasing, but smoothens the resulting map (which is desired to avoid overfitting). For further details about the loss of power in filtering algorithms such as the one in this paper we refer to Enßlin & Frommert (2011).

The fixed point power spectrum falls as a power law with index −5.5 for k> 2 × 10^-4. Our algorithm allows for spectral indices of up to −5.5. Without this limit, the power spectrum would fall²⁰ to a minimal value of q_k/ϱ_k for k ≳ 3 × 10^-4. However, introducing a hard limit speeds up the convergence, and using a slope of 5.5 made no difference toward the lower limit for the resulting maps in our tests. This hard limit can be thought of as part of the power spectrum prior.

In Fig. C.2 we show the convergence of the propagated distance variance (Eq. (20)) of a random selection of ten data points. The behavior seen in this plot is qualitatively the same for all data points we investigated. Most data points clearly reach convergence rather quickly, but there are also outliers. In this plot, the lowest line exhibits a kink after it had seemingly already converged. This behavior is unfortunately not entirely suppressible, but it appears to have very little effect on the resulting map because only a small fraction of data points shows this behavior.

In Figs. C.3 and C.4 we show the convergence of the profile functions, where we shifted the functions by a global value to line them up at β( | z | = 0) and α(r = 0). We note that the functions α and β are degenerate with respect to a global addition in their effect on the Galactic profile field and degenerate with the monopole of s as well. This is why a shift by a constant for plotting purposes is reasonable. The z-profile function β seems to reach a fixed point for | z | < 2400 pc. For higher | z | the profile function reaches no clear fixed point. However, for the Galactic profile, where the profile function is exponentiated, this makes only a small difference, since β is already three e-foldings below its values at | z | = 0. The radial profile function α seems to only correct the initial guess mildly, and it is not clear whether the result is independent of the initial guess. However, it appears that α does reach a fixed point.

Appendix D: Goodness of fit (χ²) of the reconstructions

In this section, we discuss the goodness of fit characterized by the reduced χ² value, $\begin{array}{ccc}{\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{data}}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{data}}}\frac{{\mathrm{(}\mathit{D}{\mathit{M}}_{\mathit{i}}\mathrm{-}{\begin{array}{c}\mathrm{􏽧}\\ \mathit{DM}\end{array}}_{\mathit{i}}\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(D.1)where $\begin{array}{c}{\mathrm{􏽧}}_{}\\ \mathit{DM}\end{array}$ is the dispersion data reproduced by applying the response to our reconstruction, $\begin{array}{ccc}{\begin{array}{c}\mathrm{􏽧}\\ \mathit{DM}\end{array}}_{\mathit{i}}\mathrm{=}{\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{R}\end{array}{\mathit{m}}^{\mathrm{\left(}\mathit{\rho }\mathrm{\right)}}\mathrm{)}}_{\mathit{i}}\mathit{,}& & \end{array}$(D.2)and σ_i is the crudely propagated distance uncertainty given by $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{i}}\mathrm{=}\frac{\sqrt{\mathrm{Var}\mathrm{\left[}{\mathit{d}}_{\mathit{i}}\mathrm{\right]}}}{{\mathit{d}}_{\mathit{i}}}\mathit{D}{\mathit{M}}_{\mathit{i}}\mathit{.}& & \end{array}$(D.3)

To set the χ² value into perspective, we compared them with the χ² values of the null model (m^(ρ) = 0) and the Galactic profile alone (m^(ρ) = Δ) The reduced χ² values corresponding to the maps shown in Fig. 5 are shown in Table D.1. The reconstruction of the mixed data set has a reduced χ² of 1.31. In Table D.2 we show the reduced χ² values of the maps shown in Fig. 7. For the reconstruction of the contrast-enhanced input model shown in Fig. 8 the reduced χ² value is 345 for the null model, 110 for the profile alone, and 2.6 for the full reconstruction.

It is evident from the tables that our reconstructions account for a large part of the data variance in all cases. The Galactic profile without local fluctuations also accounts for a large part of the variance, especially if the distance uncertainties are high and the fluctuation strength of the input model is weak. For our reconstructions the χ² values are close to 1 for the 25% and 15% data sets. Therefore, we assume that our inference mechanism resolved the most relevant information in the data sets and that the prior assumptions are not too restrictive for these data sets. For the 5% reconstructions the χ² values are around 3. This is a hint that the data might contain more information than the reconstruction resolves and that more elaborate prior assumptions might yield a better map. However, how to achieve this is a non-trivial question that we do not aim to answer in this work.

Appendix E: Reconstruction with real power spectrum and profile functions

The posterior map our algorithm finds depends on the prior power spectrum, the effective errors, and the profile functions, all of which are simultaneously estimated from the data. To benchmark the efficiency of this joint estimation, we investigated the case where the real power spectrum and the real profile functions are known, that is, only iterating the effective errors. The resulting map serves as an indicator whether our ansatz with unknown hyper parameters is sensible or if the problem is too unrestricted in this setting. We depict the map resulting from the real hyper parameters in Fig. E.1. The morphology of the result obviously does not change. More of the small-scale structure is resolved, and the intensity of the overdensity between the Sun and the Galactic center, which belongs to the ring the original model, is more pronounced. Consequently, this map has a better reduced χ² value of 1.95 compared to the value of 2.56 of our map with unknown hyperparameters. But considering the amount of unknowns, this is a satisfactory result. We therefore regard our estimation procedure for the hyper parameters as sensible.

Appendix F: Uncertainty map

Here, we discuss the 1σ uncertainty map of reconstructing the Galaxy model with medium fluctuations using the data set with 5000 pulsars and 5% distance uncertainty. We compare the uncertainty map for unknown profile and power spectrum and the uncertainty map for a known profile and power spectrum (see Appendix E) with the corresponding absolute errors. The density maps are shown in Fig. E.1. The uncertainty estimates σ^(ρ) (see Eq. (55)) are shown in Fig. F.1. These uncertainties are underestimated because they are calculated from the curvature of the negative log-posterior around its minimum (see Eq. (35)), not from the full distribution. By visual comparison with the absolute error²¹, | m^(ρ) − ρ |, we estimate that the uncertainty estimates are underestimated by a factor of roughly 3. However, their morphology seems to be reliable.

Appendix G: Uncertainties of the vertical fall-off

In principle, the posterior variance of α and β is the diagonal of the operator D_(ξ) (see Eq. (45)). However, this diagonal is too large because α and β are completely degenerate with respect to a constant shift (α + c and β − c yield the same profile as α and β). This degeneracy yields a large point variance, which is not instructive for quantifying the uncertainty of the vertical fall-off. Therefore, we projected out the eigenvector corresponding to the constant shift before calculating the diagonal of D_(ξ). This corrected diagonal is the squared 1σ uncertainty that we plot in Fig. 10.

All Tables

All Figures