Inferring Galactic magnetic field model parameters using IMAGINE - An Interstellar MAGnetic field INference Engine
Theo Steininger, Torsten A. Enßlin, Maksim Greiner, Tess Jaffe, Ellert van der Velden, Jiaxin Wang, Marijke Haverkorn, Jörg R. Hörandel, Jens Jasche, Jörg P. Rachen
AAstronomy & Astrophysics manuscript no. IMAGINE c (cid:13)
ESO 2018January 16, 2018
Inferring Galactic magnetic field model parameters using IMAGINE
An Interstellar MAGnetic field INference Engine
Theo Steininger , , Torsten A. Enßlin , , Maksim Greiner , Tess Ja ff e , , Ellert van der Velden , , Jiaxin Wang , ,Marijke Haverkorn , Jörg R. Hörandel , , Jens Jasche , and Jörg P. Rachen Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany Ludwig-Maximilians-Universität München, Geschwister-Scholl-Platz 1, 80539, München, Germany Excellence Cluster Universe, Technische Universität München, Boltzmannstrasse 2, 85748 Garching, Germany CRESST, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA Department of Astronomy, University of Maryland, College Park, MD 20742, USA Department of Astrophysics / IMAPP, Radboud University, Heyendaalseweg 135, 6525 AJ Nijmegen, The Netherlands Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, VIC 3122, Australia Scuola Internazionale Superiore di Studi Avanzati, Via Bonomea 265, 34136 Trieste, Italy Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, Via Bonomea 265, 34136 Trieste, Italy Nikhef, Science Park Amsterdam, 1098 XG Amsterdam, The NetherlandsReceived January 12, 2018; accepted 0000 00, 0000
ABSTRACT
Context.
The Galactic magnetic field (GMF) has a huge impact on the evolution of the Milky Way. Yet currently there exists nostandard model for it, as its structure is not fully understood. In the past many parametric GMF models of varying complexity havebeen developed that all have been fitted to an individual set of observational data complicating comparability.
Aims.
Our goal is to systematize parameter inference of GMF models. We want to enable a statistical comparison of di ff erent modelsin the future, allow for simple refitting with respect to newly available data sets and thereby increase the research area’s transparency.We aim to make state-of-the-art Bayesian methods easily available and in particular to treat the statistics related to the randomcomponents of the GMF correctly. Methods.
To achieve our goals, we built I magine , the
Interstellar Magnetic Field Inference Engine . It is a modular open sourceframework for doing inference on generic parametric models of the Galaxy. We combine highly optimized tools and technology suchas the M ulti N est sampler and the information field theory framework NIFT y in order to leverage existing expertise. Results.
We demonstrate the steps needed for robust parameter inference and model comparison. Our results show how impor-tant the combination of complementary observables like synchrotron emission and Faraday depth is while building a model andfitting its parameters to data. I magine is open-source software available under the GNU General Public License v3 (GPL-3) at:https: // gitlab.mpcdf.mpg.de / ift / IMAGINE
Key words.
Galaxy: general, Galaxy: structure, methods: numerical, methods: statistical, methods: data analysis
1. Introduction
The interstellar magnetic field in galaxies plays a key role in pro-cesses at various scales from star formation up to overall galacticevolution. Its energy density is comparable to that of the tur-bulent gas or cosmic rays (CRs), and therefore the dynamicalfeedback on the interstellar medium (ISM) must not be ignored.Galactic magnetic fields a ff ect in- and outflows of the ISM thatalready exist as well as the formation of new ones. They influ-ence the propagation of CRs, which gyrate along the field lines.Though these e ff ects are all important, it is challenging to inferthe field, since it is only accessible via indirect detection meth-ods. Additionally, since our Solar System is located within theGalactic plane, the tracers of the Galactic magnetic field (GMF)in our own Milky Way are highly degenerate as they are line-of-sight integrated quantities. This also means that the view ofthe opposite side of the Galaxy is obstructed by the interveningISM. Because of all this, the GMF is currently mainly modeledvia heuristic parametric models that have physically motivatedfeatures. The degrees of freedom in those models are morpho- logical properties, field strengths (of possibly individual spatialcomponents) of the magnetic field and the strength and charac-teristics of random contributions. Significant progress has beenmade here, which is the reason why a rather large number ofGMF models is available today. At the same time the availabledata becomes better and better. Hence, there is need for a stan-dardized platform that allows systematic parameter estimationand model comparison for a continuously expanding abundanceof models and data.
2. Bayesian Parameter Inference and ModelComparison
The GMF can naturally be thought of as a vector field withan infinite number of degrees of freedom: under the constraintof zero divergence the magnetic field can have an individualstrength and direction at every point in space. This view cor-responds to the most generic model possible, where the model’sparameters are the field’s degrees of freedom. To infer the GMFone must simplify this most generic model, for example, by dis-
Article number, page 1 of 25 a r X i v : . [ a s t r o - ph . I M ] J a n & A proofs: manuscript no. IMAGINE cretizing space. Doing so reduces the model parameters to a fi-nite but still huge number, namely twice the number of voxelsof the considered volume. However, now one can try to con-cretely infer the magnetic field voxel by voxel, a method knownas non-parametric modelling. Generally speaking, constrainingthose non-parametric models is certainly hard, because the hugenumber of degrees of freedom often are counteracted by a lim-ited amount of data. Because of this, one often builds a sim-pler model with a heavily reduced number of parameters, whichtherefore only covers a tiny slice in the full parameter space butstill represents the most important features of the modeled quan-tity. In the case of the GMF, various models have been developedthat di ff er greatly in their complexity: the number of parametersvaries between only a few and up to 40. Given a model and ob-servational data one must find an estimate for a set of the model’sparameters that explains the observed data well. However, in ad-dition to the parameter estimation of a given model, there is alsothe task of comparing the plausibility of di ff erent models. In thecase of GMF inference this is especially important since so farthere is no standard model available.In terms of Bayesian inference, parameter estimation andmodel comparison can be described by the following compo-nents: a given model m that has a set of parameters θ shall beconstrained by data d . This means, that we are interested in theposterior probability density P ( θ | d , m ). Bayes’ theorem providesus with a calculation prescription P ( θ | d , m ) = P ( d | θ, m ) P ( θ | m ) P ( d | m ) , (1)where P ( d | θ, m ) is the likelihood of the data, P ( θ | m ) is the param-eter prior, and P ( d | m ) is the model’s evidence. The latter guaran-tees the posterior’s normalization and is given by Z = P ( d | m ) = (cid:90) Ω θ P ( d | θ, m ) P ( θ | m )d θ. (2)For parameter estimation with one model, the evidence can beneglected, hence it is su ffi cient to maximize the product of thelikelihood and the prior. However, for comparing di ff erent mod-els, e.g., m and m , one needs normalized posteriors to form theratio R = P ( m | d ) P ( m | d ) = P ( d | m ) P ( m ) P ( d | m ) P ( m ) = Z P ( m ) Z P ( m ) . (3)Often there is no strong a priori reason for preferring one modelover the other which corresponds to setting the model prior ratio P ( m ) / P ( m ) to unity. In this case, the model’s evidence is theonly source of information for model selection.
3. Galactic Variance
The likelihood P ( d | θ, m ) describes the probability to measure thedata d if reality was given by θ and m . By modeling the physicalsystem this probability can be explicitly calculated for certainsets ( θ , m ). For this, one uses a forward simulation code to com-pute observables like sky-maps of Faraday rotation, synchrotronemission, and thermal dust emission. Given measured data, bymodeling the noise characteristics of the detector, a probabilitycan be assigned to the calculated maps, which is in principle astandard approach. However, when analyzing parametric mod-els of the GMF one must be careful at this step because of howthose models describe small scale structure of the magnetic field.Generally speaking, parametric models specify the large scale structure of the magnetic field explicitly by parameterizing thegeometry of its components – for example, the disk and possi-bly its arms, the halo, X-shaped components, et cetera – and thefield strength therein. Together, these components form the so-called regular field. Small scale structure, in contrast, is modeledin terms of its statistical properties rather than an explicit realiza-tion. This means, that when for a given parameter set θ a modelinstance is created, a random magnetic field is generated andadded to the regular field. Depending on the model, the randommagnetic field obeys, for example, a certain power spectrum, islocally proportional to the regular field, or shows a certain degreeof anisotropy. As a consequence, the set ( θ , m ) corresponds notonly to one, but rather infinitely many possible field realizations.For the calculation of a likelihood this means that the measuredobservables must be compared with the ensemble average, whichin practice is the simulated mean of a yet finite set of observablerealizations that result from the magnetic field realizations. Intheory one can work out the e ff ect of various types of randomfields on the used observables; for example, the total intensity ofsynchrotron emission does not depend on whether the structureof the magnetic field is ordered or completely random. Hence,one could use fudge factors to calculate the observable’s meandirectly without having to create numerous samples. However,to do a proper uncertainty quantification one must not neglectthe so-called Galactic variance , a term introduced in Ja ff e et al.(2010). This variance measures how strong the influence of therandom magnetic field on the individual pixels of an observable’ssky-map is. Regions where the influence is high, that is wherethe observable’s variance is high, must be down-weighted whenbeing compared to measured data, in contrast to regions werethe randomness of the magnetic field has little influence on theobservable’s randomness. This makes it again necessary to cal-culate instances of ( θ , m ) to be able to construct an estimate forthe Galactic variance. See Sec. 4.6 for details.
4. The I magine
Framework
As mentioned in Sec. 1, the number of available GMF modelsand the abundance and quality of observational data are continu-ously increasing. The goal of I magine is to provide scientists witha standardized framework to analyze the probability distributionsof model parameters based on physical observables. In doing so,Bayesian statistics is used to judge the mismatch between mea-sured data and model prediction. It is important to note that I mag - ine ’s inference is not limited to magnetic field models. Rather,I magine creates an instance of the Milky Way based on a set ofparameters. It is irrelevant for the framework whether the param-eters are controlling the appearance of the GMF or, for example,the properties of the free electron density or the dust density.Nevertheless, for the time being, we focus on the GMF and keepall other components fixed.It is desirable to have a flexible and open framework avail-able when doing parameter inference. The magnetic field in par-ticular must be analyzed indirectly via observables like syn-chrotron emission, Faraday rotation, dust absorption, or thermaldust emission since there is no direct detection method. This im-plies that the inference depends on the assumptions that weremade regarding further constituents of the Milky Way, for ex-ample the free electron density, the population of cosmic rays,or the dust density. Hence, it is very likely that once the self-consistent analysis of a magnetic field model is finished, newinsights regarding one or more other components make it neces-sary to redo the calculations with the new set-up. An example forthis is the NE2001 model for the Galaxy’s free electron density Article number, page 2 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE
IMAGINE
NIFTy Hammurabi PyMultiNestMultiNestHEALPixNumPy FFTW3D2O
Fig. 1: The building blocks of the I magine framework.
SamplerPrior RepositoryPipeline SampleGalaxy-GeneratorGalaxy-Instance Observable-Generator ObservablesLikelihood
Fig. 2: The structure of the I magine data processing and interpretation.(Cordes 2004). Today, updated versions like the YMW16 model(Yao et al. 2017) are available, and it would be very interestingto update parameter estimates from the past. In practice, eitherthis does not happen at all or only with a huge time delay; in-ference pipelines are usually not made public and the originatormay not have the necessary resources anymore. A standardizedand open inference framework can help here to speed up scien-tific progress and make scientific results more transparent.I magine is built on the programming language P ython to en-sure flexibility, and several external libraries for numerical ef-ficiency, cf. Fig. 1. Here, P ython is primarily used as glue toconnect individual components and external libraries. A strictlyobject-oriented design makes it easy to extend its functionalityfrom existing base-classes. The configuration of the inferenceruns is also done in P ython . No configuration files are used asthe needs for future derived custom classes can not be foreseentoday. Instead, the scientist instantiates the individual compo-nents in the main P ython script which are ultimately embracedby the I magine -pipeline.
The structure of I magine is shown in Fig. 2 and discussed here.The
Pipeline object plays the key-role as it embraces all otherobjects and orchestrates their function calls. Its partner is the
Sampler , with a functional interface for likelihood evaluations. The pipeline hides physical units and scales from the sampler.This means that the former exposes the latter N variables rang-ing from 0 to 1, each. In this way, the sampler can operate verygenerically on this unit cube [0 . . . N without the need to knowany internal details on the Galaxy models.The likelihood evaluation inside the pipeline consists of thefollowing steps. The Sampler yields a point from [0 . . . N .Hence, first, the Galaxy-Generator maps those variables tophysical parameters. Note that N does not need to be the fullnumber of all parameters a model has. All parameters that arenot marked as active in the Pipeline are set to their individ-ually configurable default value. The
Galaxy-Generator thenuses these parameters to generate a certain
Galaxy model re-alization. This means to set up all constituents of the abstractGalaxy model including, e.g., the regular and the random mag-netic field, the thermal electron density field, the dust-densityfield, et cetera. Next, the
Observable-Generator , for exam-ple H ammurabi (Waelkens et al. 2009), processes the Galaxy in-stance and computes physical
Observables , like sky-maps ofthe Faraday depth, synchrotron emission, or thermal dust emis-sion. Those simulated quantities are then compared with mea-sured data by the
Likelihood , which in turn consists of sub-likelihoods for the individual observables. Together, parame-ters, Galaxy model, observables and likelihood values form a
Sample . Finally, the pipeline can be configured to store those
Samples in a repository for post-processing and caching before
Article number, page 3 of 25 & A proofs: manuscript no. IMAGINE the likelihood value is returned to the sampler. Together with theprior, the sampler can then determine which variable configura-tion should be evaluated next.As described in Sec. 3 the GMF models may consist of a ran-dom field component to model the small scale structure stochas-tically. To deal with the resulting Galactic variance, instead ofa single simulation, a set of realizations is created for a certainparameterization. The members of that set are processed in par-allel by the
Observable-Generator such that horizontal scal-ing, i.e. using multiple computers as a cluster, can be exploitedto compensate for the massively increased computational costsone has compared to approaches which ignore the Galactic vari-ance. For this purpose, the I magine framework uses the softwarepackages NIFT y ffi cient data par-allelization, respectively. D2O is based on the Message Pass-ing Interface standard (MPI) (Message Passing Interface Forum1994, 1998) and in particular on mpi4py (Dalcín et al. 2005). Incombination with OpenMP threading (Dagum & Menon 1998)of the
Observable-Generator and the accompanying verticalscaling, I magine e ffi ciently exploits the parallel architecture of amodern high performance computing cluster as a whole as wellas its nodes. The goal of the I magine framework is to provide deep proba-bilistic insights into Galaxy models given observational data. Be-cause of the complexity of the problem, it is not su ffi cient to cal-culate point estimates like a maximum a-posteriori approxima-tion. We expect very counter intuitive interdependencies amongthe model parameters and hence need a thorough uncertaintyquantification in order to correctly interpret the observations.To achieve this, I magine uses Markov Chain Monte Carlo(MCMC) methods as described by Gelman et al. (2014). As de-picted in Sec. 2, we seek to perform parameter estimation for agiven model as well as Bayesian hypothesis testing when com-paring models. Because of its modularity, the I magine frameworkcan easily make use of the full arsenal of the Bayesian method-ology, since it is straightforward to plug in di ff erent MCMC li-braries and to write interfaces for new ones.Over the years, various sampling methods based on MCMChave been created. In the following, we briefly discuss the con-cepts of Metropolis-Hastings , Hamiltonian Monte Carlo and
Nested sampling . The Metropolis-Hastings (MH) algorithm (Metropolis et al.1953; Hastings 1970) creates a biased random walk through theparameter space. If the random walk is ergodic and its tran-sition probabilities obey detailed balance, P ( x → x (cid:48) ) P ( x ) = P ( x (cid:48) → x ) P ( x (cid:48) ), the samples generated by the random walk fol-low the probability distribution P ( x ). Typically, this is achievedby combining a suggestion step with symmetric transition prob-abilities between any pair of locations from an unbiased ran-dom walk with a rejection step that ensures detailed balance, P accept = min (cid:110) , P ( x proposed ) / P ( x old ) (cid:111) .During the walk the samples in the chain must decorrelatefrom the starting position. Hence, the e ffi ciency of an MCMCalgorithm is crucial. Choosing a small step length for that pur-pose indeed means a lower rejection ratio. However, because ofthe small steps the chain does not move. In contrast, a large step length yields a high rejection ratio and therefore a chain that doesnot move, either. This relationship gets worse with higher dimen-sions. An approach to achieve high acceptance rates is Hamilto-nian Monte Carlo sampling. Hamiltonian Monte Carlo (HMC) sampling (also known as Hy-brid Monte Carlo sampling) is a unique MCMC algorithm thatintroduces an auxiliary Gaussian random variable p of the samedimensionality as the original parameters x , cf. Brooks et al.(2011); Betancourt (2017).The auxiliary variable plays the role of a momentum, theoriginal parameters the role of a position in equations of mo-tion from Hamiltonian mechanics. The negative log-probabilitycorresponds to an energy. A new position in parameter space ofposition and momentum is generated by integrating the Hamil-tonian equations of motion in time. This new position is thentreated as the result of a proposal step in the sense of the MHalgorithm. Since the Newtonian equations of motions conserveenergy the proposed parameters should be accepted 100% of thetime, while at the same time being far away from the initial pa-rameters to ensure decorrelation of x . This makes HMC sam-pling much more e ffi cient in exploring the parameter space thanMH sampling.Although this makes an HMC sampler move much fasterthan an ordinary MH sampler it has a downside: it requires thegradient field of the desired probability density function (PDF).Especially when dealing with a high number of dimensions,this can pose a problem if finite di ff erencing must be used forgradient computation. Furthermore, some GMF models exhibitdiscontinuities that result in non-smooth likelihood landscapes,which makes gradients even more problematic. Hence, the I mag - ine pipeline primarily uses nested sampling which does not re-quire gradient information and allows for model comparison, cf.Sec. 2, too. Nested sampling is an MCMC method developed by Skilling(2006), that is capable of directly estimating the relation betweenthe likelihood function and the prior mass. It is unique in the factthat nested sampling is specifically made for usage in Bayesianproblems, giving the evidence as its primary result instead of theposterior probability.Nested sampling works with a set of live-points . In each iter-ation, the point that has the lowest likelihood value gets replacedby a new one with a higher likelihood value. As this methodprogresses, the new points sample a smaller and smaller priorvolume. The algorithm thus traverses through nested shells ofthe likelihood.
There are many parametric field models in the literature,from relatively simple axisymmetric spirals to complex multi-component models. In addition to defining the parametrizedstructure of a magnetic field model, estimates for the values ofthose parameters must be made. Usually, the term model is usedfor both the analytical structure of the magnetic field and for acertain parameter fit. Note that in the context of I magine , model refers to the analytical structure only, since the goal is to inves-tigate its parameter space. It is more straightforward to denote Article number, page 4 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE two samples from the same parameter space as belonging to thesame model instead of constituting distinct models themselves,especially when doing Bayesian model comparison.In addition to the models’ intrinsic complexities, the analysesin the literature also vary with respect to how many observablesand datasets were used in the optimization. An example for arather simple magnetic field model that was fitted to only oneobservable is the WMAP logarithmic-spiral-arm (LSA) model(Page et al. 2007). In a Galacto-centric cylindrical frame thisregular GMF model is given as B ( r , φ, z ) = B (cid:104) sin( ψ ) cos( χ ) ˆ r + cos( ψ ) cos( χ ) ˆ φ + sin( χ ) ˆ z (cid:105) , (4) ψ = ψ + ψ ln (cid:32) rR (cid:33) ,χ = χ tanh (cid:32) zz (cid:33) , where ψ represents the pitch angle of the magnetic field spiralarm which varies according to ψ and a logarithmic dependencyon the radial distance r . R is the distance between the Galacticcenter and the Sun, and ψ defines the local regular field orien-tation. The parameter χ corresponds to the o ff -disk tilting of theGalactic field, and z characterizes the vertical scale height ofthe poloidal field strength modulation. This simple LSA modelfor the coherent field was fitted to synchrotron polarization dataat 23 GHz by Page et al. (2007). Since the observable intensityof the synchrotron radiation depends on both B and the cosmicray electron (CRE) density in a degenerate way, only the otherthree parameters were fitted.At the more complicated end is the Jansson & Farrar (2012)model (JF12 hereafter) with dozens of parameters describing in-dependent spiral arm segments for regular and random fields andthin and thick disks, an X-shaped halo, and more. JF12 was op-timized against both Faraday rotation measures (RM) and syn-chrotron total and polarized intensity. The model of Ja ff e et al.(2013) (and references therein, Ja ff e13 hereafter) is in betweenin terms of number of parameters, with fewer fitted parameterscompared to JF12 though originally optimized against the sameobservables.Some analyses in the literature include only a coherent fieldcomponent, while some additionally study the random compo-nent from the turbulent ISM in a variety of ways. The JF12model includes an analytic expression for the average amountof each observable that would result from the given turbulencemodel. Ja ff e13 is notable in that it uniquely includes the e ff ect ofthe Galactic variance described in Sec. 3 explicitly in the likeli-hood. That analysis used a set of numerical realizations of eachmodel to quantify not only the average amount of emission butalso its variations for a given point in parameter space, which isa necessary step for an unbiased likelihood analysis as describedin Sec. 4.6.A further complication to this sort of analysis is how to treatthe anisotropy in the random component. As described in Ja ff eet al. (2010), from an observational point of view, the GMF canbe divided into three components: coherent, isotropic random,and a third variously called the ordered random, the anisotropicrandom, or the striated component. This third component is ex-pected to arise in the turbulent ISM due to both shocks andshears on large scales. The JF12 model includes a scalar fudge-factor to adjust the synchrotron polarization amplitude from thecoherent field to estimate this striated component. In contrast,Ja ff e13 explicitly models it by projecting the numerically sim-ulated isotropic random component onto the coherent compo-nent to generate an additional anisotropic component. These are complementary methods to model phenomenologically the ef-fect of anisotropic, turbulent, magnetohydrodynamical processesthat are computationally expensive to model physically.On an abstract level, the regular and random components of amagnetic field model are independent. Because of this, I magine distinguishes them such that the user can combine any regularwith any random field model. This is made possible not leastthrough recent developments related to I magine ’s primary ob-servable generator H ammurabi . The H ammurabi code (Waelkens et al. 2009) was built for sim-ulating Galactic polarized foreground emission, absorption, andpolarization rotation. Its core functionality is to produce 2D ob-servables in terms of HEALP ix maps (Górski et al. 2005) basedon 3D physical field configurations in the Galaxy, e.g., the mag-netic, cosmic ray and free electron fields. To analyze various dif-ferent models, H ammurabi is able to construct physical fieldsboth analytically and numerically. Both regular and randomfields covering Galactic scales can be generated with built-infield generators. The observables are produced through line-of-sight integration, including synchrotron and polarized dust emis-sion, Faraday depth, and dispersion measure. In the course of theintegration, radiative transfer and polarization rotation are evalu-ated by accumulating absorption and rotation e ff ects backwardsfrom the observer to the emitter. Technically speaking, the line-of-sight integration is conducted on a set of nested HEALP ix shells. Given R as the maximum simulation radius, the n th shellout of N total shells covers the radial distance from 2 ( n − N − R to 2 ( n − N ) R , except for the first shell which starts at the observer.The angular resolution in each shell is set by HEALP ix ’s N side parameter. The n th shell is by default set up with N side = ( n − M ,where M represents the lowest simulation resolution at the firstshell. Accumulation of observables among shells is carried outby standard HEALP ix interpolation. Within each shell, physicalquantities are estimated from inside out on discrete radial bins,where the radial bin number is proportional to the radial thick-ness of the corresponding shell. Since the observables and thephysical fields are constructed and evaluated in di ff erent coordi-nate frames, a trilinear interpolation method is used to retrieveinformation from the physical fields during the line-of-sight in-tegration. While exploring a magnetic field model’s parameter space, thelikelihood must be evaluated very often. Hence, H ammurabi and especially its random field generator must be swift to pre-serve computational feasibility. To accomplish I magine ’s scien-tific goals, H ammurabi was recently redesigned; the new versionis called H ammurabi X hereafter.In addition to numerous small to medium sized improve-ments, H ammurabi X provides two novel solutions for randommagnetic field configurations on global, i.e. Galactic, and lo-cal, i.e. Solar neighborhood scales, respectively. In the case ofglobal field generation, the focus lies on computational e ffi -ciency. Hence, a triple Fourier transform approach is used to doanisotropy enforcement, field strength rescaling and divergencecleaning. For a given power spectrum, P ( k ), a random magneticfield, ˜ B ( k ), is created in the harmonic Fourier base. The first http: // healpix.sourceforge.net https: // bitbucket.org / hammurabicode / hamxArticle number, page 5 of 25 & A proofs: manuscript no. IMAGINE
Fourier transform translates ˜ B ( k ) into the spatial domain B ( x ).There, anisotropy that may depend on the alignment of the reg-ular magnetic field is introduced. Additionally, a template fieldstrength scaling can be included in terms of a function S ( x ) as B ( x ) → B ( x ) (cid:112) S ( x ). (5)An example for such a scaling function is S ( x ) = S ( r , φ, z ) = exp (cid:32) − rh r − | z | h z (cid:33) , (6)where h r and h z are the characteristic scales of the radial and ver-tical profiles, respectively. The second Fourier transform trans-lates the re-profiled field B ( x ) back into harmonic space, wherea Gram-Schmidt procedure is used to clean up the divergence: ˜ B → ˜ B − ( k · ˜ B ) k / k (cid:12)(cid:12)(cid:12) ˜ B − ( k · ˜ B ) k / k (cid:12)(cid:12)(cid:12) | ˜ B | . (7)Finally, a last Fourier transform is applied to retrieve the desired B ( x ). Hence, the anisotropic random magnetic field is drawnfrom a one-dimensional power spectrum which in contrast cor-responds to statistical homogeneity and isotropy. Breaking theisotropy with subsequent divergence cleaning results in a fieldthat does not precisely obey the original power spectrum P ( k )anymore.In contrast to the global method, for local scale simulations astrict method including vector decomposition of the power spec-trum tensor is available in H ammurabi X. This method is notprone to the inaccuracies described above. Details with respectto the local field generator are beyond the scope of this paper butare available in the release publication of H ammurabi
X (Wanget al., in prep.).
Magnetic fields cannot be measured directly. Instead, their prop-erties need to be inferred indirectly via observables (also re-ferred to as tracers ). The most commonly used observables in-clude Faraday rotation, synchrotron radiation, dust absorptionand emission to probe properties of the GMF, as well as dis-persion measure to probe the thermal electron density. These ob-servables are briefly described below.
Faraday rotation can be described as a double refraction e ff ectwhen linearly polarized light travels through a magnetized, ion-ized medium. The polarization angle of the Faraday rotation isgiven by θ = θ + Φ λ , (8)with θ being the observed polarization angle, θ the original po-larization angle, Φ the Faraday depth and λ the wavelength of thelight ray. The Faraday depth is given by a line-of-sight integralover a distance l to an observer, Φ rad m − = . (cid:90) l n e ( l )cm − B (cid:107) ( l ) µ G d l pc , (9)with n e ( l ) and B (cid:107) ( l ) being the thermal electron density andstrength of the parallel magnetic field, respectively, at distance l away from the observer. Φ is positive (negative) when the mag-netic field is pointing towards (away from) the observer by con-vention. Assuming the emitted polarization angle θ is constantfor a specific source, the Faraday depth gives information aboutthe average strength of the line-of-sight (i.e., parallel) compo-nent of the magnetic field. The synchrotron radiation that is used for the GMF inferenceis caused by the acceleration of relativistic electrons within thisvery magnetic field. This linearly polarized electromagnetic ra-diation is emitted radially to the acceleration. Its intensity isgiven by I s ∝ N ( E ) B x ⊥ , (10)with N ( E ) being the density of relativistic electrons in the rele-vant energy range, E . The index x depends on the energy spec-trum of these electrons, typically x ≈ .
8. Even though the inten-sity of synchrotron radiation is degenerate with other emissioncomponents, like free-free and spinning dust in the microwaveband, Stokes Q and U still provide information regarding themagnetic field. The other components are assumed to be un-polarized. The random components of the GMF depolarize thesynchrotron radiation; see the classic paper by Burn (1966). Thestrength of this depolarization depends on the degree of orderingin the field, which can be written as B ⊥ , r / B ⊥ with B ⊥ , r being theregular part of B ⊥ . Using the Stokes I, Q, and U together, wecan calculate the strength of the magnetic field perpendicular tothe line-of-sight B ⊥ (using the intensity I ) and the fraction of thetotal magnetic field that is regular B ⊥ , r / B ⊥ (using the polarizedintensity PI ). This makes it a useful tool for studying the ran-dom component of magnetic fields. In addition, the lines-of-sightfor an extended source with a per se constant polarization angletraverse space with a di ff erent field configuration each. This re-sults in varying polarization angles within the instrument beam,known as Faraday beam depolarization which provides furtherinformation. Starlight polarization is caused by rotating dust grains absorbingcertain polarizations of light. In a magnetic field, a dust graintends to align its long axis perpendicular to the direction of thelocal magnetic field (see Davis & Greenstein (1951) and refer-ences therein). If the field is perpendicular to the line-of-sight,certain polarizations of the light-ray get blocked, viz. dust ab-sorption of background starlight. The resulting observed light-ray is thus polarized, which gives information about the directionof the magnetic field perpendicular to the line-of-sight betweenthe observer and the star.The approach above works well for low-density dust clouds.In high-density dust clouds, the probability that a light-ray getscompletely absorbed along the way is fairly high. However, dustheats up if it absorbs a lot of radiation, which in return will be re-emitted in the infrared. This emitted infrared light is also polar-ized according to the dust grain’s geometry, viz. polarized ther-mal dust emission. Since as already mentioned the dust grainsare aligned in the magnetic field, the polarized dust emissionprovides complementary information about the direction of B ⊥ . Article number, page 6 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE
When a neutron star forms in the course of a supernova collapsethe preserved angular momentum causes the neutron star to ro-tate rapidly. Along the neutron star’s magnetic axis, a highly fo-cused beam of radiation is emitted, and the rotational and mag-netic axes are not necessarily the same. Since the beam is highlyfocused, from an observer’s point of view this may result in ablinking pattern, which is why those stars are called pulsars. Thegroup and phase velocity of the emitted radiation are not thesame in the interstellar medium because of its ionized compo-nents, mainly free electrons. Because of this, higher frequenciesarrive earlier than lower ones. This extra time delay added at afrequency ν is given by t ( ν ) = e π m e c DM ν , (11)with DM being the so-called dispersion measure. The DM itselfis given by the line-of-sight integral,DM = (cid:90) l n e ( l ) d l . (12)If one has information on the thermal electron density, the DMsolely depends on the distance l between the source and theobserver.The dispersion measure, although it does not give any in-formation on magnetic field properties, is still a very importantobservable. With DM, the thermal electron density can be in-ferred, which in turn is needed for the inference of Faraday ro-tation, as described in Ekers et al. (1969). Using a combinationof Faraday rotation, synchrotron radiation, starlight polarizationand dispersion measure data is key for inferring the constituentsof the Galaxy. The likelihood is the probability P ( d | θ, m ) to obtain the data d from a measurement under the assumption that reality is givenby the model m that in turn is configured by the parameters θ .It is the key element to rate the probability of a stochastic sam-ple. Assuming the generic case of a measurement with linearresponse function R of a signal s which involves additive noise n , the corresponding equation for the data d reads d = R ( s ) + n . (13)If the measurement device is assumed to exhibit Gaussian noisecharacteristics with a covariance matrix N , i.e. n ← (cid:45) G ( n , N ) = | π N | / exp (cid:32) − n † N − n (cid:33) (14)the log-likelihood for a simulated signal that is the result of theevaluation of a model m with parameters θ , i.e. s (cid:48) = m ( θ ), to haveproduced the measured data d is L ( d | s (cid:48) ) = − (cid:0) d − R ( s (cid:48) ) (cid:1) † N − (cid:0) d − R ( s (cid:48) ) (cid:1) −
12 ln ( | N | ) . (15)In the context of I magine , as discussed in Sec. 3, the GMFmodels posses random components that are described by ( m , θ )only stochastically. Marginalizing over those random degrees offreedom results in a modification of the e ff ective covariance termin Eq. (15), namely that the Galactic variance must be addedto the data’s noise covariance. During the further discussion weconsider the following quantities: – The individual GMF samples within an ensemble of size N ens are named B i , with i ∈ [1 , N ens ]. – The process of creating observables from B i is encoded inthe response R . – The simulated observables are denoted by c i = R ( B i ). – The measured observable’s data is named d .Denoting furthermore the data’s noise covariance by A , theGalactic covariance by C , and the dimensionality of observablesby N dim the log-likelihood reads L ( d | c ) = −
12 ( d − ¯ c ) † ( A + C ) − ( d − ¯ c ) −
12 ln ( | A + C | ) (16)with the ensemble mean of c ¯ c = N ens N ens (cid:88) i = c i . (17)As discussed in Sec. 3 the Galactic covariance C reflects the factthat the observables posses an intrinsic variance because of therandom parts of the GMF. For example, the higher the intrin-sic variance, the more the likelihood will be flattened by the( A + C ) − term. This means that the likelihood is less respon-sive to deviations from the ensemble mean for regions of highvariance. Hence, there is the risk of overestimating random fieldcontributions, since they are favored by the likelihood. How-ever, this is compensated by the second summand in Eq. (16):the covariance matrix’ log-determinant ln ( | A + C | ). In Eq. (15)the covariance matrix and thus its determinant are constant andtherefore can be neglected as we are not interested in the abso-lute scales of the likelihood. In contrast, for Eq. (16) we haveto consider it as this determinant varies from point to point inparameter space.The Galactic covariance C is not known, hence, we must es-timate it. A classic approach for C is to evaluate the dyadic prod-uct of the samples’ deviations from their mean: C cl = N dim N ens N ens (cid:88) i = ( c i − ¯ c )( c i − ¯ c ) † = N ens N ens (cid:88) i = u i u i † (18)with u i = (cid:112) N dim (cid:16) c i − ¯ c (cid:17) . (19)Since the number of samples in an ensemble is much smallerthan the number of dimensions this classical estimator for the co-variance matrix is insu ffi cient. Most of its eigenvalues are zero,making an operator-inversion impossible. Hence, it is better touse a sophisticated estimator using a shrinkage target (e.g., adiagonal matrix) and a shrinkage factor. Here, we use the Or-acle Approximating Shrinkage (OAS) estimator by Chen et al.(2011): C = µρ + (1 − ρ ) C cl . (20)The specific quantities needed to compute the OAS estimator are µ = N dim tr ( C cl ) = N dim N ens N ens (cid:88) i = u i † u i (21) a = tr (cid:16) C † cl C cl (cid:17) = N N ens (cid:88) i = N ens (cid:88) j = (cid:18) u i † u j (cid:19) (22) r = min , (1 − / N dim ) a + N µ ( N ens + − / N dim ) (cid:0) a − N dim µ (cid:1) . (23) Article number, page 7 of 25 & A proofs: manuscript no. IMAGINE
In the likelihood one needs to apply the inverse of the sum of A and C , ( A + C ) − . Since we do not know a basis in which A + C isdiagonal, the inversion of this operator is a nontrivial task. How-ever, because of its structure, we can use the Sherman-Morrison-Woodbury matrix identity (Sherman & Morrison 1950; Wood-bury 1950) by re-sorting A + C = ( A + µ r ) + (1 − r ) C cl = B + VV † (24)with B = A + µ r and V = (cid:114) − rN ens U . (25)Namely,( B + VV † ) − = B − − B − V ( + V † B − V ) − V † B − . (26)With this formula only a matrix of size N instead of N mustbe inverted.For computing the log-determinant ln ( | A + C | ) one could usethe result of the OAS estimator and apply the generalized form ofthe matrix determinant lemma (Harville 2008) to it. Its structureis closely related to the Sherman-Morrison-Woodburry matrixidentity: it turns the problem into the calculation of the deter-minant of a matrix of size N instead of N . For our case itreads: | A + C | = (cid:12)(cid:12)(cid:12) B + VV † (cid:12)(cid:12)(cid:12) = | B | · (cid:12)(cid:12)(cid:12) + V † B − V (cid:12)(cid:12)(cid:12) (27)However, the OAS estimator has been designed for and is good atapproximating covariance matrices in terms of quadratic forms;using it for determinant estimation yields rather poor results.And in fact, it can be shown that it is not possible to constructa general purpose estimator from covariance matrix samples ifthe number of samples is lower than the number of dimensions(Cai et al. 2015). Nevertheless, heuristic as well as Bayesian es-timators have been developed trying to cover special cases, as forexample the case of sparse or diagonally dominated covariancematrices (Fitzsimons et al. 2017; Hu et al. 2017). For the timebeing we approximate the determinant | A + C | by its diagonal:ln ( | A + C | ) ≈ N dim tr ln A + N ens N ens (cid:88) i = (cid:16) c i − ¯ c (cid:17) . (28)This approximation serves the purpose of regularizing the ran-dom magnetic field strength. Future improvements could includethe usage of one of the widely used shrinkage estimators as dis-cussed in Hu et al. (2017). They work similarly to the OAS esti-mator, though exhibiting shrinkage coe ffi cients and targets tailormade for covariance determinant approximation. For those, thenEq. (27) can be used for e ffi cient computation. In either case,the inversion of the covariance matrix as well as the calcula-tion of its determinant can be done explicitly, if approximately,which therefore allows us to evaluate the ensemble likelihood inEq. (16) e ffi ciently.
5. Application
In the following we discuss possible usage scenarios of the I mag - ine pipeline. Regardless of parameter estimation or model com-parison, first, a Galaxy model must be set up. Below we will useI magine to analyze the following scenario. Our Galaxy modelconsists of the WMAP logarithmic-spiral-arm (LSA) magneticfield model (Page et al. 2007) in combination with an isotropic Gaussian random field as described in Sec. 4.4.1. In H am - murabi X, the random field’s normalization is chosen such thatits RMS field strength at the Sun’s position is given by τ . We de-note the spectral index of the random field’s power spectrum as α . Furthermore, we choose the YMW16 model (Yao et al. 2017)for the thermal electron density. Here, our goal is to infer theparameters of the magnetic field model, so the thermal electrondensity we assume to be fixed. For the input data, we considerpolarized synchrotron emission at 1 .
41 GHz (Stokes Q and U )following Wolleben et al. (2006), 408 MHz (Stokes I ) and at30 GHz (Stokes Q and U ) following (Planck Collaboration et al.2016a), and the Faraday depth map following Oppermann et al.(2012). It is advisable, before starting a large likelihood exploration,to check if the chosen observables (tracers) are sensitive to themodel parameters that are about to be inferred. In principle, allobservables used here are sensitive to the GMF configuration es-pecially near the Solar neighborhood. In terms of the WMAPLSA model, the influence of ψ on polarized synchrotron emis-sion is expected to be the most noticeable feature, cf. Fig. 3. Bydefinition of the model, ψ has greater influence than ψ and χ on the field’s configuration when r < R / e . We therefore expectthe observables to be more sensitive to ψ at low Galactic lat-itudes where line-of-sight integration accumulates informationthrough the Galactic center. However, Faraday depolarization atlow Galactic latitudes and low frequencies diminishes constrain-ing power of polarized synchrotron emission on ψ . During the development of the IMAGINE framework, initialmock data tests were performed on the base of the JF12 modeland are described in van der Velden (2017). One result of thatwork was an increased appreciation for the di ffi culty workingwith such a complex model. The first ( ω ) is solely a regularWMAP field, while the second ( ω ) additionally possesses a ran-dom component as described in Sec. 4.4.1. To test the pipelinewith a parameter set that is as realistic as possible, we used thebest fit estimates for the WMAP LSA model given in Page et al.(2007), except for ψ which would be 0 . ◦ . To conduct propertests it is helpful if the mock data generating parameter valuesare not located at the boundaries of parameter space, so we set ψ to 7 . ◦ . Since B is not given in Page et al. (2007) we useBeck & Krause (2005) as a reference and set it to 6 µ G. Further-more, for ω we set the random magnetic field’s strength aroundthe Sun τ to 2 µ G. The spectral index is set to α = . ≈ / B = . ∈ [0 . , . µ G (29) χ = . ∈ [1 . , . ◦ (30) ψ = . ∈ [6 . , . ◦ (31) ψ = . ∈ [0 , . ◦ (32) τ = . ∈ [0 . , . µ G (33) α = . ∈ [0 . , .
2] (34)After processing the mock magnetic fields with H am - murabi X, we add individual random noise samples with thevariances given in Wolleben et al. (2006); Planck Collaborationet al. (2016a); Oppermann et al. (2012) to the calculated observ-
Article number, page 8 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE
Fig. 3: Simulated synchrotron emission di ff erence maps (including 408 MHz total intensity T tot at northern hemisphere, 30 GHzStokes Q and polarized intensity T pol in mK) with di ff erent ψ or χ settings. ψ has influence mainly along Galactic longitude while χ a ff ects more the latitude direction. T tot ( = 10 o ) T tot ( = 40 o ) -500 500 Q ( = 10 o ) Q ( = 40 o ) -0.005 0.005 T pol ( = 10 o ) T pol ( = 40 o ) T tot ( = 10 o ) T tot ( = 40 o ) -500 500 Q ( = 10 o ) Q ( = 40 o ) -0.005 0.005 T pol ( = 10 o ) T pol ( = 40 o ) ables. For the Oppermann Faraday depth map there is an un-certainty map available which is based on a Bayesian Wienerfilter reconstruction. Since the pixel-wise noise is uncorrelatedon small scales, we downscale the uncertainty map to N side = N side =
32 as well, no further adaptionof this noise map is necessary. For the Planck and Wolleben syn-chrotron (Stokes Q and U in each case) we take a constant statis-tical uncertainty of 2 . µ K (Planck Collaboration et al. 2016b,Tab. 10) and 12 mK (Wolleben et al. 2006, section 5.2). For theStokes I map at 408 MHz an uncertainty map is given. We down-grade all four data sets to our simulation resolution of N side = N ens =
64. Our tests showed that for the resolution N side =
32 this is theensemble size where the classical covariance term in the ensem-ble likelihood becomes dominant over the shrinkage target, i.e. r falls below 0 .
5, cf. Sec. 4.6. To make likelihood maximizationand sampling possible, it is also necessary to stabilize the likeli-hood by fixing the ensemble member’s random seed. This intro-duces a bias, which we found, however, to be already negligiblein the case of N ens =
64 compared to the emerging Galactic vari-ance. In the future, one could try to enhance existing samplingtechniques already including simulated annealing (Kirkpatricket al. 1983) to become capable of treating the noisy likelihoodsurface directly.
First, we consider the first mock data set ω that does not con-tain random field components. For this data set we perform one-dimensional likelihood scans through the parameter space, asthis is a systematic way to check the observables’ sensitivity withrespect to the model parameters. In doing so, we vary only oneparameter at a time while keeping all others fixed to the mockdata’s generating values.Fig. 4 shows how well the di ff erent observables yield peaksin the likelihood. Since there is no random magnetic field, theensemble likelihood simplifies to a standard χ likelihood. Sev-eral comments are in order. First, one sees that the total log- Table 1: Log-likelihood maximizing parameter values for mockdata ω inferred with a Nelder-Mead optimizer; showing the firstsignificant digit of deviation. B [ µ G] ψ [ ◦ ] ψ [ ◦ ] χ [ ◦ ]Mock values 6 . . .
95 25 . .
999 26 .
99 7 .
943 25 . true mock datavalues for all four WMAP LSA parameters. Second, as expected, B shows the strongest dependence, followed by ψ and χ ; ψ a ff ects the observables as well but much more weakly than theother three parameters. Third, it is remarkable that for all pa-rameters the total log-likelihood is dominated by synchrotronemission Stokes Q & U at 30 GHz and Stokes I at 408 MHz.Faraday rotation also adds some information, but synchrotrondata at 1 .
41 GHz yields four to six orders of magnitude weakersignals in the log-likelihood. This is because of the signal-to-noise ration which is better for the Planck than for the Wollebendata set. Furthermore, due to the depolarization e ff ects that havea huge impact on low-frequency polarized synchrotron data, wesee sharp peaks for 1 .
41 GHz synchrotron data, as the morphol-ogy of the observable map tremendously changes when varyingthe magnetic field. If the GMF were regular, this would allow usto constrain the GMF parameters very precisely. However, thepresence of random magnetic fields and Faraday depolarizatione ff ects render this frequency uninformative for this analysis. Wetherefore exclude the 1 .
41 GHz data from the subsequent analy-sis. After examining the one-dimensional parameter scans, wethen check whether it is possible to infer the input parametersfrom the mock data set with simple minimization. Tab. 1 showsthe values a Nelder-Mead minimizer (Nelder & Mead 1965)yields when operating with the mock data set ω . As mentionedabove, only Faraday depth and synchrotron data at 408 MHz and30 GHz were used according to the insights we drew from theparameter scans. Article number, page 9 of 25 & A proofs: manuscript no. IMAGINE
The minimizer is able to reliably find the correct parame-ter values, which suggests that the likelihood surface is well-behaved throughout the parameter space volume and not onlyalong the optimum-intersecting axes. Note that in general itis advisable to use a gradient-free minimization scheme likeNelder-Mead due to possible non-smooth transitions that are par-ticularly part of more complex magnetic field models.Finally, we use P y M ulti N est (Buchner et al. 2014) to ex-plore the likelihood surface of the mock data ω . Fig. 5 showsthe marginalized probability density functions as well as pair-wise correlation plots. As expected, B is inferred with the high-est precision; followed by ψ and χ , and finally ψ . In the courseof this, the addition of mock noise causes the inferred parametermeans to be shifted with respect to the true values. Di ff erent ran-dom seeds for the noise yield varying o ff sets. The likelihood isinsensitive to these deviations, however, since we knew the truenoise covariance matrix and take it into account. With the highsignal-to-noise ratios, the 2 σ intervals are narrow and cover the ω ’s true parameter values. This means that the likelihood is con-sistent with the process of mock data creation and mock noisegeneration. In the following we repeat the steps from the previous sectionfor mock data set ω : scanning the parameter space, finding op-timal parameter values with Nelder-Mead minimization and do-ing a full sampling with P y M ulti N est . Fig. 6 shows that includ-ing a random magnetic field reduces the sensitivity of the en-semble likelihood considerably. In contrast to Fig. 4, now thelog-likelihood values vary over one to three instead over eightorders of magnitude. As before, the signal for B is strongest,followed by ψ and χ , and finally ψ . With respect to the param-eters of the random magnetic field component we see that τ , theparameter for the random magnetic field’s strength, and α , therandom field’s spectral index, exhibit a slight peak at their truevalues. However, as foreseen in Sec. 4.6, τ ’s likelihood flattenssignificantly for large values. The fact that for α the likelihoodhas its maximum near the true mock data value is the incidentalresult of combining contrarily biased Faraday rotation and syn-chrotron radiation likelihoods. It should be noted that such shiftsare not unexpected, since the mock data include a single real-ization of the Galactic and noise variance that can cause suchchance alignment with slightly shifted parameters. Finally, wenote that Faraday rotation data would not be able to constrain τ and α reasonably. The total likelihood’s shape around the truemock data value is rather flat for τ and α . Their influence on thelikelihood is comparably small in this mock scenario, but wouldincrease with the strength of the random field component; herethe setting is B = µ G vs. τ = µ G. ψ , τ , and α get traced bythe observables – at least slightly – which is why we keep themfor the further inference. At this point the importance of this sen-sitivity analysis becomes evident, as we can draw the followingconclusions: If we find a parameter which has completely negli-gible or even misleading influence on the likelihood it should beexcluded from inference. It would solely increase the dimension-ality of the problem and with respect to minimizers and samplersbehave in the best case as a noisy contribution and therefore dis-turb convergence.For completeness, we visualize the importance of the regu-larizing determinant in Eq. (16). Fig. 7 shows that without the de-terminant the ensemble likelihood favors too high random fieldstrengths and spectral indices. As in Sec. 5.1.2, we continue by inferring the parameter val-ues of ω using a Nelder-Mead minimizer. Tab. 2 shows the re-sults of the optimization whereby we see that as expected theaccuracy is significantly lower than without a random magneticfield component, cf. Tab. 1. One can also see the trend that τ getsoverestimated, which already became apparent in the parameterscan, cf. Fig. 6.Finally, Fig. 8 shows the marginal plots based on a P y M ulti -N est run on the mock data set ω . First, we recognize that theGalactic variance caused rather broad uncertainties. Neverthe-less, the uncertainty intervals are highly reasonable: for exam-ple, although the sample mean value for B lies rather preciselyat 6 µ G, the maximum likelihood value is significantly shiftedto the right. Furthermore, one sees that the Galactic variancewashes out almost all predictive power on ψ . The fact that thesample mean matches the mock data’s generating parameter ismainly due to the fact that the true value is at the center of theprior volume. As seen before, τ gets overestimated, while 2 µ Gstill lies within the 2 σ interval. Interestingly enough, looking atthe joint probability density plot for τ and α , one sees that forlarger τ also larger α become more likely. This is on the onehand an indicator for an unsurprising degeneracy between thetotal strength and the spectral index. On the other hand, we ex-pect that the predictive power on one of the parameters can beincreased by fixing the other by the use of strong prior informa-tion.Note, that a naive χ likelihood which, unlike the ensem-ble likelihood, does not reflect the Galactic variance massivelyunderestimates the uncertainties introduced by the random mag-netic field. Fig. 9 shows the result of P y M ulti N est maximizinga χ likelihood on the ω mock data set. The spectral index ispushed to a small value of α = .
202 making the random mag-netic field rather white. As a consequence, in the ensemble meanthe influence of the random magnetic field maximally cancelsout as the set of samples in the ensemble is finite. The otherparameters then heavily over-fit the variations in the mock-datawhich come from its specific random magnetic field realization.This illustrates the importance of taking the Galactic varianceinto account when doing model parameter inference.
One strength of sampling methods like M ulti N est is that theyproduce an estimate for the evidence . As discussed in Sec. 2,the evidence is crucial for model selection. To illustrate the pro-cedure, we set up the following scenario: Given the prevailingmock data set ω , we compare two models that are both trivialversions of the WMAP LSA model. The only free parameter isnow B . For model M the values for the hidden parameters areequal to those of the mock data, while for model M they arefixed to ψ = . ◦ , ψ = . ◦ , and χ = . ◦ . Tab. 3 shows thatthe log-evidence for M is significantly higher than for M , cor-responding to a massive Bayes factor of R = . · . Butbesides the quality of fit the evidence also takes the model’scomplexity into account. The fewer parameters a model has,the smaller is its total parameter space volume. Hence, even ifa rather complicated model has a better best-fit estimate thana simpler one, if over-fitting occurs its evidence value will beworse. Tab. 3 also shows the log-evidence for the full four-parameter WMAP LSA model ( M ). The log-evidence for M lies in between those of M and M ; the Bayes factor between M and M is R = .
18, which means that there is substantial ev-idence that M is more likely (Je ff reys 1998). Thus, one sees thepenalty coming from M ’s larger parameter space volume com- Article number, page 10 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE
Table 2: Log-likelihood maximizing parameter values for mock data ω inferred with a Nelder-Mead optimizer; showing twosignificant digits of deviation. B [ µ G] ψ [ ◦ ] ψ [ ◦ ] χ [ ◦ ] τ [ µ G] α Mock values 6 . . .
95 25 . . . .
063 26 .
86 8 .
38 23 . .
11 1 . M . However, the improvements of a better parameterfit may compensate for this penalty as the comparison with M shows. In Sec. 5.1, we verified that the I magine pipeline produces self-consistent results for the WMAP LSA model in combinationwith the chosen observables. Now we analyze the likelihoodstructure for the real synchrotron data at 408 MHz and 30 GHz(Planck Collaboration et al. 2016a), and the Faraday depth data(Oppermann et al. 2012) Generally, it is advisable to thoroughlyprepare the input data by masking regions in the sky obviouslyperturbed by local phenomena, for example supernova remnants.However, for this paper this is beyond the scope, as the goalis to illustrate the concepts behind I magine rather than produc-ing high-precision estimates. First, we use the synchrotron datato constrain the parameters of the purely ordered WMAP LSAmodel; so far no random fields are included in the magnetic fieldnor in the likelihood. The result is shown in Fig. 10.The resulting uncertainties are quantitatively consistent withthe mock-data results, cf. Fig. 5. Qualitatively speaking, B isdetermined with the highest accuracy, followed by ψ and χ ,and finally ψ . Also the inferred magnetic field strength B = . µ G is of a reasonable order of magnitude (Ruiz-Granadoset al. 2010; Han 2006). ψ = . ◦ lies within the wide range ofestimates one can find in literature, e.g., 8 ◦ (Beck 2001; Han2006) and 35 ◦ (Page et al. 2007). However, ψ = . ◦ and χ = . ◦ are far o ff from the best-fit values given in Pageet al. (2007), namely ψ = . ◦ and χ = ◦ . This, in combina-tion with the very small uncertainties, indicates that an inferenceneglecting the influence of random components in the magneticfield as well as the likelihood is making matters too easy.When using Faraday depth instead of synchrotron data (Op-permann et al. 2012) for a parameter fit, cf. Fig. 11, the lim-ited capabilities of the WMAP LSA model become clear. Eventhough the WMAP LSA model was designed for fitting syn-chrotron radiation but not Faraday depth data, it is neverthelessremarkable how incompatible they are. Not only are the esti-mates for ψ = ◦ and ψ = ◦ far o ff their reference val-ues, the magnetic field strength is pushed to values near zero.The latter indicates a general incompatibility between the modeland the data. Furthermore, the best fit value for B is negative,which in our case means that the direction of the magnetic fieldis reversed compared to Page et al. (2007). Fig. 12 illustrates theissue. In the data, one can locate a dipole as well as a quadrupolemoment both being aligned with the Galactic plane. Because ofits simple structure, the WMAP LSA model cannot account forthe double anti-axisymmetric quadrupole structure. That is ex-pected, and if such a feature is needed one can use more complex For using I magine in production it is advisable to use the raw datacompiled by Oppermann et al. (2012) as this ensures that there is noalteration of the noise information by a Wiener filter. models like JF12 (Jansson & Farrar 2012) or Ja ff e13 (Ja ff e et al.2013). But, beyond this, Fig. 12 in combination with Fig. 13 re-veals that the likelihood peak at ψ = . ◦ corresponds to aconfiguration where the model exhibits field reversals to fit thestructure in the Galactic plane. Although the Faraday rotationmap compares well to the data, such a parameter configurationis the result of a simple model fitted to a complicated datasetand is not necessarily the most physically realistic solution. Thisdemonstrates another possible pitfall and also how important itis to incorporate physical priors for the model parameters whendoing a real-life analysis. I magine provides the structure for com-prehensive studies that are not only built on powerful algorithmssuch as M ulti N est that will find parameter estimates in any casebut also regularizes them and points out problems in the recon-struction. Furthermore, irrespective of the quadrupole, for thereference parameter values also the dipole does not fit; it has thewrong sign. This means that the overall field orientation itself inthe WMAP LSA model cannot be correct. This is a fact that doesnot become apparent when solely using synchrotron data, sinceeven though synchrotron emission is sensitive to the magneticfield’s direction, it is not to its orientation. Using the I magine pipeline for parameter estimation, it is economic to include var-ious data sets from di ff erent observables, since the I magine datarepository is open and will grow through collaborative contribu-tion. Such obvious contradictions can then be avoided by a morecomprehensive approach.Going a step further, we try to find parameter estimates forthe WMAP LSA plus random magnetic field model that we pre-viously used for the mock data tests in Sec. 5.1. The results areshown in Fig. 14. With respect to the random magnetic field, welimit ourselves to the inference of τ , the strength of the randommagnetic field. For the sampling we used a wide prior volume,especially for the angular parameters ψ , ψ and χ ∈ [0 , ◦ ,each. Note that only ψ is a truly circular parameter, cf. Eq. (4),and thus was setup as such in P y M ulti N est . Comparing the re-sults of a fit based purely on synchrotron emission, given inFig. 14, to the scenario in Fig. 10 with only an ordered fieldprovides several insights. First of all, we see how approachesthat neglect the Galactic variance tremendously underestimateuncertainties when doing parameter estimation. The estimate for B is smeared out the least, but the predictive power for ψ , ψ and χ disappears when taking the Galactic variance into accountcorrectly. In the light of the above, it is noteworthy how clear theprediction for the strength of the random magnetic field turns outto be. All in all, despite their weak predictive power, the resultsshown in Fig. 14 are consistent with those in Fig. 10. For the for-mer, the regular magnetic field strength B is smaller comparedto the latter as now the random magnetic field also contains mag-netic field power in τ . Also note the reasonable anti-correlationbetween B and τ . The overall order of magnitude of τ is com-patible with Han (2006). Since ψ is a circular parameter it isnecessary to consider circular definitions of mean and standard-deviation (Watson 1983), which yield ψ = . ± . ◦ . Hence, Article number, page 11 of 25 & A proofs: manuscript no. IMAGINE
Table 3: Sample mean and log-evidence values from P y M ulti N est for di ff erent WMAP LSA plus random magnetic field modelsbased on the mock data set ω . At M all four WMAP parameters are flexible; at M and M only B is adjustable. An asterisk ( ∗ )indicates a fixed value. In any case, the random magnetic field’s parameters were kept at their mock data’s default value τ = µ G,and α = . B [ µ G] ψ [ ◦ ] ψ [ ◦ ] χ [ ◦ ]Mock values 6 . . .
95 25 . M . ± .
20 6 . ± .
253 26 . ± .
42 8 . ± .
39 23 . ± . M . ± .
18 6 . ± .
212 27 . ∗ . ∗ . ∗ M − . ± .
15 6 . ± .
232 3 . ∗ . ∗ . ∗ the estimate for ψ points towards the same order of magnitudeas ψ = . ◦ shown in Fig. 10. Furthermore, for χ we seea peak around 80 ◦ which can be interpreted to correspond tothe previous best fit value χ = . ◦ . The second peak around χ = ◦ is less clear and is likely to be a morphological de-generacy. As synchrotron emission is sensitive to the magneticfield’s direction but not to its orientation, considering Eq. (4),we expect a di ff use degeneracy in χ , which gets disturbed bythe factor tanh ( z / z ).Repeating this analysis for Faraday rotation data from Op-permann et al. (2012), results shown in Fig. 15, underlines whathas been seen in Fig. 11. The WMAP LSA model is inherentlyincompatible to Faraday rotation observations: B is pushed tovalues near zero and an increasing τ solely broadens B ’s like-lihood as discussed in Sec. 4.6 but does not add anything to theintrinsic quality of fit. The likelihoods for ψ , ψ and χ don’tposses any clear peaks nor pairwise correlations.All in all, this simple example of the WMAP LSA modelaugmented with a random magnetic field already illustrates howchallenging it is to create models of the constituents of theGalaxy with consistent geometry and to find reliable estimatesfor their parameters. While the example in this section was ratheracademic due to the model’s simplicity, the presented steps sim-ilarly apply when analyzing more complex models such as JF12and Ja ff e13.
6. Conclusion & Outlook
In this paper we presented I magine , a framework for GMF modelparameter inference. We have discussed the motivation behindBayesian parameter inference and model comparison as well asthe importance of the
Galactic variance . We then described themodular structure and extensibility of the I magine framework. Itsmost important building blocks are: – state-of-the-art parametric GMF models, – a varied set of complementary observables, – the new and improved H ammurabi X simulator, and – the di ff erent sampling algorithms that can be used withinI magine .In Sec. 5, we showed with mock data that the pipeline worksself-consistently, we illustrated the concept of Bayesian modelcomparison, and we then applied the pipeline to real data. In thecourse of this, we showed the importance of multi-observablebased parameter fitting. This analysis was, however, a simpleproof-of-concept to demonstrate the capabilities of the I magine pipeline. Now, more sophisticated analyses are in order to gainas much scientific insight from existing data sets and GMF mod-els as possible. Since I magine is uniquely suited to handle therandom component of the GMF and its uncertainties correctly, itcan be adapted into a powerful tool to study the turbulent ISM by, e.g., adding a structure function analysis to the likelihood inorder to constrain the turbulent spectral index. All those insightsshould be used to build improved models and to keep the models’best-fit parameter estimates up-to-date with respect to the everimproving data. Within this paper we solely inferred the parame-ters of the GMF while keeping the thermal electron density fixed.With an extended list of observables (e.g., the DM), more infor-mative datasets, and better models, we can extend this work to ajoint inference of the magnetic field, the thermal electron density,the cosmic ray population, and even the dust model parameters.The I magine pipeline is ready to help tackle this challenge and isavailable at: https: // gitlab.mpcdf.mpg.de / ift / IMAGINE
Acknowledgements.
We thank François Boulanger, Martin Reinecke, Luiz F. S.Rodrigues, and Anvar Shukurov for fruitful discussions and valuable sugges-tions. Part of this work was supported by the
Studienstiftung des deutschenVolkes . The original concept of I magine arose from two International Team meet-ings hosted by the International Space Science Institute in Bern. We also ac-knowledge support and hospitality of the Lorentz Center in Leiden, where theI magine project was further discussed and refined. We acknowledge the supportby the DFG Cluster of Excellence ”Origin and Structure of the Universe”. Thecomputations have been carried out on the computing facilities of the Computa-tional Center for Particle and Astrophysics (C2PAP) and the Radboud University,Nijmegen, respectively. This research has been partly supported by the DFG Re-search Unit 1254 and has made use of the NASA / IPAC Infrared Science Archive,which is operated by the Jet Propulsion Laboratory, California Institute of Tech-nology, under contract with the National Aeronautics and Space Administration.Some of the results in this paper have been derived using the HEALP ix Górskiet al. (2005) package. The corner plots where made using the corner P ython package (Foreman-Mackey 2016). References
Beck, R. 2001, Space Science Reviews, 99, 243Beck, R. & Krause, M. 2005, Astron. Nachr., 326, 414Betancourt, M. 2017, ArXiv e-prints [ arXiv:1701.02434 ]Brooks, S., Gelman, A., Jones, G., & Meng, X. 2011, Handbook of MarkovChain Monte Carlo, Chapman & Hall / CRC Handbooks of Modern StatisticalMethods (CRC Press)Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125Burn, B. J. 1966, MNRAS, 133, 67Cai, T. T., Liang, T., & Zhou, H. H. 2015, Journal of Multivariate Analysis, 137,161Chen, Y., Wiesel, A., & Hero, A. O. 2011, IEEE Transactions on Signal Process-ing, 59, 4097Cordes, J. M. 2004, in Astronomical Society of the Pacific Conference Series,Vol. 317, Milky Way Surveys: The Structure and Evolution of our Galaxy, ed.D. Clemens, R. Shah, & T. Brainerd, 211Dagum, L. & Menon, R. 1998, Computational Science & Engineering, IEEE, 5,46Dalcín, L., Paz, R., & Storti, M. 2005, Journal of Parallel and Distributed Com-puting, 65, 1108Davis, Jr., L. & Greenstein, J. L. 1951, ApJ, 114, 206Ekers, R. D., Lequeux, J., Mo ff et, A. T., & Seielstad, G. A. 1969, ApJ, 156, L21Fitzsimons, J., Cutajar, K., Osborne, M., Roberts, S., & Filippone, M. 2017,ArXiv e-prints [ arXiv:1704.01445 ] http: // / teams / bayesianmodel / http: // / lc / web / / / info.php3?wsid = Foreman-Mackey, D. 2016, The Journal of Open Source Software, 24Gelman, A., Carlin, J. B., Stern, H. S., et al. 2014, Bayesian data analysis, Vol. 2(CRC press Boca Raton, FL)Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759Han, J. L. 2006, Chinese Journal of Astronomy and Astrophysics Supplement,6, 211Harville, D. 2008, Matrix Algebra From a Statistician’s Perspective (SpringerNew York)Hastings, W. K. 1970, Biometrika, 57, 97Hu, Z., Dong, K., Dai, W., & Tong, T. 2017, The international journal of bio-statistics, 13Ja ff e, T. R., Ferrière, K. M., Banday, A. J., et al. 2013, MNRAS, 431, 683Ja ff e, T. R., Leahy, J. P., Banday, A. J., et al. 2010, MNRAS, 401, 1013Jansson, R. & Farrar, G. R. 2012, ApJL, 761, L11Je ff reys, H. 1998, The Theory of Probability, Oxford Classic Texts in the Physi-cal Sciences (OUP Oxford)Kirkpatrick, S., Gelatt, C. D., & Vecchi, M. P. 1983, Science, 220, 671Message Passing Interface Forum. 1994, International Journal of SupercomputerApplications, 8, 159Message Passing Interface Forum. 1998, High Performance Computing Appli-cations, 12, 1Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller,E. 1953, J. Chem. Phys., 21, 1087Nelder, J. A. & Mead, R. 1965, Comput. J., 7, 308Oppermann, N., Junklewitz, H., Robbers, G., et al. 2012, A&A, 542, A93Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJS, 170, 335Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016a, A&A, 594, A1Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016b, A&A, 594, A10Ruiz-Granados, B., Rubiño-Martín, J. A., & Battaner, E. 2010, A&A, 522, A73Sherman, J. & Morrison, W. J. 1950, Ann. Math. Statist., 21, 124Skilling, J. 2006, Bayesian Analysis, 1, 833Steininger, T., Dixit, J., Frank, P., et al. 2017, ArXiv e-prints[ arXiv:1708.01073 ]Steininger, T., Greiner, M., Beaujean, F., & Enßlin, T. 2016, J. Big Data, 3, 17van der Velden, E. 2017, Master’s thesis, Radboud University, Nijmegen, TheNetherlandsWaelkens, A., Ja ff e, T., Reinecke, M., Kitaura, F. S., & Enßlin, T. A. 2009, A&A,495, 697Watson, G. 1983, Statistics on spheres, University of Arkansas lecture notes inthe mathematical sciences (Wiley)Wolleben, M., Landecker, T. L., Reich, W., & Wielebinski, R. 2006, A&A, 448,411Woodbury, M. A. 1950, Statistical Research Group, Princeton University, Prince-ton, N. J., Memo. Rep., 4ppYao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29 Article number, page 13 of 25 & A proofs: manuscript no. IMAGINE
Fig. 4:
Mock data ω : Scans through the parameter space of the WMAP LSA regular magnetic field model, including simulatedadditive measurement noise according to Wolleben et al. (2006), (Planck Collaboration et al. 2016a), Oppermann et al. (2012).The mock data input parameter values are at the very center of each abscissa and indicated by the vertical line. Since the Plancksynchrotron data dominates the overall likelihood, we show the total likelihood in the leftmost plot. Note that the log-likelihoodvaries over several orders of magnitude.0 2 4 6 8 10 12 − − − − · B [ µ G] L og - L i k e li hood − − . · B [ µ G] 0 2 4 6 8 − − − − B [ µ G]10 20 30 40 50 − − − · ψ [ ◦ ] L og - L i k e li hood
10 20 30 40 50 − − . · ψ [ ◦ ] 10 20 30 40 50 − − ψ [ ◦ ]0 5 10 15 − − − · ψ [ ◦ ] L og - L i k e li hood − . − − . · ψ [ ◦ ] 0 5 10 15 − − − ψ [ ◦ ]0 10 20 30 40 50 − − − − · χ [ ◦ ] L og - L i k e li hood − − · χ [ ◦ ] 0 10 20 30 40 50 − − χ [ ◦ ]Synchr. 408 MHz ( I ) Synchr. 30 GHz ( Q ) Synchr. 30 GHz ( U ) TotalFaraday depth Synchr. 1 .
41 GHz ( Q ) Synchr. 1 .
41 GHz ( U ) Article number, page 14 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE
Fig. 5:
Mock data ω : Marginalized posterior plots and projected pairwise correlation plots from applying P y M ulti N est to mockdata ω . The dashed lines represent the 16%, 50% and 84% quantiles, respectively. The parameters for the mock data, indicated bythe solid lines, were set to B = . µ G, ψ = . ◦ , ψ = . ◦ , and χ = ◦ . B = 5.9998 +0.00010.0001 . . . . +2.7e1 = 26.9960 +0.00240.0020 . . . . = 7.9539 +0.00820.0067 . . . . B . . . . +2.5e1 . . . . +2.7e1 . . . . . . . . +2.5e1 = 24.9977 +0.00280.0038 Article number, page 15 of 25 & A proofs: manuscript no. IMAGINE
Fig. 6:
Mock data ω : Scans through the parameter space of the WMAP LSA regular plus isotropic random magnetic field model,including simulated additive measurement noise according to (Planck Collaboration et al. 2016a), Oppermann et al. (2012). Themock data input parameter values are indicated by the vertical line. Note that the plot of the total likelihood does not include thesynchrotron data at 1 .
41 GHz.2 4 6 8 10 12 − − B [ µ G] L og - L i k e li hood − − B [ µ G] 4 6 8 10 − − B [ µ G]10 20 30 40 50 − − − ψ [ ◦ ] L og - L i k e li hood
10 20 30 40 50 − − − ψ [ ◦ ] 10 20 30 40 50 − − ψ [ ◦ ]0 5 10 15 − − · − ψ [ ◦ ] L og - L i k e li hood − − · − ψ [ ◦ ] 0 5 10 15 − . − . − . − . ψ [ ◦ ]0 10 20 30 40 50 − − − χ [ ◦ ] L og - L i k e li hood − . − . − . χ [ ◦ ] 0 10 20 30 40 50 − − − − χ [ ◦ ]Synchr. 408 MHz ( I ) Synchr. 30 GHz ( Q ) Synchr. 30 GHz ( U )Faraday depth Total Article number, page 16 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE − − . − − . τ [ µ G] L og - L i k e li hood − − τ [ µ G] 1 2 3 4 − − − − τ [ µ G]0 1 2 3 − · − α L og - L i k e li hood − . − . α − . − . α Synchr. 408 MHz ( I ) Synchr. 30 GHz ( Q ) Synchr. 30 GHz ( U )Faraday depth Total Article number, page 17 of 25 & A proofs: manuscript no. IMAGINE
Fig. 7:
Mock data ω , without determinant term: Scans through the parameter space of the WMAP LSA regular plus isotropicrandom magnetic field model, including simulated additive measurement noise according to Wolleben et al. (2006), (Planck Collab-oration et al. 2016a), Oppermann et al. (2012). For these plots the ensemble likelihood was evaluated without the determinant term.The mock data input parameter values are indicated by the vertical line. Note that the plot of the total likelihood does not includethe synchrotron data at 1 .
41 GHz.1 2 3 4 − − τ [ µ G] L og - L i k e li hood − − τ [ µ G] 2 3 4 − − τ [ µ G]0 1 2 3 − . α L og - L i k e li hood − . − . − . . α − . − − . . α Synchr. 408 MHz ( I ) Synchr. 30 GHz ( Q ) Synchr. 30 GHz ( U )Faraday depth Total Article number, page 18 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE
Fig. 8:
Mock data ω : Marginalized posterior plots and projected pairwise correlation plots from applying P y M ulti N est to mockdata ω . The dashed lines represent the 16%, 50% and 84% quantiles, respectively. The parameters for the mock data, indicated bythe solid lines, were set to B = . µ G, ψ = . ◦ , ψ = . ◦ , χ = ◦ , τ = µ G, and α = . B = 6.0255 +0.28570.2986 = 26.8177 +4.40704.4468 = 8.1895 +4.86985.2569 = 23.2077 +8.54458.2916 . . . . = 3.0017 +0.51090.6159 . . . . . B . . . . . . . . . . . . . . = 2.0444 +0.76541.1499 Article number, page 19 of 25 & A proofs: manuscript no. IMAGINE
Fig. 9:
Mock data ω in combination with a χ likelihood: Marginalized posterior plots and projected pairwise correlation plotsfrom applying P y M ulti N est to mock data ω using a simple χ likelihood that not reflects the influence of the Galactic variance.The dashed lines represent the 16%, 50% and 84% quantiles, respectively. The parameters for the mock data, indicated by the solidlines, were set to B = . µ G, ψ = . ◦ , ψ = . ◦ , χ = ◦ , τ = µ G, and α = . B = 6.1402 +0.00000.0001 . . . . = 26.8561 +0.01250.0004 . . . . . = 7.0040 +0.00150.0395 .
42 3 .
22 4 .
02 4 .
82 5 . = 23.2948 +0.00140.0575 . . . . = 2.1122 +0.00310.0001 . . . . . B . . . . . . . . . . . . . . . . . . . . . . . . . . . . = 0.2018 +0.00000.0007 Article number, page 20 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE
Fig. 10:
Synchrotron data, purely ordered WMAP LSA magnetic field
Marginalized posterior plots and projected pairwisecorrelation plots from applying P y M ulti N est to 408 MHz and 30 GHz synchrotron data from Planck Collaboration et al. (2016a).The dashed lines represent the 16%, 50% and 84% quantiles, respectively. B = 4.4392 +0.00020.0002 . . . . = 6.6900 +0.00740.0076 .
23 4 .
43 4 .
63 4 . = 34.4352 +0.01910.0193 . . . . B . . . . . . . . . . . . . . . . = 78.5798 +0.00790.0081 Article number, page 21 of 25 & A proofs: manuscript no. IMAGINE
Fig. 11:
Faraday rotation data, purely ordered WMAP LSA magnetic field
Marginalized posterior plots and projected pairwisecorrelation plots from applying P y M ulti N est to Faraday rotation data from Oppermann et al. (2012). The dashed lines representthe 16%, 50% and 84% quantiles, respectively. B = 0.5110 +0.00730.0069 = 280.4312 +0.67690.6401 = 330.2488 +3.77553.6653 . . . . B = 30.2961 +0.89340.9208 Article number, page 22 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE
Fig. 12: Comparison of Faraday depth maps in rad / m . -500 500 (a) Map of the Galactic Faraday depth given by Oppermann et al. (2012) in rad / m . -500 500 (b) Map of the Galactic Faraday depth produced with H am - murabi X based on the WMAP LSA model in rad / m . Themodel’s parameters were set to B = . µ G, ψ = . ◦ , ψ = . ◦ , and χ = . ◦ , following Page et al. (2007). -500 500 (c) Map of the Galactic Faraday depth produced with H am - murabi X based on the WMAP LSA model in rad / m . Themodel’s parameters were set to B = − . µ G, ψ = ◦ , ψ = ◦ , and χ = . ◦ . Fig. 13: Streamplot of the WMAP LSA model for B = − . µ G, ψ = ◦ , ψ = ◦ , and χ = . ◦ .
20 15 10 5 0 5 10 15 20201510505101520 y [ k p c ] x-y plane at z=0.00
20 15 10 5 0 5 10 15 20x [kpc] 42024 z [ k p c ] x-z plane at y=0.00
20 15 10 5 0 5 10 15 20y [kpc] 42024 z [ k p c ] y-z plane at x=0.00WMAP LSA Article number, page 23 of 25 & A proofs: manuscript no. IMAGINE
Fig. 14:
Synchrotron radiation data, WMAP LSA plus random magnetic field
Marginalized posterior plots and projectedpairwise correlation plots from applying P y M ulti N est to 408 MHz and 30 GHz synchrotron data from Planck Collaboration et al.(2016a). The dashed lines represent the 16%, 50% and 84% quantiles, respectively. B = 2.1766 +1.61421.4709 = 175.4581 +132.7069131.8077 = 192.8189 +103.5407129.0273 = 202.7614 +99.8222131.5343 . . . . B = 8.9100 +1.83941.2987 Article number, page 24 of 25heo Steininger et al.: Inferring Galactic magnetic field model parameters using IMAGINE
Fig. 15:
Faraday rotation data, WMAP LSA plus random magnetic field
Marginalized posterior plots and projected pairwisecorrelation plots from applying P y M ulti N est to Faraday rotation data from Oppermann et al. (2012). The dashed lines representthe 16%, 50% and 84% quantiles, respectively. B = 0.0741 +2.61992.6590 = 189.2054 +116.6778129.0451 = 195.0638 +109.8785128.3471 = 184.1044 +115.9966120.4723 . . . . . B
481 21 6 8 0 1 6 0 2 4 0 3 2 0 = 9.1467 +4.39344.2122+4.39344.2122